frcfcmnl.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2015   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include "frcfcmnl.h"
#include "gausspoint.h"
#include "mathfem.h"
#include "classfactory.h"
#include "contextioerr.h"
#include "datastream.h"
#include "node.h"

#include "nonlocalmaterialext.h"

#include "../sm/Elements/structuralelement.h"

#include "floatmatrix.h"
#include "floatarray.h"

#include "dynamicinputrecord.h"
#include "unknownnumberingscheme.h"
#include "stressvector.h"
#include "strainvector.h"


namespace oofem {
REGISTER_Material(FRCFCMNL);

FRCFCMNL :: FRCFCMNL(int n, Domain *d) : FRCFCM(n, d), StructuralNonlocalMaterialExtensionInterface(d)
{}


IRResultType
FRCFCMNL :: initializeFrom(InputRecord *ir)
{
    IRResultType result;                // Required by IR_GIVE_FIELD macro

    result = FRCFCM :: initializeFrom(ir);
    if ( result != IRRT_OK ) {
        return result;
    }

    result = StructuralNonlocalMaterialExtensionInterface :: initializeFrom(ir);
    if ( result != IRRT_OK ) {
        return result;
    }

    //    IR_GIVE_FIELD(ir, participAngle, _IFT_FRCFCMNL_participAngle);
    participAngle = 90.;

    return IRRT_OK;
}


void
FRCFCMNL :: giveRealStressVector(FloatArray &answer, GaussPoint *gp,
                                 const FloatArray &totalStrain,
                                 TimeStep *tStep)
{
    // first try conventional way as if no sourrounding cracks existed
    FCMMaterial :: giveRealStressVector(answer, gp, totalStrain, tStep);

    FRCFCMNLStatus *status = static_cast< FRCFCMNLStatus * >( this->giveStatus(gp) );
    int numberOfActiveCracks = status->giveNumberOfTempCracks();

    double sigma_f_nonlocal = 0.;

    FloatArray sigma, crackVec;
    FloatMatrix sigmaL2G, sigmaG2L, crackDirs;
    double sigma_f_local = 0.;
    double crackStrain;

    if ( numberOfActiveCracks == 0 ) {
        sigma_f_nonlocal = this->computeNonlocalStressInFibersInUncracked(gp, tStep);
        sigma_f_nonlocal = max( sigma_f_nonlocal, status->giveFiberStressNL(1) );
        status->setTempFiberStressNL(1, sigma_f_nonlocal);

        return;
    } else {
        sigma = answer; // global stress


        sigmaL2G = status->giveL2GStressVectorTransformationMtrx();


        sigmaG2L = status->giveG2LStressVectorTransformationMtrx();

        sigma.rotatedWith(sigmaG2L, 'n'); // global stress transformed to local stress

        sigma_f_local = 0.;

        for ( int iCrack = 1; iCrack <= ( status->giveNumberOfTempCracks() ); iCrack++ ) {
            // get cracking strain for i-th crack
            crackStrain = status->giveTempCrackStrain(iCrack);

            if ( crackStrain > 0. ) {
                // get local fiber stress
                sigma_f_local = this->computeStressInFibersInCracked(gp, crackStrain, iCrack);
                // set local fiber stress to status
                status->setTempFiberStressLoc(iCrack, sigma_f_local);
            } else {
                // set zero fiber stress
                status->setTempFiberStressLoc(iCrack, 0.);
            }

            // and compare it to stress from surroundings sigma_f_nonlocal

            crackDirs = status->giveCrackDirs();

            crackVec.resize( status->giveMaxNumberOfCracks(gp) );
            crackVec.zero();

            for ( int i = 1; i <= status->giveMaxNumberOfCracks(gp); i++ ) {
                crackVec.at(i) = crackDirs.at(i, iCrack);
            }

            // compute nonlocal stress
            sigma_f_nonlocal = this->computeNonlocalStressInFibers(crackVec, gp, tStep);

            //test !!!!
            sigma_f_nonlocal = max( sigma_f_nonlocal, status->giveFiberStressNL(iCrack) );


            // set nonlocal fiber stress to status
            status->setTempFiberStressNL(iCrack, sigma_f_nonlocal);

            sigma.at(iCrack) += sigma_f_nonlocal;
        } // loop over crack directions

        sigma.rotatedWith(sigmaL2G, 'n'); // modified local stress transformed to global stress

        answer = sigma;
        status->letTempStressVectorBe(sigma);
    }

    return;
}




double
FRCFCMNL :: computeDebondedLength(double delta) {
    double a = 0.;

    if ( this->fiberType == FT_CAF ) { // continuous aligned fibers
        a = sqrt( ( this->Ef * this->Df * delta ) / ( 2. * this->tau_0 * ( 1. + this->eta ) ) );
    } else if (  this->fiberType == FT_SAF ) { // short aligned fibers
        // no snubbing here - debonding length is simply related to pull-out displacement and not to force...

        a = sqrt( ( this->Ef * this->Df * delta ) / ( 2. * this->tau_0 * ( 1. + this->eta ) ) );
        // threshold is related to the fibre length
        a = min( a, this->Lf / ( 2. * ( 1. + this->eta ) ) );
    } else if (  this->fiberType == FT_SRF ) { // short random fibers
        a = sqrt( ( this->Ef * this->Df * delta ) / ( 2. * this->tau_0 * ( 1. + this->eta ) ) );
        // threshold is related to the fibre length
        a = min( a, this->Lf / ( 2. * ( 1. + this->eta ) ) );
    } else {
        OOFEM_ERROR("Unknown fiber type");
    }

    return a;
}


double
FRCFCMNL :: computeDecreaseInFibreStress(double distance, double delta, double targetDebondedLength) {
    // compute decrease in fiber stress due to stress transfer to matrix
    // delta_sigma_f = 4 * tau_eff * Vf * Xca'/ ( Df * ( 1 - Vf))
    // Xca' = ( 4 * x^2 - 4 * x * Lf ) / ( -2 * pi * Lf * lambda )
    // lambda = 4 / ( pi * g )
    // where x is the distance from the crack, in this case the distance between the gauss-points


    double delta_sigma_f = 0.;
    double tau, x;

    if ( this->fiberType == FT_CAF ) { // continuous aligned fibers
        delta_sigma_f = this->Vf * 4. * this->tau_0 * distance / this->Df;
    } else {
        if ( delta < this->w_star / 2. ) {
            tau = this->tau_0;
        } else {
            // maybe delta should be changed to delta_max?!?
            tau = this->computeFiberBond(2. * delta);
        }

        x = min(distance, targetDebondedLength);

        if (  this->fiberType == FT_SAF ) { // short aligned fibers
            // delta_sigma_f is computed here in direction of the fibers, therefore no snubbing and no reduction in fibre volume
            // this expression is derived for x < a (we are in the debonding zone) and for fibers perpendicular to the crack
            delta_sigma_f = this->Vf * 4. * tau * ( this->Lf * x - x * x ) / ( this->Df * this->Lf );
        } else {
            // this expression uses the same hypothesis as in Lu & Leung 2016 (taken from Aveston 1974) that the stress transfer by friction
            // does not continue after x (section where we want to now the stress) = xd (length of the debonded zone)
            delta_sigma_f = this->Vf * 4. * tau * ( this->Lf * x - x * x ) / ( 3 * this->Df * this->Lf );
        }
    }

    return delta_sigma_f;
}


void
FRCFCMNL :: computeElementCentroid(FloatArray &answer, GaussPoint *gp) {
    double A, cx, cy, x_i, x_ii, y_i, y_ii;
    Element *elem;
    elem = gp->giveElement();
    int nnode = elem->giveNumberOfNodes();

    answer.resize(2);
    answer.zero();

    A = 0.;
    cx = cy = 0.;

    if ( ( nnode == 3 ) || ( nnode == 4 ) ) {
        for ( int inod = 1; inod < nnode; inod++ ) {
            x_i = elem->giveNode(inod)->giveCoordinate(1);
            x_ii = elem->giveNode(inod + 1)->giveCoordinate(1);

            y_i = elem->giveNode(inod)->giveCoordinate(2);
            y_ii = elem->giveNode(inod + 1)->giveCoordinate(2);

            A += 0.5 * ( x_i * y_ii - x_ii * y_i );

            cx += ( x_i + x_ii ) * ( x_i * y_ii - x_ii * y_i );
            cy += ( y_i + y_ii ) * ( x_i * y_ii - x_ii * y_i );
        }

        x_i = elem->giveNode(nnode)->giveCoordinate(1);
        x_ii = elem->giveNode(1)->giveCoordinate(1);

        y_i = elem->giveNode(nnode)->giveCoordinate(2);
        y_ii = elem->giveNode(1)->giveCoordinate(2);

        A += 0.5 * ( x_i * y_ii - x_ii * y_i );

        cx += ( x_i + x_ii ) * ( x_i * y_ii - x_ii * y_i );
        cy += ( y_i + y_ii ) * ( x_i * y_ii - x_ii * y_i );
    }

    cx /= ( 6. * A );
    cy /= ( 6. * A );

    answer.at(1) = cx;
    answer.at(2) = cy;
}


bool
FRCFCMNL :: isInElementProjection(GaussPoint *homeGp, GaussPoint *nearGp, int iNlCrack)
{
    FRCFCMNLStatus *nonlocStatus;
    nonlocStatus = static_cast< FRCFCMNLStatus * >( nearGp->giveMaterialStatus() );

    Element *nearElem;
    FloatArray examinedCoords;
    int nnode;
    double x, y, b_min, b_max, b_home, slope, x_min, x_max;

    nearElem = nearGp->giveElement();

    if ( ( fiberType == FT_CAF ) || ( fiberType == FT_SAF ) ) {
        slope = this->orientationVector.at(2) / this->orientationVector.at(1);
    } else {
        slope = nonlocStatus->giveCrackDirs().at(2, iNlCrack) / nonlocStatus->giveCrackDirs().at(1, iNlCrack);
    }

    this->computeElementCentroid(examinedCoords, homeGp);

    nnode = nearElem->giveNumberOfNodes();

    if ( fabs(slope) < 1.e6 ) {
        b_min = b_max = b_home = 0.;

        for ( int inod = 1; inod <= nnode; inod++ ) {
            x = nearElem->giveNode(inod)->giveCoordinate(1);
            y = nearElem->giveNode(inod)->giveCoordinate(2);

            // y = slope * x + b

            if ( inod == 1 ) {
                b_min = y - slope * x;
                b_max = y - slope * x;
            } else {
                b_min = min(b_min, y - slope * x);
                b_max = max(b_max, y - slope * x);
            }
        }

        b_home = examinedCoords.at(2) - slope *examinedCoords.at(1);

        if ( ( b_min < b_home ) && ( b_home < b_max ) ) {
            return true;
        } else {
            return false;
        }
    } else {
        x_min = x_max = 0.;

        for ( int inod = 1; inod < nnode; inod++ ) {
            x = nearElem->giveNode(inod)->giveCoordinate(1);

            if ( inod == 1 ) {
                x_min = x;
                x_max = x;
            } else {
                x_min = min(x_min, x);
                x_max = max(x_max, x);
            }
        }

        if ( ( x_min < examinedCoords.at(1) ) && ( examinedCoords.at(1) < x_max ) ) {
            return true;
        } else {
            return false;
        }
    }

    // happy compiler
    return true;
}



double
FRCFCMNL :: computeNonlocalStressInFibers(const FloatArray &crackVectorHome, GaussPoint *gp, TimeStep *tStep) {
    // using "non-local" approach

    // computes stress in fibers in the direction of "crackVectorHome"
    // this can be used either to check the stress-strength criterion or when computing the stress from strain

    // they key is to find the maximum stress in fibers which
    // - diminishes with the distance from crack
    // - diminishes with increasing angle wrt to its normal


    FRCFCMNLStatus *nonlocStatus;

    this->buildNonlocalPointTable(gp);

    auto *list = this->giveIPIntegrationList(gp); // !

    // two fiber stresses - from left and right part of the crack - rewrite if it turns out that one stress is fully sufficientaveraged
    double sigma_f_left = 0.;
    double sigma_f_right = 0.;

    FloatArray coordsHome, coordsTarget;

    double alpha, k_alpha_Home;

    FloatArray crackVectorTarget;
    FloatArray pointVector;

    bool side;

    double distance, targetCrackOpening, delta, a, targetDebondedLength, theta;
    double sigma_f_iCr, delta_sigma_f;


    this->computeElementCentroid(coordsHome, gp);

    crackVectorTarget.resize( gp->giveMaterial()->giveDomain()->giveNumberOfSpatialDimensions() );
    pointVector = crackVectorTarget;

    for ( auto &lir : *list ) {
        GaussPoint *neargp = lir.nearGp;

        if ( neargp->giveElement()->giveNumber() != gp->giveElement()->giveNumber() ) { // it must not be the same element
            if ( neargp->giveMaterial()->giveNumber() == gp->giveMaterial()->giveNumber() ) { // it must be the same material
                nonlocStatus = static_cast< FRCFCMNLStatus * >( neargp->giveMaterialStatus() );

                // crack(s) exist in the near gp
                if ( nonlocStatus->giveNumberOfTempCracks() > 0 ) {
                    // do just once for all cracks
                    // get coordinates

                    this->computeElementCentroid(coordsTarget, neargp);

                    // compute the distance between evaluated element centers
                    distance = 0.;
                    for ( int i = 1; i <= coordsHome.giveSize(); i++ ) {
                        distance += ( coordsHome.at(i) - coordsTarget.at(i) ) * ( coordsHome.at(i) - coordsTarget.at(i) );
                    }

                    distance = sqrt(distance);

                    // zero the orientation vector between two points
                    pointVector.zero();


                    for ( int iNlCrack = 1; iNlCrack <= nonlocStatus->giveNumberOfTempCracks(); iNlCrack++ ) {
                        // length of the tunnel (zone where the fibers are debonded from the matrix and where friction acts) at the target gauss-point
                        targetCrackOpening = this->computeNormalCrackOpening(neargp, iNlCrack);


                        // pull-out displacement
                        delta = ( targetCrackOpening - fibreActivationOpening ) / 2.;

                        if ( delta > 0. ) {
                            a = this->computeDebondedLength(delta);

                            targetDebondedLength = a;

                            if ( distance < targetDebondedLength ) {
                                if ( this->isInElementProjection(gp, neargp, iNlCrack) ) {
                                    // get stress in fibers bridging iNlCrack - local stress
                                    sigma_f_iCr = nonlocStatus->giveTempFiberStressLoc(iNlCrack);


                                    // get vector from points' coordinates - if has not been done before for the previous crack
                                    if ( pointVector.containsOnlyZeroes() ) {
                                        for ( int i = 1; i <= coordsHome.giveSize(); i++ ) {
                                            pointVector.at(i) = ( coordsHome.at(i) - coordsTarget.at(i) ) / distance;
                                        }
                                    }

                                    // the change in the bridging stress due to pulley and snubbing effects
                                    if ( ( this->fiberType == FT_CAF ) ||  ( this->fiberType == FT_SAF ) ) {
                                        theta = this->computeCrackFibreAngle(neargp, iNlCrack);
                                        sigma_f_iCr /= ( fabs( cos(theta) ) * exp(fabs(theta) * this->f) );
                                    } else if ( this->fiberType == FT_SRF ) {
                                        sigma_f_iCr *= 2. / ( 3. * this->g );
                                    } else {
                                        OOFEM_ERROR("Unknown fiber type");
                                    }

                                    // change in fiber stress due to bond friction
                                    delta_sigma_f = this->computeDecreaseInFibreStress(distance, delta, targetDebondedLength);

                                    sigma_f_iCr -= delta_sigma_f;

                                    // only positive values - can even happen to be negative?
                                    sigma_f_iCr = max(sigma_f_iCr, 0.);

                                    // decide the side
                                    alpha = this->computeAngleBetweenVectors(pointVector, crackVectorHome);

                                    if ( alpha < M_PI / 2. ) {
                                        side = true;
                                    } else {
                                        side = false;
                                    }

                                    if ( this->fiberType == FT_CAF ) {
                                        alpha = this->computeAngleBetweenVectors(crackVectorHome, this->orientationVector);
                                    } else if ( this->fiberType == FT_SAF ) {
                                        alpha = this->computeAngleBetweenVectors(crackVectorHome, this->orientationVector);
                                    } else if ( this->fiberType == FT_SRF ) {
                                        // angle factor - point vector vs. target crack
                                        crackVectorTarget.zero();
                                        for ( int i = 1; i <= coordsHome.giveSize(); i++ ) {
                                            crackVectorTarget.at(i) = nonlocStatus->giveCrackDirs().at(i, iNlCrack);
                                        }

                                        alpha = this->computeAngleBetweenVectors(crackVectorHome, crackVectorTarget);
                                    } else {
                                        OOFEM_ERROR("Unknown fiber type");
                                    }

                                    // necessary correction
                                    if ( alpha > M_PI / 2. ) {
                                        alpha = M_PI - alpha;
                                    }

                                    // add some snubbing effect?
                                    if ( alpha >  participAngle * M_PI / 180. ) {
                                        k_alpha_Home = 0.;
                                    } else {
                                        k_alpha_Home = cos( alpha * (  90. / participAngle ) );
                                    }


                                    // get fiber stress
                                    if ( side ) {
                                        sigma_f_left = max(k_alpha_Home * sigma_f_iCr, sigma_f_left);
                                    } else {
                                        sigma_f_right = max(k_alpha_Home * sigma_f_iCr, sigma_f_right);
                                    }
                                } // is in element projection
                            } // is in the debonded zone
                        } // positive delta
                    } // loop over target crack directions
                } // condition for at least crack
            } // the same material material
        } // not the same element
    } // loop over all gauss-points


    return ( max(sigma_f_left, sigma_f_right) );
    //  return ( sigma_f_left + sigma_f_right );
}


double
FRCFCMNL :: computeNonlocalStressInFibersInUncracked(GaussPoint *gp, TimeStep *tStep)
{
    // using "non-local" approach

    FRCFCMNLStatus *nonlocStatus;

    this->buildNonlocalPointTable(gp);

    auto *list = this->giveIPIntegrationList(gp); // !

    FloatArray coordsHome, coordsTarget;

    double distance, targetCrackOpening, delta, a, targetDebondedLength, theta, sigma_f_iCr, delta_sigma_f;
    double sigma_f = 0.;

    this->computeElementCentroid(coordsHome, gp);


    for ( auto &lir : *list ) {
        GaussPoint *neargp = lir.nearGp;

        if ( neargp->giveElement()->giveNumber() != gp->giveElement()->giveNumber() ) { // it must not be the same element
            if ( neargp->giveMaterial()->giveNumber() == gp->giveMaterial()->giveNumber() ) { // it must be the same material
                nonlocStatus = static_cast< FRCFCMNLStatus * >( neargp->giveMaterialStatus() );

                // crack(s) exist in the near gp
                if ( nonlocStatus->giveNumberOfTempCracks() > 0 ) {
                    // do just once for all cracks
                    // get coordinates

                    this->computeElementCentroid(coordsTarget, neargp);

                    // compute the distance between evaluated element centers
                    distance = 0.;
                    for ( int i = 1; i <= coordsHome.giveSize(); i++ ) {
                        distance += ( coordsHome.at(i) - coordsTarget.at(i) ) * ( coordsHome.at(i) - coordsTarget.at(i) );
                    }

                    distance = sqrt(distance);


                    for ( int iNlCrack = 1; iNlCrack <= nonlocStatus->giveNumberOfTempCracks(); iNlCrack++ ) {
                        // length of the tunnel (zone where the fibers are debonded from the matrix and where friction acts) at the target gauss-point
                        targetCrackOpening = this->computeNormalCrackOpening(neargp, iNlCrack);


                        // pull-out displacement
                        delta = ( targetCrackOpening - fibreActivationOpening ) / 2.;

                        if ( delta > 0. ) {
                            a = this->computeDebondedLength(delta);

                            targetDebondedLength = a;

                            if ( distance < targetDebondedLength ) {
                                if ( this->isInElementProjection(gp, neargp, iNlCrack) ) {
                                    // get stress in fibers bridging iNlCrack - local stress
                                    sigma_f_iCr = nonlocStatus->giveTempFiberStressLoc(iNlCrack);

                                    // the change in the bridging stress due to pulley and snubbing effects
                                    if ( ( this->fiberType == FT_CAF ) ||  ( this->fiberType == FT_SAF ) ) {
                                        theta = this->computeCrackFibreAngle(neargp, iNlCrack);
                                        sigma_f_iCr /= ( fabs( cos(theta) ) * exp(fabs(theta) * this->f) );
                                    } else if ( this->fiberType == FT_SRF ) {
                                        sigma_f_iCr /= this->g;
                                    } else {
                                        OOFEM_ERROR("Unknown fiber type");
                                    }

                                    // change in fiber stress due to bond friction
                                    delta_sigma_f = this->computeDecreaseInFibreStress(distance, delta, targetDebondedLength);

                                    sigma_f_iCr -= delta_sigma_f;

                                    // only positive values - can even happen to be negative?
                                    sigma_f_iCr = max(sigma_f_iCr, 0.);

                                    // get fiber stress
                                    sigma_f = max(sigma_f_iCr, sigma_f);
                                } // is in element projection
                            } // is in the debonded zone
                        } // positive delta
                    } // loop over target crack directions
                } // condition for at least crack
            } // the same material material
        } // not the same element
    } // loop over all gauss-points


    return sigma_f;
}





bool
FRCFCMNL :: isStrengthExceeded(const FloatMatrix &base, GaussPoint *gp, TimeStep *tStep, int iCrack, double trialStress) {
    // evaluates cheaper function than is the nonlocal approach
    if ( !FRCFCM :: isStrengthExceeded(base, gp, tStep, iCrack, trialStress) ) {
        return false;
    }

    FRCFCMNLStatus *status = static_cast< FRCFCMNLStatus * >( this->giveStatus(gp) );

    double sigma_f, sigma_m;

    FloatArray crackVec;

    crackVec.resize( status->giveMaxNumberOfCracks(gp) );
    crackVec.zero();

    for ( int i = 1; i <= status->giveMaxNumberOfCracks(gp); i++ ) {
        crackVec.at(i) = base.at(i, iCrack);
    }

    sigma_f = this->computeNonlocalStressInFibers(crackVec, gp, tStep);

    //test !!!!
    sigma_f = max( sigma_f, status->giveFiberStressNL(iCrack) );

    status->setTempFiberStressNL(iCrack, sigma_f);

    // get neat stress in matrix
    sigma_m = ( trialStress - sigma_f ) / ( 1. - this->Vf );

    // compare matrix strength with computed stress in the matrix

    if ( FCMMaterial :: isStrengthExceeded(base, gp, tStep, iCrack, sigma_m) ) {
        return true;
    } else {
        return false;
    }
}



void
FRCFCMNL :: giveMaterialStiffnessMatrix(FloatMatrix &answer,
                                        MatResponseMode rMode,
                                        GaussPoint *gp,
                                        TimeStep *tStep)
//
// returns effective material stiffness matrix in full form
// for gp stress strain mode
//
{
    if (  rMode == ElasticStiffness ) {
        FCMMaterial :: giveMaterialStiffnessMatrix(answer, rMode, gp, tStep);
    } else if ( rMode == SecantStiffness ) {
        /*FRCFCMNLStatus *status = static_cast< FRCFCMNLStatus * >( this->giveStatus(gp) );
         * int numberOfActiveCracks = status->giveNumberOfTempCracks();
         * double localSigmaF;
         * double NLsigmaF;
         *
         * for (int iCrack = 1; iCrack <= numberOfActiveCracks; iCrack++) {
         *
         * localSigmaF = status->giveTempFiberStressLoc(iCrack);
         * NLsigmaF =  status->giveTempFiberStressNL(iCrack);
         *
         * //      if ( NLsigmaF > localSigmaF ) {
         * if ( NLsigmaF > 0. ) {
         *  FCMMaterial :: giveMaterialStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
         *  return;
         * }
         * }*/

        FCMMaterial :: giveMaterialStiffnessMatrix(answer, rMode, gp, tStep);
    } else { // tangent stiffness
        FCMMaterial :: giveMaterialStiffnessMatrix(answer, rMode, gp, tStep);
    }

    return;
}


double
FRCFCMNL :: computeAngleBetweenVectors(const FloatArray &vec1, const FloatArray &vec2)
{
    // compute angle between two vectors

    if ( vec1.giveSize() != vec2.giveSize() ) {
        OOFEM_ERROR("the length of vec1 and vec2 is not of the same");
    }


    //return 0.;


    double length1 = 0., length2 = 0.;
    double theta = 0.;

    for ( int i = 1; i <= min( vec1.giveSize(), vec2.giveSize() ); i++ ) {
        length1 += pow(vec1.at(i), 2);
        length2 += pow(vec2.at(i), 2);
        theta += vec1.at(i) * vec2.at(i);
    }
    length1 = sqrt(length1);
    length2 = sqrt(length2);

    // normalize
    theta /= ( length1 * length2 );


    // can be exceeded due to truncation error
    if ( theta > 1 ) {
        return 0.;
    } else if ( theta < -1. ) {
        return M_PI;
    }


    return acos(theta);
}


Interface *
FRCFCMNL :: giveInterface(InterfaceType type)
{
    if ( type == NonlocalMaterialExtensionInterfaceType ) {
        return static_cast< StructuralNonlocalMaterialExtensionInterface * >(this);
    } else {
        return NULL;
    }
}



int
FRCFCMNL :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    FRCFCMNLStatus *status = static_cast< FRCFCMNLStatus * >( this->giveStatus(gp) );

    double local = 0.;
    double nl = 0.;

    for ( int i = 1; i <= status->giveMaxNumberOfCracks(gp); i++ ) {
        local = max(status->giveFiberStressLoc(i), local);
        nl = max(status->giveFiberStressNL(i), nl);
    }


    // stress in fibers (per continuum)
    if  ( type == IST_FiberStressLocal ) {
        answer.resize(1);
        answer.at(1) = local;
        return 1;

        // stress in fibers (per continuum)
    } else if  ( type == IST_FiberStressNL ) {
        answer.resize(1);
        answer.at(1) = nl;
        return 1;
    } else {
        return FRCFCM :: giveIPValue(answer, gp, type, tStep);
    }
}



//                      FRC FCM NL STATUS                       ///


FRCFCMNLStatus :: FRCFCMNLStatus(int n, Domain *d, GaussPoint *gp) :
    FRCFCMStatus(n, d, gp), StructuralNonlocalMaterialStatusExtensionInterface(),
    fiberStressLoc(), tempFiberStressLoc(), fiberStressNL(), tempFiberStressNL()
{
    fiberStressLoc.resize(this->nMaxCracks);
    fiberStressLoc.zero();
    tempFiberStressLoc = fiberStressLoc;
    fiberStressNL = fiberStressLoc;
    tempFiberStressNL = fiberStressLoc;
}


FRCFCMNLStatus :: ~FRCFCMNLStatus()
{}



void
FRCFCMNLStatus :: printOutputAt(FILE *file, TimeStep *tStep)
{
    FRCFCMStatus :: printOutputAt(file, tStep);

    fprintf(file, "maxFiberStressLocal: {");
    for ( int i = 1; i <= this->giveMaxNumberOfCracks(gp); i++ ) {
        fprintf( file, " %f", this->giveFiberStressLoc(i) );
    }
    fprintf(file, "}\n");

    fprintf(file, "maxFiberStressNL: {");
    for ( int i = 1; i <= this->giveMaxNumberOfCracks(gp); i++ ) {
        fprintf( file, " %f", this->giveFiberStressNL(i) );
    }
    fprintf(file, "}\n");
}


void
FRCFCMNLStatus :: initTempStatus()
//
// initializes temp variables according to variables form previous equlibrium state.
// builds new crackMap
//
{
    FRCFCMStatus :: initTempStatus();

    this->tempFiberStressLoc = this->fiberStressLoc;
    this->tempFiberStressNL = this->fiberStressNL;
}



void
FRCFCMNLStatus :: updateYourself(TimeStep *tStep)
//
// updates variables (nonTemp variables describing situation at previous equilibrium state)
// after a new equilibrium state has been reached
// temporary variables correspond to newly reched equilibrium.
//
{
    FRCFCMStatus :: updateYourself(tStep);

    this->fiberStressLoc = this->tempFiberStressLoc;
    this->fiberStressNL = this->tempFiberStressNL;
}



contextIOResultType
FRCFCMNLStatus :: saveContext(DataStream &stream, ContextMode mode, void *obj)
//
// saves full information stored in this Status
// no temp variables stored
//
{
    contextIOResultType iores;

    // save parent class status
    if ( ( iores = FRCFCMStatus :: saveContext(stream, mode, obj) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    if ( ( iores = fiberStressLoc.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    if ( ( iores = fiberStressNL.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    return CIO_OK;
}

contextIOResultType
FRCFCMNLStatus :: restoreContext(DataStream &stream, ContextMode mode, void *obj)
//
// restores full information stored in stream to this Status
//
{
    contextIOResultType iores;

    // read parent class status
    if ( ( iores = FRCFCMStatus :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }


    if ( ( iores = fiberStressLoc.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    if ( ( iores = fiberStressNL.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    return CIO_OK; // return succes
}


Interface *
FRCFCMNLStatus :: giveInterface(InterfaceType type)
{
    if ( type == NonlocalMaterialStatusExtensionInterfaceType ) {
        return static_cast< StructuralNonlocalMaterialStatusExtensionInterface * >(this);
    } else {
        return ConcreteFCMStatus :: giveInterface(type);
    }
}
} // end namespace oofem