C++ StructuralMaterial::giveStatus方法代码示例

void
SimpleCrossSection :: giveGeneralizedStress_Shell(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
  /**Note: (by bp): This assumes that the behaviour is elastic
     there exist a nuumber of nonlinear integral material models for beams/plates/shells
     defined directly in terms of integral forces and moments and corresponding 
     deformations and curvatures. This would require to implement support at material model level.
     Mikael: See earlier response to comment
  */
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    FloatArray elasticStrain, et, e0;
    FloatMatrix tangent;
    elasticStrain = strain;
    this->giveTemperatureVector(et, gp, tStep);
    if ( et.giveSize() ) {
        double thick = this->give(CS_Thickness, gp);
        mat->giveThermalDilatationVector(e0, gp, tStep);
        elasticStrain.at(1) -= e0.at(1) * ( et.at(1) - mat->giveReferenceTemperature() );
        elasticStrain.at(2) -= e0.at(2) * ( et.at(1) - mat->giveReferenceTemperature() );
        if ( et.giveSize() > 1 ) {
            elasticStrain.at(4) -= e0.at(1) * et.at(2) / thick;     // kappa_x
            elasticStrain.at(5) -= e0.at(2) * et.at(2) / thick;     // kappa_y
        }
    }
    this->give3dShellStiffMtrx(tangent, ElasticStiffness, gp, tStep);
    answer.beProductOf(tangent, elasticStrain);
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}

void
SimpleCrossSection :: giveGeneralizedStress_Beam2d(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
  /**Note: (by bp): This assumes that the behaviour is elastic
     there exist a nuumber of nonlinear integral material models for beams defined directly in terms of integral forces and moments and corresponding 
     deformations and curvatures. This would require to implement support at material model level.
     Mikael: That would not be a continuum material model, but it would highly depend on the cross-section shape, thus, it should be a special cross-section model instead.
     This cross-section assumes you can split the response into inertia moments and pure material response. This is only possible for a constant, elastic response (i.e. elastic).
     Therefore, this cross-section may only be allowed to give the elastic response.
  */
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    FloatArray elasticStrain, et, e0;
    FloatMatrix tangent;
    elasticStrain = strain;
    this->giveTemperatureVector(et, gp, tStep);
    if ( et.giveSize() > 0 ) {
        mat->giveThermalDilatationVector(e0, gp, tStep);
        double thick = this->give(CS_Thickness, gp);
        elasticStrain.at(1) -= e0.at(1) * ( et.at(1) - mat->giveReferenceTemperature() );
        if ( et.giveSize() > 1 ) {
            elasticStrain.at(2) -= e0.at(1) * et.at(2) / thick;     // kappa_x
        }
    }
    this->give2dBeamStiffMtrx(tangent, ElasticStiffness, gp, tStep);
    answer.beProductOf(tangent, elasticStrain);
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}

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

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