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

int
CBS :: checkConsistency()
{
    // check internal consistency
    // if success returns nonzero
    Domain *domain = this->giveDomain(1);

    // check for proper element type
    for ( auto &elem : domain->giveElements() ) {
        if ( !dynamic_cast< CBSElement * >( elem.get() ) ) {
            OOFEM_WARNING("Element %d has no CBS base", elem->giveLabel());
            return 0;
        }
    }

    EngngModel :: checkConsistency();


    // scale boundary and initial conditions
    if ( equationScalingFlag ) {
        for ( auto &bc: domain->giveBcs() ) {
            if ( bc->giveBCValType() == VelocityBVT ) {
                bc->scale(1. / uscale);
            } else if ( bc->giveBCValType() == PressureBVT ) {
                bc->scale( 1. / this->giveVariableScale(VST_Pressure) );
            } else if ( bc->giveBCValType() == ForceLoadBVT ) {
                bc->scale( 1. / this->giveVariableScale(VST_Force) );
            } else {
                OOFEM_WARNING("unknown bc/ic type");
                return 0;
            }
        }

        for ( auto &ic : domain->giveIcs() ) {
            if ( ic->giveICValType() == VelocityBVT ) {
                ic->scale(VM_Total, 1. / uscale);
            } else if ( ic->giveICValType() == PressureBVT ) {
                ic->scale( VM_Total, 1. / this->giveVariableScale(VST_Pressure) );
            } else {
                OOFEM_WARNING("unknown bc/ic type");
                return 0;
            }
        }
    }

    return 1;
}

bool
StaticStructural :: requiresEquationRenumbering(TimeStep *tStep)
{
    if ( tStep->isTheFirstStep() ) {
        return true;
    }
    // Check if Dirichlet b.c.s has changed.
    Domain *d = this->giveDomain(1);
    for ( auto &gbc : d->giveBcs() ) {
        ActiveBoundaryCondition *active_bc = dynamic_cast< ActiveBoundaryCondition * >(gbc.get());
        BoundaryCondition *bc = dynamic_cast< BoundaryCondition * >(gbc.get());
        // We only need to consider Dirichlet b.c.s
        if ( bc || ( active_bc && ( active_bc->requiresActiveDofs() || active_bc->giveNumberOfInternalDofManagers() ) ) ) {
            // Check of the dirichlet b.c. has changed in the last step (if so we need to renumber)
            if ( gbc->isImposed(tStep) != gbc->isImposed(tStep->givePreviousStep()) ) {
                return true;
            }
        }
    }
    return false;
}

bool
TransientTransportProblem :: requiresEquationRenumbering(TimeStep *tStep)
{
    ///@todo This method should be set as the default behavior instead of relying on a user specified flag. Then this function should be removed.
    if ( tStep->isTheFirstStep() ) {
        return true;
    }
    // Check if Dirichlet b.c.s has changed.
    Domain *d = this->giveDomain(1);
    for ( auto &gbc : d->giveBcs() ) {
        ActiveBoundaryCondition *active_bc = dynamic_cast< ActiveBoundaryCondition * >(gbc.get());
        BoundaryCondition *bc = dynamic_cast< BoundaryCondition * >(gbc.get());
        // We only need to consider Dirichlet b.c.s
        if ( bc || ( active_bc && ( active_bc->requiresActiveDofs() || active_bc->giveNumberOfInternalDofManagers() ) ) ) {
            // Check of the dirichlet b.c. has changed in the last step (if so we need to renumber)
            if ( gbc->isImposed(tStep) != gbc->isImposed(tStep->givePreviousStep()) ) {
                return true;
            }
        }
    }
    return false;
}

void
MatlabExportModule :: doOutputSpecials(TimeStep *tStep,    FILE *FID)
{
//    FloatArray v_hat, GradPTemp, v_hatTemp;

    Domain *domain  = emodel->giveDomain(1);
/*
    v_hat.clear();

        ///@todo Sort out this hack in a nicer (modular) way / Mikael
#if 0
    for ( auto &elem : domain->giveElements() ) {
#ifdef __FM_MODULE

        if ( Tr21Stokes *Tr = dynamic_cast< Tr21Stokes * >( elem.get() ) ) {
            Tr->giveIntegratedVelocity(v_hatTemp, tStep);
            v_hat.add(v_hatTemp);
        } else if ( Tet21Stokes *Tet = dynamic_cast< Tet21Stokes * >( elem.get() ) ) {
            Tet->giveIntegratedVelocity(v_hatTemp, tStep);
            v_hat.add(v_hatTemp);
        }

#endif
    }
#endif

    // Compute intrinsic area/volume
    double intrinsicSize = 1.0;

    std :: vector <double> V;

    for ( int i = 0; i < (int)smax.size(); i++ ) {
        intrinsicSize *= ( smax.at(i) - smin.at(i) );
    }

    for ( double vh: v_hat ) {
        V.push_back(vh);
    }

    fprintf(FID, "\tspecials.velocitymean=[");
    if (V.size()>0) {
        for (int i=0; i<ndim; i++) {
            fprintf(FID, "%e", V.at(i));
            if (i!=(ndim-1)) fprintf (FID, ", ");
        }
        fprintf(FID, "];\n");
    } else {
        fprintf(FID, "]; %% No velocities\n");
    }

    */

    // Output weak periodic boundary conditions
    unsigned int wpbccount = 1, sbsfcount = 1, mcount = 1;

    for ( auto &gbc : domain->giveBcs() ) {
        WeakPeriodicBoundaryCondition *wpbc = dynamic_cast< WeakPeriodicBoundaryCondition * >( gbc.get() );
        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 ( Dof *dof: *wpbc->giveInternalDofManager(j) ) {
                    double X = dof->giveUnknown(VM_Total, tStep);
                    fprintf(FID, "%e\t", X);
                }

                fprintf(FID, "];\n");
                wpbccount++;
            }
        }
        SolutionbasedShapeFunction *sbsf = dynamic_cast< SolutionbasedShapeFunction *>( gbc.get());
        if (sbsf) {
            fprintf(FID, "\tspecials.solutionbasedsf{%u}.values=[", sbsfcount);
            for ( Dof *dof: *sbsf->giveInternalDofManager(1) ) {                  // Only one internal dof manager
                double X = dof->giveUnknown(VM_Total, tStep);
                fprintf(FID, "%e\t", X);
            }
            fprintf(FID, "];\n");
            sbsfcount++;
        }
        PrescribedMean *m = dynamic_cast<PrescribedMean *> ( gbc.get() );
        if (m) {
            fprintf(FID, "\tspecials.prescribedmean{%u}.value=[", mcount);
            for ( Dof *dof: *m->giveInternalDofManager(1)) {
                double X = dof->giveUnknown(VM_Total, tStep);
                fprintf(FID, "%e\t", X);
            }
            fprintf(FID, "];\n");
            mcount++;
        }
    }
}

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

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

int Skyline :: buildInternalStructure(EngngModel *eModel, int di, const UnknownNumberingScheme &s)
{
    // first create array of
    // maximal column height for assembled characteristics matrix
    //

    int maxle;
    int ac1;
    int neq;
    if ( s.isDefault() ) {
        neq = eModel->giveNumberOfDomainEquations(di, s);
    } else {
        neq = s.giveRequiredNumberOfDomainEquation();
    }
    if ( neq == 0 ) {
        if ( mtrx ) {
            delete mtrx;
        }
        mtrx = NULL;
        adr.clear();
        return true;
    }

    IntArray loc;
    IntArray mht(neq);
    Domain *domain = eModel->giveDomain(di);

    for ( int j = 1; j <= neq; j++ ) {
        mht.at(j) = j; // initialize column height, maximum is line number (since it only stores upper triangular)
    }

    // loop over elements code numbers
    for ( auto &elem : domain->giveElements() ) {
        elem->giveLocationArray(loc, s);
        maxle = INT_MAX;
        for ( int ieq : loc ) {
            if ( ieq != 0 ) {
                maxle = min(maxle, ieq);
            }
        }

        for ( int ieq : loc ) {
            if ( ieq != 0 ) {
                mht.at(ieq) = min( maxle, mht.at(ieq) );
            }
        }
    }

    // loop over active boundary conditions (e.g. relative kinematic constraints)
    std :: vector< IntArray >r_locs;
    std :: vector< IntArray >c_locs;

    for ( auto &gbc : domain->giveBcs() ) {
        ActiveBoundaryCondition *bc = dynamic_cast< ActiveBoundaryCondition * >( gbc.get() );
        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 ];
                maxle = INT_MAX;
                for ( int ii : krloc ) {
                    if ( ii > 0 ) {
                        maxle = min(maxle, ii);
                    }
                }
                for ( int jj : kcloc ) {
                    if ( jj > 0 ) {
                        mht.at(jj) = min( maxle, mht.at(jj) );
                    }
                }
            }
        }
    }

    
    if ( domain->hasContactManager() ) {
        ContactManager *cMan = domain->giveContactManager();
            
        for ( int i =1; i <= cMan->giveNumberOfContactDefinitions(); i++ ) {
            ContactDefinition *cDef = cMan->giveContactDefinition(i);
            for ( int k = 1; k <= cDef->giveNumbertOfContactElements(); k++ ) {
                ContactElement *cEl = cDef->giveContactElement(k);
                cEl->giveLocationArray(loc, s);
            
                maxle = INT_MAX;
                for ( int ieq : loc ) {
                    if ( ieq != 0 ) {
                        maxle = min(maxle, ieq);
                    }
                }

                for ( int ieq : loc ) {
                    if ( ieq != 0 ) {
                        mht.at(ieq) = min( maxle, mht.at(ieq) );
                    }
                  
                }
            }
        }
    }
//.........这里部分代码省略.........

int DSSMatrix :: buildInternalStructure(EngngModel *eModel, int di, const UnknownNumberingScheme &s)
{
    IntArray loc;
    Domain *domain = eModel->giveDomain(di);
    int neq = eModel->giveNumberOfDomainEquations(di, s);
    unsigned long indx;
    // allocation map
    std :: vector< std :: set< int > >columns(neq);

    unsigned long nz_ = 0;

    for ( auto &elem : domain->giveElements() ) {
        elem->giveLocationArray(loc, s);

        for ( int ii : loc ) {
            if ( ii > 0 ) {
                for ( int jj : loc ) {
                    if ( jj > 0 ) {
                        columns [ jj - 1 ].insert(ii - 1);
                    }
                }
            }
        }
    }

    // loop over active boundary conditions
    std::vector<IntArray> r_locs;
    std::vector<IntArray> c_locs;
    
    for ( auto &gbc : domain->giveBcs() ) {
        ActiveBoundaryCondition *bc = dynamic_cast< ActiveBoundaryCondition * >( gbc.get() );
        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 ii : krloc ) {
                    if ( ii > 0 ) {
                        for ( int jj : kcloc ) {
                            if ( jj > 0 ) {
                                columns [ jj - 1 ].insert(ii - 1);
                            }
                        }
                    }
                }
            }
        }
    }
    
    for ( int 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("free store exhausted, exiting");
    }

    indx = 0;

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

    colptr_ [ neq ] = indx;

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

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