C++ Domain::giveBc方法代码示例

void
StructuralEngngModel :: updateInternalState(TimeStep *tStep)
{
    int nnodes;
    Domain *domain;

    for ( int idomain = 1; idomain <= this->giveNumberOfDomains(); idomain++ ) {
        domain = this->giveDomain(idomain);

        nnodes = domain->giveNumberOfDofManagers();
        if ( requiresUnknownsDictionaryUpdate() ) {
            for ( int j = 1; j <= nnodes; j++ ) {
                this->updateDofUnknownsDictionary(domain->giveDofManager(j), tStep);
            }
        }

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

            if ( ( abc = dynamic_cast< ActiveBoundaryCondition * >(bc) ) ) {
                int ndman = abc->giveNumberOfInternalDofManagers();
                for ( int j = 1; j <= ndman; j++ ) {
                    this->updateDofUnknownsDictionary(abc->giveInternalDofManager(j), tStep);
                }
            }
        }

        if ( internalVarUpdateStamp != tStep->giveSolutionStateCounter() ) {
            int nelem = domain->giveNumberOfElements();
            for ( int j = 1; j <= nelem; j++ ) {
                domain->giveElement(j)->updateInternalState(tStep);
            }

            internalVarUpdateStamp = tStep->giveSolutionStateCounter();
        }
    }
}

void GnuplotExportModule::outputGradient(int bc, Domain &d, FloatArray &grad, TimeStep *tStep)
{
    // Homogenized strain
    double timeFactor = d.giveBc(bc)->giveTimeFunction()->evaluateAtTime(tStep->giveTargetTime());
    printf("timeFactor: %e\n", timeFactor );
    grad.times(timeFactor);
    printf("Mean grad computed in Gnuplot export module: "); grad.printYourself();

    double time = tStep->giveTargetTime();


    std :: stringstream strMeanGrad;
    strMeanGrad << "PrescribedGradientGnuplotMeanGrad" << bc << "Time" << time << ".dat";
    std :: string nameMeanGrad = strMeanGrad.str();
    std::vector<double> componentArrayGrad, gradArray;

    for(int i = 1; i <= grad.giveSize(); i++) {
        componentArrayGrad.push_back(i);
        gradArray.push_back(grad.at(i));
    }

    XFEMDebugTools::WriteArrayToGnuplot(nameMeanGrad, componentArrayGrad, gradArray);
}

//.........这里部分代码省略.........
     */
    int neq = eModel->giveNumberOfDomainEquations(di, s);

#ifndef DynCompCol_USE_STL_SETS
    IntArray loc;
    Domain *domain = eModel->giveDomain(di);
    int nelem = domain->giveNumberOfElements();
    int i, ii, j, jj, n;
    Element *elem;

    nColumns = nRows = neq;

    if ( rowind_ ) {
        for ( i = 0; i < nColumns; i++ ) {
            delete this->rowind_ [ i ];
        }

        delete this->rowind_;
    }

    rowind_ = ( IntArray ** ) new IntArray * [ neq ];
    for ( j = 0; j < neq; j++ ) {
        rowind_ [ j ] = new IntArray();
    }

    // allocate value array
    if ( columns_ ) {
        for ( i = 0; i < nColumns; i++ ) {
            delete this->columns_ [ i ];
        }

        delete this->columns_;
    }

    columns_ = ( FloatArray ** ) new FloatArray * [ neq ];
    for ( j = 0; j < neq; j++ ) {
        columns_ [ j ] = new FloatArray();
    }

    for ( n = 1; n <= nelem; n++ ) {
        elem = domain->giveElement(n);
        elem->giveLocationArray(loc, s);

        for ( i = 1; i <= loc.giveSize(); i++ ) {
            if ( ( ii = loc.at(i) ) ) {
                for ( j = 1; j <= loc.giveSize(); j++ ) {
                    if ( ( jj = loc.at(j) ) ) {
                        this->insertRowInColumn(ii - 1, jj - 1);
                    }
                }
            }
        }
    }

    // loop over active boundary conditions
    int nbc = domain->giveNumberOfBoundaryConditions();
    std :: vector< IntArray >r_locs;
    std :: vector< IntArray >c_locs;

    for ( int i = 1; i <= nbc; ++i ) {
        ActiveBoundaryCondition *bc = dynamic_cast< ActiveBoundaryCondition * >( domain->giveBc(i) );
        if ( bc != NULL ) {
            bc->giveLocationArrays(r_locs, c_locs, UnknownCharType, s, s);
            for ( std :: size_t k = 0; k < r_locs.size(); k++ ) {
                IntArray &krloc = r_locs [ k ];
                IntArray &kcloc = c_locs [ k ];
                for ( int i = 1; i <= krloc.giveSize(); i++ ) {
                    if ( ( ii = krloc.at(i) ) ) {
                        for ( int j = 1; j <= kcloc.giveSize(); j++ ) {
                            if ( ( jj = kcloc.at(j) ) ) {
                                this->insertRowInColumn(jj - 1, ii - 1);
                            }
                        }
                    }
                }
            }
        }
    }

    int nz_ = 0;
    for ( j = 0; j < neq; j++ ) {
        nz_ += this->rowind_ [ j ]->giveSize();
    }

    OOFEM_LOG_DEBUG("DynCompCol info: neq is %d, nelem is %d\n", neq, nz_);
#else
    nColumns = nRows = neq;

    columns.resize( neq );
    for ( auto &col: columns ) {
        col.clear();
    }

#endif

    // increment version
    this->version++;

    return true;
}

void IncrementalLinearStatic :: solveYourselfAt(TimeStep *tStep)
{
    Domain *d = this->giveDomain(1);
    // Creates system of governing eq's and solves them at given time step


    // >>> beginning PH
    // The following piece of code updates assignment of boundary conditions to dofs
    // (this allows to have multiple boundary conditions assigned to one dof
    // which can be arbitrarily turned on and off in time)
    // Almost the entire section has been copied from domain.C
    std :: vector< std :: map< int, int > > dof_bc( d->giveNumberOfDofManagers() );

    for ( int i = 1; i <= d->giveNumberOfBoundaryConditions(); ++i ) {
        GeneralBoundaryCondition *gbc = d->giveBc(i);

        if ( gbc->isImposed(tStep) ) {

            if ( gbc->giveSetNumber() > 0 ) { ///@todo This will eventually not be optional.
                // Loop over nodes in set and store the bc number in each dof.
                Set *set = d->giveSet( gbc->giveSetNumber() );
                ActiveBoundaryCondition *active_bc = dynamic_cast< ActiveBoundaryCondition * >(gbc);
                BoundaryCondition *bc = dynamic_cast< BoundaryCondition * >(gbc);
                if ( bc || ( active_bc && active_bc->requiresActiveDofs() ) ) {
                    const IntArray &appliedDofs = gbc->giveDofIDs();
                    const IntArray &nodes = set->giveNodeList();
                    for ( int inode = 1; inode <= nodes.giveSize(); ++inode ) {
                        for ( int idof = 1; idof <= appliedDofs.giveSize(); ++idof ) {

                            if  ( dof_bc [ nodes.at(inode) - 1 ].find( appliedDofs.at(idof) ) == dof_bc [ nodes.at(inode) - 1 ].end() ) {
                                // is empty
                                dof_bc [ nodes.at(inode) - 1 ] [ appliedDofs.at(idof) ] = i;

                                DofManager * dofman = d->giveDofManager( nodes.at(inode) );
                                Dof * dof = dofman->giveDofWithID( appliedDofs.at(idof) );

                                dof->setBcId(i);

                            } else {
                                // another bc has been already prescribed at this time step to this dof
                                OOFEM_WARNING("More than one boundary condition assigned at time %f to node %d dof %d. Considering boundary condition %d", tStep->giveTargetTime(),  nodes.at(inode), appliedDofs.at(idof), dof_bc [ nodes.at(inode) - 1 ] [appliedDofs.at(idof)] );


                            }
                        }
                    }
                }
            }
        }
    }

    // to get proper number of equations
    this->forceEquationNumbering();
    // <<< end PH



    // Initiates the total displacement to zero.
    if ( tStep->isTheFirstStep() ) {
        for ( auto &dofman : d->giveDofManagers() ) {
            for ( Dof *dof: *dofman ) {
                dof->updateUnknownsDictionary(tStep->givePreviousStep(), VM_Total, 0.);
                dof->updateUnknownsDictionary(tStep, VM_Total, 0.);
            }
        }

        for ( auto &bc : d->giveBcs() ) {
            ActiveBoundaryCondition *abc;

            if ( ( abc = dynamic_cast< ActiveBoundaryCondition * >(bc.get()) ) ) {
                int ndman = abc->giveNumberOfInternalDofManagers();
                for ( int i = 1; i <= ndman; i++ ) {
                    DofManager *dofman = abc->giveInternalDofManager(i);
                    for ( Dof *dof: *dofman ) {
                        dof->updateUnknownsDictionary(tStep->givePreviousStep(), VM_Total, 0.);
                        dof->updateUnknownsDictionary(tStep, VM_Total, 0.);
                    }
                }
            }
        }
    }

    // Apply dirichlet b.c's on total values
    for ( auto &dofman : d->giveDofManagers() ) {
        for ( Dof *dof: *dofman ) {
            double tot = dof->giveUnknown( VM_Total, tStep->givePreviousStep() );
            if ( dof->hasBc(tStep) ) {
                tot += dof->giveBcValue(VM_Incremental, tStep);
            }

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

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

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

//.........这里部分代码省略.........

/////////////////////////////////////////////////
// Help functions
void GnuplotExportModule::outputReactionForces(TimeStep *tStep)
{
    // Add sum of reaction forces to arrays
    // Compute sum of reaction forces for each BC number
    Domain *domain = emodel->giveDomain(1);
    StructuralEngngModel *seMod = dynamic_cast<StructuralEngngModel* >(emodel);
    if(seMod == NULL) {
        OOFEM_ERROR("failed to cast to StructuralEngngModel.");
    }

    IntArray ielemDofMask;
    FloatArray reactions;
    IntArray dofManMap, dofidMap, eqnMap;

    // test if solution step output is active
    if ( !testTimeStepOutput(tStep) ) {
        return;
    }

    // map contains corresponding dofmanager and dofs numbers corresponding to prescribed equations
    // sorted according to dofmanger number and as a minor crit. according to dof number
    // this is necessary for extractor, since the sorted output is expected
    seMod->buildReactionTable(dofManMap, dofidMap, eqnMap, tStep, 1);

    // compute reaction forces
    seMod->computeReaction(reactions, tStep, 1);

    // Find highest index of prescribed dofs
    int maxIndPresDof = 0;
    for ( int i = 1; i <= dofManMap.giveSize(); i++ ) {
        maxIndPresDof = std::max(maxIndPresDof, dofidMap.at(i));
    }

    int numBC = domain->giveNumberOfBoundaryConditions();

    while ( mReactionForceHistory.size() < size_t(numBC) ) {
        std::vector<FloatArray> emptyArray;
        mReactionForceHistory.push_back( emptyArray );
    }

    maxIndPresDof = domain->giveNumberOfSpatialDimensions();

    while ( mDispHist.size() < size_t(numBC) ) {
        std::vector<double> emptyArray;
        mDispHist.push_back( emptyArray );
    }

    for(int bcInd = 0; bcInd < numBC; bcInd++) {
        FloatArray fR(maxIndPresDof), disp(numBC);
        fR.zero();


        for ( int i = 1; i <= dofManMap.giveSize(); i++ ) {
            DofManager *dMan = domain->giveDofManager( dofManMap.at(i) );
            Dof *dof = dMan->giveDofWithID( dofidMap.at(i) );

            if ( dof->giveBcId() == bcInd+1 ) {
                fR.at( dofidMap.at(i) ) += reactions.at( eqnMap.at(i) );

                // Slightly dirty
                BoundaryCondition *bc = dynamic_cast<BoundaryCondition*> (domain->giveBc(bcInd+1));
                if ( bc != NULL ) {
                    disp.at(bcInd+1) = bc->give(dof, VM_Total, tStep->giveTargetTime());
                }
                ///@todo This function should be using the primaryfield instead of asking BCs directly. / Mikael
            }
        }

        mDispHist[bcInd].push_back(disp.at(bcInd+1));
        mReactionForceHistory[bcInd].push_back(fR);



        // X
        FILE * pFileX;
        char fileNameX[100];
        sprintf(fileNameX, "ReactionForceGnuplotBC%dX.dat", bcInd+1);
        pFileX = fopen ( fileNameX , "wb" );

        fprintf(pFileX, "#u Fx\n");
        for ( size_t j = 0; j < mDispHist[bcInd].size(); j++ ) {
            fprintf(pFileX, "%e %e\n", mDispHist[bcInd][j], mReactionForceHistory[bcInd][j].at(1) );
        }

        fclose(pFileX);

        // Y
        FILE * pFileY;
        char fileNameY[100];
        sprintf(fileNameY, "ReactionForceGnuplotBC%dY.dat", bcInd+1);
        pFileY = fopen ( fileNameY , "wb" );

        fprintf(pFileY, "#u Fx\n");
        for ( size_t j = 0; j < mDispHist[bcInd].size(); j++ ) {
            if( mReactionForceHistory[bcInd][j].giveSize() >= 2 ) {
                fprintf(pFileY, "%e %e\n", mDispHist[bcInd][j], mReactionForceHistory[bcInd][j].at(2) );
            }
        }
//.........这里部分代码省略.........

void GnuplotExportModule::doOutput(TimeStep *tStep, bool forcedOutput)
{
    if (!(testTimeStepOutput(tStep) || forcedOutput)) {
        return;
    }

    // Export the sum of reaction forces for each Dirichlet BC
    if(mExportReactionForces) {
        outputReactionForces(tStep);
    }

    Domain *domain = emodel->giveDomain(1);

    // Export output from boundary conditions
    if(mExportBoundaryConditions) {
        int numBC = domain->giveNumberOfBoundaryConditions();

        for(int i = 1; i <= numBC; i++) {

            PrescribedGradient *presGradBC = dynamic_cast<PrescribedGradient*>( domain->giveBc(i) );
            if(presGradBC != NULL) {
                outputBoundaryCondition(*presGradBC, tStep);
            }


            PrescribedGradientBCNeumann *presGradBCNeumann = dynamic_cast<PrescribedGradientBCNeumann*>( domain->giveBc(i) );
            if(presGradBCNeumann != NULL) {
                outputBoundaryCondition(*presGradBCNeumann, tStep);
            }

            PrescribedGradientBCWeak *presGradBCWeak = dynamic_cast<PrescribedGradientBCWeak*>( domain->giveBc(i) );
            if(presGradBCWeak != NULL) {
                outputBoundaryCondition(*presGradBCWeak, tStep);
            }

        }
    }

    mTimeHist.push_back( tStep->giveTargetTime() );

    if(mExportXFEM) {
        if(domain->hasXfemManager()) {
            XfemManager *xMan = domain->giveXfemManager();

            int numEI = xMan->giveNumberOfEnrichmentItems();

            std::vector< std::vector<FloatArray> > points;

            for(int i = 1; i <= numEI; i++) {
                EnrichmentItem *ei = xMan->giveEnrichmentItem(i);
                ei->callGnuplotExportModule(*this, tStep);

                GeometryBasedEI *geoEI = dynamic_cast<GeometryBasedEI*>(ei);
                if(geoEI != NULL) {
                    std::vector<FloatArray> eiPoints;
                    geoEI->giveSubPolygon(eiPoints, 0.0, 1.0);
                    points.push_back(eiPoints);
                }
            }

            outputXFEMGeometry(points);
        }
    }

    if(mExportMesh) {
        outputMesh(*domain);
    }

    if(mMonitorNodeIndex != -1) {
        DofManager *dMan = domain->giveDofManager(mMonitorNodeIndex);
        outputNodeDisp(*dMan, tStep);
    }
}

//.........这里部分代码省略.........
     *
     * // increment version
     * this->version++;
     *
     * return true;
     */
    IntArray loc;
    Domain *domain = eModel->giveDomain(di);
    int neq = eModel->giveNumberOfDomainEquations(di, s);
    int nelem = domain->giveNumberOfElements();
    int i, ii, j, jj, n, indx;
    Element *elem;
    // allocation map
    std :: vector< std :: set< int > > columns(neq);
    /*
     * std::set<int> **columns = new std::set<int>*[neq];
     * for (j=0; j<neq; j++) {
     * columns[j] = new std::set<int>;
     * }
     */

    this->nz_ = 0;

    for ( n = 1; n <= nelem; n++ ) {
        elem = domain->giveElement(n);
        elem->giveLocationArray(loc, s);

        for ( i = 1; i <= loc.giveSize(); i++ ) {
            if ( ( ii = loc.at(i) ) ) {
                for ( j = 1; j <= loc.giveSize(); j++ ) {
                    if ( ( jj = loc.at(j) ) ) {
                        columns [ jj - 1 ].insert(ii - 1);
                    }
                }
            }
        }
    }


    // loop over active boundary conditions
    int nbc = domain->giveNumberOfBoundaryConditions();
    std :: vector< IntArray >r_locs;
    std :: vector< IntArray >c_locs;

    for ( int i = 1; i <= nbc; ++i ) {
        ActiveBoundaryCondition *bc = dynamic_cast< ActiveBoundaryCondition * >( domain->giveBc(i) );
        if ( bc != NULL ) {
            bc->giveLocationArrays(r_locs, c_locs, UnknownCharType, s, s);
            for ( std :: size_t k = 0; k < r_locs.size(); k++ ) {
                IntArray &krloc = r_locs [ k ];
                IntArray &kcloc = c_locs [ k ];
                for ( int i = 1; i <= krloc.giveSize(); i++ ) {
                    if ( ( ii = krloc.at(i) ) ) {
                        for ( int j = 1; j <= kcloc.giveSize(); j++ ) {
                            if ( ( jj = kcloc.at(j) ) ) {
                                columns [ jj - 1 ].insert(ii - 1);
                            }
                        }
                    }
                }
            }
        }
    }

    for ( i = 0; i < neq; i++ ) {
        this->nz_ += columns [ i ].size();
    }

    rowind_.resize(nz_);
    colptr_.resize(neq + 1);
    indx = 0;

    for ( j = 0; j < neq; j++ ) { // column loop
        colptr_(j) = indx;
        for ( int row: columns [  j ] ) { // row loop
            rowind_(indx++) = row;
        }
    }

    colptr_(neq) = indx;

    /*
     * // delete map
     * for (i=0; i< neq; i++) {columns[i]->clear(); delete columns[i];}
     * delete columns;
     */

    // allocate value array
    val_.resize(nz_);
    val_.zero();

    OOFEM_LOG_DEBUG("CompCol info: neq is %d, nwk is %d\n", neq, nz_);

    dim_ [ 0 ] = dim_ [ 1 ] = nColumns = nRows = neq;

    // increment version
    this->version++;

    return true;
}

void
MatlabExportModule :: doOutputSpecials(TimeStep *tStep,    FILE *FID)
{
    FloatMatrix GradP, v_hat, GradPTemp, v_hatTemp;
    double dPdx, dPdy, Vu, Vv;

    Domain *domain  = emodel->giveDomain(1);

    GradP.resize(2, 1);
    GradP.zero();

    v_hat.resize(2, 1);
    v_hat.zero();

    for ( int i = 1; i <= domain->giveNumberOfElements(); i++ ) {
#ifdef __FM_MODULE
        Tr21Stokes *T = dynamic_cast< Tr21Stokes * >( domain->giveElement(i) );

        if ( T ) {
            T->giveGradP(GradPTemp, tStep);
            T->giveIntegratedVelocity(v_hatTemp, tStep);
            GradP.add(GradPTemp);
            v_hat.add(v_hatTemp);
        }

#endif
    }

    dPdx = GradP.at(1, 1);
    dPdy = GradP.at(2, 1);

    dPdx = dPdx / Area;
    dPdy = dPdy / Area;

    Vu = v_hat.at(1, 1) / ( ( xmax - xmin ) * ( ymax - ymin ) );
    Vv = v_hat.at(2, 1) / ( ( xmax - xmin ) * ( ymax - ymin ) );

    fprintf(FID, "\tspecials.gradpmean=[%f, %f];\n", dPdx, dPdy);
    fprintf(FID, "\tspecials.velocitymean=[%e, %e];\n", Vu, Vv);

    // Output weak periodic boundary conditions
    int wpbccount = 1;

    for ( int i = 1; i <= domain->giveNumberOfBoundaryConditions(); i++ ) {
        WeakPeriodicbc *wpbc = dynamic_cast< WeakPeriodicbc * >( domain->giveBc(i) );
        if ( wpbc ) {
            for ( int j = 1; j <= wpbc->giveNumberOfInternalDofManagers(); j++ ) {
                fprintf( FID, "\tspecials.weakperiodic{%u}.descType=%u;\n", wpbccount, wpbc->giveBasisType() );
                fprintf(FID, "\tspecials.weakperiodic{%u}.coefficients=[", wpbccount);
                for ( int k = 1; k <= wpbc->giveInternalDofManager(j)->giveNumberOfDofs(); k++ ) {
                    FloatArray unknowns;
                    IntArray DofMask;
                    double X = wpbc->giveInternalDofManager(j)->giveDof(k)->giveUnknown(EID_MomentumBalance_ConservationEquation, VM_Total, tStep);
                    fprintf(FID, "%e\t", X);
                }

                fprintf(FID, "];\n");
                wpbccount++;
            }
        }
    }
}

//.........这里部分代码省略.........
        MatCreate(PETSC_COMM_WORLD, & mtrx);
        MatSetSizes(mtrx, leqs, leqs, geqs, geqs);
        MatSetType(mtrx, MATMPIAIJ);
        MatSetFromOptions(mtrx);
        MatSetUp(mtrx);
        MatMPIAIJSetPreallocation(mtrx, 0, d_nnz.givePointer(), 0, o_nnz.givePointer());
        //MatMPIBAIJSetPreallocation( mtrx, PETSC_DECIDE, 0, d_nnz.givePointer(), onz, onnz );
        //MatMPISBAIJSetPreallocation( mtrx, PETSC_DECIDE, 0, d_nnz_sym.givePointer(), onz, onnz );

        MatSetOption(mtrx, MAT_ROW_ORIENTED, PETSC_FALSE); // To allow the insertion of values using MatSetValues in column major order
        MatSetOption(mtrx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

 #ifdef __VERBOSE_PARALLEL
        VERBOSEPARALLEL_PRINT("PetscSparseMtrx:: buildInternalStructure", "done", rank);
 #endif
    } else {
#endif

    leqs = geqs = eModel->giveNumberOfEquations(ut);
    int i, ii, j, jj, n;
    Element *elem;
    // allocation map
    std :: vector< std :: set< int > >rows(leqs);
    std :: vector< std :: set< int > >rows_sym(leqs);

    nelem = domain->giveNumberOfElements();
    for ( n = 1; n <= nelem; n++ ) {
        elem = domain->giveElement(n);
        elem->giveLocationArray(loc, ut, s);
        for ( i = 1; i <= loc.giveSize(); i++ ) {
            if ( ( ii = loc.at(i) ) ) {
                for ( j = 1; j <= loc.giveSize(); j++ ) {
                    jj = loc.at(j);
                    if ( jj ) {
                        rows [ ii - 1 ].insert(jj - 1);
                        if ( jj >= ii ) {
                            rows_sym [ ii - 1 ].insert(jj - 1);
                        }
                    }
                }
            }
        }
    }

    // Structure from active boundary conditions.
    AList<IntArray> locs, temp;
    for ( n = 1; n <= domain->giveNumberOfBoundaryConditions(); n++ ) {
        ActiveBoundaryCondition *activebc = dynamic_cast<ActiveBoundaryCondition*>(domain->giveBc(n));
        if (activebc) {
            ///@todo Deal with the CharType here.
            activebc->giveLocationArrays(locs, temp, ut, TangentStiffnessMatrix, s, s, domain);
            for (int k = 1; k < locs.giveSize(); k++) {
                IntArray *kloc = locs.at(k);
                for ( i = 1; i <= kloc->giveSize(); i++ ) {
                    if ( ( ii = kloc->at(i) ) ) {
                        for ( j = 1; j <= kloc->giveSize(); j++ ) {
                            jj = kloc->at(j);
                            if ( jj ) {
                                rows [ ii - 1 ].insert(jj - 1);
                                if ( jj >= ii ) {
                                    rows_sym [ ii - 1 ].insert(jj - 1);
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    IntArray d_nnz(leqs);
    IntArray d_nnz_sym(leqs);
    for ( i = 0; i < leqs; i++ ) {
        d_nnz(i) = rows [ i ].size();
        d_nnz_sym(i) = rows_sym [ i ].size();
    }

    MatCreate(PETSC_COMM_SELF, & mtrx);
    MatSetSizes(mtrx, leqs, leqs, geqs, geqs);
    MatSetType(mtrx, MATSEQAIJ);
    //MatSetType(mtrx, MATSBAIJ);
    //MatSetType(mtrx, MATDENSE);
    MatSetFromOptions(mtrx);

    MatSetUp(mtrx);
    MatSeqAIJSetPreallocation( mtrx, 0, d_nnz.givePointer() );
    MatSeqBAIJSetPreallocation( mtrx, PETSC_DECIDE, 0, d_nnz.givePointer() );
    MatSeqSBAIJSetPreallocation( mtrx, PETSC_DECIDE, 0, d_nnz_sym.givePointer() );

    MatSetOption(mtrx, MAT_ROW_ORIENTED, PETSC_FALSE); // To allow the insertion of values using MatSetValues in column major order
    MatSetOption(mtrx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

#ifdef __PARALLEL_MODE
    }
#endif

    nRows = nColumns = geqs;
    this->newValues = true;
    return true;
}

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
//.........这里部分代码省略.........

int DSSMatrix :: buildInternalStructure(EngngModel *eModel, int di, EquationID ut, const UnknownNumberingScheme &s)
{
    IntArray loc;
    Domain *domain = eModel->giveDomain(di);
    int neq = eModel->giveNumberOfDomainEquations(di, s);
    int nelem = domain->giveNumberOfElements();
    int i, ii, j, jj, n;
    unsigned long indx;
    Element *elem;
    // allocation map
    std :: vector< std :: set< int > >columns(neq);

    unsigned long nz_ = 0;

    for ( n = 1; n <= nelem; n++ ) {
        elem = domain->giveElement(n);
        elem->giveLocationArray(loc, ut, s);

        for ( i = 1; i <= loc.giveSize(); i++ ) {
            if ( ( ii = loc.at(i) ) ) {
                for ( j = 1; j <= loc.giveSize(); j++ ) {
                    if ( ( jj = loc.at(j) ) ) {
                        columns [ jj - 1 ].insert(ii - 1);
                    }
                }
            }
        }
    }

    // loop over active boundary conditions
    int nbc = domain->giveNumberOfBoundaryConditions();
    std::vector<IntArray> r_locs;
    std::vector<IntArray> c_locs;
    
    for ( i = 1; i <= nbc; ++i ) {
        ActiveBoundaryCondition *bc = dynamic_cast< ActiveBoundaryCondition * >( domain->giveBc(i) );
        if ( bc != NULL ) {
            bc->giveLocationArrays(r_locs, c_locs, ut, UnknownCharType, s, s);
	    for (std::size_t k = 0; k < r_locs.size(); k++) {
	      IntArray &krloc = r_locs[k];
	      IntArray &kcloc = c_locs[k];
	      for ( int ri = 1; ri <= krloc.giveSize(); ri++ ) {
		if ( ( ii = krloc.at(ri) ) ) {
		  for ( j = 1; j <= kcloc.giveSize(); j++ ) {
		    if ( (jj = kcloc.at(j) ) ) {
		      columns [ jj - 1 ].insert(ii - 1);
		    }
		  }
		}
	      }
	    }
	}
    }
    
    for ( i = 0; i < neq; i++ ) {
        nz_ += columns [ i ].size();
    }

    unsigned long *rowind_ = new unsigned long [ nz_ ];
    unsigned long *colptr_ = new unsigned long [ neq + 1 ];
    if ( ( rowind_ == NULL ) || ( colptr_ == NULL ) ) {
        OOFEM_ERROR("DSSMatrix::buildInternalStructure: free store exhausted, exiting");
    }

    indx = 0;

    std :: set< int > :: iterator pos;
    for ( j = 0; j < neq; j++ ) { // column loop
        colptr_ [ j ] = indx;
        for ( pos = columns [ j ].begin(); pos != columns [ j ].end(); ++pos ) { // row loop
            rowind_ [ indx++ ] = * pos;
        }
    }

    colptr_ [ neq ] = indx;

    if ( _sm ) {
        delete _sm;
    }

    if ( ( _sm = new SparseMatrixF(neq, NULL, rowind_, colptr_, 0, 0, true) ) == NULL ) {
        OOFEM_ERROR("DSSMatrix::buildInternalStructure: free store exhausted, exiting");
    }

    int bsize = eModel->giveDomain(1)->giveDefaultNodeDofIDArry().giveSize();
    /*
     *  Assemble block to equation mapping information
     */

    bool _succ = true;
    int _ndofs, _neq, ndofmans = domain->giveNumberOfDofManagers();
    int ndofmansbc = 0;
    // count number of internal dofmans on active bc
    for (n=1; n<=nbc; n++) {
      ndofmansbc+=domain->giveBc(n)->giveNumberOfInternalDofManagers();
    }

    long *mcn = new long [ (ndofmans+ndofmansbc) * bsize ];
    long _c = 0;
    DofManager *dman;
//.........这里部分代码省略.........