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

void
MatlabExportModule :: doOutputMesh(TimeStep *tStep, FILE *FID)
{
    Domain *domain  = emodel->giveDomain(1);

    fprintf(FID, "\tmesh.p=[");
    for ( auto &dman : domain->giveDofManagers() ) {
        for ( int j = 1; j <= domain->giveNumberOfSpatialDimensions(); j++) {
            double c = dman->giveCoordinate(j);
            fprintf(FID, "%f, ", c);
        }
        fprintf(FID, "; ");
    }

    fprintf(FID, "]';\n");

    int numberOfDofMans=domain->giveElement(1)->giveNumberOfDofManagers();

    fprintf(FID, "\tmesh.t=[");
    for ( auto &elem : domain->giveElements() ) {
        if ( elem->giveNumberOfDofManagers() == numberOfDofMans ) {
            for ( int j = 1; j <= elem->giveNumberOfDofManagers(); j++ ) {
                fprintf( FID, "%d,", elem->giveDofManagerNumber(j) );
            }
        }
        fprintf(FID, ";");
    }

    fprintf(FID, "]';\n");
}

void
NLTransientTransportProblem :: copyUnknownsInDictionary(ValueModeType mode, TimeStep *fromTime, TimeStep *toTime)
{
    Domain *domain = this->giveDomain(1);

    for ( auto &node : domain->giveDofManagers() ) {
        for ( Dof *dof: *node ) {
            double val = dof->giveUnknown(mode, fromTime);
            dof->updateUnknownsDictionary(toTime, mode, val);
        }
    }
}

void 
NonLinearDynamic :: initializeYourself(TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    int neq = this->giveNumberOfDomainEquations(1, EModelDefaultEquationNumbering());

    if ( tStep->isTheFirstStep() && initFlag ) {
        // Initialization
        incrementOfDisplacement.resize(neq);
        incrementOfDisplacement.zero();
        totalDisplacement.resize(neq);
        totalDisplacement.zero();
        velocityVector.resize(neq);
        velocityVector.zero();
        accelerationVector.resize(neq);
        accelerationVector.zero();
        internalForces.resize(neq);
        internalForces.zero();
        previousIncrementOfDisplacement.resize(neq);
        previousIncrementOfDisplacement.zero();
        previousTotalDisplacement.resize(neq);
        previousTotalDisplacement.zero();
        previousVelocityVector.resize(neq);
        previousVelocityVector.zero();
        previousAccelerationVector.resize(neq);
        previousAccelerationVector.zero();
        previousInternalForces.resize(neq);
        previousInternalForces.zero();

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

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

                int jj = dof->__giveEquationNumber();
                if ( jj ) {
                    totalDisplacement.at(jj)  = dof->giveUnknown(VM_Total, stepWhenIcApply);
                    velocityVector.at(jj)     = dof->giveUnknown(VM_Velocity, stepWhenIcApply);
                    accelerationVector.at(jj) = dof->giveUnknown(VM_Acceleration, stepWhenIcApply);
                }
            }
        }
        this->giveInternalForces(internalForces, true, 1, tStep);
    }
}

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

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

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

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

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

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

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

void
MatlabExportModule :: doOutputData(TimeStep *tStep, FILE *FID)
{
    Domain *domain  = emodel->giveDomain(1);
    std :: vector< int >DofIDList;
    std :: vector< int > :: iterator it;
    std :: vector< std :: vector< double > >valuesList;

    for ( auto &dman : domain->giveDofManagers() ) {
        for ( Dof *thisDof: *dman ) {
            it = std :: find( DofIDList.begin(), DofIDList.end(), thisDof->giveDofID() );

            double value = thisDof->giveUnknown(VM_Total, tStep);
            if ( it == DofIDList.end() ) {
                DofIDList.push_back( thisDof->giveDofID() );
                valuesList.push_back({value});
            } else {
                std::size_t pos = it - DofIDList.begin();
                valuesList[pos].push_back(value);
            }
        }
    }

    fprintf(FID, "\tdata.DofIDs=[");
    for ( auto &dofid : DofIDList ) {
        fprintf( FID, "%d, ", dofid );
    }

    fprintf(FID, "];\n");

    for ( size_t i = 0; i < valuesList.size(); i++ ) {
        fprintf(FID, "\tdata.a{%lu}=[", static_cast< long unsigned int >(i) + 1);
        for ( double val: valuesList[i] ) {
            fprintf( FID, "%f,", val );
        }

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

}

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

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

void
NonStationaryTransportProblem :: applyIC(TimeStep *stepWhenIcApply)
{
    Domain *domain = this->giveDomain(1);
    int neq =  this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );
    FloatArray *solutionVector;
    double val;

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

    UnknownsField->advanceSolution(stepWhenIcApply);
    solutionVector = UnknownsField->giveSolutionVector(stepWhenIcApply);
    solutionVector->resize(neq);
    solutionVector->zero();

    for ( auto &node : domain->giveDofManagers() ) {

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

            int jj = dof->__giveEquationNumber();
            if ( jj ) {
                val = dof->giveUnknown(VM_Total, stepWhenIcApply);
                solutionVector->at(jj) = val;
                //update in dictionary, if the problem is growing/decreasing
                if ( this->changingProblemSize ) {
                    dof->updateUnknownsDictionary(stepWhenIcApply, VM_Total, val);
                }
            }
        }
    }


    //project initial temperature to integration points

    //     for ( int j = 1; j <= nelem; j++ ) {
    //         domain->giveElement(j)->updateInternalState(stepWhenIcApply);
    //     }

#ifdef __CEMHYD_MODULE
    // Not relevant in linear case, but needed for CemhydMat for temperature averaging before solving balance equations
    // Update element state according to given ic
    for ( auto &elem : domain->giveElements() ) {
        TransportElement *element = static_cast< TransportElement * >( elem.get() );
        CemhydMat *cem = dynamic_cast< CemhydMat * >( element->giveMaterial() );
        //assign status to each integration point on each element
        if ( cem ) {
            cem->initMaterial(element); //create microstructures and statuses on specific GPs
            element->updateInternalState(stepWhenIcApply);   //store temporary unequilibrated temperature
            element->updateYourself(stepWhenIcApply);   //store equilibrated temperature
            cem->clearWeightTemperatureProductVolume(element);
            cem->storeWeightTemperatureProductVolume(element, stepWhenIcApply);
        }
    }

    //perform averaging on each material instance of CemhydMatClass
    for ( auto &mat : domain->giveMaterials() ) {
        CemhydMat *cem = dynamic_cast< CemhydMat * >( mat.get() );
        if ( cem ) {
            cem->averageTemperature();
        }
    }

#endif //__CEMHYD_MODULE
}

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

    _sm.reset( new SparseMatrixF(neq, NULL, rowind_, colptr_, 0, 0, true) ); 
    if ( !_sm ) {
        OOFEM_FATAL("free store exhausted, exiting");
    }


    /*
     *  Assemble block to equation mapping information
     */

    bool _succ = true;
    int _ndofs, _neq, ndofmans = domain->giveNumberOfDofManagers();
    int ndofmansbc = 0;

    ///@todo This still misses element internal dofs.
    // count number of internal dofmans on active bc
    for ( auto &bc : domain->giveBcs() ) {
        ndofmansbc += bc->giveNumberOfInternalDofManagers();
    }

    int bsize = 0;
    if ( ndofmans > 0 ) {
        bsize = domain->giveDofManager(1)->giveNumberOfDofs();
    }

    long *mcn = new long [ (ndofmans+ndofmansbc) * bsize ];
    long _c = 0;

    if ( mcn == NULL ) {
        OOFEM_FATAL("free store exhausted, exiting");
    }

    for ( auto &dman : domain->giveDofManagers() ) {
        _ndofs = dman->giveNumberOfDofs();
        if ( _ndofs > bsize ) {
            _succ = false;
            break;
        }

        for ( Dof *dof: *dman ) {
            if ( dof->isPrimaryDof() ) {
                _neq = dof->giveEquationNumber(s);
                if ( _neq > 0 ) {
                    mcn [ _c++ ] = _neq - 1;
                } else {
                    mcn [ _c++ ] = -1; // no corresponding row in sparse mtrx structure
                }
            } else {
                    mcn [ _c++ ] = -1; // no corresponding row in sparse mtrx structure
            }
        }

        for ( int i = _ndofs + 1; i <= bsize; i++ ) {
            mcn [ _c++ ] = -1;                         // no corresponding row in sparse mtrx structure
        }
    }

    // loop over internal dofmans of active bc
    for ( auto &bc : domain->giveBcs() ) {
      int ndman = bc->giveNumberOfInternalDofManagers();
      for (int idman = 1; idman <= ndman; idman ++) {
            DofManager *dman = bc->giveInternalDofManager(idman);
            _ndofs = dman->giveNumberOfDofs();
            if ( _ndofs > bsize ) {
                _succ = false;

void TransientTransportProblem :: solveYourselfAt(TimeStep *tStep)
{
    Domain *d = this->giveDomain(1);
    int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );

    if ( tStep->isTheFirstStep() ) {
        this->applyIC();
    }

    field->advanceSolution(tStep);

#if 1
    // This is what advanceSolution should be doing, but it can't be there yet 
    // (backwards compatibility issues due to inconsistencies in other solvers).
    TimeStep *prev = tStep->givePreviousStep();
    for ( auto &dman : d->giveDofManagers() ) {
        static_cast< DofDistributedPrimaryField* >(field.get())->setInitialGuess(*dman, tStep, prev);
    }

    for ( auto &elem : d->giveElements() ) {
        int ndman = elem->giveNumberOfInternalDofManagers();
        for ( int i = 1; i <= ndman; i++ ) {
            static_cast< DofDistributedPrimaryField* >(field.get())->setInitialGuess(*elem->giveInternalDofManager(i), tStep, prev);
        }
    }

    for ( auto &bc : d->giveBcs() ) {
        int ndman = bc->giveNumberOfInternalDofManagers();
        for ( int i = 1; i <= ndman; i++ ) {
            static_cast< DofDistributedPrimaryField* >(field.get())->setInitialGuess(*bc->giveInternalDofManager(i), tStep, prev);
        }
    }
#endif

    field->applyBoundaryCondition(tStep);
    field->initialize(VM_Total, tStep, solution, EModelDefaultEquationNumbering());


    if ( !effectiveMatrix ) {
        effectiveMatrix.reset( classFactory.createSparseMtrx(sparseMtrxType) );
        effectiveMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );

        if ( lumped ) {
            capacityDiag.resize(neq);
            this->assembleVector( capacityDiag, tStep, LumpedMassVectorAssembler(), VM_Total, EModelDefaultEquationNumbering(), d );
        } else {
            capacityMatrix.reset( classFactory.createSparseMtrx(sparseMtrxType) );
            capacityMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );
            this->assemble( *capacityMatrix, tStep, MassMatrixAssembler(), EModelDefaultEquationNumbering(), d );
        }
        
        if ( this->keepTangent ) {
            this->assemble( *effectiveMatrix, tStep, TangentAssembler(TangentStiffness),
                           EModelDefaultEquationNumbering(), d );
            effectiveMatrix->times(alpha);
            if ( lumped ) {
                effectiveMatrix->addDiagonal(1./tStep->giveTimeIncrement(), capacityDiag);
            } else {
                effectiveMatrix->add(1./tStep->giveTimeIncrement(), *capacityMatrix);
            }
        }
    }


    OOFEM_LOG_INFO("Assembling external forces\n");
    FloatArray externalForces(neq);
    externalForces.zero();
    this->assembleVector( externalForces, tStep, ExternalForceAssembler(), VM_Total, EModelDefaultEquationNumbering(), d );
    this->updateSharedDofManagers(externalForces, EModelDefaultEquationNumbering(), LoadExchangeTag);

    // set-up numerical method
    this->giveNumericalMethod( this->giveCurrentMetaStep() );
    OOFEM_LOG_INFO("Solving for %d unknowns...\n", neq);

    internalForces.resize(neq);

    FloatArray incrementOfSolution;
    double loadLevel;
    int currentIterations;
    this->nMethod->solve(*this->effectiveMatrix,
                         externalForces,
                         NULL, // ignore
                         NULL,
                         this->solution,
                         incrementOfSolution,
                         this->internalForces,
                         this->eNorm,
                         loadLevel, // ignore
                         SparseNonLinearSystemNM :: rlm_total, // ignore
                         currentIterations, // ignore
                         tStep);
}