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 ) {
        ///@todo This was almost definitely wrong before. It never seems to be used. Is this code even relevant?
        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);
}

void
SUPG :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    // update element stabilization
    for ( auto &elem : d->giveElements() ) {
        static_cast< FMElement * >( elem.get() )->updateStabilizationCoeffs(tStep);
    }

    if ( cmpn == InternalRhs ) {
        this->internalForces.zero();
        this->assembleVector(this->internalForces, tStep, SUPGInternalForceAssembler(lscale, dscale, uscale), VM_Total,
                             EModelDefaultEquationNumbering(), d, & this->eNorm);
        this->updateSharedDofManagers(this->internalForces, EModelDefaultEquationNumbering(), InternalForcesExchangeTag);
        return;
    } else if ( cmpn == NonLinearLhs ) {
        this->lhs->zero();
        //if ( 1 ) { //if ((nite > 5)) // && (rnorm < 1.e4))
        this->assemble( *lhs, tStep, SUPGTangentAssembler(TangentStiffness, lscale, dscale, uscale, alpha),
                        EModelDefaultEquationNumbering(), d );
       // } else {
       //     this->assemble( lhs, tStep, SUPGTangentAssembler(SecantStiffness),
       //                     EModelDefaultEquationNumbering(), d );
       // }
        return;
    } else {
        OOFEM_ERROR("Unknown component");
    }
}

void DarcyFlow :: DumpMatricesToFile(FloatMatrix *LHS, FloatArray *RHS, FloatArray *SolutionVector)
{
    FILE *rhsFile = fopen("RHS.txt", "w");
    // rhs.printYourself();

    for ( int i = 1; i <= RHS->giveSize(); i++ ) {
        fprintf( rhsFile, "%0.15e\n", RHS->at(i) );
    }
    fclose(rhsFile);

    FILE *lhsFile = fopen("LHS.txt", "w");

    for ( int i = 1; i <= this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ); i++ ) {
        for ( int j = 1; j <= this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ); j++ ) {
            fprintf( lhsFile, "%0.15e\t", LHS->at(i, j) );
        }
        fprintf(lhsFile, "\n");
    }

    fclose(lhsFile);

    if ( SolutionVector == NULL ) {
        return;
    }

    FILE *SolutionFile = fopen("SolutionVector.txt", "w");
    for ( int i = 1; i <= SolutionVector->giveSize(); i++ ) {
        fprintf( SolutionFile, "%0.15e\n", SolutionVector->at(i) );
    }
    fclose(SolutionFile);
}

void
TransientTransportProblem :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    // F(T) + C*dT/dt = Q
    // Linearized:
    // F(T^(k)) + K*a*dT_1 = Q - C * dT/dt^(k) - C/dt * dT_1
    // Rearranged
    // (a*K + C/dt) * dT_1 = Q - (F(T^(k)) + C * dT/dt^(k))
    // K_eff        * dT_1 = Q - F_eff
    // Update:
    // T_1 += dT_1
    
    ///@todo NRSolver should report when the solution changes instead of doing it this way.
    this->field->update(VM_Total, tStep, solution, EModelDefaultEquationNumbering());
    ///@todo Need to reset the boundary conditions properly since some "update" is doing strange 
    /// things such as applying the (wrong) boundary conditions. This call will be removed when that code can be removed.
    this->field->applyBoundaryCondition(tStep);

    if ( cmpn == InternalRhs ) {
        // F_eff = F(T^(k)) + C * dT/dt^(k)
        this->internalForces.zero();
        this->assembleVector(this->internalForces, tStep, InternalForceAssembler(), VM_Total,
                             EModelDefaultEquationNumbering(), this->giveDomain(1), & this->eNorm);
        this->updateSharedDofManagers(this->internalForces, EModelDefaultEquationNumbering(), InternalForcesExchangeTag);

        if ( lumped ) {
            // Note, inertia contribution cannot be computed on element level when lumped mass matrices are used.
            FloatArray oldSolution, vel;
            this->field->initialize(VM_Total, tStep->givePreviousStep(), oldSolution, EModelDefaultEquationNumbering());
            vel.beDifferenceOf(solution, oldSolution);
            vel.times( 1./tStep->giveTimeIncrement() );
            for ( int i = 0; i < vel.giveSize(); ++i ) {
                this->internalForces[i] += this->capacityDiag[i] * vel[i];
            }
        } else {
            FloatArray tmp;
            this->assembleVector(this->internalForces, tStep, InertiaForceAssembler(), VM_Total,
                                EModelDefaultEquationNumbering(), this->giveDomain(1), & tmp);
            this->eNorm.add(tmp); ///@todo Fix this, assembleVector shouldn't zero eNorm inside the functions. / Mikael
        }

    } else if ( cmpn == NonLinearLhs ) {
        // K_eff = (a*K + C/dt)
        if ( !this->keepTangent ) {
            this->effectiveMatrix->zero();
            this->assemble( *effectiveMatrix, tStep, TangentAssembler(TangentStiffness),
                           EModelDefaultEquationNumbering(), this->giveDomain(1) );
            effectiveMatrix->times(alpha);
            if ( lumped ) {
                effectiveMatrix->addDiagonal(1./tStep->giveTimeIncrement(), capacityDiag);
            } else {
                effectiveMatrix->add(1./tStep->giveTimeIncrement(), *capacityMatrix);
            }
        }
    } else {
        OOFEM_ERROR("Unknown component");
    }
}

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() ) ) ) {
#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->giveInternalForces(internalForces, true, d->giveNumber(), tStep);
        break;

    default:
        OOFEM_ERROR("Unknown Type of component.");
    }
}

void EigenValueDynamic :: solveYourselfAt(TimeStep *tStep)
{
    //
    // creates system of governing eq's and solves them at given time step
    //
    // first assemble problem at current time step

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling stiffness and mass matrices\n");
#endif

    if ( tStep->giveNumber() == 1 ) {
        //
        // first step  assemble stiffness Matrix
        //

        stiffnessMatrix = classFactory.createSparseMtrx(sparseMtrxType);
        stiffnessMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );

        massMatrix = classFactory.createSparseMtrx(sparseMtrxType);
        massMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );

        this->assemble( stiffnessMatrix, tStep, StiffnessMatrix,
                       EModelDefaultEquationNumbering(), this->giveDomain(1) );
        this->assemble( massMatrix, tStep, MassMatrix,
                       EModelDefaultEquationNumbering(), this->giveDomain(1) );
        //
        // create resulting objects eigVec and eigVal
        //
        eigVec.resize(this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ), numberOfRequiredEigenValues);
        eigVec.zero();
        eigVal.resize(numberOfRequiredEigenValues);
        eigVal.zero();
    }

    //
    // set-up numerical model
    //
    this->giveNumericalMethod( this->giveMetaStep( tStep->giveMetaStepNumber() ) );

    //
    // call numerical model to solve arised problem
    //
#ifdef VERBOSE
    OOFEM_LOG_INFO("Solving ...\n");
#endif

    nMethod->solve(stiffnessMatrix, massMatrix, & eigVal, & eigVec, rtolv, numberOfRequiredEigenValues);

    delete stiffnessMatrix;
    delete massMatrix;
    stiffnessMatrix = massMatrix = NULL;
}

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;
}

int
EIPrimaryUnknownMapper :: mapAndUpdate(FloatArray &answer, ValueModeType mode,
                                       Domain *oldd, Domain *newd,  TimeStep *tStep)
{
    int inode, nd_nnodes = newd->giveNumberOfDofManagers();
    int nsize = newd->giveEngngModel()->giveNumberOfDomainEquations( newd->giveNumber(), EModelDefaultEquationNumbering() );
    FloatArray unknownValues;
    IntArray dofidMask, locationArray;
    IntArray reglist;
#ifdef OOFEM_MAPPING_CHECK_REGIONS
    ConnectivityTable *conTable = newd->giveConnectivityTable();
    const IntArray *nodeConnectivity;
#endif

    answer.resize(nsize);
    answer.zero();

    for ( inode = 1; inode <= nd_nnodes; inode++ ) {
        DofManager *node = newd->giveNode(inode);
        /* HUHU CHEATING */
#ifdef __PARALLEL_MODE
        if ( ( node->giveParallelMode() == DofManager_null ) ||
            ( node->giveParallelMode() == DofManager_remote ) ) {
            continue;
        }

#endif

#ifdef OOFEM_MAPPING_CHECK_REGIONS
        // build up region list for node
        nodeConnectivity = conTable->giveDofManConnectivityArray(inode);
        reglist.resize( nodeConnectivity->giveSize() );
        reglist.clear();
        for ( int indx = 1; indx <= nodeConnectivity->giveSize(); indx++ ) {
            reglist.insertSortedOnce( newd->giveElement( nodeConnectivity->at(indx) )->giveRegionNumber() );
        }

#endif
        ///@todo Shouldn't we pass a primary field or something to this function?
        if ( this->evaluateAt(unknownValues, dofidMask, mode, oldd, * node->giveCoordinates(), reglist, tStep) ) {
            ///@todo This doesn't respect local coordinate systems in nodes. Supporting that would require major reworking.
            for ( int ii = 1; ii <= dofidMask.giveSize(); ii++ ) {
                // exclude slaves; they are determined from masters
                auto it = node->findDofWithDofId((DofIDItem)dofidMask.at(ii));
                if ( it != node->end() ) {
                    Dof *dof = *it;
                    if ( dof->isPrimaryDof() ) {
                        int eq = dof->giveEquationNumber(EModelDefaultEquationNumbering());
                        answer.at( eq ) += unknownValues.at(ii);
                    }
                }
            }
        } else {
            OOFEM_ERROR("evaluateAt service failed for node %d", inode);
        }
    }

    return 1;
}

void
NlDEIDynamic :: computeLoadVector(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    answer.resize( this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ) );
    answer.zero();

    //
    // Assemble the nodal part of load vector.
    //
    this->assembleVector( answer, tStep, ExternalForceAssembler(), mode,
                         EModelDefaultEquationNumbering(), this->giveDomain(1) );

    //
    // Exchange contributions.
    //
    this->updateSharedDofManagers(answer, EModelDefaultEquationNumbering(), LoadExchangeTag);
}

void
ParallelContext :: init(int newDi)
{
    di = newDi;
    ///@todo Should we even do this here? The user of the requested ParallelContext will just set this manually instead.
#ifdef __PARALLEL_MODE
    if ( emodel->isParallel() ) {
        ///@todo This shouldn't be hardcoded to just the default numbering schemes. In fact, this shouldn't even have "prescribed" and "free", just use the given numbering.
        n2g.init( emodel, di, EModelDefaultEquationNumbering() );
        n2l.init( emodel, di, EModelDefaultEquationNumbering() );

 #ifdef __VERBOSE_PARALLEL
        fprintf( stderr, "[%d] ParallelContext :: init - leq:%d, neq:%d, geq:%d\n", emodel->giveRank(), giveNumberOfLocalEqs(), giveNumberOfNaturalEqs(), giveNumberOfGlobalEqs() );
 #endif
    }

#endif
}

void
StokesFlowVelocityHomogenization :: computeTangent(FloatMatrix &answer, TimeStep *tStep)
{
    IntArray loc, col;

    Domain *domain = this->giveDomain(1);
    int nsd = domain->giveNumberOfSpatialDimensions();
    int ndof = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );

    // Build F matrix
    IntArray dofs(nsd);
    for ( int i = 0; i < nsd; ++i ) {
        dofs[i] = V_u + i; ///@todo This is a hack. We should have these as user input instead.
    }
    FloatMatrix F(ndof, nsd), Fe, N;
    col.enumerate(nsd);

    for ( auto &elem : domain->giveElements() ) {

        this->integrateNMatrix(N, *elem, tStep);
        
        elem->giveLocationArray( loc, dofs, EModelDefaultEquationNumbering() );
        Fe.beTranspositionOf(N);
        F.assemble(Fe, loc, col);
    }

    FloatMatrix H;

    std :: unique_ptr< SparseLinearSystemNM > solver( classFactory.createSparseLinSolver(solverType, this->giveDomain(1), this) );

    H.resize( F.giveNumberOfRows(), F.giveNumberOfColumns() );
    H.zero();

    // For correct homogenization, the tangent at the converged values should be used
    // (the one from the Newton iterations from solveYourselfAt is not updated to contain the latest values).
    SparseMtrxType stype = solver->giveRecommendedMatrix(true);
    std :: unique_ptr< SparseMtrx > Kff( classFactory.createSparseMtrx( stype ) );
    Kff->buildInternalStructure(this, domain->giveNumber(), EModelDefaultEquationNumbering() );
    this->assemble(*Kff, tStep, TangentStiffnessMatrix, EModelDefaultEquationNumbering(), domain);
    solver->solve(*Kff, F, H);

    answer.beTProductOf(H, F);
    answer.times( 1. / this->giveAreaOfRVE() );
}

void
StructuralEngngModel :: giveInternalForces(FloatArray &answer, bool normFlag, int di, TimeStep *tStep)
{
    // Simply assembles contributions from each element in domain
    Domain *domain = this->giveDomain(di);
    // Update solution state counter
    tStep->incrementStateCounter();

    answer.resize( this->giveNumberOfDomainEquations( di, EModelDefaultEquationNumbering() ) );
    answer.zero();
    this->assembleVector(answer, tStep, InternalForcesVector, VM_Total,
                         EModelDefaultEquationNumbering(), domain, normFlag ? & this->internalForcesEBENorm : NULL);

    // Redistributes answer so that every process have the full values on all shared equations
    this->updateSharedDofManagers(answer, EModelDefaultEquationNumbering(), InternalForcesExchangeTag);

    // Remember last internal vars update time stamp.
    internalVarUpdateStamp = tStep->giveSolutionStateCounter();
}

void StaticStructural :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    if ( cmpn == InternalRhs ) {
        // Updates the solution in case it has changed 
        ///@todo NRSolver should report when the solution changes instead of doing it this way.
        this->field->update(VM_Total, tStep, this->solution, EModelDefaultEquationNumbering());
        this->field->applyBoundaryCondition(tStep);///@todo Temporary hack to override the incorrect vavues that is set by "update" above. Remove this when that is fixed.

        this->internalForces.zero();
        this->assembleVector(this->internalForces, tStep, InternalForceAssembler(), VM_Total,
                             EModelDefaultEquationNumbering(), d, & this->eNorm);
        this->updateSharedDofManagers(this->internalForces, EModelDefaultEquationNumbering(), InternalForcesExchangeTag);

        internalVarUpdateStamp = tStep->giveSolutionStateCounter(); // Hack for linearstatic
    } else if ( cmpn == NonLinearLhs ) {
        this->stiffnessMatrix->zero();
        this->assemble(*this->stiffnessMatrix, tStep, TangentAssembler(TangentStiffness), EModelDefaultEquationNumbering(), d);
    } else {
        OOFEM_ERROR("Unknown component");
    }
}

void DarcyFlow :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    switch ( cmpn ) {
    case InternalRhs:
        this->internalForces.zero();
        this->assembleVector(this->internalForces, tStep,  InternalForcesVector, VM_Total,
                             EModelDefaultEquationNumbering(), d, & this->ebeNorm);
        this->updateSharedDofManagers(this->externalForces, EModelDefaultEquationNumbering(), InternalForcesExchangeTag);
        break;

    case NonLinearLhs:

        this->stiffnessMatrix->zero();
        this->assemble( *this->stiffnessMatrix, tStep, TangentStiffnessMatrix,
                        EModelDefaultEquationNumbering(), this->giveDomain(1) );
        break;

    default:
        OOFEM_ERROR("Unknown component id (%d)", ( int ) cmpn);
    }
}

int
ParallelContext :: giveNumberOfGlobalEqs()
{
#ifdef __PARALLEL_MODE
    if ( emodel->isParallel() ) {
        return n2g.giveNumberOfGlobalEqs();
    } else {
#endif
    return emodel->giveNumberOfDomainEquations( di, EModelDefaultEquationNumbering() );

#ifdef __PARALLEL_MODE
}
#endif
}

void
StationaryTransportProblem :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    if ( cmpn == InternalRhs ) {
        this->internalForces.zero();
        this->assembleVector( this->internalForces, tStep, EID_ConservationEquation, InternalForcesVector, VM_Total,
                              EModelDefaultEquationNumbering(), this->giveDomain(1), & this->eNorm );
        return;
    } else if ( cmpn == NonLinearLhs ) {
        if ( !this->keepTangent ) {
            // Optimization for linear problems, we can keep the old matrix (which could save its factorization)
            this->conductivityMatrix->zero();
            ///@todo We should use some problem-neutral names instead of "ConductivityMatrix" (and something nicer for LHSBCMatrix)
            this->assemble( conductivityMatrix, tStep, EID_ConservationEquation, ConductivityMatrix,
                        EModelDefaultEquationNumbering(), this->giveDomain(1) );
            this->assemble( conductivityMatrix, tStep, EID_ConservationEquation, LHSBCMatrix,
                        EModelDefaultEquationNumbering(), this->giveDomain(1) );
        }
        return;
    } else {
        OOFEM_ERROR("StationaryTransportProblem::updateComponent - Unknown component");
    }
}

void FreeWarping :: 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

        this->initializeCommMaps();
    }

    StructuralEngngModel :: solveYourself();
}

void StokesFlow :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    velocityPressureField->update(VM_Total, tStep, solutionVector, EModelDefaultEquationNumbering());

    // update element stabilization
    for ( auto &elem : d->giveElements() ) {
        static_cast< FMElement * >( elem.get() )->updateStabilizationCoeffs(tStep);
    }

    if ( cmpn == InternalRhs ) {
        this->internalForces.zero();
        this->assembleVector(this->internalForces, tStep, InternalForceAssembler(), VM_Total,
                             EModelDefaultEquationNumbering(), d, & this->eNorm);
        this->updateSharedDofManagers(this->internalForces, EModelDefaultEquationNumbering(), InternalForcesExchangeTag);
        return;
    } else if ( cmpn == NonLinearLhs ) {
        this->stiffnessMatrix->zero();
        this->assemble(*stiffnessMatrix, tStep, TangentAssembler(TangentStiffness),
                       EModelDefaultEquationNumbering(), d);
        return;
    } else {
        OOFEM_ERROR("Unknown component");
    }
}

void StokesFlow :: updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
{
    // update element stabilization
    int nelem = d->giveNumberOfElements();
    for ( int i = 1; i <= nelem; ++i ) {
        static_cast< FMElement* >( d->giveElement(i) )->updateStabilizationCoeffs(tStep);
    }

    if ( cmpn == InternalRhs ) {
        this->internalForces.zero();
        this->assembleVector( this->internalForces, tStep, EID_MomentumBalance_ConservationEquation, InternalForcesVector, VM_Total,
                              EModelDefaultEquationNumbering(), this->giveDomain(1), &this->eNorm );
        return;

    } else if ( cmpn == NonLinearLhs ) {
        this->stiffnessMatrix->zero();
        this->assemble(this->stiffnessMatrix, tStep, EID_MomentumBalance_ConservationEquation, StiffnessMatrix,
                EModelDefaultEquationNumbering(), d);
        return;

    } else {
        OOFEM_ERROR("StokesFlow::updateComponent - Unknown component");
    }
}

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();
}