mplasticmaterial.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   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 "mplasticmaterial.h"
#include "gausspoint.h"
#include "floatmatrix.h"
#include "floatarray.h"
#include "intarray.h"
#include "datastream.h"
#include "contextioerr.h"
 
namespace oofem {
#define YIELD_TOL 1.e-4
#define RES_TOL   1.e-4
#define PLASTIC_MATERIAL_MAX_ITERATIONS 90
 
MPlasticMaterial :: MPlasticMaterial(int n, Domain *d)  : StructuralMaterial(n, d)
{
}
 
 
MPlasticMaterial :: ~MPlasticMaterial()
{
    delete linearElasticMaterial;
}
 
 
bool
MPlasticMaterial :: hasMaterialModeCapability(MaterialMode mode) const
//
// returns whether receiver supports given mode
//
{
    return mode == _3dMat ||
           mode == _1dMat ||
           //<RESTRICTED_SECTION>
           mode == _PlaneStress  ||
           //</RESTRICTED_SECTION>
           mode == _PlaneStrain;
}
 
 
std::unique_ptr<MaterialStatus> 
MPlasticMaterial :: CreateStatus(GaussPoint *gp) const
{
    return std::make_unique<MPlasticMaterialStatus>(gp, this->giveSizeOfReducedHardeningVarsVector(gp));
}
 
 
void
MPlasticMaterial :: giveRealStressVector(FloatArray &answer,
                                         GaussPoint *gp,
                                         const FloatArray &totalStrain,
                                         TimeStep *tStep) const
//
// returns real stress vector in 3d stress space of receiver according to
// previous level of stress and current
// strain increment, the only way, how to correctly update gp records
//
// completely formulated in Reduced stress-strain space
{
    FloatArray strainSpaceHardeningVariables;
    FloatArray fullStressVector;
    FloatArray strainVectorR, plasticStrainVectorR;
    FloatArray helpVec;
    IntArray activeConditionMap(this->nsurf);
    FloatArray gamma;
 
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
 
    this->initTempStatus(gp);
 
    // subtract stress independent part
    // note: eigenStrains (temperature) is not contained in mechanical strain stored in gp
    // therefore it is necessary to subtract always the total eigen strain value
    this->giveStressDependentPartOfStrainVector(strainVectorR, gp, totalStrain,
                                                tStep, VM_Total);
 
    /*
     * stress return algorithm
     */
 
    //if (0) this->cuttingPlaneReturn (fullStressVector, activeConditionMap, gamma, form, gp, totalStrain, tStep);
    //else this->closestPointReturn (fullStressVector, activeConditionMap, gamma, form, gp, totalStrain, tStep);
    if ( rmType == mpm_ClosestPoint ) {
        this->closestPointReturn(fullStressVector, activeConditionMap, gamma, gp, strainVectorR,
                                 plasticStrainVectorR, strainSpaceHardeningVariables, tStep);
    } else {
        this->cuttingPlaneReturn(fullStressVector, activeConditionMap, gamma, gp, strainVectorR,
                                 plasticStrainVectorR, strainSpaceHardeningVariables, tStep);
    }
 
 
    status->letTempStrainVectorBe(totalStrain);
    StructuralMaterial :: giveReducedSymVectorForm( helpVec, fullStressVector, gp->giveMaterialMode() );
    status->letTempStressVectorBe(helpVec);
 
    status->letTempPlasticStrainVectorBe(plasticStrainVectorR);
    status->letTempStrainSpaceHardeningVarsVectorBe(strainSpaceHardeningVariables);
 
    status->setTempGamma(gamma);
    status->setTempActiveConditionMap(activeConditionMap);
 
 
    // update state flag
    int newState, state = status->giveStateFlag();
    bool yieldFlag = false;
    for ( int i = 1; i <= nsurf; i++ ) {
        if ( gamma.at(i) > 0. ) {
            yieldFlag = true;
        }
    }
 
    if ( yieldFlag ) {
        newState = MPlasticMaterialStatus :: PM_Yielding;         // test if plastic loading occur
    }
    // no plastic loading - check for unloading
    else if ( ( state == MPlasticMaterialStatus :: PM_Yielding ) || ( state == MPlasticMaterialStatus :: PM_Unloading ) ) {
        newState = MPlasticMaterialStatus :: PM_Unloading;
    } else {
        newState = MPlasticMaterialStatus :: PM_Elastic;
    }
 
    status->letTempStateFlagBe(newState);
 
    answer = status->giveTempStressVector();
}
 
 
void
MPlasticMaterial :: closestPointReturn(FloatArray &answer,
                                       IntArray &activeConditionMap,
                                       FloatArray &gamma,
                                       GaussPoint *gp,
                                       const FloatArray &totalStrain,
                                       FloatArray &plasticStrainVectorR,
                                       FloatArray &strainSpaceHardeningVariables,
                                       TimeStep *tStep) const
{
    FloatArray fullStressVector;
    FloatArray elasticStrainVectorR;
    FloatArray fullStressSpaceHardeningVars, residualVectorR;
 
    FloatArray helpVector, helpVector2;
    FloatArray dgamma, tempGamma;
    FloatMatrix elasticModuli, hardeningModuli, consistentModuli;
    FloatMatrix elasticModuliInverse, hardeningModuliInverse;
    FloatMatrix helpMtrx;
    FloatMatrix gmat;
    std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf), * yieldGradVecPtr, * loadGradVecPtr;
    FloatArray rhs;
    double yieldValue;
    int nIterations = 0;
    int i, j, strSize, totSize;
    int elastic, restart, actSurf, indx;
    bool yieldConsistency;
 
 
    if ( this->plType == associatedPT ) {
        yieldGradVecPtr = loadGradVecPtr = & yieldGradVec;
    } else {
        yieldGradVecPtr = & yieldGradVec;
        loadGradVecPtr  = & loadGradVec;
    }
 
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
 
    plasticStrainVectorR = status->givePlasticStrainVector();
    strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
 
    dgamma.resize(nsurf);
    dgamma.zero();
    gamma.resize(nsurf);
    gamma.zero();
 
    strSize = totalStrain.giveSize(); // size of reducedStrain Vector
    totSize = strSize + strainSpaceHardeningVariables.giveSize();
 
    // compute elastic moduli and its inverse
    this->computeReducedElasticModuli(elasticModuli, gp, tStep);
    elasticModuliInverse.beInverseOf(elasticModuli);
    // delete elasticModuli;
 
    //
    // Elastic Predictor
    //
 
    elasticStrainVectorR = totalStrain;
    elasticStrainVectorR.subtract(plasticStrainVectorR);
    // stress vector in full form due to computational convinience
    //if (fullStressVector) delete fullStressVector;
    this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
    this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
 
    //
    // check for plastic process
    //
    elastic = 1;
    actSurf = 0;
    activeConditionMap.zero();
    for ( i = 1; i <= this->nsurf; i++ ) {
        yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
        if ( yieldValue > YIELD_TOL ) {
            elastic = 0;
            actSurf += 1;
            activeConditionMap.at(i) = actSurf;
        }
    }
 
    if ( !elastic ) {
        /*
         * plastic multisurface closest-point corrector
         */
 
 
        do {
            do { // restart loop
                elasticStrainVectorR = totalStrain;
                elasticStrainVectorR.subtract(plasticStrainVectorR);
                // stress vector in full form due to computational convinience
                //if (fullStressVector) delete fullStressVector;
                this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
                this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
 
                // compute gradients
                for ( i = 1; i <= nsurf; i++ ) {
                    this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
                }
 
                if ( this->plType == nonassociatedPT ) {
                    for ( i = 1; i <= nsurf; i++ ) {
                        this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
                    }
                }
 
 
                //this->computeGradientVector (gradientVectorR, gp, &fullStressVector, fullStressSpaceHardeningVars);
                this->computeResidualVector(residualVectorR, gp, gamma, activeConditionMap, plasticStrainVectorR,
                                            strainSpaceHardeningVariables, * loadGradVecPtr);
 
                // check convergence
                yieldConsistency = true;
                for ( i = 1; i <= nsurf; i++ ) {
                    if ( activeConditionMap.at(i) ) {
                        yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
                        if ( ( yieldValue >= YIELD_TOL ) ) {
                            yieldConsistency = false;
                        }
                    }
                }
 
                if ( yieldConsistency && ( residualVectorR.computeNorm() < RES_TOL ) ) {
                    answer = fullStressVector;
                    printf(" (%d iterations)", nIterations);
                    return;
                }
 
                // compute  consistent tangent moduli
                this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
                if ( hardeningModuli.giveNumberOfRows() ) {
                    hardeningModuliInverse.beInverseOf(hardeningModuli);
                    // delete hardeningModuli;
                } else {
                    hardeningModuliInverse.clear();
                }
 
                this->computeAlgorithmicModuli(consistentModuli, gp, elasticModuliInverse,
                                               hardeningModuliInverse, gamma, activeConditionMap,
                                               fullStressVector, fullStressSpaceHardeningVars);
 
                rhs.resize(actSurf);
                gmat.resize(actSurf, actSurf);
                for ( i = 1; i <= nsurf; i++ ) {
                    if ( activeConditionMap.at(i) ) {
                        yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
                        helpVector.beTProductOf(consistentModuli, ( * yieldGradVecPtr ) [ i - 1 ]);
 
                        for ( j = 1; j <= this->nsurf; j++ ) {
                            if ( activeConditionMap.at(j) ) {
                                gmat.at(i, j) = ( * loadGradVecPtr ) [ j - 1 ].dotProduct(helpVector);
                            }
                        }
 
                        rhs.at(i) = yieldValue - residualVectorR.dotProduct(helpVector);
                    }
                }
 
                // obtain increment to consistency parameter
                gmat.solveForRhs(rhs, helpVector);
 
                // assign gammas
                for ( i = 1; i <= this->nsurf; i++ ) {
                    if ( ( indx = activeConditionMap.at(i) ) ) {
                        dgamma.at(i) = helpVector.at(indx);
                    } else {
                        dgamma.at(i) = 0.0;
                    }
                }
 
                tempGamma = gamma;
                tempGamma.add(dgamma);
 
                //
                // check active constraint
                //
 
                restart = 0;
                for ( i = 1; i <= this->nsurf; i++ ) {
                    if ( tempGamma.at(i) < 0.0 ) {
                        restart = 1;
                    }
                }
 
                if ( restart ) {
                    // build up new activeConditionMap and restart the iterartion from begining
                    actSurf = 0;
                    activeConditionMap.zero();
                    for ( i = 1; i <= this->nsurf; i++ ) {
                        if ( tempGamma.at(i) > 0.0 ) {
                            actSurf += 1;
                            activeConditionMap.at(i) = actSurf;
                        }
                    }
 
                    continue;
                } else {
                    break;
                }
            } while ( 1 ); // end restart loop
 
            // obtain incremental plastic strains and internal variables
            // residualVectorR vector and gradVec array are changed
            // but they are computed once again when new iteration starts
            for ( i = 1; i <= nsurf; i++ ) {
                if ( activeConditionMap.at(i) ) {
                    ( * loadGradVecPtr ) [ i - 1 ].times( dgamma.at(i) );
                    residualVectorR.add( ( * loadGradVecPtr ) [ i - 1 ] );
                }
            }
 
            helpVector.beProductOf(consistentModuli, residualVectorR);
            this->computeDiagModuli(helpMtrx, gp, elasticModuliInverse, hardeningModuliInverse);
            helpVector2.beProductOf(helpMtrx, helpVector);
 
            // Update state variables and consistency parameter
            for ( i = 1; i <= strSize; i++ ) {
                plasticStrainVectorR.at(i) += helpVector2.at(i);
            }
 
            for ( i = strSize + 1; i <= totSize; i++ ) {
                strainSpaceHardeningVariables.at(i - strSize) += helpVector2.at(i);
            }
 
            gamma.add(dgamma);
 
            // increment iteration count
            nIterations++;
 
            if ( nIterations > PLASTIC_MATERIAL_MAX_ITERATIONS ) {
                OOFEM_WARNING("local equlibrium not reached in %d iterations\nElement %d, gp %d, continuing",
                          PLASTIC_MATERIAL_MAX_ITERATIONS, gp->giveElement()->giveNumber(), gp->giveNumber() );
                answer = fullStressVector;
                break;
            }
        } while ( 1 );
    } else {
        answer = fullStressVector;
    }
}
 
 
void
MPlasticMaterial :: cuttingPlaneReturn(FloatArray &answer,
                                       IntArray &activeConditionMap,
                                       FloatArray &gamma,
                                       GaussPoint *gp,
                                       const FloatArray &totalStrain,
                                       FloatArray &plasticStrainVectorR,
                                       FloatArray &strainSpaceHardeningVariables,
                                       TimeStep *tStep) const
{
    FloatArray elasticStrainVectorR;
    FloatArray fullStressVector, fullStressSpaceHardeningVars;
    FloatArray helpVector, trialStressIncrement, rhs;
    FloatMatrix elasticModuli, hardeningModuli, dmat, gmatInv;
    std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf), * yieldGradVecPtr, * loadGradVecPtr;
    FloatArray dgamma(this->nsurf);
    double yieldValue;
    int nIterations = 0;
    int size, sizeR, i, j, elastic, restart, actSurf, indx, iindx, jindx;
 
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
 
    if ( this->plType == associatedPT ) {
        yieldGradVecPtr = loadGradVecPtr = & yieldGradVec;
    } else {
        yieldGradVecPtr = & yieldGradVec;
        loadGradVecPtr  = & loadGradVec;
    }
 
    // huhu:
    plasticStrainVectorR = status->givePlasticStrainVector();
    strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
 
    dgamma.resize(nsurf);
    dgamma.zero();
    gamma.resize(nsurf);
    gamma.zero();
    // compute elastic moduli
    this->computeReducedElasticModuli(elasticModuli, gp, tStep);
    // compute  hardening moduli and its inverse
    this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
 
    // initialize activeConditionMap
    elasticStrainVectorR = totalStrain;
    elasticStrainVectorR.subtract(plasticStrainVectorR);
    // stress vector in full form due to computational convinience
    this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
    StructuralMaterial :: giveReducedSymVectorForm( trialStressIncrement, fullStressVector, gp->giveMaterialMode() );
    trialStressIncrement.subtract( status->giveStressVector() );
    this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
 
    elastic = 1;
    actSurf = 0;
    activeConditionMap.zero();
    for ( i = 1; i <= this->nsurf; i++ ) {
        yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
        if ( yieldValue > YIELD_TOL ) {
            elastic = 0;
            actSurf += 1;
            activeConditionMap.at(i) = actSurf;
        }
    }
 
    if ( !elastic ) {
        //
        // compute consistent moduli
        //
        sizeR = elasticModuli.giveNumberOfRows();
        size = sizeR + hardeningModuli.giveNumberOfRows();
 
        dmat.resize(size, size);
        dmat.zero();
 
        for ( i = 1; i <= sizeR; i++ ) {
            for ( j = 1; j <= sizeR; j++ ) {
                dmat.at(i, j) = elasticModuli.at(i, j);
            }
        }
 
        for ( i = sizeR + 1; i <= size; i++ ) {
            for ( j = sizeR + 1; j <= size; j++ ) {
                dmat.at(i, j) = hardeningModuli.at(i - sizeR, j - sizeR);
            }
        }
 
        do {            // local equilibrium loop
            // update  hardening moduli and its inverse
            this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
            for ( i = sizeR + 1; i <= size; i++ ) {
                for ( j = sizeR + 1; j <= size; j++ ) {
                    dmat.at(i, j) = hardeningModuli.at(i - sizeR, j - sizeR);
                }
            }
 
            do {              // loop for restart
                gmatInv.resize(actSurf, actSurf);
                gmatInv.zero();
                rhs.resize(actSurf);
                rhs.zero();
 
                // compute gradients
                for ( i = 1; i <= nsurf; i++ ) {
                    this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
                }
 
                if ( this->plType == nonassociatedPT ) {
                    for ( i = 1; i <= nsurf; i++ ) {
                        this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, fullStressSpaceHardeningVars);
                    }
                }
 
 
 
                for ( i = 1; i <= nsurf; i++ ) {
                    if ( ( iindx = activeConditionMap.at(i) ) ) {
                        rhs.at(iindx) = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
                        helpVector.beTProductOf(dmat, ( * yieldGradVecPtr ) [ i - 1 ]);
 
                        for ( j = 1; j <= this->nsurf; j++ ) {
                            if ( ( jindx = activeConditionMap.at(j) ) ) {
                                gmatInv.at(iindx, jindx) = ( * loadGradVecPtr ) [ j - 1 ].dotProduct(helpVector);
                            }
                        }
                    }
                }
 
                /*
                 *                              // compute plastic multipliers for active conds.
                 *      //computee gmatInv
                 *      gradMat.resize(actSurf, size);
                 *      for (i=1; i<=nsurf; i++) {
                 *        if ((iindx = activeConditionMap.at(i))) {
                 *          rhs.at (iindx) = this-> computeYieldValueAt (gp, i, fullStressVector, fullStressSpaceHardeningVars);
                 *          for (j=1; j<=size; j++) {
                 *            gradMat.at(iindx, j) = (*yieldGradVecPtr)[i-1].at(j);
                 *          }
                 *        }
                 *      }
                 *      helpMtrx.resize(actSurf, size);
                 *      helpMtrx.beProductOf (gradMat, dmat);
                 *      gmatInv.beProductTOf (helpMtrx, gradMat);
                 */
                // solve for plastic multiplier increments
                gmatInv.solveForRhs(rhs, helpVector);
                // assign gammas
                for ( i = 1; i <= this->nsurf; i++ ) {
                    if ( ( indx = activeConditionMap.at(i) ) ) {
                        dgamma.at(i) = helpVector.at(indx);
                    } else {
                        dgamma.at(i) = 0.0;
                    }
                }
 
 
 
                restart = 0;
                // check for active conditions
                for ( i = 1; i <= this->nsurf; i++ ) {
                    if ( ( gamma.at(i) + dgamma.at(i) ) < 0.0 ) {
                        gamma.at(i) = 0.0;
                        restart = 1;
                    }
                }
 
                if ( restart ) {
                    // build up new activeConditionMap and restart the iterartion from begining
                    actSurf = 0;
                    activeConditionMap.zero();
                    for ( i = 1; i <= this->nsurf; i++ ) {
                        if ( ( gamma.at(i) + dgamma.at(i) ) > 0.0 ) {
                            actSurf += 1;
                            activeConditionMap.at(i) = actSurf;
                        }
                    }
 
                    continue;
                } else {
                    break;
                }
            } while ( 1 );           // end restart loop
 
            // compute plastic strain vector
            helpVector.resize(sizeR);
            helpVector.zero();
            for ( i = 1; i <= this->nsurf; i++ ) {
                if ( activeConditionMap.at(i) ) {
                    ( * loadGradVecPtr ) [ i - 1 ].times( dgamma.at(i) );
                    for ( j = 1; j <= sizeR; j++ ) {
                        helpVector.at(j) += ( * loadGradVecPtr ) [ i - 1 ].at(j);
                    }
                }
            }
 
            plasticStrainVectorR.add(helpVector);
            // update strain space hardening vector
            helpVector.resize(size - sizeR);
            helpVector.zero();
            for ( i = 1; i <= this->nsurf; i++ ) {
                if ( activeConditionMap.at(i) ) {
                    //gradVec[i-1].times(gamma.at(i));
                    for ( j = sizeR + 1; j <= size; j++ ) {
                        helpVector.at(j - sizeR) += ( * loadGradVecPtr ) [ i - 1 ].at(j);
                    }
                }
            }
 
            strainSpaceHardeningVariables.add(helpVector);
            // update plastic multipliers
            gamma.add(dgamma);
 
            elasticStrainVectorR = totalStrain;
            elasticStrainVectorR.subtract(plasticStrainVectorR);
            // stress vector in full form due to computational convinience
            this->computeTrialStressIncrement(fullStressVector, gp, elasticStrainVectorR, tStep);
            this->computeStressSpaceHardeningVars(fullStressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
 
            elastic = 1;
            for ( i = 1; i <= this->nsurf; i++ ) {
                yieldValue = this->computeYieldValueAt(gp, i, fullStressVector, fullStressSpaceHardeningVars);
                // check for end of iteration
                if ( yieldValue > YIELD_TOL ) {
                    elastic = 0;
                }
            }
 
            if ( elastic ) {
                break;
            }
 
            // increment iteration count
            nIterations++;
 
            if ( nIterations > PLASTIC_MATERIAL_MAX_ITERATIONS ) {
                char buff [ 200 ], buff1 [ 150 ];
                sprintf(buff, "GiveRealStressVector: local equlibrium not reached in %d iterations", PLASTIC_MATERIAL_MAX_ITERATIONS);
                sprintf( buff1, "Element %d, gp %d, continuing", gp->giveElement()->giveNumber(), gp->giveNumber() );
                OOFEM_WARNING(buff, buff1);
                //nIterations = 0; goto huhu;
                break;
            }
        } while ( 1 );
    }
 
    printf(" (%d iterations)", nIterations);
    answer = fullStressVector;
}
 
 
void
MPlasticMaterial :: computeGradientVector(FloatArray &answer, functType ftype, int isurf,
                                          GaussPoint *gp,
                                          const FloatArray &fullStressVector,
                                          const FloatArray &fullStressSpaceHardeningVars) const
{
    /*
     * Computes gradient vector R in reduced form.
     * R = {\der{\sigma}{f_{n+1}}, \der{q}{f_{n+1}}}^T
     * Gradient vector contains partial derivatives of yield function
     * with respect to stresses and with respect to
     * strain space hardening variables
     *
     * Note: variables with R postfix are in reduced stress-space.
     */
    FloatArray stressGradient, stressGradientR;
    FloatArray stressSpaceHardVarGradient;
    int isize, size;
 
    this->computeStressGradientVector(stressGradient, ftype, isurf, gp, fullStressVector,
                                      fullStressSpaceHardeningVars);
 
    StructuralMaterial :: giveReducedSymVectorForm( stressGradientR, stressGradient, gp->giveMaterialMode() );
 
    this->computeStressSpaceHardeningVarsReducedGradient(stressSpaceHardVarGradient, ftype, isurf, gp,
                                                         fullStressVector, fullStressSpaceHardeningVars);
 
    isize = stressGradientR.giveSize();
    if ( stressSpaceHardVarGradient.isNotEmpty() ) {
        size = isize + stressSpaceHardVarGradient.giveSize();
    } else {
        size = isize;
    }
 
    answer.resize(size);
    for ( int i = 1; i <= isize; i++ ) {
        answer.at(i) = stressGradientR.at(i);
    }
 
    for ( int i = isize + 1; i <= size; i++ ) {
        answer.at(i) = stressSpaceHardVarGradient.at(i - isize);
    }
}
 
 
void
MPlasticMaterial :: computeResidualVector(FloatArray &answer, GaussPoint *gp, const FloatArray &gamma,
                                          const IntArray &activeConditionMap, const FloatArray &plasticStrainVectorR,
                                          const FloatArray &strainSpaceHardeningVariables,
                                          std :: vector< FloatArray > &gradientVectorR) const
{
    /* Computes Residual vector for closest point projection algorithm */
 
    FloatArray oldPlasticStrainVectorR, oldStrainSpaceHardeningVariables;
    int isize, size;
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
 
    isize = plasticStrainVectorR.giveSize();
    size =  gradientVectorR [ 0 ].giveSize();
 
    answer.resize(size);
    oldPlasticStrainVectorR = status->givePlasticStrainVector();
    oldStrainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
 
    for ( int i = 1; i <= isize; i++ ) {
        answer.at(i) = oldPlasticStrainVectorR.at(i) - plasticStrainVectorR.at(i);
    }
 
    for ( int i = isize + 1; i <= size; i++ ) {
        answer.at(i) = oldStrainSpaceHardeningVariables.at(i - isize) - strainSpaceHardeningVariables.at(i - isize);
    }
 
    for ( int i = 1; i <= this->nsurf; i++ ) {
        if ( activeConditionMap.at(i) ) {
            for ( int j = 1; j <= size; j++ ) {
                answer.at(j) += gamma.at(i) * gradientVectorR [ i - 1 ].at(j);
            }
        }
    }
}
 
 
void
MPlasticMaterial :: computeTrialStressIncrement(FloatArray &answer, GaussPoint *gp,
                                                const FloatArray &elasticStrainVectorR,
                                                TimeStep *tStep) const
{
    /* Computes the full trial elastic stress vector */
 
    FloatMatrix de;
    FloatArray reducedAnswer;
 
    this->computeReducedElasticModuli(de, gp, tStep);
    //this->giveLinearElasticMaterial()->giveCharacteristicMatrix(de, TangentStiffness,                                                                                                                                                                gp, tStep);
    reducedAnswer.beProductOf(de, elasticStrainVectorR);
    StructuralMaterial :: giveFullSymVectorForm( answer, reducedAnswer, gp->giveMaterialMode() );
}
 
 
void
MPlasticMaterial :: computeAlgorithmicModuli(FloatMatrix &answer,
                                             GaussPoint *gp,
                                             const FloatMatrix &elasticModuliInverse,
                                             const FloatMatrix &hardeningModuliInverse,
                                             const FloatArray &gamma,
                                             const IntArray &activeConditionMap,
                                             const FloatArray &fullStressVector,
                                             const FloatArray &fullStressSpaceHardeningVars) const
{
    /* returns consistent moduli in reduced form.
     * Note: elasticModuli and hardeningModuli will be inverted
     *
     * stressVector in full form
     */
    FloatMatrix gradientMatrix;
    FloatMatrix helpInverse;
    int i, j, isize, size;
 
    isize = elasticModuliInverse.giveNumberOfRows();
    if ( this->hasHardening() ) {
        size = isize + hardeningModuliInverse.giveNumberOfRows();
    } else {
        size = isize;
    }
 
    // assemble consistent moduli
    helpInverse.resize(size, size);
    helpInverse.zero();
    for ( i = 1; i <= nsurf; i++ ) {
        if ( activeConditionMap.at(i) ) {
            // ask gradientMatrix
            this->computeReducedGradientMatrix(gradientMatrix,  i, gp, fullStressVector,
                                               fullStressSpaceHardeningVars);
 
            gradientMatrix.times( gamma.at(i) );
            helpInverse.add(gradientMatrix);
        }
    }
 
    for ( i = 1; i <= isize; i++ ) {
        for ( j = 1; j <= isize; j++ ) {
            helpInverse.at(i, j) += elasticModuliInverse.at(i, j);
        }
    }
 
    if ( hardeningModuliInverse.giveNumberOfRows() > 0 ) {
        for ( i = isize + 1; i <= size; i++ ) {
            for ( j = isize + 1; j <= size; j++ ) {
                helpInverse.at(i, j) += hardeningModuliInverse.at(i - isize, j - isize);
            }
        }
    }
 
    if ( hasHardening() && ( hardeningModuliInverse.giveNumberOfRows() == 0 ) ) {
        FloatMatrix help;
        help.beInverseOf(helpInverse);
        answer.resize( isize + fullStressSpaceHardeningVars.giveSize(), isize + fullStressSpaceHardeningVars.giveSize() );
        answer.zero();
        for ( i = 1; i <= size; i++ ) {
            for ( j = 1; j <= size; j++ ) {
                answer.at(i, j) = help.at(i, j);
            }
        }
    } else {
        answer.beInverseOf(helpInverse);
    }
}
 
// ----------------------------------------------------------------------------//
 
 
FloatMatrix
MPlasticMaterial :: giveConsistentStiffnessMatrix(MatResponseMode mode,
                                                  GaussPoint *gp,
                                                  TimeStep *tStep) const
{
    //
    // returns receiver material matrix for given reached state
    // described by temp variables.
    //
 
    FloatMatrix consistentModuli, elasticModuli, hardeningModuli;
    FloatMatrix elasticModuliInverse, hardeningModuliInverse;
    FloatMatrix gmatInv, gmat, gradMat, helpMtrx, helpMtrx2, nmat, sbm1;
    FloatArray stressVector, fullStressVector;
    FloatArray stressSpaceHardeningVars;
    FloatArray strainSpaceHardeningVariables;
    std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf), * yieldGradVecPtr, * loadGradVecPtr;
 
    IntArray activeConditionMap;
    FloatArray gamma;
    int size, sizeR, i, j, iindx, actSurf = 0;
 
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
 
    if ( this->plType == associatedPT ) {
        yieldGradVecPtr = loadGradVecPtr = & yieldGradVec;
    } else {
        yieldGradVecPtr = & yieldGradVec;
        loadGradVecPtr  = & loadGradVec;
    }
 
    // ask for plastic consistency parameter
    gamma = status->giveTempGamma();
    activeConditionMap = status->giveTempActiveConditionMap();
    //
    // check for elastic cases
    //
    if ( ( status->giveTempStateFlag() == MPlasticMaterialStatus :: PM_Elastic ) || ( status->giveTempStateFlag() == MPlasticMaterialStatus :: PM_Unloading ) ) {
        FloatMatrix answer;
        this->linearElasticMaterial->giveStiffnessMatrix(answer, TangentStiffness, gp, tStep);
        return answer;
    }
 
    //
    // plastic case
    //
    // determine number of active surfaces
    for ( i = 1; i <= nsurf; i++ ) {
        if ( activeConditionMap.at(i) ) {
            actSurf++;
        }
    }
 
    // compute elastic moduli and its inverse
    this->computeReducedElasticModuli(elasticModuli, gp, tStep);
    elasticModuliInverse.beInverseOf(elasticModuli);
    sizeR = elasticModuliInverse.giveNumberOfRows();
 
    // compute  hardening moduli and its inverse
    this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
    if ( hardeningModuli.giveNumberOfRows() ) {
        hardeningModuliInverse.beInverseOf(hardeningModuli);
    } else {
        hardeningModuliInverse.clear();
    }
 
    stressVector = status->giveStressVector();
    StructuralMaterial :: giveFullSymVectorForm( fullStressVector, stressVector, gp->giveMaterialMode() );
    strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
    this->computeStressSpaceHardeningVars(stressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
 
 
    //
    // compute consistent moduli
    //
    this->computeAlgorithmicModuli(consistentModuli, gp, elasticModuliInverse, hardeningModuliInverse, gamma,
                                   activeConditionMap, fullStressVector, stressSpaceHardeningVars);
 
    //computee gmatInv
    for ( i = 1; i <= nsurf; i++ ) {
        this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
    }
 
    if ( this->plType == nonassociatedPT ) {
        for ( i = 1; i <= nsurf; i++ ) {
            this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
        }
    }
 
 
    size = consistentModuli.giveNumberOfColumns();
    gradMat.resize(actSurf, size);
    for ( i = 1; i <= nsurf; i++ ) {
        if ( ( iindx = activeConditionMap.at(i) ) ) {
            for ( j = 1; j <= size; j++ ) {
                gradMat.at(iindx, j) = ( * yieldGradVecPtr ) [ i - 1 ].at(j);
            }
        }
    }
 
    if ( this->plType == associatedPT ) {
        helpMtrx.resize(actSurf, size);
        helpMtrx.beProductTOf(consistentModuli, gradMat);
        gmatInv.beProductOf(gradMat, helpMtrx);
        gmat.beInverseOf(gmatInv);
 
        nmat.resize(sizeR, size);
        sbm1.beSubMatrixOf(consistentModuli, 1, sizeR, 1, size);
        nmat.beProductTOf(sbm1, gradMat);
        helpMtrx.beProductOf(nmat, gmat);
        sbm1.beSubMatrixOf(consistentModuli, 1, size, 1, sizeR);
        nmat.beProductOf(gradMat, sbm1);
 
        helpMtrx2.beProductOf(helpMtrx, nmat);
 
        FloatMatrix answer;
        answer.beSubMatrixOf(consistentModuli, 1, sizeR, 1, sizeR);
        answer.subtract(helpMtrx2);
        return answer;
    } else {
        FloatMatrix lgradMat(actSurf, size);
        for ( i = 1; i <= nsurf; i++ ) {
            if ( ( iindx = activeConditionMap.at(i) ) ) {
                for ( j = 1; j <= size; j++ ) {
                    lgradMat.at(iindx, j) = ( * loadGradVecPtr ) [ i - 1 ].at(j);
                }
            }
        }
 
        helpMtrx.resize(actSurf, size);
        helpMtrx.beProductTOf(consistentModuli, lgradMat);
        gmatInv.beProductOf(gradMat, helpMtrx);
        gmat.beInverseOf(gmatInv);
 
        nmat.resize(sizeR, size);
        sbm1.beSubMatrixOf(consistentModuli, 1, sizeR, 1, size);
        nmat.beProductTOf(sbm1, lgradMat);
        helpMtrx.beProductOf(nmat, gmat);
        sbm1.beSubMatrixOf(consistentModuli, 1, size, 1, sizeR);
        nmat.beProductOf(gradMat, sbm1);
        helpMtrx2.beProductOf(helpMtrx, nmat);
 
        FloatMatrix answer;
        answer.beSubMatrixOf(consistentModuli, 1, sizeR, 1, sizeR);
        answer.subtract(helpMtrx2);
        return answer;
    }
}
 
 
void
MPlasticMaterial :: giveElastoPlasticStiffnessMatrix(FloatMatrix &answer,
                                                     MatResponseMode mode,
                                                     GaussPoint *gp,
                                                     TimeStep *tStep) const
{
    //
    // returns receiver material matrix for given reached state
    // described by temp variables.
    //
 
    FloatMatrix dmat, elasticModuli, hardeningModuli;
    FloatMatrix gmatInv, gmat, gradMat, helpMtrx, helpMtrx2;
    FloatArray stressVector, fullStressVector;
    FloatArray stressSpaceHardeningVars;
    FloatArray strainSpaceHardeningVariables;
    std :: vector< FloatArray > yieldGradVec(this->nsurf), loadGradVec(this->nsurf);
 
    IntArray activeConditionMap;
    int size, sizeR, i, j, iindx, actSurf = 0;
 
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
 
    // ask for plastic consistency parameter
    activeConditionMap = status->giveTempActiveConditionMap();
    //
    // check for elastic cases
    //
    if ( ( status->giveTempStateFlag() == MPlasticMaterialStatus :: PM_Elastic ) || ( status->giveTempStateFlag() == MPlasticMaterialStatus :: PM_Unloading ) ) {
        this->linearElasticMaterial->giveStiffnessMatrix(answer, TangentStiffness, gp, tStep);
        return;
    }
 
    //
    // plastic case
    //
    // determine number of active surfaces
    for ( i = 1; i <= nsurf; i++ ) {
        if ( activeConditionMap.at(i) ) {
            actSurf++;
        }
    }
 
    // compute elastic moduli and its inverse
    this->computeReducedElasticModuli(elasticModuli, gp, tStep);
    sizeR = elasticModuli.giveNumberOfRows();
 
    // compute  hardening moduli and its inverse
    this->computeHardeningReducedModuli(hardeningModuli, gp, strainSpaceHardeningVariables, tStep);
 
    stressVector = status->giveStressVector();
    StructuralMaterial :: giveFullSymVectorForm( fullStressVector, stressVector, gp->giveMaterialMode() );
    strainSpaceHardeningVariables = status->giveStrainSpaceHardeningVars();
    this->computeStressSpaceHardeningVars(stressSpaceHardeningVars, gp, strainSpaceHardeningVariables);
 
 
    //
    // compute consistent moduli
    //
    size = elasticModuli.giveNumberOfRows() + hardeningModuli.giveNumberOfRows();
    dmat.resize(size, size);
    dmat.zero();
 
    for ( i = 1; i <= sizeR; i++ ) {
        for ( j = 1; j <= sizeR; j++ ) {
            dmat.at(i, j) = elasticModuli.at(i, j);
        }
    }
 
    for ( i = sizeR + 1; i <= size; i++ ) {
        for ( j = sizeR + 1; j <= size; j++ ) {
            dmat.at(i, j) = hardeningModuli.at(i - sizeR, j - sizeR);
        }
    }
 
 
    //computee gmatInv
    for ( i = 1; i <= nsurf; i++ ) {
        this->computeGradientVector(yieldGradVec [ i - 1 ], yieldFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
    }
 
    if ( this->plType == nonassociatedPT ) {
        for ( i = 1; i <= nsurf; i++ ) {
            this->computeGradientVector(loadGradVec [ i - 1 ], loadFunction, i, gp, fullStressVector, stressSpaceHardeningVars);
        }
    }
 
    gradMat.resize(actSurf, size);
    for ( i = 1; i <= nsurf; i++ ) {
        if ( ( iindx = activeConditionMap.at(i) ) ) {
            for ( j = 1; j <= size; j++ ) {
                gradMat.at(iindx, j) = yieldGradVec [ i - 1 ].at(j);
            }
        }
    }
 
    if ( this->plType == associatedPT ) {
        helpMtrx.resize(actSurf, size);
        helpMtrx.beProductOf(gradMat, dmat);
        gmatInv.beProductTOf(helpMtrx, gradMat);
        gmat.beInverseOf(gmatInv);
 
        helpMtrx.beProductOf(gradMat, elasticModuli);
        helpMtrx2.beProductOf(gmat, helpMtrx);
        answer.beTProductOf(helpMtrx, helpMtrx2);
        answer.negated();
        answer.add(elasticModuli);
    } else {
        FloatMatrix lgradMat(actSurf, size);
        for ( i = 1; i <= nsurf; i++ ) {
            if ( ( iindx = activeConditionMap.at(i) ) ) {
                for ( j = 1; j <= size; j++ ) {
                    lgradMat.at(iindx, j) = loadGradVec [ i - 1 ].at(j);
                }
            }
        }
 
        helpMtrx.resize(actSurf, size);
        helpMtrx.beProductOf(gradMat, dmat);
        gmatInv.beProductTOf(helpMtrx, lgradMat);
        gmat.beInverseOf(gmatInv);
 
        helpMtrx.beProductOf(gradMat, elasticModuli);
        helpMtrx2.beProductOf(gmat, helpMtrx);
        helpMtrx.beProductTOf(elasticModuli, lgradMat);
        answer.beTProductOf(helpMtrx, helpMtrx2);
        answer.negated();
        answer.add(elasticModuli);
    }
}
 
 
void
MPlasticMaterial :: computeDiagModuli(FloatMatrix &answer,
                                      GaussPoint *gp, FloatMatrix &elasticModuliInverse,
                                      FloatMatrix &hardeningModuliInverse) const
{
    //
    // assembles diagonal moduli from elasticModuliInverse and hardeningModuliInverse
    //
    int size1, size2;
 
    size1 = elasticModuliInverse.giveNumberOfRows();
    if ( hardeningModuliInverse.giveNumberOfRows() ) {
        size2 = size1 + hardeningModuliInverse.giveNumberOfRows();
    } else {
        size2 = size1;
    }
 
    answer.resize(size2, size2);
    answer.zero();
 
    for ( int i = 1; i <= size1; i++ ) {
        for ( int j = 1; j <= size1; j++ ) {
            answer.at(i, j) = elasticModuliInverse.at(i, j);
        }
    }
 
    for ( int i = size1 + 1; i <= size2; i++ ) {
        for ( int j = size1 + 1; j <= size2; j++ ) {
            answer.at(i, j) = hardeningModuliInverse.at(i - size1, j - size1);
        }
    }
}
 
 
void
MPlasticMaterial :: computeReducedElasticModuli(FloatMatrix &answer,
                                                GaussPoint *gp,
                                                TimeStep *tStep) const
{  /* Returns elastic moduli in reduced stress-strain space*/
    this->linearElasticMaterial->giveStiffnessMatrix(answer, ElasticStiffness, gp, tStep);
}
 
 
// overloaded from structural material
 
FloatMatrixF<6,6>
MPlasticMaterial :: give3dMaterialStiffnessMatrix(MatResponseMode mode,
                                                  GaussPoint *gp,
                                                  TimeStep *tStep) const
//
//
//
// computes full 3d constitutive matrix for case of 3d stress-strain state.
// it returns elasto-plastic stiffness material matrix.
// if strainIncrement == NULL a loading is assumed
// for detailed description see (W.F.Chen: Plasticity in Reinforced Concrete, McGraw-Hill, 1982,
// chapter 6.)
//
// if derived material would like to implement failure behaviour
// it must redefine basic Give3dMaterialStiffnessMatrix function
// in order to take possible failure (tension cracking) into account
//
//
//
{
    // we can force 3d response, and we obtain correct 3d tangent matrix,
    // but in fact, stress integration algorithm will not work
    // because in stress integration algorithm we are unable to recognize
    // which reduction from 3d case should be performed to obtain correct result.
    // so for new stressStrain state, instead of programming 3d reduction,
    // you should enhance imposeConstraints functions for ne state, and
    // then programming simple inteface function for you stressstrain state
    // calling GiveMaterailStiffenssMatrix, which imposes constrains correctly.
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->give3dMaterialStiffnessMatrix(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
FloatMatrixF<3,3>
MPlasticMaterial :: givePlaneStressStiffMtrx(MatResponseMode mode,
                                             GaussPoint *gp,
                                             TimeStep *tStep) const
 
//
// returns receiver's 2dPlaneStressMtrx
// (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
//
// standard method from Material Class overloaded, because no inversion is needed.
// the reduction from 3d case will not work
// this implementation should be faster.
{
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->givePlaneStressStiffMtrx(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
FloatMatrixF<4,4>
MPlasticMaterial :: givePlaneStrainStiffMtrx(MatResponseMode mode,
                                             GaussPoint *gp,
                                             TimeStep *tStep) const
 
//
// return receiver's 2dPlaneStrainMtrx constructed from
// general 3dMatrialStiffnessMatrix
// (2dPlaneStrain ==> eps_z = gamma_xz = gamma_yz = 0.)
//
{
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->givePlaneStrainStiffMtrx(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        this->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
FloatMatrixF<1,1>
MPlasticMaterial :: give1dStressStiffMtrx(MatResponseMode mode,
                                          GaussPoint *gp,
                                          TimeStep *tStep) const
 
//
// returns receiver's 1dMaterialStiffnessMAtrix
// (1d case ==> sigma_y = sigma_z = tau_yz = tau_zx = tau_xy  = 0.)
{
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->give1dStressStiffMtrx(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        const_cast<MPlasticMaterial*>(this)->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
FloatMatrixF<2,2>
MPlasticMaterial :: give2dBeamLayerStiffMtrx(MatResponseMode mode,
                                             GaussPoint *gp,
                                             TimeStep *tStep) const
//
// returns receiver's 2dBeamLayerStiffMtrx.
// (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
//
// standard method from Material Class overloaded, because no inversion is needed.
// the reduction from 3d case will not work
// this implementation should be faster.
{
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->give2dBeamLayerStiffMtrx(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        const_cast<MPlasticMaterial*>(this)->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
FloatMatrixF<5,5>
MPlasticMaterial :: givePlateLayerStiffMtrx(MatResponseMode mode,
                                            GaussPoint *gp,
                                            TimeStep *tStep) const
//
// returns receiver's 2dPlateLayerMtrx
// (2dPlaneStres ==> sigma_z = tau_xz = tau_yz = 0.)
//
// standard method from Material Class overloaded, because no inversion is needed.
// the reduction from 3d case will not work
// this implementation should be faster.
{
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->givePlateLayerStiffMtrx(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        const_cast<MPlasticMaterial*>(this)->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
FloatMatrixF<3,3>
MPlasticMaterial :: giveFiberStiffMtrx(MatResponseMode mode,
                                       GaussPoint *gp,
                                       TimeStep *tStep) const
//
// returns receiver's Fiber
// (1dFiber ==> sigma_y = sigma_z = tau_yz = 0.)
//
// standard method from Material Class overloaded, because no inversion is needed.
// the reduction from 3d case will not work
// this implementation should be faster.
{
    if ( mode == ElasticStiffness ) {
        return this->linearElasticMaterial->giveFiberStiffMtrx(mode, gp, tStep);
    } else if ( rmType == mpm_ClosestPoint ) {
        return this->giveConsistentStiffnessMatrix(mode, gp, tStep);
    } else {
        FloatMatrix answer;
        const_cast<MPlasticMaterial*>(this)->giveElastoPlasticStiffnessMatrix(answer, mode, gp, tStep);
        return answer;
    }
}
 
 
int
MPlasticMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    MPlasticMaterialStatus *status = static_cast< MPlasticMaterialStatus * >( this->giveStatus(gp) );
    if ( type == IST_PlasticStrainTensor ) {
        const FloatArray &ep = status->givePlasticStrainVector();
        StructuralMaterial :: giveFullSymVectorForm( answer, ep, gp->giveMaterialMode() );
        return 1;
    } else if ( type == IST_PrincipalPlasticStrainTensor ) {
        FloatArray st;
 
        const FloatArray &s = status->givePlasticStrainVector();
        StructuralMaterial :: giveFullSymVectorForm( st, s, gp->giveMaterialMode() );
 
        this->computePrincipalValues(answer, st, principal_strain);
        return 1;
    } else {
        return StructuralMaterial :: giveIPValue(answer, gp, type, tStep);
    }
}
 
 
MPlasticMaterialStatus :: MPlasticMaterialStatus(GaussPoint *g, int statusSize) :
    StructuralMaterialStatus(g),
    strainSpaceHardeningVarsVector(statusSize), tempStrainSpaceHardeningVarsVector(statusSize)
{}
 
 
void
MPlasticMaterialStatus :: printOutputAt(FILE *file, TimeStep *tStep) const
{
    StructuralMaterialStatus :: printOutputAt(file, tStep);
    fprintf(file, "status { ");
    if ( ( state_flag == MPlasticMaterialStatus :: PM_Yielding ) || ( state_flag == MPlasticMaterialStatus :: PM_Unloading ) ) {
        if ( state_flag == MPlasticMaterialStatus :: PM_Yielding ) {
            fprintf(file, " Yielding, ");
        } else {
            fprintf(file, " Unloading, ");
        }
 
        fprintf(file, " plastic strains ");
        for ( auto &val : plasticStrainVector ) {
            fprintf( file, " %.4e", val );
        }
 
        if ( strainSpaceHardeningVarsVector.giveSize() ) {
            fprintf(file, ", strain space hardening vars ");
            for ( auto &val : strainSpaceHardeningVarsVector ) {
                fprintf( file, " %.4e", val );
            }
        }
    }
 
    fprintf(file, "}\n");
}
 
 
void MPlasticMaterialStatus :: initTempStatus()
//
// initializes temp variables according to variables form previous equlibrium state.
//
{
    StructuralMaterialStatus :: initTempStatus();
 
    if ( plasticStrainVector.giveSize() == 0 ) {
        plasticStrainVector.resize( StructuralMaterial :: giveSizeOfVoigtSymVector( gp->giveMaterialMode() ) );
        plasticStrainVector.zero();
    }
 
    tempPlasticStrainVector = plasticStrainVector;
 
    tempStrainSpaceHardeningVarsVector = strainSpaceHardeningVarsVector;
 
 
    tempGamma = gamma;
    tempActiveConditionMap = activeConditionMap;
}
 
 
void
MPlasticMaterialStatus :: updateYourself(TimeStep *tStep)
//
// updates variables (nonTemp variables describing situation at previous equilibrium state)
// after a new equilibrium state has been reached
// temporary variables are having values corresponding to newly reched equilibrium.
//
{
    StructuralMaterialStatus :: updateYourself(tStep);
 
    plasticStrainVector = tempPlasticStrainVector;
    strainSpaceHardeningVarsVector = tempStrainSpaceHardeningVarsVector;
 
    state_flag = temp_state_flag;
    gamma = tempGamma;
    activeConditionMap = tempActiveConditionMap;
}
 
 
void
MPlasticMaterialStatus :: saveContext(DataStream &stream, ContextMode mode)
{
    StructuralMaterialStatus :: saveContext(stream, mode);
 
    contextIOResultType iores;
    if ( ( iores = plasticStrainVector.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = strainSpaceHardeningVarsVector.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( !stream.write(state_flag) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( ( iores = gamma.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = activeConditionMap.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
}
 
 
void
MPlasticMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode)
{
    StructuralMaterialStatus :: restoreContext(stream, mode);
 
    contextIOResultType iores;
    if ( ( iores = plasticStrainVector.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = strainSpaceHardeningVarsVector.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( !stream.read(state_flag) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( ( iores = gamma.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = activeConditionMap.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
}
} // end namespace oofem