C++ StructuralMaterialStatus类代码示例

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 SimpleVitrificationMaterial :: giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp,
                                                            const FloatArray &reducedStrain, TimeStep *tStep)
{
    FloatArray strainVector;
    FloatMatrix d;
    FloatArray deltaStrain;

    StructuralMaterialStatus *status = dynamic_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );

    this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);

    deltaStrain.beDifferenceOf( strainVector, status->giveStrainVector() );

    this->give3dMaterialStiffnessMatrix(d, TangentStiffness, gp, tStep);

    FloatArray deltaStress;
    deltaStress.beProductOf(d, deltaStrain);

    answer = status->giveStressVector();
    answer.add(deltaStress);

    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
}

void
PhaseFieldElement :: computeStiffnessMatrix_ud(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    FloatMatrix B, DB, N_d, DN, S(3,1);
    FloatArray stress;
    NLStructuralElement *el = this->giveElement();
    StructuralCrossSection *cs = dynamic_cast<StructuralCrossSection* > (el->giveCrossSection() );

    answer.clear();

    for ( auto &gp: el->giveIntegrationRule(0) ) {
        StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

        double dV = el->computeVolumeAround(gp);
        // compute int_V ( B^t * D_B * B )dV

        this->computeNd_matrixAt(*gp->giveNaturalCoordinates(), N_d);

        // stress   
        FloatArray reducedStrain, a_u;
        FloatMatrix B_u;
        el->computeBmatrixAt( gp, B_u );
        stress = matStat->giveTempStressVector();
        stress.times( this->computeGPrim( gp, VM_Total, tStep ) );
        
        S.setColumn(stress,1);
        DN.beProductOf(S, N_d);
        answer.plusProductUnsym(B_u, DN, dV);
    }

}

void
HyperElasticMaterial :: giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &totalStrain, TimeStep *tStep)
{
    double J2;
    FloatMatrix C(3, 3);
    FloatMatrix invC;
    FloatArray strainVector;

    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    this->giveStressDependentPartOfStrainVector_3d(strainVector, gp,
                                                totalStrain,
                                                tStep, VM_Total);

    C.at(1, 1) = 1. + 2. * strainVector.at(1);
    C.at(2, 2) = 1. + 2. * strainVector.at(2);
    C.at(3, 3) = 1. + 2. * strainVector.at(3);
    C.at(1, 2) = C.at(2, 1) = strainVector.at(6);
    C.at(1, 3) = C.at(3, 1) = strainVector.at(5);
    C.at(2, 3) = C.at(3, 2) = strainVector.at(4);
    invC.beInverseOf(C);
    J2 = C.giveDeterminant();

    answer.resize(6);
    double aux = ( K - 2. / 3. * G ) * ( J2 - 1. ) / 2. - G;
    answer.at(1) = aux * invC.at(1, 1) + G;
    answer.at(2) = aux * invC.at(2, 2) + G;
    answer.at(3) = aux * invC.at(3, 3) + G;
    answer.at(4) = aux * invC.at(2, 3);
    answer.at(5) = aux * invC.at(1, 3);
    answer.at(6) = aux * invC.at(1, 2);

    // update gp
    status->letTempStrainVectorBe(totalStrain);
    status->letTempStressVectorBe(answer);
}

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

double
PhaseFieldElement :: computeFreeEnergy(GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );
    FloatArray strain, stress;
    stress = matStat->giveTempStressVector();
    strain = matStat->giveTempStrainVector();
    return 0.5 * stress.dotProduct( strain );
}

void
CCTPlate3d :: giveCharacteristicTensor(FloatMatrix &answer, CharTensor type, GaussPoint *gp, TimeStep *tStep)
// returns characteristic tensor of the receiver at given gp and tStep
// strain vector = (Kappa_x, Kappa_y, Kappa_xy, Gamma_zx, Gamma_zy)
{
    FloatArray charVect;
    Material *mat = this->giveMaterial();
    StructuralMaterialStatus *ms = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );

    answer.resize(3, 3);
    answer.zero();

    if ( ( type == LocalForceTensor ) || ( type == GlobalForceTensor ) ) {
        //this->computeStressVector(charVect, gp, tStep);
        charVect = ms->giveStressVector();

        answer.at(1, 3) = charVect.at(4);
        answer.at(3, 1) = charVect.at(4);
        answer.at(2, 3) = charVect.at(5);
        answer.at(3, 2) = charVect.at(5);
    } else if ( ( type == LocalMomentumTensor ) || ( type == GlobalMomentumTensor ) ) {
        //this->computeStressVector(charVect, gp, tStep);
        charVect = ms->giveStressVector();

        answer.at(1, 1) = charVect.at(1);
        answer.at(2, 2) = charVect.at(2);
        answer.at(1, 2) = charVect.at(3);
        answer.at(2, 1) = charVect.at(3);
    } else if ( ( type == LocalStrainTensor ) || ( type == GlobalStrainTensor ) ) {
        //this->computeStrainVector(charVect, gp, tStep);
        charVect = ms->giveStrainVector();

        answer.at(1, 3) = charVect.at(4) / 2.;
        answer.at(3, 1) = charVect.at(4) / 2.;
        answer.at(2, 3) = charVect.at(5) / 2.;
        answer.at(3, 2) = charVect.at(5) / 2.;
    } else if ( ( type == LocalCurvatureTensor ) || ( type == GlobalCurvatureTensor ) ) {
        //this->computeStrainVector(charVect, gp, tStep);
        charVect = ms->giveStrainVector();

        answer.at(1, 1) = charVect.at(1);
        answer.at(2, 2) = charVect.at(2);
        answer.at(1, 2) = charVect.at(3) / 2.;
        answer.at(2, 1) = charVect.at(3) / 2.;
    } else {
        _error("GiveCharacteristicTensor: unsupported tensor mode");
        exit(1);
    }

    if ( ( type == GlobalForceTensor  ) || ( type == GlobalMomentumTensor  ) ||
        ( type == GlobalStrainTensor ) || ( type == GlobalCurvatureTensor ) ) {
        this->computeGtoLRotationMatrix();
        answer.rotatedWith(* GtoLRotationMatrix);
    }
}

void 
WinklerPasternakMaterial::giveRealStressVector_2dPlateSubSoil(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedE, TimeStep *tStep)
{
  FloatMatrix tangent;
  this->give2dPlateSubSoilStiffMtrx(tangent, ElasticStiffness, gp, tStep);
  answer.beProductOf(tangent, reducedE);

  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
  status->letTempStrainVectorBe(reducedE);
  status->letTempStressVectorBe(answer);
}

void
HyperElasticMaterial :: give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode, GaussPoint *gp, TimeStep *tStep)
{
    double J2, c11, c22, c33, c12, c13, c23, A, B;
    FloatMatrix C(3, 3);
    FloatMatrix invC;

    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );

    C.at(1, 1) = 1. + 2. * status->giveTempStrainVector().at(1);
    C.at(2, 2) = 1. + 2. * status->giveTempStrainVector().at(2);
    C.at(3, 3) = 1. + 2. * status->giveTempStrainVector().at(3);
    C.at(2, 3) = C.at(3, 2) = status->giveTempStrainVector().at(4);
    C.at(1, 3) = C.at(3, 1) = status->giveTempStrainVector().at(5);
    C.at(1, 2) = C.at(2, 1) = status->giveTempStrainVector().at(6);

    invC.beInverseOf(C);
    J2 = C.giveDeterminant();

    c11 = invC.at(1, 1);
    c22 = invC.at(2, 2);
    c33 = invC.at(3, 3);
    c12 = invC.at(1, 2);
    c13 = invC.at(1, 3);
    c23 = invC.at(2, 3);

    A = ( K - 2. / 3. * G ) * J2;
    B = -( K - 2. / 3. * G ) * ( J2 - 1. ) + 2. * G;

    answer.resize(6, 6);

    answer.at(1, 1) = ( A + B ) * c11 * c11;
    answer.at(2, 2) = ( A + B ) * c22 * c22;
    answer.at(3, 3) = ( A + B ) * c33 * c33;
    answer.at(4, 4) = A * c23 * c23 + B / 2. * ( c22 * c33 + c23 * c23 );
    answer.at(5, 5) = A * c13 * c13 + B / 2. * ( c11 * c33 + c13 * c13 );
    answer.at(6, 6) = A * c12 * c12 + B / 2. * ( c11 * c22 + c12 * c12 );
    answer.at(1, 2) = answer.at(2, 1) = A * c11 * c22 + B * c12 * c12;
    answer.at(1, 3) = answer.at(3, 1) = A * c11 * c33 + B * c13 * c13;
    answer.at(1, 4) = answer.at(4, 1) = A * c11 * c23 + B * c12 * c13;
    answer.at(1, 5) = answer.at(5, 1) = A * c11 * c13 + B * c11 * c13;
    answer.at(1, 6) = answer.at(6, 1) = A * c11 * c12 + B * c11 * c12;
    answer.at(2, 3) = answer.at(3, 2) = A * c22 * c33 + B * c23 * c23;
    answer.at(2, 4) = answer.at(4, 2) = A * c22 * c23 + B * c22 * c23;
    answer.at(2, 5) = answer.at(5, 2) = A * c22 * c13 + B * c12 * c23;
    answer.at(2, 6) = answer.at(6, 2) = A * c22 * c12 + B * c22 * c12;
    answer.at(3, 4) = answer.at(4, 3) = A * c33 * c23 + B * c33 * c23;
    answer.at(3, 5) = answer.at(5, 3) = A * c33 * c13 + B * c33 * c13;
    answer.at(3, 6) = answer.at(6, 3) = A * c33 * c12 + B * c13 * c23;
    answer.at(4, 5) = answer.at(5, 4) = A * c23 * c13 + B / 2. * ( c12 * c33 + c13 * c23 );
    answer.at(4, 6) = answer.at(6, 4) = A * c23 * c12 + B / 2. * ( c12 * c23 + c22 * c13 );
    answer.at(5, 6) = answer.at(6, 5) = A * c13 * c12 + B / 2. * ( c11 * c23 + c12 * c13 );
}

void
FiberedCrossSection :: giveGeneralizedStress_Beam3d(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
    double fiberThick, fiberWidth, fiberZCoord, fiberYCoord;
    FloatArray fiberStrain, reducedFiberStress;
    StructuralElement *element = static_cast< StructuralElement * >( gp->giveElement() );
    FiberedCrossSectionInterface *interface;

    if ( ( interface = static_cast< FiberedCrossSectionInterface * >( element->giveInterface(FiberedCrossSectionInterfaceType) ) ) == NULL ) {
        OOFEM_ERROR("element with no fiber support encountered");
    }

    answer.resize(6);
    answer.zero();

    for ( int i = 1; i <= numberOfFibers; i++ ) {
        GaussPoint *fiberGp = this->giveSlaveGaussPoint(gp, i - 1);
        StructuralMaterial *fiberMat = static_cast< StructuralMaterial * >( domain->giveMaterial( fiberMaterials.at(i) ) );
        // the question is whether this function should exist ?
        // if yes the element details will be hidden.
        // good idea also should be existence of element::GiveBmatrixOfLayer
        // and computing strains here - but first idea looks better
        // but treating of geometric non-linearities may become more complicated
        // another approach - use several functions with assumed kinematic constraints

        // resolve current layer z-coordinate
        fiberThick  = this->fiberThicks.at(i);
        fiberWidth  = this->fiberWidths.at(i);
        fiberYCoord = fiberGp->giveNaturalCoordinate(1);
        fiberZCoord = fiberGp->giveNaturalCoordinate(2);

        interface->FiberedCrossSectionInterface_computeStrainVectorInFiber(fiberStrain, strain, fiberGp, tStep);

        fiberMat->giveRealStressVector_Fiber(reducedFiberStress, fiberGp, fiberStrain, tStep);

        // perform integration
        // 1) membrane terms N, Qz, Qy
        answer.at(1) += reducedFiberStress.at(1) * fiberWidth * fiberThick;
        answer.at(2) += reducedFiberStress.at(2) * fiberWidth * fiberThick;
        answer.at(3) += reducedFiberStress.at(3) * fiberWidth * fiberThick;
        // 2) bending terms mx, my, mxy
        answer.at(4) += ( reducedFiberStress.at(2) * fiberWidth * fiberThick * fiberYCoord -
                          reducedFiberStress.at(3) * fiberWidth * fiberThick * fiberZCoord );
        answer.at(5) += reducedFiberStress.at(1) * fiberWidth * fiberThick * fiberZCoord;
        answer.at(6) -= reducedFiberStress.at(1) * fiberWidth * fiberThick * fiberYCoord;
    }

    // now we must update master gp ///@ todo simply chosen the first fiber material as master material /JB
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >
                                       ( domain->giveMaterial( fiberMaterials.at(1) )->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}

void
StructuralMaterialSettable :: giveRealStressVector_3d(FloatArray &answer,
                                                  GaussPoint *gp,
                                                  const FloatArray &totalStrain,
                                                  TimeStep *atTime)
{
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
    const FloatArray& stressVector = status->giveStressVector();

    status->letTempStrainVectorBe(totalStrain);
    status->letTempStressVectorBe(stressVector);
    answer = stressVector;
}

int
FiberedCrossSection :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    Material *mat = this->giveDomain()->giveMaterial( fiberMaterials.at(1) ); ///@todo For now, create material status according to the first fiber material
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );

    if ( type == IST_BeamForceMomentTensor ) {
        answer = status->giveStressVector();
        return 1;
    } else if ( type == IST_BeamStrainCurvatureTensor ) {
        answer = status->giveStrainVector();
        return 1;
    }
    return CrossSection :: giveIPValue(answer, gp, type, tStep);
}

void
SimpleCrossSection :: giveGeneralizedStress_MembraneRot(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
    FloatMatrix tangent;
    this->giveMembraneRotStiffMtrx(tangent, ElasticStiffness, gp, tStep);
    answer.beProductOf(tangent, strain);

    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveMaterial(gp)->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);


    ///@todo We should support nonlinear behavior for the membrane part. In fact, should be even bundle the rotation part with the membrane?
    /// We gain nothing from this design anyway as the rotation field is always separate. Separate manual integration by the element would be an option.
}

void
LIBeam3dNL :: computeTempCurv(FloatArray &answer, TimeStep *tStep)
{
    IntegrationRule *iRule = this->giveDefaultIntegrationRulePtr();
    GaussPoint *gp = iRule->getIntegrationPoint(0);

    FloatArray ui(3), ac(3), PrevEpsilon;
    FloatMatrix sc(3, 3), tmid(3, 3);

    // update curvature at midpoint
    // first, compute Tmid
    // ask increments
    this->computeVectorOf(VM_Incremental, tStep, ui);

    ac.at(1) = 0.5 * ( ui.at(10) - ui.at(4) );
    ac.at(2) = 0.5 * ( ui.at(11) - ui.at(5) );
    ac.at(3) = 0.5 * ( ui.at(12) - ui.at(6) );
    this->computeSMtrx(sc, ac);
    sc.times(1. / 2.);
    // compute I+sc
    sc.at(1, 1) += 1.0;
    sc.at(2, 2) += 1.0;
    sc.at(3, 3) += 1.0;
    tmid.beProductOf(sc, this->tc);

    // update curvature at centre
    ac.at(1) = ( ui.at(10) - ui.at(4) );
    ac.at(2) = ( ui.at(11) - ui.at(5) );
    ac.at(3) = ( ui.at(12) - ui.at(6) );

    answer.beTProductOf(tmid, ac);
    answer.times(1 / this->l0);

    // ask for previous kappa
    StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );
    if ( matStat ) {
        PrevEpsilon = matStat->giveStrainVector();
    }
    if ( PrevEpsilon.giveSize() ) {
        answer.at(1) += PrevEpsilon.at(4);
        answer.at(2) += PrevEpsilon.at(5);
        answer.at(3) += PrevEpsilon.at(6);
    }
}

void
LinQuad3DPlaneStress :: giveCharacteristicTensor(FloatMatrix &answer, CharTensor type, GaussPoint *gp, TimeStep *tStep)
// returns characteristic tensor of the receiver at given gp and tStep
// strain vector = (Eps_X, Eps_y, Gamma_xy, Kappa_z)
{
    FloatArray charVect;
    StructuralMaterialStatus *ms = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

    answer.resize(3, 3);
    answer.zero();

    if ( ( type == LocalForceTensor ) || ( type == GlobalForceTensor ) ) {
        //this->computeStressVector(charVect, gp, tStep);
        charVect = ms->giveStressVector();

        answer.at(1, 1) = charVect.at(1);
        answer.at(2, 2) = charVect.at(2);
        answer.at(1, 2) = charVect.at(3);
        answer.at(2, 1) = charVect.at(3);
    } else if ( ( type == LocalMomentTensor ) || ( type == GlobalMomentTensor ) ) {
    } else if ( ( type == LocalStrainTensor ) || ( type == GlobalStrainTensor ) ) {
        //this->computeStrainVector(charVect, gp, tStep);
        charVect = ms->giveStrainVector();

        answer.at(1, 1) = charVect.at(1);
        answer.at(2, 2) = charVect.at(2);
        answer.at(1, 2) = charVect.at(3) / 2.;
        answer.at(2, 1) = charVect.at(3) / 2.;
    } else if ( ( type == LocalCurvatureTensor ) || ( type == GlobalCurvatureTensor ) ) {
    } else {
        OOFEM_ERROR("unsupported tensor mode");
        exit(1);
    }

    if ( ( type == GlobalForceTensor  ) || ( type == GlobalMomentTensor  ) ||
        ( type == GlobalStrainTensor ) || ( type == GlobalCurvatureTensor ) ) {
        this->computeGtoLRotationMatrix();
        answer.rotatedWith(* GtoLRotationMatrix);
    }
}