自学教程：C++ EModelDefaultEquationNumbering函数代码示例

voidNonLinearStatic :: 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);}

voidSUPG :: 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);}

voidTransientTransportProblem :: 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");    }}

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

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

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

voidNlDEIDynamic :: 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);}

voidParallelContext :: 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}

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

voidStructuralEngngModel :: 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);    }}

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

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

voidNlDEIDynamic :: computeLoadVector(FloatArray &answer, ValueModeType mode, TimeStep *stepN){    answer.resize( this->giveNumberOfEquations(EID_MomentumBalance) );    answer.zero();    //    // Assemble the nodal part of load vector.    //    this->assembleVector( answer, stepN, EID_MomentumBalance, ExternalForcesVector, mode,                          EModelDefaultEquationNumbering(), this->giveDomain(1));    //    // Exchange contributions.    //#ifdef __PARALLEL_MODE    this->updateSharedDofManagers( answer, LoadExchangeTag );#endif}

intTrabBoneNL3D :: giveLocalNonlocalStiffnessContribution(GaussPoint *gp, IntArray &loc, const UnknownNumberingScheme &s,                                                       FloatArray &lcontrib, TimeStep *tStep){    TrabBoneNL3DStatus *nlStatus = static_cast< TrabBoneNL3DStatus * >( this->giveStatus(gp) );    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );    int nrows, nsize;    double sum, nlKappa, dDamFunc, dam, tempDam;    FloatArray localNu;    FloatMatrix b;    this->computeCumPlastStrain(nlKappa, gp, tStep);    dam = nlStatus->giveDam();    tempDam = nlStatus->giveTempDam();    if ( ( tempDam - dam ) > 0.0 ) {        elem->giveLocationArray(loc, s);        localNu = nlStatus->giveTempEffectiveStress();        elem->giveLocationArray(loc, EModelDefaultEquationNumbering() );        elem->computeBmatrixAt(gp, b);        dDamFunc = expDam * critDam * exp(-expDam * nlKappa);        nrows = b.giveNumberOfColumns();        nsize = localNu.giveSize();        lcontrib.resize(nrows);        for ( int i = 1; i <= nrows; i++ ) {            sum = 0.0;            for ( int j = 1; j <= nsize; j++ ) {                sum += b.at(j, i) * localNu.at(j);            }            lcontrib.at(i) = mParam * dDamFunc * sum;        }        return 1;    } else {        loc.clear();        return 0;    }}

void XFEMStatic :: setValsFromDofMap(FloatArray &oArray, const FloatArray &iArray){    int neq = 0;    for ( int domainIndex = 1; domainIndex <= this->giveNumberOfDomains(); domainIndex++ ) {        neq += this->giveNumberOfDomainEquations( domainIndex, EModelDefaultEquationNumbering() );    }    int numEqOld = iArray.giveSize();    printf("Setting values from dof map. neq: %d numEqOld: %d/n", neq, numEqOld);    oArray.resize(neq);    oArray.zero();    for ( int domainIndex = 1; domainIndex <= this->giveNumberOfDomains(); domainIndex++ ) {        Domain *domain = this->giveDomain(domainIndex);        for ( int dManIndex = 1; dManIndex <= domain->giveNumberOfDofManagers(); dManIndex++ ) {            DofManager *dMan = domain->giveDofManager(dManIndex);            for ( Dof *dof: *dMan ) {                int eqNumNew = dof->giveEqn();                if ( eqNumNew > 0 ) {                    std :: vector< int > key(3);                    key [ 0 ] = domainIndex;                    key [ 1 ] = dManIndex;                    key [ 2 ] = dof->giveDofID();                    if ( mDofEqnNumMap.find(key) != mDofEqnNumMap.end() ) {                        int eqNumOld = mDofEqnNumMap [ key ];                        //                      printf("eqNumNew: %d eqNumOld: %d/n", eqNumNew, eqNumOld);                        if ( eqNumOld > 0 && eqNumOld <= numEqOld  ) {                            oArray.at(eqNumNew) = iArray.at(eqNumOld);                        }                    }                }            }        }    }}

voidNonStationaryTransportProblem :: updateInternalState(TimeStep *tStep){    for ( auto &domain: domainList ) {        if ( requiresUnknownsDictionaryUpdate() ) {            //update temperature vector            UnknownsField->update( VM_Total, tStep, * ( this->UnknownsField->giveSolutionVector(tStep) ), EModelDefaultEquationNumbering() );            //update Rhs vector            UnknownsField->update(VM_RhsTotal, tStep, bcRhs, EModelDefaultEquationNumbering());        }        if ( internalVarUpdateStamp != tStep->giveSolutionStateCounter() ) {            for ( auto &elem : domain->giveElements() ) {                elem->updateInternalState(tStep);            }            internalVarUpdateStamp = tStep->giveSolutionStateCounter();        }    }}

void NonStationaryTransportProblem :: solveYourselfAt(TimeStep *tStep){    // Creates system of governing eq's and solves them at given tStep    // The solution is stored in UnknownsField. If the problem is growing/decreasing, the UnknownsField is projected on DoFs when needed.    // If equations are not renumbered, the algorithm is efficient without projecting unknowns to DoFs (nodes).    //Right hand side    FloatArray rhs;    int neq = this->giveNumberOfEquations(EID_ConservationEquation);#ifdef VERBOSE    OOFEM_LOG_RELEVANT( "Solving [step number %8d, time %15e]/n", tStep->giveNumber(), tStep->giveTargetTime() );#endif    //Solution at the first time step needs history. Therefore, return back one time increment and create it.    if ( tStep->isTheFirstStep() ) {        this->giveSolutionStepWhenIcApply();        bcRhs.resize(neq); //rhs vector from solution step i-1        bcRhs.zero();        this->applyIC(stepWhenIcApply);        //project initial conditions to have temorary temperature in integration points                //edge or surface load on elements        this->assembleVectorFromElements( bcRhs, stepWhenIcApply, EID_ConservationEquation, ElementBCTransportVector,                                         VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );        //add prescribed value, such as temperature, on nodes        this->assembleDirichletBcRhsVector( bcRhs, stepWhenIcApply, EID_ConservationEquation, VM_Total,                                           NSTP_MidpointLhs, EModelDefaultEquationNumbering(), this->giveDomain(1) );        //add internal source vector on elements        this->assembleVectorFromElements( bcRhs, stepWhenIcApply, EID_ConservationEquation, ElementInternalSourceVector,                                         VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );        //add nodal load        this->assembleVectorFromDofManagers( bcRhs, stepWhenIcApply, EID_ConservationEquation, ExternalForcesVector,                                            VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );    }    //Create a new lhs matrix if necessary    if ( tStep->isTheFirstStep() || this->changingProblemSize ) {        if ( conductivityMatrix ) {            delete conductivityMatrix;        }        conductivityMatrix = CreateUsrDefSparseMtrx(sparseMtrxType);        if ( conductivityMatrix == NULL ) {            _error("solveYourselfAt: sparse matrix creation failed");        }        conductivityMatrix->buildInternalStructure( this, 1, EID_ConservationEquation, EModelDefaultEquationNumbering() );#ifdef VERBOSE        OOFEM_LOG_INFO("Assembling conductivity and capacity matrices/n");#endif        this->assemble( conductivityMatrix, stepWhenIcApply, EID_ConservationEquation, LHSBCMatrix,                       EModelDefaultEquationNumbering(), this->giveDomain(1) );        conductivityMatrix->times(alpha);        this->assemble( conductivityMatrix, stepWhenIcApply, EID_ConservationEquation, NSTP_MidpointLhs,                       EModelDefaultEquationNumbering(), this->giveDomain(1) );    }    //obtain the last Rhs vector from DoFs directly    if ( !tStep->isTheFirstStep() && this->changingProblemSize ) {        UnknownsField->initialize(VM_RhsTotal, tStep, bcRhs);    }    //prepare position in UnknownsField to store the results    UnknownsField->advanceSolution(tStep);    FloatArray *solutionVector = UnknownsField->giveSolutionVector(tStep);    solutionVector->resize(neq);    solutionVector->zero();    #ifdef VERBOSE    OOFEM_LOG_INFO("Assembling rhs/n");#endif    // assembling load from elements    rhs = bcRhs;    rhs.times(1. - alpha);    bcRhs.zero();    this->assembleVectorFromElements( bcRhs, tStep, EID_ConservationEquation, ElementBCTransportVector,                                     VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );    this->assembleDirichletBcRhsVector( bcRhs, tStep, EID_ConservationEquation, VM_Total, NSTP_MidpointLhs,                                       EModelDefaultEquationNumbering(), this->giveDomain(1) );    this->assembleVectorFromElements( bcRhs, tStep, EID_ConservationEquation, ElementInternalSourceVector,                                     VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );    // assembling load from nodes    this->assembleVectorFromDofManagers( bcRhs, tStep, EID_ConservationEquation, InternalForcesVector, VM_Total,                                        EModelDefaultEquationNumbering(), this->giveDomain(1) );    for ( int i = 1; i <= neq; i++ ) {        rhs.at(i) += bcRhs.at(i) * alpha;    }    // add the rhs part depending on previous solution    assembleAlgorithmicPartOfRhs( rhs, EID_ConservationEquation,                                 EModelDefaultEquationNumbering(), tStep->givePreviousStep() );    // set-up numerical model//.........这里部分代码省略.........

void FreeWarping :: solveYourselfAt(TimeStep *tStep){    //    // creates system of governing eq's and solves them at given time step    //    // first assemble problem at current time step    if ( initFlag ) {#ifdef VERBOSE        OOFEM_LOG_DEBUG("Assembling stiffness matrix/n");#endif        //        // first step  assemble stiffness Matrix        //        stiffnessMatrix = classFactory.createSparseMtrx(sparseMtrxType);        if ( stiffnessMatrix == NULL ) {            OOFEM_ERROR("sparse matrix creation failed");        }        stiffnessMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );        this->assemble( *stiffnessMatrix, tStep, TangentAssembler(TangentStiffness),                       EModelDefaultEquationNumbering(), this->giveDomain(1) );        initFlag = 0;    }#ifdef VERBOSE    OOFEM_LOG_DEBUG("Assembling load/n");#endif    //    // allocate space for displacementVector    //    displacementVector.resize( this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ) );    displacementVector.zero();    //    // assembling the load vector    //    loadVector.resize( this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ) );    loadVector.zero();    this->assembleVector( loadVector, tStep, ExternalForceAssembler(), VM_Total,                         EModelDefaultEquationNumbering(), this->giveDomain(1) );    //    // internal forces (from Dirichlet b.c's, or thermal expansion, etc.)    //    // no such forces exist for the free warping problem    /*    FloatArray internalForces( this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ) );    internalForces.zero();    this->assembleVector( internalForces, tStep, InternalForceAssembler(), VM_Total,                         EModelDefaultEquationNumbering(), this->giveDomain(1) );    loadVector.subtract(internalForces);    */    this->updateSharedDofManagers(loadVector, EModelDefaultEquationNumbering(), ReactionExchangeTag);    //    // set-up numerical model    //    this->giveNumericalMethod( this->giveMetaStep( tStep->giveMetaStepNumber() ) );    //    // call numerical model to solve arose problem    //#ifdef VERBOSE    OOFEM_LOG_INFO("/n/nSolving .../n/n");#endif    NM_Status s = nMethod->solve(*stiffnessMatrix, loadVector, displacementVector);    if ( !( s & NM_Success ) ) {        OOFEM_ERROR("No success in solving system.");    }    tStep->incrementStateCounter();            // update solution state counter}

void StokesFlow :: solveYourselfAt(TimeStep *tStep){    Domain *d = this->giveDomain(1);    FloatArray externalForces;    FloatArray incrementOfSolution;    if ( d->giveNumberOfElements() == 0 && d->giveTopology() ) {        d->giveTopology()->replaceFEMesh();        this->meshqualityee.reset( new MeshQualityErrorEstimator( 1, d ) );    }    if ( d->giveTopology() && this->meshqualityee ) {        // Check the quality of the deformed mesh.        double meshdeformation = this->meshqualityee->giveValue(globalErrorEEV, tStep);        if ( this->giveProblemScale() == macroScale ) {            OOFEM_LOG_INFO("StokesFlow :: solveYourselfAt - Mesh deformation at %e/n", meshdeformation);        }        if ( this->ts == TS_NeedsRemeshing || meshdeformation > this->maxdef ) {            d->giveTopology()->replaceFEMesh();            OOFEM_LOG_INFO( "StokesFlow :: solveYourselfAt - New mesh created (%d elements)./n", d->giveNumberOfElements() );            /*meshdeformation =*/ this->meshqualityee->giveValue(globalErrorEEV, tStep);            this->giveExportModuleManager()->initialize();        }    }    int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );    // Move solution space to current time step    velocityPressureField->advanceSolution(tStep);    // Point pointer SolutionVector to current solution in velocityPressureField    velocityPressureField->initialize(VM_Total, tStep, solutionVector, EModelDefaultEquationNumbering() );    // Create "stiffness matrix"    if ( !this->stiffnessMatrix ) {        this->stiffnessMatrix.reset( classFactory.createSparseMtrx(sparseMtrxType) );        if ( !this->stiffnessMatrix ) {            OOFEM_ERROR("Couldn't create requested sparse matrix of type %d", sparseMtrxType);        }        this->stiffnessMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );    }    incrementOfSolution.resize(neq);    this->internalForces.resize(neq);    // Build initial/external load (LoadVector)    externalForces.resize(neq);    externalForces.zero();    this->assembleVector( externalForces, tStep, ExternalForceAssembler(), VM_Total,                         EModelDefaultEquationNumbering(), d );    this->updateSharedDofManagers(externalForces, EModelDefaultEquationNumbering(), LoadExchangeTag);    if ( this->giveProblemScale() == macroScale ) {        OOFEM_LOG_INFO("StokesFlow :: solveYourselfAt - Solving step %d, metastep %d, (neq = %d)/n", tStep->giveNumber(), tStep->giveMetaStepNumber(), neq);    }    this->giveNumericalMethod( this->giveCurrentMetaStep() );    this->initMetaStepAttributes( this->giveCurrentMetaStep() );#if 1    double loadLevel;    int currentIterations;    this->updateComponent( tStep, InternalRhs, d );    NM_Status status = this->nMethod->solve(*this->stiffnessMatrix,                                            externalForces,                                            NULL,                                            solutionVector,                                            incrementOfSolution,                                            this->internalForces,                                            this->eNorm,                                            loadLevel, // Only relevant for incrementalBCLoadVector?                                            SparseNonLinearSystemNM :: rlm_total,                                            currentIterations,                                            tStep);#else    this->updateComponent( tStep, InternalRhs, d );    this->updateComponent( tStep, NonLinearLhs, d );    this->internalForces.negated();    this->internalForces.add(externalForces);    NM_Status status = this->nMethod->giveLinearSolver()->solve(this->stiffnessMatrix.get(), & ( this->internalForces ), & solutionVector);    this->updateComponent( tStep, NonLinearLhs, d );#endif    if ( !( status & NM_Success ) ) {        OOFEM_ERROR("No success in solving problem at time step", tStep->giveNumber());    }}

void DEIDynamic :: solveYourselfAt(TimeStep *tStep){    //    // creates system of governing eq's and solves them at given time step    //    // this is an explicit problem: we assemble governing equating at time t    // and solution is obtained for time t+dt    //    // first assemble problem at current time step to obtain results in following    // time step.    // and then print results for this step also.    // for first time step we need special start code    Domain *domain = this->giveDomain(1);    int nelem = domain->giveNumberOfElements();    int nman = domain->giveNumberOfDofManagers();    IntArray loc;    Element *element;    DofManager *node;    Dof *iDof;    int nDofs, neq;    int i, k, n, j, jj, kk, init = 0;    double coeff, maxDt, maxOmi, maxOm = 0., maxOmEl, c1, c2, c3;    FloatMatrix charMtrx, charMtrx2;    FloatArray previousDisplacementVector;    neq = this->giveNumberOfEquations(EID_MomentumBalance);    if ( tStep->giveNumber() == giveNumberOfFirstStep() ) {        init = 1;#ifdef VERBOSE        OOFEM_LOG_INFO("Assembling mass matrix/n");#endif        //        // first step  assemble mass Matrix        //        massMatrix.resize(neq);        massMatrix.zero();        EModelDefaultEquationNumbering dn;        for ( i = 1; i <= nelem; i++ ) {            element = domain->giveElement(i);            element->giveLocationArray(loc, EID_MomentumBalance, dn);            element->giveCharacteristicMatrix(charMtrx,  LumpedMassMatrix, tStep);            // charMtrx.beLumpedOf(fullCharMtrx);            element->giveCharacteristicMatrix(charMtrx2, StiffnessMatrix, tStep);            //            // assemble it manually            //#ifdef DEBUG            if ( ( n = loc.giveSize() ) != charMtrx.giveNumberOfRows() ) {                _error("solveYourselfAt : dimension mismatch");            }#endif            n = loc.giveSize();            maxOmEl = 0.;            for ( j = 1; j <= n; j++ ) {                if ( charMtrx.at(j, j) > ZERO_MASS ) {                    maxOmi =  charMtrx2.at(j, j) / charMtrx.at(j, j);                    if ( init ) {                        maxOmEl = ( maxOmEl > maxOmi ) ? ( maxOmEl ) : ( maxOmi );                    }                }            }            maxOm = ( maxOm > maxOmEl ) ? ( maxOm ) : ( maxOmEl );            for ( j = 1; j <= n; j++ ) {                jj = loc.at(j);                if ( ( jj ) && ( charMtrx.at(j, j) <= ZERO_MASS ) ) {                    charMtrx.at(j, j) = charMtrx2.at(j, j) / maxOmEl;                }            }            for ( j = 1; j <= n; j++ ) {                jj = loc.at(j);                if ( jj ) {                    massMatrix.at(jj) += charMtrx.at(j, j);                }            }        }        // if init - try to determine the best deltaT        if ( init ) {            maxDt = 2 / sqrt(maxOm);            if ( deltaT > maxDt ) {                OOFEM_LOG_RELEVANT("DEIDynamic: deltaT reduced to %e/n", maxDt);                deltaT = maxDt;                tStep->setTimeIncrement(deltaT);            }        }        //        // special init step - compute displacements at tstep 0        ////.........这里部分代码省略.........

void NLTransientTransportProblem :: solveYourselfAt(TimeStep *tStep){    // creates system of governing eq's and solves them at given time step    // first assemble problem at current time step    // Right hand side    FloatArray rhs;    double solutionErr, incrementErr;    int neq =  this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );#ifdef VERBOSE    OOFEM_LOG_RELEVANT( "Solving [step number %8d, time %15e]/n", tStep->giveNumber(), tStep->giveTargetTime() );#endif    //Delete lhs matrix and create a new one. This is necessary due to growing/decreasing number of equations.    if ( tStep->isTheFirstStep() || this->changingProblemSize ) {        if ( conductivityMatrix ) {            delete conductivityMatrix;        }        conductivityMatrix = classFactory.createSparseMtrx(sparseMtrxType);        if ( conductivityMatrix == NULL ) {            OOFEM_ERROR("sparse matrix creation failed");        }        conductivityMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );#ifdef VERBOSE        OOFEM_LOG_INFO("Assembling conductivity and capacity matrices/n");#endif    }    //create previous solution from IC or from previous tStep    if ( tStep->isTheFirstStep() ) {        if ( !stepWhenIcApply ) {            stepWhenIcApply.reset( new TimeStep( *tStep->givePreviousStep() ) );        }        this->applyIC(stepWhenIcApply.get()); //insert solution to hash=1(previous), if changes in equation numbering    }    double dTTau = tStep->giveTimeIncrement();    double Tau = tStep->giveTargetTime() - ( 1. - alpha ) * tStep->giveTimeIncrement();    //Time step in which material laws are taken into account    TimeStep TauStep(tStep->giveNumber(), this, tStep->giveMetaStepNumber(), Tau, dTTau, tStep->giveSolutionStateCounter() + 1);    //Predictor    FloatArray *solutionVector;    UnknownsField->advanceSolution(tStep);    solutionVector = UnknownsField->giveSolutionVector(tStep);    //Initialize and give solutionVector from previous solution    if ( changingProblemSize ) {        if ( !tStep->isTheFirstStep() ) {            //copy recent solution to previous position, copy from hash=0 to hash=1(previous)            copyUnknownsInDictionary( VM_Total, tStep, tStep->givePreviousStep() );        }        UnknownsField->initialize( VM_Total, tStep->givePreviousStep(), *solutionVector, EModelDefaultEquationNumbering() );    } else {        //copy previous solution vector to actual        *solutionVector = *UnknownsField->giveSolutionVector( tStep->givePreviousStep() );    }    this->updateInternalState(& TauStep); //insert to hash=0(current), if changes in equation numbering    FloatArray solutionVectorIncrement(neq);    int nite = 0;    OOFEM_LOG_INFO("Time            Iter       ResidNorm       IncrNorm/n__________________________________________________________/n");    do {        nite++;        // Corrector#ifdef VERBOSE        // printf("/nAssembling conductivity and capacity matrices");#endif        if ( ( nite == 1 ) || ( NR_Mode == nrsolverFullNRM ) || ( ( NR_Mode == nrsolverAccelNRM ) && ( nite % MANRMSteps == 0 ) ) ) {            conductivityMatrix->zero();            //Assembling left hand side - start with conductivity matrix            this->assemble( *conductivityMatrix, & TauStep, IntSourceLHSMatrix,                           EModelDefaultEquationNumbering(), this->giveDomain(1) );            conductivityMatrix->times(alpha);            //Add capacity matrix            this->assemble( *conductivityMatrix, & TauStep, NSTP_MidpointLhs,                           EModelDefaultEquationNumbering(), this->giveDomain(1) );        }        rhs.resize(neq);        rhs.zero();        //edge or surface load on element        this->assembleVectorFromElements( rhs, & TauStep, ElementBCTransportVector, VM_Total,                                         EModelDefaultEquationNumbering(), this->giveDomain(1) );        //add internal source vector on elements        this->assembleVectorFromElements( rhs, & TauStep, ElementInternalSourceVector, VM_Total,                                         EModelDefaultEquationNumbering(), this->giveDomain(1) );        //add nodal load        this->assembleVectorFromDofManagers( rhs, & TauStep, ExternalForcesVector, VM_Total,                                            EModelDefaultEquationNumbering(), this->giveDomain(1) );//.........这里部分代码省略.........

voidSUPG :: solveYourselfAt(TimeStep *tStep){    int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );    FloatArray *solutionVector = NULL, *prevSolutionVector = NULL;    FloatArray externalForces(neq);    this->internalForces.resize(neq);    if ( tStep->isTheFirstStep() ) {        TimeStep *stepWhenIcApply = tStep->givePreviousStep();        if ( materialInterface ) {            materialInterface->initialize();        }        this->applyIC(stepWhenIcApply);        //if (this->fsflag) this->updateDofManActivityMap(tStep);    }    if ( !requiresUnknownsDictionaryUpdate() ) {        VelocityPressureField->advanceSolution(tStep);        solutionVector = VelocityPressureField->giveSolutionVector(tStep);        prevSolutionVector = VelocityPressureField->giveSolutionVector( tStep->givePreviousStep() );    }    if ( initFlag ) {        if ( !requiresUnknownsDictionaryUpdate() ) {            //previousAccelerationVector.resize(neq);            accelerationVector.resize(neq);            solutionVector->resize(neq);            prevSolutionVector->resize(neq);        }        incrementalSolutionVector.resize(neq);        lhs.reset( classFactory.createSparseMtrx(sparseMtrxType) );        if ( !lhs ) {            OOFEM_ERROR("sparse matrix creation failed");        }        lhs->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );        if ( materialInterface ) {            this->updateElementsForNewInterfacePosition(tStep);        }        initFlag = 0;    } else if ( requiresUnknownsDictionaryUpdate() ) {        // rebuild lhs structure and resize solution vector        incrementalSolutionVector.resize(neq);        lhs->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );    }    if ( !requiresUnknownsDictionaryUpdate() ) {        *solutionVector = *prevSolutionVector;    }    //previousAccelerationVector=accelerationVector;    // evaluate element supg and sppg stabilization coeffs    this->evaluateElementStabilizationCoeffs(tStep);    //    // predictor    //    if ( requiresUnknownsDictionaryUpdate() ) {        this->updateDofUnknownsDictionary_predictor(tStep);    } else {        this->updateSolutionVectors_predictor(* solutionVector, accelerationVector, tStep);    }    if ( tStep->giveNumber() != 1 ) {        if ( materialInterface ) {            //if (this->fsflag) updateDofManVals(tStep);#ifdef SUPG_IMPLICIT_INTERFACE #ifdef TIME_REPORT            Timer timer;            timer.startTimer(); #endif            materialInterface->updatePosition( this->giveCurrentStep() );            updateElementsForNewInterfacePosition(tStep); #ifdef TIME_REPORT            timer.stopTimer();            OOFEM_LOG_INFO("SUPG info: user time consumed by updating interfaces: %.2fs/n", timer.getUtime();                           ); #endif#else            //updateElementsForNewInterfacePosition (tStep);#endif            //if (this->fsflag) this->updateDofManActivityMap(tStep);        }    }