nlinearstatic.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 "sm/EngineeringModels/nlinearstatic.h"
#include "sm/Elements/structuralelement.h"
#include "nummet.h"
#include "timestep.h"
#include "metastep.h"
#include "error.h"
#include "verbose.h"
#include "sparsenonlinsystemnm.h"
#include "nrsolver.h"
#include "calmls.h"
#include "outputmanager.h"
#include "datastream.h"
#include "classfactory.h"
#include "timer.h"
#include "contextioerr.h"
#include "sparsemtrx.h"
#include "errorestimator.h"
#include "mathfem.h"
#include "dofmanager.h"
#include "dof.h"
#include "unknownnumberingscheme.h"
#include "function.h"
 
#ifdef __MPI_PARALLEL_MODE
 #include "problemcomm.h"
 #include "communicator.h"
 #include "loadbalancer.h"
#endif
 
namespace oofem {
REGISTER_EngngModel(NonLinearStatic);
 
NonLinearStatic :: NonLinearStatic(int i, EngngModel *_master) : LinearStatic(i, _master),
    totalDisplacement(), incrementOfDisplacement(), internalForces(), initialLoadVector(), incrementalLoadVector(),
    initialLoadVectorOfPrescribed(), incrementalLoadVectorOfPrescribed()
{
    //
    // constructor
    //
    prevStepLength = 0.;
    currentStepLength = 0.;
    loadLevel = cumulatedLoadLevel = 0.;
    mstepCumulateLoadLevelFlag = 0;
    numMetStatus = CR_UNKNOWN;
    stiffMode = nls_tangentStiffness; // default
    internalVarUpdateStamp = 0;
    initFlag = loadInitFlag = 1;
    controlMode = nls_indirectControl;
    refLoadInputMode = SparseNonLinearSystemNM :: rlm_total;
    nMethod = NULL;
    initialGuessType = IG_None;
    deltaT = 1.0;
    dtFunction = 0;
    currentIterations = 0;
}
 
 
NonLinearStatic :: ~NonLinearStatic()
{
    delete nMethod;
}
 
 
NumericalMethod *NonLinearStatic :: giveNumericalMethod(MetaStep *mStep)
{
    if ( mStep == NULL ) {
        OOFEM_ERROR("undefined meta step");
    }
 
    int _val = 0;
    IR_GIVE_OPTIONAL_FIELD( ( mStep->giveAttributesRecord() ), _val, _IFT_NonLinearStatic_controlmode );
    IR_GIVE_OPTIONAL_FIELD( ( mStep->giveAttributesRecord() ), _val, "controllmode" ); 
    NonLinearStatic_controlType mode = ( NonLinearStatic_controlType ) _val;
 
    if ( mode == nls_indirectControl ) {
        if ( nMethod ) {
            if ( dynamic_cast< CylindricalALM * >(nMethod) ) {
                return nMethod;
            } else {
                delete nMethod;
            }
        }
 
        this->nMethod = new CylindricalALM(this->giveDomain(1), this);
    } else if ( mode == nls_directControl ) {
        if ( nMethod ) {
            if ( dynamic_cast< NRSolver * >(nMethod) ) {
                return nMethod;
            } else {
                delete nMethod;
            }
        }
 
        this->nMethod = new NRSolver(this->giveDomain(1), this);
    } else {
        OOFEM_ERROR("unsupported controlMode");
    }
 
    return this->nMethod;
}
 
 
void
NonLinearStatic :: updateAttributes(MetaStep *mStep)
{
    MetaStep *mStep1 = this->giveMetaStep( mStep->giveNumber() ); //this line ensures correct input file in staggered problem
    auto &ir = mStep1->giveAttributesRecord();
 
    LinearStatic :: updateAttributes(mStep1);
 
    /*
     * if ((mstep->giveFirstStepNumber() == tStep->giveNumber()) && hasString(initString, "fixload")) {
     * double factor;
     *
     * printf ("NonLinearStatic: fixed load level");
     * if (initialLoadVector.isEmpty()) initialLoadVector.resize(loadVector.giveSize());
     * if ((controlMode == nls_directControl) || (controlMode == nls_directControl2)) factor = 1.0;
     * else factor = loadLevel;
     * loadVector.times (factor);
     * initialLoadVector.add(loadVector);
     * loadVector.zero();
     * this->loadInitFlag = 1;
     * this->loadLevel = 0.0;
     * }
     */
    int _val = nls_indirectControl;
    IR_GIVE_OPTIONAL_FIELD(ir, _val, _IFT_NonLinearStatic_controlmode);
    IR_GIVE_OPTIONAL_FIELD(ir, _val, "controllmode"); 
    this->controlMode = ( NonLinearStatic_controlType ) _val;
 
    _val = IG_None;
    IR_GIVE_OPTIONAL_FIELD(ir, _val, _IFT_EngngModel_initialGuess);
    this->initialGuessType = ( InitialGuess ) _val;
 
    deltaT = 1.0;
    IR_GIVE_OPTIONAL_FIELD(ir, deltaT, _IFT_NonLinearStatic_deltat);
    if ( deltaT < 0. ) {
        OOFEM_ERROR("deltaT < 0");
    }
 
    dtFunction = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, dtFunction, _IFT_NonLinearStatic_deltatfunction);
 
    _val = nls_tangentStiffness;
    IR_GIVE_OPTIONAL_FIELD(ir, _val, _IFT_NonLinearStatic_stiffmode);
    this->stiffMode = ( NonLinearStatic_stiffnessMode ) _val;
 
    _val = SparseNonLinearSystemNM :: rlm_total;
    IR_GIVE_OPTIONAL_FIELD(ir, _val, _IFT_NonLinearStatic_refloadmode);
    this->refLoadInputMode = ( SparseNonLinearSystemNM :: referenceLoadInputModeType ) _val;
 
    mstepCumulateLoadLevelFlag = ir.hasField(_IFT_NonLinearStatic_keepll);
 
    // called just to mark field as recognized, used later
   ir.hasField(_IFT_NonLinearStatic_donotfixload);
}
 
 
void
NonLinearStatic :: initializeFrom(InputRecord &ir)
{
    LinearStatic :: initializeFrom(ir);
 
    nonlocalStiffnessFlag = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, nonlocalStiffnessFlag, _IFT_NonLinearStatic_nonlocstiff);
 
    updateElasticStiffnessFlag = false;
    if ( ir.hasField(_IFT_NonLinearStatic_updateElasticStiffnessFlag) ) {
        updateElasticStiffnessFlag = true;
    }
 
#ifdef __MPI_PARALLEL_MODE
    if ( isParallel() ) {
        //commBuff = new CommunicatorBuff (this->giveNumberOfProcesses(), CBT_dynamic);
        commBuff = new CommunicatorBuff(this->giveNumberOfProcesses(), CBT_static);
        communicator = new NodeCommunicator( this, commBuff, this->giveRank(),
                                             this->giveNumberOfProcesses() );
 
        if ( ir.hasField(_IFT_NonLinearStatic_nonlocalext) ) {
            nonlocalExt = 1;
            nonlocCommunicator = new ElementCommunicator( this, commBuff, this->giveRank(),
                                                          this->giveNumberOfProcesses() );
        }
    }
#endif
}
 
 
double NonLinearStatic :: giveUnknownComponent(ValueModeType mode, TimeStep *tStep, Domain *d, Dof *dof)
// returns unknown quantity like displacement, velocity of equation eq
// This function translates this request to numerical method language
{
    if ( tStep != this->giveCurrentStep() ) {
        OOFEM_ERROR("unknown time step encountered");
    }
 
    int eq = dof->__giveEquationNumber();
    if (mode == VM_Residual) {
        // evaluate the residual of momentum balance for specific unknown
        // evaluate the residual of momentum balance for specific unknown
        
        if (eq && internalForces.isNotEmpty()) {
            double ans = loadLevel*incrementalLoadVector.at(eq)-internalForces.at(eq);
            if (initialLoadVector.isNotEmpty()) {
                ans += initialLoadVector.at(eq);
            }
            return ans;
        } else {
            return 0.;
        }
    } else if (mode == VM_Total) {
#ifdef DEBUG
        if ( eq == 0 ) OOFEM_ERROR("invalid equation number");
#endif
        if ( totalDisplacement.isNotEmpty() ) {
            return totalDisplacement.at(eq);
        } else {
            return 0.;
        }
    } else if (mode == VM_Incremental) {
#ifdef DEBUG
        if ( eq == 0 ) OOFEM_ERROR("invalid equation number");
#endif
        // return incrementOfDisplacement -> at(eq);
        // return nMethod-> giveUnknownComponent(IncrementOfSolution, eq);
        if ( incrementOfDisplacement.isNotEmpty() ) {
            return incrementOfDisplacement.at(eq);
        } else {
            return 0.;
        }
    } else if (mode == VM_Velocity) {
#ifdef DEBUG
        if ( eq == 0 ) OOFEM_ERROR("invalid equation number");
#endif
        if ( incrementOfDisplacement.isNotEmpty() ) {
            return incrementOfDisplacement.at(eq) / tStep->giveTimeIncrement();
        } else {
            return 0.;
        }
    } else {
       OOFEM_ERROR("Unknown is of undefined ValueModeType for this problem"); 
    }
}
 
TimeStep *NonLinearStatic :: giveSolutionStepWhenIcApply(bool force)
{
    if ( master && ( !force ) ) {
        return master->giveSolutionStepWhenIcApply();
    } else {
        if ( !stepWhenIcApply ) {
            int inin = giveNumberOfTimeStepWhenIcApply();
            //        int nFirst = giveNumberOfFirstStep();
            stepWhenIcApply = std :: make_unique< TimeStep >(inin, this, 0, -deltaT, deltaT, 0);
        }
 
        return stepWhenIcApply.get();
    }
}
 
 
TimeStep *NonLinearStatic :: giveNextStep()
{
    int istep = giveNumberOfFirstStep();
    int mStepNum = 1;
    double totalTime = 0.0;
    StateCounterType counter = 1;
    // double deltaTtmp = deltaT;
    double deltaTtmp = this->giveDeltaT(istep);
 
    if ( currentStep ) {
        istep = currentStep->giveNumber() + 1;
    deltaTtmp = this->giveDeltaT(istep);
    totalTime = currentStep->giveTargetTime() + deltaTtmp;
        counter = currentStep->giveSolutionStateCounter() + 1;
        mStepNum = currentStep->giveMetaStepNumber();
    /*
        if ( !this->giveMetaStep(mStepNum)->isStepValid(istep) ) {
            mStepNum++;
            if ( mStepNum > nMetaSteps ) {
                OOFEM_ERROR("no next step available, mStepNum=%d > nMetaSteps=%d", mStepNum, nMetaSteps);
            }
        }
    */
    } else {
        // first step -> generate initial step
        TimeStep *newStep = giveSolutionStepWhenIcApply();
        currentStep = std :: make_unique< TimeStep >(* newStep);
    }
 
    previousStep = std :: move(currentStep);
    currentStep = std :: make_unique< TimeStep >(istep, this, mStepNum, totalTime, deltaTtmp, counter);
    // dt variable are set eq to 0 for statics - has no meaning
    // *Wrong* It has meaning for viscoelastic materials.
 
    return currentStep.get();
}
 
Function *
NonLinearStatic :: giveDtFunction()
// Returns the load-time function of the receiver.
{
    if ( !dtFunction ) {
        return NULL;
    }
 
    return giveDomain(1)->giveFunction(dtFunction);
}
 
 
double
NonLinearStatic :: giveDeltaT(int n)
{
 
    //do not increase deltaT on microproblem
    if ( pScale == microScale ) {
      return 0.;
    }
  
    if ( giveDtFunction() ) {
      return giveDtFunction()->evaluateAtTime(n);
    }
 
    return deltaT;
}
 
 
void NonLinearStatic :: solveYourself()
{
    if ( this->isParallel() ) {
#ifdef __VERBOSE_PARALLEL
        // force equation numbering before setting up comm maps
        int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );
        OOFEM_LOG_INFO("[process rank %d] neq is %d\n", this->giveRank(), neq);
#endif
 
        // set up communication patterns
        this->initializeCommMaps();
        // init remote dofman list
        // this->initRemoteDofManList ();
    }
    StructuralEngngModel :: solveYourself();
}
 
 
void
NonLinearStatic :: solveYourselfAt(TimeStep *tStep)
{
    proceedStep(1, tStep);
}
 
 
void
NonLinearStatic :: terminate(TimeStep *tStep)
{
    this->doStepOutput(tStep);
    this->updateLoadVectors(tStep);
    this->saveStepContext(tStep, CM_State | CM_Definition);
}
 
 
void
NonLinearStatic :: updateLoadVectors(TimeStep *tStep)
{
    MetaStep *mstep = this->giveMetaStep( tStep->giveMetaStepNumber() );
    double msFinalTime = mstep->giveFinalTime() - this->giveInitialTime();
    bool isLastMetaStep = (fabs(tStep->giveTargetTime() - msFinalTime) <= 1.e-6 * mstep->giveMinDeltaT());
    if ( controlMode == nls_indirectControl ) {
        //if ((tStep->giveNumber() == mstep->giveLastStepNumber()) &&ir.hasField("fixload")) {
        if ( isLastMetaStep ) {
            if ( !mstep->giveAttributesRecord().hasField(_IFT_NonLinearStatic_donotfixload) ) {
                OOFEM_LOG_INFO("Fixed load level\n");
 
                //update initialLoadVector
                initialLoadVector.add(loadLevel, incrementalLoadVector);
 
                initialLoadVectorOfPrescribed.add(loadLevel, incrementalLoadVectorOfPrescribed);
 
                incrementalLoadVector.zero();
                incrementalLoadVectorOfPrescribed.zero();
 
                this->loadInitFlag = 1;
            }
 
            //if (!mstep->giveAttributesRecord()->hasField("keepll")) this->loadLevelInitFlag = 1;
        }
    } else { // direct control
        //update initialLoadVector after each step of direct control
        //(here the loading is not proportional)
        OOFEM_LOG_DEBUG("Fixed load level\n");
 
        initialLoadVector.add(loadLevel, incrementalLoadVector);
 
        initialLoadVectorOfPrescribed.add(loadLevel, incrementalLoadVectorOfPrescribed);
 
        incrementalLoadVector.zero();
        incrementalLoadVectorOfPrescribed.zero();
 
        this->loadInitFlag = 1;
    }
 
 
    // if (isLastMetaStep) {
    if ( isLastMetaStep && !mstep->giveAttributesRecord().hasField(_IFT_NonLinearStatic_donotfixload) ) {
#ifdef VERBOSE
        OOFEM_LOG_INFO("Reseting load level\n");
#endif
        if ( mstepCumulateLoadLevelFlag ) {
            cumulatedLoadLevel += loadLevel;
        } else {
            cumulatedLoadLevel = 0.0;
        }
 
        this->loadLevel = 0.0;
    }
}
 
 
void
NonLinearStatic :: proceedStep(int di, TimeStep *tStep)
{
    //
    // creates system of governing eq's and solves them at given time step
    //
    // first assemble problem at current time step
    //
 
    if ( initFlag ) {
        //
        // first step  create space for stiffness Matrix
        //
        if ( !stiffnessMatrix ) {
            stiffnessMatrix = classFactory.createSparseMtrx(sparseMtrxType);
        }
 
        if ( !stiffnessMatrix ) {
            OOFEM_ERROR("sparse matrix creation failed");
        }
 
        if ( nonlocalStiffnessFlag ) {
            if ( !stiffnessMatrix->isAsymmetric() ) {
                OOFEM_ERROR("stiffnessMatrix does not support asymmetric storage");
            }
        }
 
        stiffnessMatrix->buildInternalStructure( this, di, EModelDefaultEquationNumbering() );
    }
 
#if 0
    if ( ( mstep->giveFirstStepNumber() == tStep->giveNumber() ) ) {
 #ifdef VERBOSE
        OOFEM_LOG_INFO("Resetting load level\n");
 #endif
        if ( mstepCumulateLoadLevelFlag ) {
            cumulatedLoadLevel += loadLevel;
        } else {
            cumulatedLoadLevel = 0.0;
        }
        this->loadLevel = 0.0;
    }
#endif
 
    if ( loadInitFlag || controlMode == nls_directControl ) {
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling reference load\n");
#endif
        //
        // assemble the incremental reference load vector
        //
        this->assembleIncrementalReferenceLoadVectors(incrementalLoadVector, incrementalLoadVectorOfPrescribed,
                                                      refLoadInputMode, this->giveDomain(di), tStep);
 
        loadInitFlag = 0;
    }
 
    if ( tStep->giveNumber() == 1 ) {
        int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );
        totalDisplacement.resize(neq);
        totalDisplacement.zero();
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
    }
 
    //
    //    ->   BEGINNING OF LOAD (OR DISPLACEMENT) STEP  <-
    //
 
    //
    // set-up numerical model
    //
    this->giveNumericalMethod( this->giveMetaStep( tStep->giveMetaStepNumber() ) );
    //
    // call numerical model to solve arise problem
    //
#ifdef VERBOSE
    OOFEM_LOG_RELEVANT( "\n\nSolving       [step number %5d.%d, time = %e]\n\n", tStep->giveNumber(), tStep->giveVersion(), tStep->giveIntrinsicTime() );
#endif
 
    if ( this->initialGuessType == IG_Tangent ) {
#ifdef VERBOSE
        OOFEM_LOG_RELEVANT("Computing initial guess\n");
#endif
        FloatArray extrapolatedForces;
        this->assemblePrescribedExtrapolatedForces( extrapolatedForces, tStep, TangentStiffnessMatrix, this->giveDomain(di) );
        extrapolatedForces.negated();
        this->updateMatrix( * stiffnessMatrix, tStep, this->giveDomain(di) );
        SparseLinearSystemNM *linSolver = nMethod->giveLinearSolver();
        OOFEM_LOG_RELEVANT("solving for increment\n");
        linSolver->solve(* stiffnessMatrix, extrapolatedForces, incrementOfDisplacement);
        OOFEM_LOG_RELEVANT("initial guess found\n");
        totalDisplacement.add(incrementOfDisplacement);
    } else if ( this->initialGuessType != IG_None ) {
        OOFEM_ERROR("Initial guess type: %d not supported", initialGuessType);
    } else {
        incrementOfDisplacement.zero();
    }
 
    //totalDisplacement.printYourself();
    if ( initialLoadVector.isNotEmpty() ) {
        numMetStatus = nMethod->solve(* stiffnessMatrix, incrementalLoadVector, & initialLoadVector,
                                      totalDisplacement, incrementOfDisplacement, internalForces,
                                      internalForcesEBENorm, loadLevel, refLoadInputMode, currentIterations, tStep);
    } else {
        numMetStatus = nMethod->solve(* stiffnessMatrix, incrementalLoadVector, NULL,
                                      totalDisplacement, incrementOfDisplacement, internalForces,
                                      internalForcesEBENorm, loadLevel, refLoadInputMode, currentIterations, tStep);
    }
    OOFEM_LOG_RELEVANT("Equilibrium reached at load level = %f in %d iterations\n", cumulatedLoadLevel + loadLevel, currentIterations);
    prevStepLength =  currentStepLength;
    tStep->numberOfIterations = currentIterations;
    tStep->convergedReason = numMetStatus;
}
 
 
void
NonLinearStatic :: updateSolution(FloatArray &solutionVector, TimeStep *tStep, Domain *d)
{
    // No-op: This can't really be supported in any nice way in nlinearstatic.
}
 
 
void
NonLinearStatic :: updateMatrix(SparseMtrx &mat, TimeStep *tStep, Domain *d)
{
    if ( stiffMode == nls_tangentStiffness ) {
        mat.zero(); // zero stiffness matrix
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling tangent stiffness matrix\n");
#endif
        this->assemble(mat, tStep, TangentAssembler(TangentStiffness), EModelDefaultEquationNumbering(), d);
    } else if ( ( stiffMode == nls_secantStiffness ) || ( stiffMode == nls_secantInitialStiffness && initFlag ) ) {
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling secant stiffness matrix\n");
#endif
        mat.zero(); // zero stiffness matrix
        this->assemble(mat, tStep, TangentAssembler(SecantStiffness), EModelDefaultEquationNumbering(), d);
        initFlag = 0;
    } else if ( ( stiffMode == nls_elasticStiffness ) && ( initFlag ||
                                                           ( this->giveMetaStep( tStep->giveMetaStepNumber() )->giveFirstStepNumber() == tStep->giveNumber() ) || ( updateElasticStiffnessFlag ) ) ) {
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling elastic stiffness matrix\n");
#endif
        mat.zero(); // zero stiffness matrix
        this->assemble(mat, tStep, TangentAssembler(ElasticStiffness),
                       EModelDefaultEquationNumbering(), d);
        initFlag = 0;
    } else {
        // currently no action , this method is mainly intended to
        // assemble new tangent stiffness after each iteration
        // when secantStiffMode is on, we use the same stiffness
        // during iteration process
    }
}
 
 
void
NonLinearStatic :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
//
// updates some component, which is used by numerical method
// to newly reached state. used mainly by numerical method
// when new tangent stiffness is needed during finding
// of new equilibrium stage.
//
{
    switch ( cmpn ) {
    case NonLinearLhs:
        if ( stiffMode == nls_tangentStiffness ) {
            stiffnessMatrix->zero(); // zero stiffness matrix
#ifdef VERBOSE
            OOFEM_LOG_DEBUG("Assembling tangent stiffness matrix\n");
#endif
            this->assemble(* stiffnessMatrix, tStep, TangentAssembler(TangentStiffness),
                           EModelDefaultEquationNumbering(), d);
        } else if ( ( stiffMode == nls_secantStiffness ) || ( stiffMode == nls_secantInitialStiffness && initFlag ) ) {
#ifdef VERBOSE
            OOFEM_LOG_DEBUG("Assembling secant stiffness matrix\n");
#endif
            stiffnessMatrix->zero(); // zero stiffness matrix
            this->assemble(* stiffnessMatrix, tStep, TangentAssembler(SecantStiffness),
                           EModelDefaultEquationNumbering(), d);
            initFlag = 0;
        } else if ( ( stiffMode == nls_elasticStiffness ) && ( initFlag ||
                                                               ( this->giveMetaStep( tStep->giveMetaStepNumber() )->giveFirstStepNumber() == tStep->giveNumber() ) || ( updateElasticStiffnessFlag ) ) ) {
#ifdef VERBOSE
            OOFEM_LOG_DEBUG("Assembling elastic stiffness matrix\n");
#endif
            stiffnessMatrix->zero(); // zero stiffness matrix
            this->assemble(* stiffnessMatrix, tStep, TangentAssembler(ElasticStiffness),
                           EModelDefaultEquationNumbering(), d);
            initFlag = 0;
        } else {
            // currently no action , this method is mainly intended to
            // assemble new tangent stiffness after each iteration
            // when secantStiffMode is on, we use the same stiffness
            // during iteration process
        }
 
        break;
    case InternalRhs:
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Updating internal forces\n");
#endif
        // update internalForces and internalForcesEBENorm concurrently
        this->updateInternalRHS(internalForces, tStep, d, & this->internalForcesEBENorm);
        break;
 
    case ExternalRhs:
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Updating external forces\n");
#endif
        this->assembleIncrementalReferenceLoadVectors(incrementalLoadVector, incrementalLoadVectorOfPrescribed, this->refLoadInputMode, d, tStep);
        break;
 
    default:
        OOFEM_ERROR("Unknown Type of component.");
    }
}
 
 
void
NonLinearStatic :: printOutputAt(FILE *file, TimeStep *tStep)
{
    if ( !this->giveDomain(1)->giveOutputManager()->testTimeStepOutput(tStep) ) {
        return;                                                                      // do not print even Solution step header
    }
 
    fprintf( file, "\n\nOutput for time %.3e, solution step number %d\n", tStep->giveTargetTime(), tStep->giveNumber() );
    fprintf(file, "Reached load level : %20.6f in %d iterations\n\n",
            cumulatedLoadLevel + loadLevel, currentIterations);
 
    nMethod->printState(file);
 
    this->giveDomain(1)->giveOutputManager()->doDofManOutput(file, tStep);
    this->giveDomain(1)->giveOutputManager()->doElementOutput(file, tStep);
    this->printReactionForces(tStep, 1, file);
}
 
 
void
NonLinearStatic :: saveContext(DataStream &stream, ContextMode mode)
{
    contextIOResultType iores;
 
    EngngModel :: saveContext(stream, mode);
 
    if ( ( iores = totalDisplacement.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = incrementOfDisplacement.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    int _cm = controlMode;
    if ( !stream.write(_cm) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.write(loadLevel) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.write(cumulatedLoadLevel) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( ( iores = initialLoadVector.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = initialLoadVectorOfPrescribed.storeYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
}
 
 
void
NonLinearStatic :: restoreContext(DataStream &stream, ContextMode mode)
{
    contextIOResultType iores;
 
    EngngModel :: restoreContext(stream, mode);
 
    //if ((iores = this->giveNumericalMethod(giveCurrentStep())->restoreContext (stream)) !=CIO_OK) THROW_CIOERR(iores);
 
    if ( ( iores = totalDisplacement.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = incrementOfDisplacement.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    int _cm;
    if ( !stream.read(_cm) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    controlMode = ( NonLinearStatic_controlType ) _cm;
    if ( !stream.read(loadLevel) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
    if ( !stream.read(cumulatedLoadLevel) ) {
        THROW_CIOERR(CIO_IOERR);
    }
 
 
    // store InitialLoadVector
    if ( ( iores = initialLoadVector.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
 
    if ( ( iores = initialLoadVectorOfPrescribed.restoreYourself(stream) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }
}
 
 
void
NonLinearStatic :: updateDomainLinks()
{
    LinearStatic :: updateDomainLinks();
 
    this->giveNumericalMethod( this->giveCurrentMetaStep() )->setDomain( this->giveDomain(1) );
#ifdef __MPI_PARALLEL_MODE
    if ( this->giveLoadBalancer() ) {
        this->giveLoadBalancer()->setDomain( this->giveDomain(1) );
    }
 
#endif
}
 
 
void
NonLinearStatic :: assemble(SparseMtrx &answer, TimeStep *tStep, const MatrixAssembler &ma,
                            const UnknownNumberingScheme &s, Domain *domain)
{
#ifdef TIME_REPORT
    Timer _timer;
    _timer.startTimer();
#endif
 
    LinearStatic :: assemble(answer, tStep, ma, s, domain);
 
    if ( ( nonlocalStiffnessFlag ) && dynamic_cast< const TangentAssembler * >(& ma) ) {
        // add nonlocal contribution
        for ( auto &elem : domain->giveElements() ) {
            static_cast< StructuralElement * >( elem.get() )->addNonlocalStiffnessContributions(answer, s, tStep);
        }
 
        // print storage statistics
        answer.printStatistics();
    }
 
#ifdef TIME_REPORT
    _timer.stopTimer();
    OOFEM_LOG_DEBUG( "NonLinearStatic: User time consumed by assembly: %.2fs\n", _timer.getUtime() );
#endif
}
 
 
#ifdef __OOFEG
void
NonLinearStatic :: showSparseMtrxStructure(int type, oofegGraphicContext &gc, TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    CharType ctype;
 
    if ( type != 1 ) {
        return;
    }
 
    if ( stiffMode == nls_tangentStiffness ) {
        ctype = TangentStiffnessMatrix;
    } else if ( stiffMode == nls_secantStiffness ) {
        ctype = SecantStiffnessMatrix;
    } else {
        ctype = SecantStiffnessMatrix;
    }
 
    for ( auto &elem : domain->giveElements() ) {
        elem->showSparseMtrxStructure(ctype, gc, tStep);
    }
 
    for ( auto &elem : domain->giveElements() ) {
        elem->showExtendedSparseMtrxStructure(ctype, gc, tStep);
    }
}
#endif
 
 
void
NonLinearStatic :: computeExternalLoadReactionContribution(FloatArray &reactions, TimeStep *tStep, int di)
{
    if ( ( di == 1 ) && ( tStep == this->giveCurrentStep() ) ) {
        reactions = initialLoadVectorOfPrescribed;
        reactions.add(loadLevel, incrementalLoadVectorOfPrescribed);
    } else {
        OOFEM_ERROR("unable to respond due to invalid solution step or domain");
    }
}
 
 
void
NonLinearStatic :: assembleIncrementalReferenceLoadVectors(FloatArray &_incrementalLoadVector,
                                                           FloatArray &_incrementalLoadVectorOfPrescribed,
                                                           SparseNonLinearSystemNM :: referenceLoadInputModeType _refMode,
                                                           Domain *sourceDomain, TimeStep *tStep)
{
    _incrementalLoadVector.resize( sourceDomain->giveEngngModel()->giveNumberOfDomainEquations( sourceDomain->giveNumber(), EModelDefaultEquationNumbering() ) );
    _incrementalLoadVector.zero();
    _incrementalLoadVectorOfPrescribed.resize( sourceDomain->giveEngngModel()->giveNumberOfDomainEquations( sourceDomain->giveNumber(), EModelDefaultPrescribedEquationNumbering() ) );
    _incrementalLoadVectorOfPrescribed.zero();
 
    if ( _refMode == SparseNonLinearSystemNM :: rlm_incremental ) {
        this->assembleVector(_incrementalLoadVector, tStep, ExternalForceAssembler(),
                             VM_Incremental, EModelDefaultEquationNumbering(), sourceDomain);
 
        this->assembleVector(_incrementalLoadVectorOfPrescribed, tStep, ExternalForceAssembler(),
                             VM_Incremental, EModelDefaultPrescribedEquationNumbering(), sourceDomain);
    } else {
        this->assembleVector(_incrementalLoadVector, tStep, ExternalForceAssembler(),
                             VM_Total, EModelDefaultEquationNumbering(), sourceDomain);
 
        this->assembleVector(_incrementalLoadVectorOfPrescribed, tStep, ExternalForceAssembler(),
                             VM_Total, EModelDefaultPrescribedEquationNumbering(), sourceDomain);
    }
 
    this->updateSharedDofManagers(_incrementalLoadVector, EModelDefaultEquationNumbering(), LoadExchangeTag);
}
 
 
int
NonLinearStatic :: estimateMaxPackSize(IntArray &commMap, DataStream &buff, int packUnpackType)
{
    int count = 0, pcount = 0;
    Domain *domain = this->giveDomain(1);
 
    if ( packUnpackType == 0 ) { 
        for ( int map : commMap ) {
            DofManager *dman = domain->giveDofManager(map);
            for ( Dof *dof : *dman ) {
                if ( dof->isPrimaryDof() && ( dof->__giveEquationNumber() ) ) {
                    count++;
                } else {
                    pcount++;
                }
            }
        }
 
        //printf ("\nestimated count is %d\n",count);
        return ( buff.givePackSizeOfDouble(1) * max(count, pcount) );
    } else if ( packUnpackType == 1 ) {
        for ( int map : commMap ) {
            count += domain->giveElement(map)->estimatePackSize(buff);
        }
 
        return count;
    }
 
    return 0;
}
 
 
#ifdef __MPI_PARALLEL_MODE
LoadBalancer *
NonLinearStatic :: giveLoadBalancer()
{
    if ( lb ) {
        return lb.get();
    }
 
    if ( loadBalancingFlag ) {
        lb = classFactory.createLoadBalancer( "parmetis", this->giveDomain(1) );
        return lb.get();
    } else {
        return nullptr;
    }
}
 
 
LoadBalancerMonitor *
NonLinearStatic :: giveLoadBalancerMonitor()
{
    if ( lbm ) {
        return lbm.get();
    }
 
    if ( loadBalancingFlag ) {
        lbm = classFactory.createLoadBalancerMonitor("wallclock", this);
        return lbm.get();
    } else {
        return nullptr;
    }
}
#endif
 
void
NonLinearStatic :: packMigratingData(TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    int ndofman = domain->giveNumberOfDofManagers();
    bool initialLoadVectorEmpty = initialLoadVector.isEmpty();
    bool initialLoadVectorOfPrescribedEmpty = initialLoadVectorOfPrescribed.isEmpty();
 
    for ( int idofman = 1; idofman <= ndofman; idofman++ ) {
        DofManager *_dm = domain->giveDofManager(idofman);
        for ( Dof *_dof : *_dm ) {
            if ( _dof->isPrimaryDof() ) {
                int _eq;
                if ( ( _eq = _dof->__giveEquationNumber() ) ) {
                    // pack values in solution vectors
                    _dof->updateUnknownsDictionary( tStep, VM_Total, totalDisplacement.at(_eq) );
                    if ( initialLoadVectorEmpty ) {
                        _dof->updateUnknownsDictionary(tStep, VM_RhsInitial, 0.0);
                    } else {
                        _dof->updateUnknownsDictionary( tStep, VM_RhsInitial, initialLoadVector.at(_eq) );
                    }
 
                    _dof->updateUnknownsDictionary( tStep, VM_RhsIncremental, incrementalLoadVector.at(_eq) );
                } else if ( ( _eq = _dof->__givePrescribedEquationNumber() ) ) {
                    // pack values in prescribed solution vectors
                    if ( initialLoadVectorOfPrescribedEmpty ) {
                        _dof->updateUnknownsDictionary(tStep, VM_RhsInitial, 0.0);
                    } else {
                        _dof->updateUnknownsDictionary( tStep, VM_RhsInitial, initialLoadVectorOfPrescribed.at(_eq) );
                    }
 
                    _dof->updateUnknownsDictionary( tStep, VM_RhsIncremental, incrementalLoadVectorOfPrescribed.at(_eq) );
                }
            } // end primary dof
        } // end dof loop
    } // end dofman loop
}
 
 
void
NonLinearStatic :: unpackMigratingData(TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    int ndofman = domain->giveNumberOfDofManagers();
    //int myrank = this->giveRank();
 
    // resize target arrays
    int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );
    totalDisplacement.resize(neq);
    incrementOfDisplacement.resize(neq);
    incrementalLoadVector.resize(neq);
    initialLoadVector.resize(neq);
    initialLoadVectorOfPrescribed.resize( giveNumberOfDomainEquations( 1, EModelDefaultPrescribedEquationNumbering() ) );
    incrementalLoadVectorOfPrescribed.resize( giveNumberOfDomainEquations( 1, EModelDefaultPrescribedEquationNumbering() ) );
 
    for ( int idofman = 1; idofman <= ndofman; idofman++ ) {
        DofManager *_dm = domain->giveDofManager(idofman);
        for ( Dof *_dof : *_dm ) {
            if ( _dof->isPrimaryDof() ) {
                int _eq;
                if ( ( _eq = _dof->__giveEquationNumber() ) ) {
                    // pack values in solution vectors
                    totalDisplacement.at(_eq) = _dof->giveUnknownsDictionaryValue(tStep, VM_Total);
                    initialLoadVector.at(_eq) = _dof->giveUnknownsDictionaryValue(tStep, VM_RhsInitial);
                    incrementalLoadVector.at(_eq) = _dof->giveUnknownsDictionaryValue(tStep, VM_RhsIncremental);
 
#if 0
                    // debug print
                    if ( _dm->giveParallelMode() == DofManager_shared ) {
                        fprintf(stderr, "[%d] Shared: %d(%d) -> %d\n", myrank, idofman, idof, _eq);
                    } else {
                        fprintf(stderr, "[%d] Local : %d(%d) -> %d\n", myrank, idofman, idof, _eq);
                    }
 
#endif
                } else if ( ( _eq = _dof->__givePrescribedEquationNumber() ) ) {
                    // pack values in prescribed solution vectors
                    initialLoadVectorOfPrescribed.at(_eq) = _dof->giveUnknownsDictionaryValue(tStep, VM_RhsInitial);
                    incrementalLoadVectorOfPrescribed.at(_eq) = _dof->giveUnknownsDictionaryValue(tStep, VM_RhsIncremental);
 
#if 0
                    // debug print
                    fprintf(stderr, "[%d] %d(%d) -> %d\n", myrank, idofman, idof, -_eq);
#endif
                }
            } // end primary dof
        } // end dof loop
    } // end dofman loop
 
    this->initializeCommMaps(true);
    nMethod->reinitialize();
    // reinitialize error estimator (if any)
    if ( this->giveDomainErrorEstimator(1) ) {
        this->giveDomainErrorEstimator(1)->reinitialize();
    }
 
    initFlag = true;
}
} // end namespace oofem