C++ DofManager::giveNumberOfDofs方法代码示例

int
StructuralEngngModel :: packDofManagers(FloatArray *src, ProcessCommunicator &processComm, bool prescribedEquations)
{
    int result = 1;
    int i, size;
    int j, ndofs, eqNum;
    Domain *domain = this->giveDomain(1);
    IntArray const *toSendMap = processComm.giveToSendMap();
    ProcessCommunicatorBuff *pcbuff = processComm.giveProcessCommunicatorBuff();
    DofManager *dman;
    Dof *jdof;

    size = toSendMap->giveSize();
    for ( i = 1; i <= size; i++ ) {
        dman = domain->giveDofManager( toSendMap->at(i) );
        ndofs = dman->giveNumberOfDofs();
        for ( j = 1; j <= ndofs; j++ ) {
            jdof = dman->giveDof(j);
            if ( prescribedEquations ) {
                eqNum = jdof->__givePrescribedEquationNumber();
            } else {
                eqNum = jdof->__giveEquationNumber();
            }
            if ( jdof->isPrimaryDof() && eqNum ) {
                    result &= pcbuff->packDouble( src->at(eqNum) );
            }
        }
    }

    return result;
}

void
LumpedMassElement :: giveDofManDofIDMask(int inode, EquationID eid, IntArray &answer) const
{
    answer.resize(0, 6);
    DofManager *dman = this->giveDofManager(inode);
    int _i, _ndof = dman->giveNumberOfDofs();
    DofIDItem _dofid;

    // simply collect all "structural" dofs of element node
    for ( _i = 1; _i <= _ndof; _i++ ) {
        _dofid = dman->giveDof(_i)->giveDofID();
        if ( ( _dofid == D_u ) || ( _dofid == D_v ) || ( _dofid == D_w ) ||
            ( _dofid == R_u ) || ( _dofid == R_v ) || ( _dofid == R_w ) ) {
            answer.followedBy(_dofid);
        }
    }
}

void
PetscNatural2LocalOrdering :: init(EngngModel *emodel, EquationID ut, int di,  EquationType et)
{
    Domain *d = emodel->giveDomain(di);
    int i, j, n_eq = 0, ndofs, ndofman = d->giveNumberOfDofManagers(), loc_eq = 1;
    bool lFlag;
    DofManager *dman;
    EModelDefaultEquationNumbering dn;
    EModelDefaultPrescribedEquationNumbering dpn;


    // determine number of local eqs + number of those shared DOFs which are numbered by receiver
    // shared dofman is numbered on partition with lovest rank number

    if ( et == et_standard ) {
        n2l.resize( emodel->giveNumberOfEquations(ut) );
    } else {
        n2l.resize( emodel->giveNumberOfPrescribedEquations(ut) );
    }

    for ( i = 1; i <= ndofman; i++ ) {
        dman = d->giveDofManager(i);
        lFlag = isLocal(dman);
        ndofs = dman->giveNumberOfDofs();
        for ( j = 1; j <= ndofs; j++ ) {
            if ( dman->giveDof(j)->isPrimaryDof() ) {
                if ( et == et_standard ) {
                    n_eq = dman->giveDof(j)->giveEquationNumber(dn);
                } else {
                    n_eq = dman->giveDof(j)->giveEquationNumber(dpn);
                }

                if ( n_eq == 0 ) {
                    continue;
                }

                if ( lFlag ) {
                    n2l.at(n_eq) = loc_eq++;
                } else {
                    n2l.at(n_eq) = 0;
                }
            }
        }
    }
}

int
LumpedMassElement :: computeNumberOfDofs(EquationID ut)
{
    DofManager *dman = this->giveDofManager(1);
    int _i, _ndof = dman->giveNumberOfDofs();
    int answer = 0;
    DofIDItem _dofid;

    // simply count all "structural" dofs of element node
    for ( _i = 1; _i <= _ndof; _i++ ) {
        _dofid = dman->giveDof(_i)->giveDofID();
        if ( ( _dofid == D_u ) || ( _dofid == D_v ) || ( _dofid == D_w ) ||
            ( _dofid == R_u ) || ( _dofid == R_v ) || ( _dofid == R_w ) ) {
            answer++;
        }
    }

    return answer;
}

void
StructuralInterfaceElementPhF :: computeLocationArrayOfDofIDs( const IntArray &dofIdArray, IntArray &answer )
{
    // Routine to compute the local ordering array an element given a dofid array.
    answer.resize( 0 );
    int k = 0;
    for(int i = 1; i <= this->giveNumberOfDofManagers(); i++) {
        DofManager *dMan = this->giveDofManager( i );
        for(int j = 1; j <= dofIdArray.giveSize( ); j++) {

            if(dMan->hasDofID( (DofIDItem) dofIdArray.at( j ) )) {
                // hack
                
                answer.followedBy( k + j );
            }
        }
        k += dMan->giveNumberOfDofs( );
    }
}

void
PhaseFieldElement :: computeLocationArrayOfDofIDs( const IntArray &dofIdArray, IntArray &answer )
{
    // Routine to extract compute the location array an element given an dofid array.
    answer.clear();
    NLStructuralElement *el = this->giveElement();
    int k = 0;
    for ( int i = 1; i <= el->giveNumberOfDofManagers(); i++ ) {
        DofManager *dMan = el->giveDofManager( i );
        for ( int j = 1; j <= dofIdArray.giveSize( ); j++ ) {

            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at( j ) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at( j ) );
                answer.followedBy( k + d->giveNumber( ) );
            }
        }
        k += dMan->giveNumberOfDofs( );
    }
}

int
StructuralEngngModel :: unpackDofManagers(FloatArray *dest, ProcessCommunicator &processComm, bool prescribedEquations)
{
    int result = 1;
    int i, size;
    int j, ndofs, eqNum;
    Domain *domain = this->giveDomain(1);
    dofManagerParallelMode dofmanmode;
    IntArray const *toRecvMap = processComm.giveToRecvMap();
    ProcessCommunicatorBuff *pcbuff = processComm.giveProcessCommunicatorBuff();
    DofManager *dman;
    Dof *jdof;
    double value;


    size = toRecvMap->giveSize();
    for ( i = 1; i <= size; i++ ) {
        dman = domain->giveDofManager( toRecvMap->at(i) );
        ndofs = dman->giveNumberOfDofs();
        dofmanmode = dman->giveParallelMode();
        for ( j = 1; j <= ndofs; j++ ) {
            jdof = dman->giveDof(j);
            if ( prescribedEquations ) {
                eqNum = jdof->__givePrescribedEquationNumber();
            } else {
                eqNum = jdof->__giveEquationNumber();
            }
            if ( jdof->isPrimaryDof() && eqNum ) {
                result &= pcbuff->unpackDouble(value);
                if ( dofmanmode == DofManager_shared ) {
                    dest->at(eqNum) += value;
                } else if ( dofmanmode == DofManager_remote ) {
                    dest->at(eqNum)  = value;
                } else {
                    _error("unpackReactions: unknown dof namager parallel mode");
                }
            }
        }
    }

    return result;
}

void
CoupledFieldsElement :: computeLocationArrayOfDofIDs(const IntArray &dofIdArray, IntArray &answer)
{
    // Routine to extract compute the location array an element given an dofid array.
    answer.resize(0);
    
    int k = 0;
    for ( int i = 1; i <= numberOfDofMans; i++ ) {
        DofManager *dMan = this->giveDofManager(i);        
        for (int j = 1; j <= dofIdArray.giveSize(); j++ ) {   
            
            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at(j) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at(j) );
                answer.followedBy( k + d->giveNumber() );
                //answer.followedBy( k + j );
            }
        }
        k += dMan->giveNumberOfDofs( );
    }
}

void
StructuralEngngModel :: buildReactionTable(IntArray &restrDofMans, IntArray &restrDofs,
                                           IntArray &eqn, TimeStep *tStep, int di)
{
    // determine number of restrained dofs
    Domain *domain = this->giveDomain(di);
    int numRestrDofs = this->giveNumberOfPrescribedDomainEquations(di, EID_MomentumBalance);
    int ndofMan = domain->giveNumberOfDofManagers();
    int i, j, indofs, rindex, count = 0;
    DofManager *inode;
    Dof *jdof;

    // initialize corresponding dofManagers and dofs for each restrained dof
    restrDofMans.resize(numRestrDofs);
    restrDofs.resize(numRestrDofs);
    eqn.resize(numRestrDofs);

    for ( i = 1; i <= ndofMan; i++ ) {
        inode = domain->giveDofManager(i);
        indofs = inode->giveNumberOfDofs();
        for ( j = 1; j <= indofs; j++ ) {
            jdof = inode->giveDof(j);
            if ( ( jdof->giveClassID() != SimpleSlaveDofClass ) && ( jdof->hasBc(tStep) ) ) { // skip slave dofs
                rindex = jdof->__givePrescribedEquationNumber();
                if ( rindex ) {
                    count++;
                    restrDofMans.at(count) = i;
                    restrDofs.at(count) = j;
                    eqn.at(count) = rindex;
                } else {
                    // NullDof has no equation number and no prescribed equation number
                    //_error ("No prescribed equation number assigned to supported DOF");
                }
            }
        }
    }
}

void
PetscNatural2GlobalOrdering :: init(EngngModel *emodel, EquationID ut, int di, EquationType et)
{
    Domain *d = emodel->giveDomain(di);
    int i, j, k, p, ndofs, ndofman = d->giveNumberOfDofManagers();
    int myrank = emodel->giveRank();
    DofManager *dman;
    // determine number of local eqs + number of those shared DOFs which are numbered by receiver
    // shared dofman is numbered on partition with lovest rank number
    EModelDefaultEquationNumbering dn;
    EModelDefaultPrescribedEquationNumbering dpn;

#ifdef __VERBOSE_PARALLEL
    VERBOSEPARALLEL_PRINT("PetscNatural2GlobalOrdering :: init", "initializing N2G ordering", myrank);
#endif

    l_neqs = 0;
    for ( i = 1; i <= ndofman; i++ ) {
        dman = d->giveDofManager(i);
        /*
         *  if (dman->giveParallelMode() == DofManager_local) { // count all dofman eqs
         *    ndofs = dman->giveNumberOfDofs ();
         *    for (j=1; j<=ndofs; j++) {
         *      if (dman->giveDof(j)->isPrimaryDof()) {
         *        if (dman->giveDof(j)->giveEquationNumber()) l_neqs++;
         *      }
         *    }
         *  } else if (dman->giveParallelMode() == DofManager_shared) {
         *    // determine if problem is the lowest one sharing the dofman; if yes the receiver is responsible to
         *    // deliver number
         *    IntArray *plist = dman->givePartitionList();
         *    int n = plist->giveSize();
         *    int minrank = myrank;
         *    for (j=1; j<=n; j++) minrank = min (minrank, plist->at(j));
         *    if (minrank == myrank) { // count eqs
         *      ndofs = dman->giveNumberOfDofs ();
         *      for (j=1; j<=ndofs; j++) {
         *        if (dman->giveDof(j)->isPrimaryDof()) {
         *          if (dman->giveDof(j)->giveEquationNumber()) l_neqs++;
         *        }
         *      }
         *    }
         *  } // end shared dman
         */
        if ( isLocal(dman) ) {
            ndofs = dman->giveNumberOfDofs();
            for ( j = 1; j <= ndofs; j++ ) {
                if ( dman->giveDof(j)->isPrimaryDof() ) {
                    if ( et == et_standard ) {
                        if ( dman->giveDof(j)->giveEquationNumber(dn) ) {
                            l_neqs++;
                        }
                    } else {
                        if ( dman->giveDof(j)->giveEquationNumber(dpn) ) {
                            l_neqs++;
                        }
                    }
                }
            }
        }
    }

    // exchange with other procs the number of eqs numbered on particular procs
    int *leqs = new int [ emodel->giveNumberOfProcesses() ];
    MPI_Allgather(& l_neqs, 1, MPI_INT, leqs, 1, MPI_INT, MPI_COMM_WORLD);
    // compute local offset
    int offset = 0;
    for ( j = 0; j < myrank; j++ ) {
        offset += leqs [ j ];
    }

    // count global number of eqs
    for ( g_neqs = 0, j = 0; j < emodel->giveNumberOfProcesses(); j++ ) {
        g_neqs += leqs [ j ];
    }

    // send numbered shared ones
    if ( et == et_standard ) {
        locGlobMap.resize( emodel->giveNumberOfEquations(ut) );
    } else {
        locGlobMap.resize( emodel->giveNumberOfPrescribedEquations(ut) );
    }

    // determine shared dofs
    int psize, nproc = emodel->giveNumberOfProcesses();
    IntArray sizeToSend(nproc), sizeToRecv(nproc), nrecToReceive(nproc);
#ifdef __VERBOSE_PARALLEL
    IntArray nrecToSend(nproc);
#endif
    const IntArray *plist;
    for ( i = 1; i <= ndofman; i++ ) {
        // if (domain->giveDofManager(i)->giveParallelMode() == DofManager_shared) {
        if ( isShared( d->giveDofManager(i) ) ) {
            int n = d->giveDofManager(i)->giveNumberOfDofs();
            plist = d->giveDofManager(i)->givePartitionList();
            psize = plist->giveSize();
            int minrank = myrank;
            for ( j = 1; j <= psize; j++ ) {
                minrank = min( minrank, plist->at(j) );
            }
//.........这里部分代码省略.........

bool
NRSolver :: checkConvergence(FloatArray &RT, FloatArray &F, FloatArray &rhs,  FloatArray &ddX, FloatArray &X,
                             double RRT, const FloatArray &internalForcesEBENorm,
                             int nite, bool &errorOutOfRange, TimeStep *tNow)
{
    double forceErr, dispErr;
    FloatArray dg_forceErr, dg_dispErr, dg_totalLoadLevel, dg_totalDisp;
    bool answer;
    EModelDefaultEquationNumbering dn;
 #ifdef __PARALLEL_MODE
  #ifdef __PETSC_MODULE
    PetscContext *parallel_context = engngModel->givePetscContext(this->domain->giveNumber());
    Natural2LocalOrdering *n2l = parallel_context->giveN2Lmap();
  #endif
 #endif

    /*
     * The force errors are (if possible) evaluated as relative errors.
     * If the norm of applied load vector is zero (one may load by temperature, etc)
     * then the norm of reaction forces is used in relative norm evaluation.
     *
     * Note: This is done only when all dofs are included (nccdg = 0). Not implemented if
     * multiple convergence criteria are used.
     *
     */

    answer = true;
    errorOutOfRange = false;

    if ( internalForcesEBENorm.giveSize() > 1 ) { // Special treatment when just one norm is given; No grouping
        int nccdg = this->domain->giveMaxDofID();
        // Keeps tracks of which dof IDs are actually in use;
        IntArray idsInUse(nccdg);
        idsInUse.zero();
        // zero error norms per group
        dg_forceErr.resize(nccdg); dg_forceErr.zero();
        dg_dispErr.resize(nccdg); dg_dispErr.zero();
        dg_totalLoadLevel.resize(nccdg); dg_totalLoadLevel.zero();
        dg_totalDisp.resize(nccdg); dg_totalDisp.zero();
        // loop over dof managers
        int ndofman = domain->giveNumberOfDofManagers();
        for ( int idofman = 1; idofman <= ndofman; idofman++ ) {
            DofManager *dofman = domain->giveDofManager(idofman);
 #if ( defined ( __PARALLEL_MODE ) && defined ( __PETSC_MODULE ) )
            if ( !parallel_context->isLocal(dofman) ) {
                continue;
            }

 #endif

            // loop over individual dofs
            int ndof = dofman->giveNumberOfDofs();
            for ( int idof = 1; idof <= ndof; idof++ ) {
                Dof *dof = dofman->giveDof(idof);
                if ( !dof->isPrimaryDof() ) continue;
                int eq = dof->giveEquationNumber(dn);
                int dofid = dof->giveDofID();
                if ( !eq ) continue;
 
                dg_forceErr.at(dofid) += rhs.at(eq) * rhs.at(eq);
                dg_dispErr.at(dofid) += ddX.at(eq) * ddX.at(eq);
                dg_totalLoadLevel.at(dofid) += RT.at(eq) * RT.at(eq);
                dg_totalDisp.at(dofid) += X.at(eq) * X.at(eq);
                idsInUse.at(dofid) = 1;
            } // end loop over DOFs
        } // end loop over dof managers

        // loop over elements and their DOFs
        int nelem = domain->giveNumberOfElements();
        for ( int ielem = 1; ielem <= nelem; ielem++ ) {
            Element *elem = domain->giveElement(ielem);
 #ifdef __PARALLEL_MODE
            if ( elem->giveParallelMode() != Element_local ) {
                continue;
            }

 #endif
            // loop over element internal Dofs
            for ( int idofman = 1; idofman <= elem->giveNumberOfInternalDofManagers(); idofman++) {
                DofManager *dofman = elem->giveInternalDofManager(idofman);
                int ndof = dofman->giveNumberOfDofs();
                // loop over individual dofs
                for ( int idof = 1; idof <= ndof; idof++ ) {
                    Dof *dof = dofman->giveDof(idof);
                    if ( !dof->isPrimaryDof() ) continue;
                    int eq = dof->giveEquationNumber(dn);
                    int dofid = dof->giveDofID();
                    
                    if ( !eq ) continue;
 #if ( defined ( __PARALLEL_MODE ) && defined ( __PETSC_MODULE ) )
                    if ( engngModel->isParallel() && !n2l->giveNewEq(eq) ) continue;
 #endif
                    dg_forceErr.at(dofid) += rhs.at(eq) * rhs.at(eq);
                    dg_dispErr.at(dofid) += ddX.at(eq) * ddX.at(eq);
                    dg_totalLoadLevel.at(dofid) += RT.at(eq) * RT.at(eq);
                    dg_totalDisp.at(dofid) += X.at(eq) * X.at(eq);
                    idsInUse.at(dofid) = 1;
                } // end loop over DOFs
            } // end loop over element internal dofmans
        } // end loop over elements
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                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
        //
        displacementVector.resize(neq);
        displacementVector.zero();
        nextDisplacementVector.resize(neq);
        nextDisplacementVector.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();


        for ( j = 1; j <= nman; j++ ) {
            node = domain->giveDofManager(j);
            nDofs = node->giveNumberOfDofs();

            for ( k = 1; k <= nDofs; k++ ) {
                // ask for initial values obtained from
                // bc (boundary conditions) and ic (initial conditions)
                // now we are setting initial cond. for step -1.
                iDof  =  node->giveDof(k);
                if ( !iDof->isPrimaryDof() ) {
                    continue;
                }

                jj = iDof->__giveEquationNumber();
                if ( jj ) {
                    nextDisplacementVector.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Total, tStep);
                    // become displacementVector after init
                    velocityVector.at(jj)     = iDof->giveUnknown(EID_MomentumBalance, VM_Velocity, tStep);
                    // accelerationVector = iDof->giveUnknown(AccelerartionVector,tStep) ;
                }
            }
        }

        for ( j = 1; j <= neq; j++ ) {
            nextDisplacementVector.at(j) -= velocityVector.at(j) * ( deltaT );
        }

        return;
    } // end of init step

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling right hand side\n");
#endif

void NlDEIDynamic :: solveYourselfAt(TimeStep *tStep)
{
    //
    // Creates system of governing eq's and solves them at given time step.
    //

    Domain *domain = this->giveDomain(1);
    int neq = this->giveNumberOfEquations(EID_MomentumBalance);
    int nman  = domain->giveNumberOfDofManagers();

    DofManager *node;
    Dof *iDof;

    int nDofs;
    int i, k, j, jj;
    double coeff, maxDt, maxOm = 0.;
    double prevIncrOfDisplacement, incrOfDisplacement;

    if ( initFlag ) {
#ifdef VERBOSE
        OOFEM_LOG_DEBUG("Assembling mass matrix\n");
#endif

        //
        // Assemble mass matrix.
        //
        this->computeMassMtrx(massMatrix, maxOm, tStep);

        if ( drFlag ) {
            // If dynamic relaxation: Assemble amplitude load vector.
            loadRefVector.resize(neq);
            loadRefVector.zero();

            this->computeLoadVector(loadRefVector, VM_Total, tStep);

#ifdef __PARALLEL_MODE
            // Compute the processor part of load vector norm pMp
            this->pMp = 0.0;
            double my_pMp = 0.0, coeff = 1.0;
            int eqNum, ndofs, ndofman = domain->giveNumberOfDofManagers();
            dofManagerParallelMode dofmanmode;
            DofManager *dman;
            Dof *jdof;
            for ( int dm = 1; dm <= ndofman; dm++ ) {
                dman = domain->giveDofManager(dm);
                ndofs = dman->giveNumberOfDofs();
                dofmanmode = dman->giveParallelMode();

                // Skip all remote and null dofmanagers
                coeff = 1.0;
                if ( ( dofmanmode == DofManager_remote ) || ( ( dofmanmode == DofManager_null ) ) ) {
                    continue;
                } else if ( dofmanmode == DofManager_shared ) {
                    coeff = 1. / dman->givePartitionsConnectivitySize();
                }

                // For shared nodes we add locally an average = 1/givePartitionsConnectivitySize()*contribution,
                for ( j = 1; j <= ndofs; j++ ) {
                    jdof = dman->giveDof(j);
                    if ( jdof->isPrimaryDof() && ( eqNum = jdof->__giveEquationNumber() ) ) {
                        my_pMp += coeff * loadRefVector.at(eqNum) * loadRefVector.at(eqNum) / massMatrix.at(eqNum);
                    }
                }
            }

            // Sum up the contributions from processors.
            MPI_Allreduce(& my_pMp, & pMp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#else
            this->pMp = 0.0;
            for ( i = 1; i <= neq; i++ ) {
                pMp += loadRefVector.at(i) * loadRefVector.at(i) / massMatrix.at(i);
            }
#endif
            // Solve for rate of loading process (parameter "c") (undamped system assumed),
            if ( dumpingCoef < 1.e-3 ) {
                c = 3.0 * this->pyEstimate / pMp / Tau / Tau;
            } else {
                c = this->pyEstimate * Tau * dumpingCoef * dumpingCoef * dumpingCoef / pMp /
                    ( -3.0 / 2.0 + dumpingCoef * Tau + 2.0 * exp(-dumpingCoef * Tau) - 0.5 * exp(-2.0 * dumpingCoef * Tau) );
            }
        }

        initFlag = 0;
    }


    if ( tStep->giveNumber() == giveNumberOfFirstStep() ) {
        //
        // Special init step - Compute displacements at tstep 0.
        //
        displacementVector.resize(neq);
        displacementVector.zero();
        previousIncrementOfDisplacementVector.resize(neq);
        previousIncrementOfDisplacementVector.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();

        for ( j = 1; j <= nman; j++ ) {
//.........这里部分代码省略.........

void
NonLinearDynamic :: proceedStep(int di, TimeStep *tStep)
{
    // creates system of governing eq's and solves them at given time step
    // first assemble problem at current time step

    int neq = this->giveNumberOfEquations(EID_MomentumBalance);

    // Time-stepping constants
    double dt2 = deltaT * deltaT;

    if ( tStep->giveTimeDiscretization() == TD_Newmark ) {
        OOFEM_LOG_DEBUG("Solving using Newmark-beta method\n");
        a0 = 1 / ( beta * dt2 );
        a1 = gamma / ( beta * deltaT );
        a2 = 1 / ( beta * deltaT );
        a3 = 1 / ( 2 *  beta ) - 1;
        a4 = ( gamma / beta ) - 1;
        a5 = deltaT / 2 * ( gamma / beta - 2 );
        a6 = 0;
    } else if ( ( tStep->giveTimeDiscretization() == TD_TwoPointBackward ) || ( tStep->giveNumber() == giveNumberOfFirstStep() ) ) {
        if ( tStep->giveTimeDiscretization() != TD_ThreePointBackward ) {
            OOFEM_LOG_DEBUG("Solving using Backward Euler method\n");
        } else {
            OOFEM_LOG_DEBUG("Solving initial step using Three-point Backward Euler method\n");
        }
        a0 = 1 / dt2;
        a1 = 1 / deltaT;
        a2 = 1 / deltaT;
        a3 = 0;
        a4 = 0;
        a5 = 0;
        a6 = 0;
    } else if ( tStep->giveTimeDiscretization() == TD_ThreePointBackward ) {
        OOFEM_LOG_DEBUG("Solving using Three-point Backward Euler method\n");
        a0 = 2 / dt2;
        a1 = 3 / ( 2 * deltaT );
        a2 = 2 / deltaT;
        a3 = 0;
        a4 = 0;
        a5 = 0;
        a6 = 1 / ( 2 * deltaT );
    } else {
        _error("NonLinearDynamic: Time-stepping scheme not found!\n")
    }

    if ( tStep->giveNumber() == giveNumberOfFirstStep() ) {
        // Initialization
        previousIncrementOfDisplacement.resize(neq);
        previousIncrementOfDisplacement.zero();
        previousTotalDisplacement.resize(neq);
        previousTotalDisplacement.zero();
        totalDisplacement.resize(neq);
        totalDisplacement.zero();
        previousInternalForces.resize(neq);
        previousInternalForces.zero();
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();

        TimeStep *stepWhenIcApply = new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 0,
                                                 -deltaT, deltaT, 0);

        int nDofs, j, k, jj;
        int nman  = this->giveDomain(di)->giveNumberOfDofManagers();
        DofManager *node;
        Dof *iDof;

        // Considering initial conditions.
        for ( j = 1; j <= nman; j++ ) {
            node = this->giveDomain(di)->giveDofManager(j);
            nDofs = node->giveNumberOfDofs();

            for ( k = 1; k <= nDofs; k++ ) {
                // Ask for initial values obtained from
                // bc (boundary conditions) and ic (initial conditions).
                iDof  =  node->giveDof(k);
                if ( !iDof->isPrimaryDof() ) {
                    continue;
                }

                jj = iDof->__giveEquationNumber();
                if ( jj ) {
                    incrementOfDisplacement.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Total, stepWhenIcApply);
                    velocityVector.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Velocity, stepWhenIcApply);
                    accelerationVector.at(jj) = iDof->giveUnknown(EID_MomentumBalance, VM_Acceleration, stepWhenIcApply);
                }
            }
        }
    } else {
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
    }

    if ( initFlag ) {
        // First assemble problem at current time step.
        // Option to take into account initial conditions.
//.........这里部分代码省略.........

void IncrementalLinearStatic :: solveYourselfAt(TimeStep *tStep)
{
    // Creates system of governing eq's and solves them at given time step

    // Initiates the total displacement to zero.
    if ( tStep->isTheFirstStep() ) {
        Domain *d = this->giveDomain(1);
        for ( int i = 1; i <= d->giveNumberOfDofManagers(); i++ ) {
            DofManager *dofman = d->giveDofManager(i);
            for ( int j = 1; j <= dofman->giveNumberOfDofs(); j++ ) {
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total_Old, 0.);
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total, 0.);
                // This is actually redundant now;
                //dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Incremental, 0.);
            }
        }

	int nbc = d->giveNumberOfBoundaryConditions();
	for ( int ibc = 1; ibc <= nbc; ++ibc ) {
	  GeneralBoundaryCondition *bc = d->giveBc(ibc);
	  ActiveBoundaryCondition *abc;

	  if ( ( abc = dynamic_cast< ActiveBoundaryCondition * >( bc ) ) ) {
	    int ndman = abc->giveNumberOfInternalDofManagers();
	    for ( int i = 1; i <= ndman; i++ ) {
	      DofManager *dofman = abc->giveInternalDofManager(i);
	      for ( int j = 1; j <= dofman->giveNumberOfDofs(); j++ ) {
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total_Old, 0.);
                dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Total, 0.);
                // This is actually redundant now;
                //dofman->giveDof(j)->updateUnknownsDictionary(tStep, VM_Incremental, 0.);
	      }
	    }
	  }
	}
    }

    // Apply dirichlet b.c's on total values
    Domain *d = this->giveDomain(1);
    for ( int i = 1; i <= d->giveNumberOfDofManagers(); i++ ) {
        DofManager *dofman = d->giveDofManager(i);
        for ( int j = 1; j <= dofman->giveNumberOfDofs(); j++ ) {
            Dof *d = dofman->giveDof(j);
            double tot = d->giveUnknown(VM_Total_Old, tStep);
            if ( d->hasBc(tStep) ) {
                tot += d->giveBcValue(VM_Incremental, tStep);
            }

            d->updateUnknownsDictionary(tStep, VM_Total, tot);
        }
    }


#ifdef VERBOSE
    OOFEM_LOG_RELEVANT( "Solving [step number %8d, time %15e]\n", tStep->giveNumber(), tStep->giveTargetTime() );
#endif

    int neq = this->giveNumberOfDomainEquations(1, EModelDefaultEquationNumbering());

    if (neq == 0) { // Allows for fully prescribed/empty problems.
        return;
    }

    incrementOfDisplacementVector.resize(neq);
    incrementOfDisplacementVector.zero();

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling load\n");
#endif
    // Assembling the element part of load vector
    internalLoadVector.resize(neq);
    internalLoadVector.zero();
    this->assembleVector( internalLoadVector, tStep, EID_MomentumBalance, InternalForcesVector,
                          VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );

    loadVector.resize(neq);
    loadVector.zero();
    this->assembleVector( loadVector, tStep, EID_MomentumBalance, ExternalForcesVector,
                          VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );

    loadVector.subtract(internalLoadVector);

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling stiffness matrix\n");
#endif
    if ( stiffnessMatrix ) {
        delete stiffnessMatrix;
    }

    stiffnessMatrix = classFactory.createSparseMtrx(sparseMtrxType);
    if ( stiffnessMatrix == NULL ) {
        _error("solveYourselfAt: sparse matrix creation failed");
    }

    stiffnessMatrix->buildInternalStructure( this, 1, EID_MomentumBalance, EModelDefaultEquationNumbering() );
    stiffnessMatrix->zero();
    this->assemble( stiffnessMatrix, tStep, EID_MomentumBalance, StiffnessMatrix,
                   EModelDefaultEquationNumbering(), this->giveDomain(1) );

#ifdef VERBOSE
//.........这里部分代码省略.........