void
NonLinearStatic :: updateLoadVectors(TimeStep *tStep)
{
    MetaStep *mstep = this->giveMetaStep( tStep->giveMetaStepNumber() );
    bool isLastMetaStep = ( tStep->giveNumber() == mstep->giveLastStepNumber() );

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

NM_Status
IMLSolver :: solve(SparseMtrx &A, FloatArray &b, FloatArray &x)
{
    int result;

    if ( x.giveSize() != b.giveSize() ) {
        OOFEM_ERROR("size mismatch");
    }


    // check preconditioner
    if ( M ) {
        if ( ( precondInit ) || ( Lhs != &A ) || ( this->lhsVersion != A.giveVersion() ) ) {
            M->init(A);
        }
    } else {
        OOFEM_ERROR("preconditioner creation error");
    }

    Lhs = &A;
    this->lhsVersion = A.giveVersion();

#ifdef TIME_REPORT
    Timer timer;
    timer.startTimer();
#endif


    if ( solverType == IML_ST_CG ) {
        int mi = this->maxite;
        double t = this->tol;
        result = CG(* Lhs, x, b, * M, mi, t);
        OOFEM_LOG_INFO("CG(%s): flag=%d, nite %d, achieved tol. %g\n", M->giveClassName(), result, mi, t);
    } else if ( solverType == IML_ST_GMRES ) {
        int mi = this->maxite, restart = 100;
        double t = this->tol;
        FloatMatrix H(restart + 1, restart); // storage for upper Hesenberg
        result = GMRES(* Lhs, x, b, * M, H, restart, mi, t);
        OOFEM_LOG_INFO("GMRES(%s): flag=%d, nite %d, achieved tol. %g\n", M->giveClassName(), result, mi, t);
    } else {
        OOFEM_ERROR("unknown lsover type");
    }

#ifdef TIME_REPORT
    timer.stopTimer();
    OOFEM_LOG_INFO( "IMLSolver info: user time consumed by solution: %.2fs\n", timer.getUtime() );
#endif


    //solved = 1;
    return NM_Success;
}

int
SloanGraph :: findBestRoot()
/* this is a very expensive method */
{
    int BestRoot = 0;
    int Diameter = 0;
    int i, nnodes = domain->giveNumberOfDofManagers();
    SloanLevelStructure *LSC;
    clock_t time_1, time_0 = :: clock();
    for ( i = 1; i <= nnodes; i++ ) {
        LSC = new SloanLevelStructure(this, i);
        int Depth = LSC->giveDepth();
        if ( Depth > Diameter ) {
            Diameter = Depth;
            BestRoot = i;
        }

        delete LSC;
        time_1 = :: clock();
        if ( ( time_1 - time_0 ) / CLOCKS_PER_SEC > SLOAN_TIME_CHUNK ) {
            OOFEM_LOG_INFO("%d roots (%5.1f per cent) checked: largest pseudo-diameter = %d\n", i, float ( 100 * i ) / nnodes, Diameter);
            //fflush(stdout);
            //if (get_yes_or_no() == NO) break;
            time_0 = time_1;
        }
    }

    return BestRoot;
}

int
Targe2Interface :: createInput(Domain *d, TimeStep *stepN)
{
    int nelem = d->giveNumberOfElements();
    FILE *outputStrem;
    Element *ielem;
    RemeshingCriteria *rc = d->giveErrorEstimator()->giveRemeshingCrit();

    outputStrem = fopen("targe2.bmf", "w");
    // print header for 2D

    for ( int i = 1; i <= nelem; i++ ) {
        ielem = d->giveElement(i);
        fprintf( outputStrem, "MC-T %e %e %e %e %e %e   %e %e %e\n",
                ielem->giveNode(1)->giveCoordinate(1), ielem->giveNode(1)->giveCoordinate(2),
                ielem->giveNode(2)->giveCoordinate(1), ielem->giveNode(2)->giveCoordinate(2),
                ielem->giveNode(3)->giveCoordinate(1), ielem->giveNode(3)->giveCoordinate(2),
                rc->giveRequiredDofManDensity(ielem->giveNode(1)->giveNumber(), stepN),
                rc->giveRequiredDofManDensity(ielem->giveNode(2)->giveNumber(), stepN),
                rc->giveRequiredDofManDensity(ielem->giveNode(3)->giveNumber(), stepN) );
    }

    fclose(outputStrem);

    OOFEM_LOG_INFO("Targe2 .bmf file created\n");
    return 1;
}

contextIOResultType LinearStability :: restoreContext(DataStream *stream, ContextMode mode, void *obj)
//
// restore state variable - displacement vector
//
{
    int activeVector, version;
    int istep = 1, iversion = 1;
    int closeFlag = 0;
    contextIOResultType iores;
    FILE *file = NULL;

    this->resolveCorrespondingStepNumber(activeVector, version, obj);
    if ( eigVal.isEmpty() ) { // not restored before
        if ( stream == NULL ) {
            if ( !this->giveContextFile(& file, istep, iversion, contextMode_read) ) {
                THROW_CIOERR(CIO_IOERR); // override
            }

            stream = new FileDataStream(file);
            closeFlag = 1;
        }

        if ( ( iores = StructuralEngngModel :: restoreContext(stream, mode, ( void * ) & istep) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( ( iores = displacementVector.restoreYourself(stream, mode) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( ( iores = eigVal.restoreYourself(stream, mode) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( ( iores = eigVec.restoreYourself(stream, mode) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( closeFlag ) {
            fclose(file);
            delete stream;
            stream = NULL;
        }                                                    // ensure consistent records

    }

    //  if (istep > numberOfRequiredEigenValues) istep = numberOfRequiredEigenValues ;
    //  printf( "Restoring - corresponding index is %d, EigenValue is %lf\n",
    //     istep,eigVal.at(istep));
    //  setActiveVector (istep);
    if ( activeVector > numberOfRequiredEigenValues ) {
        activeVector = numberOfRequiredEigenValues;
    }

    OOFEM_LOG_INFO( "Restoring - corresponding index is %d, EigenValue is %f\n",
                   activeVector, eigVal.at(activeVector) );
    this->giveCurrentStep()->setTime( ( double ) activeVector );

    return CIO_OK;
}

int
AdaptiveNonLinearStatic :: initializeAdaptive(int tStepNumber)
{
    int stepinfo [ 2 ];

    stepinfo [ 0 ] = tStepNumber;
    stepinfo [ 1 ] = 0;

    try {
        this->restoreContext(NULL, CM_State, ( void * ) stepinfo);
    } catch(ContextIOERR & c) {
        c.print();
        exit(1);
    }

    this->initStepIncrements();

    int sernum = this->giveDomain(1)->giveSerialNumber();
    OOFEM_LOG_INFO("restoring domain %d.%d\n", 1, sernum + 1);
    Domain *dNew = new Domain(2, sernum + 1, this);
    DataReader *domainDr = this->GiveDomainDataReader(1, sernum + 1, contextMode_read);
    if ( !dNew->instanciateYourself(domainDr) ) {
        OOFEM_ERROR("domain Instanciation failed");
    }

    delete domainDr;

    // remap solution to new domain
    return this->adaptiveRemap(dNew);
}

void
ErrorCheckingExportModule :: doOutput(TimeStep *tStep, bool forcedOutput)
{
#if 0
    if ( !( testTimeStepOutput(tStep) || forcedOutput ) ) {
        return;
    }
#endif

    // Error checking rules are hardcoded to domain 1 always.
    Domain *domain = emodel->giveDomain(1);

    OOFEM_LOG_INFO("Checking rules...\n");
    for ( auto &rule: this->errorCheckingRules ) {
        this->allPassed &= rule->check(domain, tStep);
    }

    if ( !tStep->isNotTheLastStep() ) {
        if ( !this->allPassed ) {
            OOFEM_ERROR("Rule not passed, exiting with error");
        }
    }

    if ( this->writeChecks ) {
        this->writeCheck(domain, tStep);
    }
}

TimeStep *
CBS :: giveNextStep()
{
    double dt = deltaT;
    if ( !currentStep ) {
        // first step -> generate initial step
        currentStep.reset( new TimeStep( *giveSolutionStepWhenIcApply() ) );
    }

    previousStep = std :: move(currentStep);

    Domain *domain = this->giveDomain(1);
    // check for critical time step
    for ( auto &elem : domain->giveElements() ) {
        dt = min( dt, static_cast< CBSElement & >( *elem ).computeCriticalTimeStep(previousStep.get()) );
    }

    dt *= 0.6;
    dt = max(dt, minDeltaT);
    dt /= this->giveVariableScale(VST_Time);

    currentStep.reset( new TimeStep(*previousStep, dt) );

    OOFEM_LOG_INFO( "SolutionStep %d : t = %e, dt = %e\n", currentStep->giveNumber(),
                    currentStep->giveTargetTime() * this->giveVariableScale(VST_Time), dt * this->giveVariableScale(VST_Time) );

    return currentStep.get();
}

double
Concrete3 :: computeStrength(GaussPoint *gp, double charLength)
//
// computes strength for given gp,
// which may be reduced according to length of "fracture process zone"
// to be energetically correct
//
{
    double Ee, Gf, Ft;

    Ee = this->give(pscm_Ee, gp);
    Gf = this->give(pscm_Gf, gp);
    Ft = this->give(pscm_Ft, gp);

    if ( this->checkSizeLimit(gp, charLength) ) {
        ;
    } else {
        // we reduce Ft and there is no softening but sudden drop
        Ft = sqrt(2. * Ee * Gf / charLength);
        //
        OOFEM_LOG_INFO("Reducing Ft to %f in element %d, gp %d, Le %f",
               Ft, gp->giveElement()->giveNumber(), gp->giveNumber(), charLength);
        //
    }

    return Ft;
}

contextIOResultType EigenValueDynamic :: restoreContext(DataStream *stream, ContextMode mode, void *obj)
//
// restore state variable - displacement vector
//
{
    int closeFlag = 0;
    int activeVector = this->resolveCorrespondingEigenStepNumber(obj);
    int istep = 1, iversion = 0;
    contextIOResultType iores;
    FILE *file;

    if ( restoreFlag == 0 ) { // not restored before
        if ( stream == NULL ) {
            if ( !this->giveContextFile(& file, istep, iversion, contextMode_read) ) {
                THROW_CIOERR(CIO_IOERR); // override
            }

            stream = new FileDataStream(file);
            closeFlag = 1;
        }

        // save element context

        if ( ( iores = EngngModel :: restoreContext(stream, mode, ( void * ) & istep) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( ( iores = eigVal.restoreYourself(stream, mode) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( ( iores = eigVec.restoreYourself(stream, mode) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }

        if ( closeFlag ) {
            fclose(file);
            delete stream;
            stream = NULL;
        } // ensure consistent records

    }

    if ( activeVector > numberOfRequiredEigenValues ) {
        activeVector = numberOfRequiredEigenValues;
    }

    OOFEM_LOG_INFO( "Restoring - corresponding index is %d, EigenValue is %f\n", activeVector, eigVal.at(activeVector) );
    this->giveCurrentStep()->setTime( ( double ) activeVector );
    this->restoreFlag = 1;

    return CIO_OK;
}

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

    // !!! should use separate field, e.g. hydramname #hydramnumber
    // Hydration>0  ->  Model starting at value, hydration<0 -> Constant at given value
    value = -2.;
    constantHydrationDegree = 1.0;
    IR_GIVE_OPTIONAL_FIELD(ir, value, _IFT_HydrationModelInterface_hydration);
    if ( value >= 0. ) {
        OOFEM_LOG_INFO("HydratingMaterial: creating HydrationModel.");
        hydrationModel.reset( new HydrationModel() );
        if ( !hydrationModel ) {
            OOFEM_WARNING("Could not create HydrationModel instance.");
            return IRRT_BAD_FORMAT;
        }

        hydrationModel->initializeFrom(ir);
    }
    // constant hydration degree
    else if ( value >= -1. ) {
        constantHydrationDegree = -value;
        OOFEM_LOG_INFO("HydratingMaterial: Hydration degree set to %.2f.", -value);
    } else {
        OOFEM_WARNING("Hydration degree input incorrect, use -1..<0 for constant hydration degree, 0..1 to set initial material hydration degree.");
        return IRRT_BAD_FORMAT;
    }

    // Material cast time - start of hydration
    // 11/3/2004 OK *unfinished in Hellmat, needs to be checked in hm_Interface->updateInternalState
    castAt = 0.;
    IR_GIVE_OPTIONAL_FIELD(ir, castAt, _IFT_HydrationModelInterface_castAt);
    if ( castAt >= 0. ) {
        OOFEM_LOG_INFO("HydratingMaterial: Hydration starts at time %.2g.", castAt);
    }

    return IRRT_OK;
}

void
NodalAveragingRecoveryModel :: initCommMaps()
{
    if ( initCommMap ) {
        EngngModel *emodel = domain->giveEngngModel();
        commBuff = new CommunicatorBuff(emodel->giveNumberOfProcesses(), CBT_dynamic);
        communicator = new NodeCommunicator(emodel, commBuff, emodel->giveRank(),
                                            emodel->giveNumberOfProcesses());
        communicator->setUpCommunicationMaps(domain->giveEngngModel(), true, true);
        OOFEM_LOG_INFO("NodalAveragingRecoveryModel :: initCommMaps: initialized comm maps\n");
        initCommMap = false;
    }
}

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

Interface *
QTRSpace :: giveInterface(InterfaceType interface)
{
    if ( interface == ZZNodalRecoveryModelInterfaceType ) {
        return static_cast< ZZNodalRecoveryModelInterface * >( this );
    } else if ( interface == SPRNodalRecoveryModelInterfaceType ) {
        return static_cast< SPRNodalRecoveryModelInterface * >( this );
    } else if ( interface == NodalAveragingRecoveryModelInterfaceType ) {
        return static_cast< NodalAveragingRecoveryModelInterface * >( this );
    }

    OOFEM_LOG_INFO("Interface on QTRSpace element not supported");
    return NULL;
}

void
CBS :: applyIC(TimeStep *stepWhenIcApply)
{
    Domain *domain = this->giveDomain(1);
    int mbneq = this->giveNumberOfDomainEquations(1, vnum);
    int pdneq = this->giveNumberOfDomainEquations(1, pnum);
    FloatArray *velocityVector, *pressureVector;

#ifdef VERBOSE
    OOFEM_LOG_INFO("Applying initial conditions\n");
#endif

    VelocityField.advanceSolution(stepWhenIcApply);
    velocityVector = VelocityField.giveSolutionVector(stepWhenIcApply);
    velocityVector->resize(mbneq);
    velocityVector->zero();

    PressureField.advanceSolution(stepWhenIcApply);
    pressureVector = PressureField.giveSolutionVector(stepWhenIcApply);
    pressureVector->resize(pdneq);
    pressureVector->zero();

    for ( auto &node : domain->giveDofManagers() ) {
        for ( Dof *iDof: *node ) {
            // ask for initial values obtained from
            // bc (boundary conditions) and ic (initial conditions)
            if ( !iDof->isPrimaryDof() ) {
                continue;
            }

            int jj = iDof->__giveEquationNumber();
            if ( jj ) {
                DofIDItem type = iDof->giveDofID();
                if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) {
                    velocityVector->at(jj) = iDof->giveUnknown(VM_Total, stepWhenIcApply);
                } else {
                    pressureVector->at(jj) = iDof->giveUnknown(VM_Total, stepWhenIcApply);
                }
            }
        }
    }

    // update element state according to given ic
    for ( auto &elem : domain->giveElements() ) {
        CBSElement *element = static_cast< CBSElement * >( elem.get() );
        element->updateInternalState(stepWhenIcApply);
        element->updateYourself(stepWhenIcApply);
    }
}

int
FreemInterface :: createInput(Domain *d, TimeStep *stepN)
{
    int i;
    int nnodes = d->giveNumberOfDofManagers(), nelem = d->giveNumberOfElements();
    double density;
    FILE *outputStrem;
    Node *inode;
    Element *ielem;

    outputStrem = fopen("freem.bmf", "w");
    // print header for 2D
    fprintf(outputStrem, "nbnodes %d nbelem %d \n", nnodes, nelem);

    /* mesh densities smoothing */
    // query nodal absolute densities
    FloatArray nodalDensities(nnodes);
    for ( i = 1; i <= nnodes; i++ ) {
        nodalDensities.at(i) = d->giveErrorEstimator()->giveRemeshingCrit()->giveRequiredDofManDensity(i, stepN);
    }

    this->smoothNodalDensities(d,  nodalDensities, stepN);
    /* end of smoothing */

    // loop over nodes
    for ( i = 1; i <= nnodes; i++ ) {
        //density = d->giveErrorEstimator ()->giveRemeshingCrit()->giveRequiredDofManDensity (i, stepN, 1);
        //density = d->giveErrorEstimator ()->giveRemeshingCrit()->giveDofManDensity (i) / nodalDensities.at(i);
        density = nodalDensities.at(i);
        inode = d->giveNode(i);
        fprintf(outputStrem, "backgroungMeshNode %d x %e y %e density %e\n", i, inode->giveCoordinate(1), inode->giveCoordinate(2), density);
    }

    for ( i = 1; i <= nelem; i++ ) {
        ielem = d->giveElement(i);
        if ( ielem->giveClassID() != PlaneStress2dClass ) {
            OOFEM_ERROR("FreemInterface::createInput : unsupported element type");
        }

        fprintf( outputStrem, "backgroundMeshElem %d  nodes 4 %d %d %d %d\n", i,
                ielem->giveNode(1)->giveNumber(), ielem->giveNode(2)->giveNumber(),
                ielem->giveNode(3)->giveNumber(), ielem->giveNode(4)->giveNumber() );
    }

    fclose(outputStrem);

    OOFEM_LOG_INFO("freem.bmf file created\n");
    return 1;
}

OOFEMTXTDataReader :: OOFEMTXTDataReader(std :: string inputfilename) : DataReader(),
    dataSourceName(std :: move(inputfilename)), recordList()
{
    std :: list< std :: pair< int, std :: string > >lines;
    // Read all the lines in the main input file:
    {
        std :: ifstream inputStream(dataSourceName);
        if ( !inputStream.is_open() ) {
            OOFEM_ERROR("Can't open input stream (%s)", dataSourceName.c_str());
        }

        int lineNumber = 0;
        std :: string line;

        this->giveRawLineFromInput(inputStream, lineNumber, outputFileName);
        this->giveRawLineFromInput(inputStream, lineNumber, description);

        while (this->giveLineFromInput(inputStream, lineNumber, line)) {
            lines.emplace_back(make_pair(lineNumber, line));
        }
    }
    // Check for included files: @include "somefile"
    for ( auto it = lines.begin(); it != lines.end(); ++it ) {
        if ( it->second.compare(0, 8, "@include") == 0 ) {
            std :: string fname = it->second.substr(10, it->second.length()-11);
            OOFEM_LOG_INFO("Reading included file: %s\n", fname.c_str());

            // Remove the include line
            lines.erase(it++);
            // Add all the included lines:
            int includedLine = 0;
            std :: string line;
            std :: ifstream includedStream(fname);
            if ( !includedStream.is_open() ) {
                OOFEM_ERROR("Can't open input stream (%s)", fname.c_str());
            }
            while (this->giveLineFromInput(includedStream, includedLine, line)) {
                lines.emplace(it, make_pair(includedLine, line));
            }
        }
    }
    ///@todo This could be parallelized, but I'm not sure it is worth it 
    /// (might make debugging faulty input files harder for users as well)
    for ( auto &line: lines ) {
        //printf("line: %s\n", line.second.c_str());
        this->recordList.emplace_back(line.first, line.second);
    }
    this->it = this->recordList.begin();
}

void
SPRNodalRecoveryModel :: initCommMaps()
{
 #ifdef __PARALLEL_MODE
    if ( initCommMap ) {
        EngngModel *emodel = domain->giveEngngModel();
        commBuff = new CommunicatorBuff(emodel->giveNumberOfProcesses(), CBT_dynamic);
        communicator = new NodeCommunicator(emodel, commBuff, emodel->giveRank(),
                                            emodel->giveNumberOfProcesses());
        communicator->setUpCommunicationMaps(domain->giveEngngModel(), true, true);
        OOFEM_LOG_INFO("SPRNodalRecoveryModel :: initCommMaps: initialized comm maps");
        initCommMap = false;
    }

 #endif
}

Interface *
QWedge_ht :: giveInterface(InterfaceType interface)
{
    if ( interface == ZZNodalRecoveryModelInterfaceType ) {
        return static_cast< ZZNodalRecoveryModelInterface * >(this);
    } else if ( interface == SPRNodalRecoveryModelInterfaceType ) {
        return static_cast< SPRNodalRecoveryModelInterface * >(this);
    } else if ( interface == NodalAveragingRecoveryModelInterfaceType ) {
        return static_cast< NodalAveragingRecoveryModelInterface * >(this);
    } else if ( interface == SpatialLocalizerInterfaceType ) {
        return static_cast< SpatialLocalizerInterface * >(this);
    }

    OOFEM_LOG_INFO("Interface on Lwedge element not supported");
    return NULL;
}