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

void
NonLinearStatic :: showSparseMtrxStructure(int type, oofegGraphicContext &gc, TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    CharType ctype;

    if ( type != 1 ) {
        return;
    }

    if ( stiffMode == nls_tangentStiffness ) {
        ctype = TangentStiffnessMatrix;
    } else if ( stiffMode == nls_secantStiffness ) {
        ctype = SecantStiffnessMatrix;
    } else {
        ctype = SecantStiffnessMatrix;
    }

    for ( auto &elem : domain->giveElements() ) {
        elem->showSparseMtrxStructure(ctype, gc, tStep);
    }

    for ( auto &elem : domain->giveElements() ) {
        elem->showExtendedSparseMtrxStructure(ctype, gc, tStep);
    }
}

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

TimeStep *
CBS :: giveNextStep()
{
    double dt = deltaT;
    if ( !currentStep ) {
        // first step -> generate initial step
        currentStep.reset( new TimeStep( *giveSolutionStepWhenIcApply() ) );
    }

    previousStep = std :: move(currentStep);

    Domain *domain = this->giveDomain(1);
    // check for critical time step
    for ( auto &elem : domain->giveElements() ) {
        dt = min( dt, static_cast< CBSElement & >( *elem ).computeCriticalTimeStep(previousStep.get()) );
    }

    dt *= 0.6;
    dt = max(dt, minDeltaT);
    dt /= this->giveVariableScale(VST_Time);

    currentStep.reset( new TimeStep(*previousStep, dt) );

    OOFEM_LOG_INFO( "SolutionStep %d : t = %e, dt = %e\n", currentStep->giveNumber(),
                    currentStep->giveTargetTime() * this->giveVariableScale(VST_Time), dt * this->giveVariableScale(VST_Time) );

    return currentStep.get();
}

// needed for CemhydMat
void
NonStationaryTransportProblem :: averageOverElements(TimeStep *tStep)
{
    ///@todo Verify this, the function is completely unused.
    Domain *domain = this->giveDomain(1);
    FloatArray vecTemperature;

    for ( auto &elem : domain->giveElements() ) {
        TransportMaterial *mat = dynamic_cast< CemhydMat * >( elem->giveMaterial() );
        if ( mat ) {
            for ( GaussPoint *gp: *elem->giveDefaultIntegrationRulePtr() ) {
                elem->giveIPValue(vecTemperature, gp, IST_Temperature, tStep);
                //mat->IP_volume += dV;
                //mat->average_temp += vecState.at(1) * dV;
            }
        }
    }

    for ( auto &mat : domain->giveMaterials() ) {
        CemhydMat *cem = dynamic_cast< CemhydMat * >( mat.get() );
        if ( cem ) {
            //cem->average_temp /= mat->IP_volume;
        }
    }
}

void
NonLinearDynamic :: showSparseMtrxStructure(int type, oofegGraphicContext &gc, TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    CharType ctype;

    if ( type != 1 ) {
        return;
    }

    ctype = TangentStiffnessMatrix;

    for ( auto &elem : domain->giveElements() ) {
        elem->showSparseMtrxStructure(ctype, gc, tStep);
    }

    for ( auto &elem : domain->giveElements() ) {
        elem->showExtendedSparseMtrxStructure(ctype, gc, tStep);
    }
}

void
StructuralEngngModel :: showSparseMtrxStructure(int type, oofegGraphicContext &gc, TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);

    if ( type != 1 ) {
        return;
    }

    for ( auto &elem : domain->giveElements() ) {
        elem->showSparseMtrxStructure(TangentStiffnessMatrix, gc, tStep);
    }
}

int StokesFlow :: checkConsistency()
{
    Domain *domain = this->giveDomain(1);

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

    return EngngModel :: checkConsistency();
}

void
NLTransientTransportProblem :: applyIC(TimeStep *stepWhenIcApply)
{
    Domain *domain = this->giveDomain(1);

    NonStationaryTransportProblem :: applyIC(stepWhenIcApply);

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

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
NonlocalMaterialExtensionInterface :: updateDomainBeforeNonlocAverage(TimeStep *tStep)
{
    Domain *d = this->giveDomain();

    if ( d->giveNonlocalUpdateStateCounter() == tStep->giveSolutionStateCounter() ) {
        return; // already updated
    }

    for ( auto &elem : d->giveElements() ) {
        elem->updateBeforeNonlocalAverage(tStep);
    }

    // mark last update counter to prevent multiple updates
    d->setNonlocalUpdateStateCounter( tStep->giveSolutionStateCounter() );
}

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

void
TransientTransportProblem :: applyIC()
{
    Domain *domain = this->giveDomain(1);
    OOFEM_LOG_INFO("Applying initial conditions\n");

    this->field->applyDefaultInitialCondition();

    ///@todo It's rather strange that the models need the initial values.
    // update element state according to given ic
    TimeStep *s = this->giveSolutionStepWhenIcApply();
    for ( auto &elem : domain->giveElements() ) {
        TransportElement *element = static_cast< TransportElement * >( elem.get() );
        element->updateInternalState(s);
        element->updateYourself(s);
    }
}

void GnuplotExportModule::outputMesh(Domain &iDomain)
{
    std::vector< std::vector<FloatArray> > pointArray;

    if(iDomain.giveNumberOfSpatialDimensions() == 2) {
        // Write all element edges to gnuplot
        for ( auto &el : iDomain.giveElements() ) {
            int numEdges = el->giveNumberOfNodes();


            for ( int edgeIndex = 1; edgeIndex <= numEdges; edgeIndex++ ) {
                std::vector<FloatArray> points;

                IntArray bNodes;
                el->giveInterpolation()->boundaryGiveNodes(bNodes, edgeIndex);

                int niLoc = bNodes.at(1);
                const FloatArray &xS = *(el->giveNode(niLoc)->giveCoordinates() );
                points.push_back(xS);

                int njLoc = bNodes.at( bNodes.giveSize() );
                const FloatArray &xE = *(el->giveNode(njLoc)->giveCoordinates() );
                points.push_back(xE);

                pointArray.push_back(points);
            }

        }


        double time = 0.0;

        TimeStep *ts = emodel->giveCurrentStep();
        if ( ts != NULL ) {
            time = ts->giveTargetTime();
        }

        std :: stringstream strMesh;
        strMesh << "MeshGnuplotTime" << time << ".dat";
        std :: string nameMesh = strMesh.str();

        WritePointsToGnuplot(nameMesh, pointArray);
    }
}

void
StokesFlowVelocityHomogenization :: computeTangent(FloatMatrix &answer, TimeStep *tStep)
{
    IntArray loc, col;

    Domain *domain = this->giveDomain(1);
    int nsd = domain->giveNumberOfSpatialDimensions();
    int ndof = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );

    // Build F matrix
    IntArray dofs(nsd);
    for ( int i = 0; i < nsd; ++i ) {
        dofs[i] = V_u + i; ///@todo This is a hack. We should have these as user input instead.
    }
    FloatMatrix F(ndof, nsd), Fe, N;
    col.enumerate(nsd);

    for ( auto &elem : domain->giveElements() ) {

        this->integrateNMatrix(N, *elem, tStep);
        
        elem->giveLocationArray( loc, dofs, EModelDefaultEquationNumbering() );
        Fe.beTranspositionOf(N);
        F.assemble(Fe, loc, col);
    }

    FloatMatrix H;

    std :: unique_ptr< SparseLinearSystemNM > solver( classFactory.createSparseLinSolver(solverType, this->giveDomain(1), this) );

    H.resize( F.giveNumberOfRows(), F.giveNumberOfColumns() );
    H.zero();

    // For correct homogenization, the tangent at the converged values should be used
    // (the one from the Newton iterations from solveYourselfAt is not updated to contain the latest values).
    SparseMtrxType stype = solver->giveRecommendedMatrix(true);
    std :: unique_ptr< SparseMtrx > Kff( classFactory.createSparseMtrx( stype ) );
    Kff->buildInternalStructure(this, domain->giveNumber(), EModelDefaultEquationNumbering() );
    this->assemble(*Kff, tStep, TangentStiffnessMatrix, EModelDefaultEquationNumbering(), domain);
    solver->solve(*Kff, F, H);

    answer.beTProductOf(H, F);
    answer.times( 1. / this->giveAreaOfRVE() );
}

int
NonStationaryTransportProblem :: 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< TransportElement * >( elem.get() ) ) {
            OOFEM_WARNING("Element %d has no TransportElement base", elem->giveLabel());
            return 0;
        }
    }

    EngngModel :: checkConsistency();

    return 1;
}