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

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

    for ( int ielem = 1; ielem <= nelem; ielem++ ) {
        TransportElement *element = static_cast< TransportElement * >( domain->giveElement(ielem) );
        TransportMaterial *mat = static_cast< CemhydMat * >( element->giveMaterial() );
        if ( mat ) {
            for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
                element->giveIPValue(vecTemperature, gp, IST_Temperature, tStep);
                //mat->IP_volume += dV;
                //mat->average_temp += vecState.at(1) * dV;
            }
        }
    }

    for ( int i = 1; i <= domain->giveNumberOfMaterialModels(); i++ ) {
        CemhydMat *mat = static_cast< CemhydMat * >( domain->giveMaterial(i) );
        if ( mat ) {
            //mat->average_temp /= mat->IP_volume;
        }
    }
}

void
NonStationaryTransportProblem :: updateYourself(TimeStep *tStep)
{
    this->updateInternalState(tStep);
    EngngModel :: updateYourself(tStep);

    ///@todo Find a cleaner way to do these cemhyd hacks
#ifdef __CEMHYD_MODULE
    for ( int idomain = 1; idomain <= ndomains; idomain++ ) {
        Domain *d = this->giveDomain(idomain);
        for ( int i = 1; i <= d->giveNumberOfElements(); ++i ) {
            TransportElement *elem = static_cast< TransportElement * >( d->giveElement(i) );
            //store temperature and associated volume on each GP before performing averaging
            CemhydMat *cem = dynamic_cast< CemhydMat * >( elem->giveMaterial() );
            if ( cem ) {
                cem->clearWeightTemperatureProductVolume(elem);
                cem->storeWeightTemperatureProductVolume(elem, tStep);
            }
        }
        //perform averaging on each material instance
        for ( int i = 1; i <= d->giveNumberOfMaterialModels(); i++ ) {
            CemhydMat *cem = dynamic_cast< CemhydMat * >( d->giveMaterial(i) );
            if ( cem ) {
                cem->averageTemperature();
            }
        }
    }
 #ifdef VERBOSE
    VERBOSE_PRINT0("Updated Materials ", 0)
 #endif
#endif
}

void
StaticFracture :: optimize(TimeStep *tStep)
{
    // Main optimization loop

    MinCompliance *objFunc = dynamic_cast< MinCompliance * >( this->objFuncList[0] ); 

    for ( int subProb = 1; subProb <= this->giveNumberOfSlaveProblems(); subProb++ ) {
        
        EngngModel *sp = this->giveSlaveProblem(subProb);
        Domain *d = sp->giveDomain(1);   

        double cost = 0.0; // should lie in obj fnc
        double dce = 0.0;
        for (int i = 1; i <= d->giveNumberOfElements(); i++) {  
            Element *el = d->giveElement(i);
            cost += objFunc->evaluateYourself(el, dce, sp->giveCurrentStep() ); // add cost for each element
        }


        // Filter sensitivities
        this->filterSensitivities(objFunc);

        // Update design variables based on some method. For now use the 'standard' optimality criteria
        this->optimalityCriteria(objFunc);
     

        // Update material parameters
        for (int i = 1; i <= d->giveNumberOfMaterialModels(); i++) {
            DynamicInputRecord ir;
            Material *mat = d->giveMaterial(i);
            mat->giveInputRecord(ir);
            double E0 = 1.0;
            double fac = pow( objFunc->designVarList.at(i), objFunc->penalty);
            ir.setField( E0 * fac, _IFT_IsotropicLinearElasticMaterial_e);
            mat->initializeFrom(&ir);
        }

        printf("\n costfunction %e & sum sensitivity %e & sum x %e \n", cost, objFunc->sensitivityList.sum(), 
            objFunc->designVarList.sum() );

    }
}

// needed for CemhydMat
void
NonStationaryTransportProblem :: averageOverElements(TimeStep *tStep)
{
    Domain *domain = this->giveDomain(1);
    int ielem, i;
    int nelem = domain->giveNumberOfElements();
    double dV;
    TransportElement *element;
    IntegrationRule *iRule;
    GaussPoint *gp;
    FloatArray vecTemperature;
    TransportMaterial *mat;




    for ( ielem = 1; ielem <= nelem; ielem++ ) {
        element = ( TransportElement * ) domain->giveElement(ielem);
        mat = ( TransportMaterial * ) element->giveMaterial();
        if ( mat->giveClassID() == CemhydMatClass ) {
            iRule = element->giveDefaultIntegrationRulePtr();
            for ( i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
                gp  = iRule->getIntegrationPoint(i);
                dV  = element->computeVolumeAround(gp);
                element->giveIPValue(vecTemperature, gp, IST_Temperature, tStep);
                //mat->IP_volume += dV;
                //mat->average_temp += vecState.at(1) * dV;
            }
        }
    }

    for ( i = 1; i <= domain->giveNumberOfMaterialModels(); i++ ) {
        mat = ( TransportMaterial * ) domain->giveMaterial(i);
        if ( mat->giveClassID() == CemhydMatClass ) {
            //mat->average_temp /= mat->IP_volume;
        }
    }
}

void StructuralMaterialEvaluator :: doStepOutput(TimeStep *tStep)
{
    FloatArray outputValue;
    Domain *d = this->giveDomain(1);
    if ( tStep->isTheFirstStep() ) {
        this->outfile << "# Time";
        for ( int var : this->vars ) {
            this->outfile << ", " << __InternalStateTypeToString( ( InternalStateType ) var );
        }

        this->outfile << '\n';
    }

    outfile << tStep->giveIntrinsicTime();
    for ( int i = 1; i <= d->giveNumberOfMaterialModels(); i++ ) {
        Material *mat = d->giveMaterial(i);
        for ( int var : this->vars ) {
            mat->giveIPValue(outputValue, gps[i-1].get(), ( InternalStateType ) var, tStep);
            outfile << " " << outputValue;
        }
    }

    outfile << std :: endl;
}

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

#ifdef VERBOSE
    OOFEM_LOG_INFO("Applying initial conditions\n");
#endif
    int nDofs, j, k, jj;
    int nman  = domain->giveNumberOfDofManagers();
    DofManager *node;
    Dof *iDof;

    UnknownsField->advanceSolution(stepWhenIcApply);
    solutionVector = UnknownsField->giveSolutionVector(stepWhenIcApply);
    solutionVector->resize(neq);
    solutionVector->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)
            iDof  =  node->giveDof(k);
            if ( !iDof->isPrimaryDof() ) {
                continue;
            }

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

    int nelem = domain->giveNumberOfElements();
    
    //project initial temperature to integration points

//     for ( 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
    TransportElement *element;
    CemhydMat *cem;
    for ( j = 1; j <= nelem; j++ ) {
        element = ( TransportElement * ) domain->giveElement(j);
        //assign status to each integration point on each element
        if ( element->giveMaterial()->giveClassID() == CemhydMatClass ) {
            element->giveMaterial()->initMaterial(element); //create microstructures and statuses on specific GPs
            element->updateInternalState(stepWhenIcApply);   //store temporary unequilibrated temperature
            element->updateYourself(stepWhenIcApply);   //store equilibrated temperature
            cem = ( CemhydMat * ) element->giveMaterial();
            cem->clearWeightTemperatureProductVolume(element);
            cem->storeWeightTemperatureProductVolume(element, stepWhenIcApply);
        }
    }

    //perform averaging on each material instance of CemhydMatClass
    int nmat = domain->giveNumberOfMaterialModels();
    for ( j = 1; j <= nmat; j++ ) {
        if ( domain->giveMaterial(j)->giveClassID() == CemhydMatClass ) {
            cem = ( CemhydMat * ) domain->giveMaterial(j);
            cem->averageTemperature();
        }
    }
#endif //__CEMHYD_MODULE
}

void StructuralMaterialEvaluator :: solveYourself()
{
    Domain *d = this->giveDomain(1);

    MaterialMode mode = _3dMat;
    FloatArray initialStrain(6);
    gps.clear();
    gps.reserve(d->giveNumberOfMaterialModels());
    for ( int i = 1; i <= d->giveNumberOfMaterialModels(); i++ ) {
        std :: unique_ptr< GaussPoint > gp = std::make_unique<GaussPoint>(nullptr, i, FloatArray(0), 1, mode);
        gps.emplace_back( std :: move(gp) );
        // Initialize the strain vector;
        StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( d->giveMaterial(i)->giveStatus( gps[i-1].get() ) );
        status->letStrainVectorBe(initialStrain);
    }

    std :: string outname = this->giveOutputBaseFileName() + ".matdata";
    this->outfile.open( outname.c_str() );

    this->timer.startTimer(EngngModelTimer :: EMTT_AnalysisTimer);

    TimeStep *tStep = giveNextStep();

    // Note, strain == strain-rate (kept as strain for brevity)
    int maxiter = 100; // User input?
    FloatArray stressC, deltaStrain, strain, stress, res;
    stressC.resize( sControl.giveSize() );
    res.resize( sControl.giveSize() );

    FloatMatrix tangent, reducedTangent;
    for ( int istep = 1; istep <= this->numberOfSteps; ++istep ) {
        this->timer.startTimer(EngngModelTimer :: EMTT_SolutionStepTimer);
        for ( int imat = 1; imat <= d->giveNumberOfMaterialModels(); ++imat ) {
            GaussPoint *gp = gps[imat-1].get();
            StructuralMaterial *mat = static_cast< StructuralMaterial * >( d->giveMaterial(imat) );
            StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );

            strain = status->giveStrainVector();
            // Update the controlled parts
            for ( int j = 1; j <= eControl.giveSize(); ++j ) {
                int p = eControl.at(j);
                strain.at(p) = d->giveFunction( cmpntFunctions.at(p) )->evaluateAtTime( tStep->giveIntrinsicTime() );
            }

            for ( int j = 1; j <= sControl.giveSize(); ++j ) {
                int p = sControl.at(j);
                stressC.at(j) = d->giveFunction( cmpntFunctions.at(p) )->evaluateAtTime( tStep->giveIntrinsicTime() );
            }

            //strain.add(-100, {6.27e-06,  6.27e-06, 6.27e-06, 0, 0, 0});
            for ( int iter = 1; iter < maxiter; iter++ ) {
#if 0
                // Debugging:
                mat->give3dMaterialStiffnessMatrix(tangent, TangentStiffness, gp, tStep);
                tangent.printYourself("# tangent");
                
                strain.zero();
                mat->giveRealStressVector_3d(stress, gp, strain, tStep);
                FloatArray strain2;
                tangent.solveForRhs(stress, strain2);
                strain2.printYourself("# thermal expansion");
                break;
#endif

                strain.printYourself("Macro strain guess");
                mat->giveRealStressVector_3d(stress, gp, strain, tStep);
                for ( int j = 1; j <= sControl.giveSize(); ++j ) {
                    res.at(j) = stressC.at(j) - stress.at( sControl.at(j) );
                }

                OOFEM_LOG_INFO("*** Time step: %d (t = %.2e), Material %d, Iteration: %d,  Residual = %e (tolerance %.2e)\n", 
                              istep, tStep->giveIntrinsicTime(), imat, iter, res.computeNorm(), tolerance);

                if ( res.computeNorm() <= tolerance ) {
                    break;
                } else {
                    if ( tangent.giveNumberOfRows() == 0 || !keepTangent ) {
                        mat->give3dMaterialStiffnessMatrix(tangent, TangentStiffness, gp, tStep);
                    }

                    // Pick out the stress-controlled part;
                    reducedTangent.beSubMatrixOf(tangent, sControl, sControl);

                    // Update stress-controlled part of the strain
                    reducedTangent.solveForRhs(res, deltaStrain);
                    //deltaStrain.printYourself("deltaStrain");
                    for ( int j = 1; j <= sControl.giveSize(); ++j ) {
                        strain.at( sControl.at(j) ) += deltaStrain.at(j);
                    }
                }
            }

            if ( res.computeNorm() > tolerance ) {
                OOFEM_WARNING("Residual did not converge!");
            }

            // This material model has converged, so we update it and go on to the next.
            gp->updateYourself(tStep);
        }

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