C++ StructuralMaterialStatus::giveStressVector方法代码示例

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

//.........这里部分代码省略.........
                            globalCoord.zero();

                            for ( int i = 1; i <= gpEl->giveNumberOfDofManagers(); i++ ) {
                                DofManager *dMan = gpEl->giveDofManager(i);
                                globalCoord.at(1) += N.at(i) * dMan->giveCoordinate(1);
                                globalCoord.at(2) += N.at(i) * dMan->giveCoordinate(2);
                            }


                            ////////////////////////////////////////
                            // Compute weight of kernel function

                            FloatArray tipToGP;
                            tipToGP.beDifferenceOf(globalCoord, xT);
                            bool inFrontOfCrack = true;
                            if ( tipToGP.dotProduct(t) < 0.0 ) {
                                inFrontOfCrack = false;
                            }

                            double r = circPoints [ pointIndex ].distance(globalCoord);

                            if ( r < l && inFrontOfCrack ) {
                                double w = ( ( l - r ) / ( pow(2.0 * M_PI, 1.5) * pow(l, 3) ) ) * exp( -0.5 * pow(r, 2) / pow(l, 2) );

                                // Compute gp volume
                                double V = gpEl->computeVolumeAround(gp_i);

                                // Get stress
                                StructuralMaterialStatus *ms = dynamic_cast< StructuralMaterialStatus * >( gp_i->giveMaterialStatus() );
                                if ( ms == NULL ) {
                                    OOFEM_ERROR("failed to fetch MaterialStatus.");
                                }

                                FloatArray stressVecGP = ms->giveStressVector();

                                if ( sumQiWiVi.giveSize() != stressVecGP.giveSize() ) {
                                    sumQiWiVi.resize( stressVecGP.giveSize() );
                                    sumQiWiVi.zero();
                                }

                                // Add to numerator
                                sumQiWiVi.add(w * V, stressVecGP);

                                // Add to denominator
                                sumWiVi += w * V;
                            }
                        }
                    }


                    if ( fabs(sumWiVi) > 1.0e-12 ) {
                        stressVec.beScaled(1.0 / sumWiVi, sumQiWiVi);
                    } else {
                        // Take stress from closest Gauss point
                        int region = 1;
                        bool useCZGP = false;
                        GaussPoint &gp = * ( localizer->giveClosestIP(circPoints [ pointIndex ], region, useCZGP) );


                        // Compute stresses
                        StructuralMaterialStatus *ms = dynamic_cast< StructuralMaterialStatus * >( gp.giveMaterialStatus() );
                        if ( ms == NULL ) {
                            OOFEM_ERROR("failed to fetch MaterialStatus.");
                        }

                        stressVec = ms->giveStressVector();

void
NLStructuralElement :: giveInternalForcesVector_withIRulesAsSubcells(FloatArray &answer,
                                                                     TimeStep *tStep, int useUpdatedGpRecord)
{
    /*
     * Returns nodal representation of real internal forces computed from first Piola-Kirchoff stress
     * if useGpRecord == 1 then stresses stored in the gp are used, otherwise stresses are computed
     * this must be done if you want internal forces after element->updateYourself() has been called
     * for the same time step.
     * The integration procedure uses an integrationRulesArray for numerical integration.
     * Each integration rule is considered to represent a separate sub-cell/element. Typically this would be used when
     * integration of the element domain needs special treatment, e.g. when using the XFEM.
     */

    FloatMatrix B;
    FloatArray vStress, vStrain;

    IntArray irlocnum;
    FloatArray *m = & answer, temp;
    if ( this->giveInterpolation() && this->giveInterpolation()->hasSubPatchFormulation() ) {
        m = & temp;
    }

    // zero answer will resize accordingly when adding first contribution
    answer.clear();


    // loop over individual integration rules
    for ( auto &iRule: integrationRulesArray ) {
        for ( GaussPoint *gp: *iRule ) {
            StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

            if ( nlGeometry == 0 ) {
                this->computeBmatrixAt(gp, B);
                if ( useUpdatedGpRecord == 1 ) {
                    vStress = matStat->giveStressVector();
                } else {
                    this->computeStrainVector(vStrain, gp, tStep);
                    this->computeStressVector(vStress, vStrain, gp, tStep);
                }
            } else if ( nlGeometry == 1 ) {
                if ( this->domain->giveEngngModel()->giveFormulation() == AL ) { // Cauchy stress
                    if ( useUpdatedGpRecord == 1 ) {
                        vStress = matStat->giveCVector();
                    } else {
                        this->computeCauchyStressVector(vStress, gp, tStep);
                    }

                    this->computeBmatrixAt(gp, B);
                } else { // First Piola-Kirchhoff stress
                    if ( useUpdatedGpRecord == 1 ) {
                        vStress = matStat->givePVector();
                    } else {
                        this->computeFirstPKStressVector(vStress, gp, tStep);
                    }

                    this->computeBHmatrixAt(gp, B);
                }
            }

            if ( vStress.giveSize() == 0 ) { //@todo is this really necessary?
                break;
            }

            // compute nodal representation of internal forces at nodes as f = B^T*stress dV
            double dV = this->computeVolumeAround(gp);
            m->plusProduct(B, vStress, dV);

            // localize irule contribution into element matrix
            if ( this->giveIntegrationRuleLocalCodeNumbers(irlocnum, *iRule) ) {
                answer.assemble(* m, irlocnum);
                m->clear();
            }
        }
    }

    // if inactive: update fields but do not give any contribution to the structure
    if ( !this->isActivated(tStep) ) {
        answer.zero();
        return;
    }
}

void
NLStructuralElement :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    FloatMatrix B;
    FloatArray vStress, vStrain, u;

    // This function can be quite costly to do inside the loops when one has many slave dofs.
    this->computeVectorOf(VM_Total, tStep, u);
    // subtract initial displacements, if defined
    if ( initialDisplacements ) {
        u.subtract(* initialDisplacements);
    }

    // zero answer will resize accordingly when adding first contribution
    answer.clear();

    for ( auto &gp: *this->giveDefaultIntegrationRulePtr() ) {
        StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

        // Engineering (small strain) stress
        if ( nlGeometry == 0 ) {
            this->computeBmatrixAt(gp, B);
            if ( useUpdatedGpRecord == 1 ) {
                vStress = matStat->giveStressVector();
            } else {
                ///@todo Is this really what we should do for inactive elements?
                if ( !this->isActivated(tStep) ) {
                    vStrain.resize( StructuralMaterial :: giveSizeOfVoigtSymVector( gp->giveMaterialMode() ) );
                    vStrain.zero();
                }
                vStrain.beProductOf(B, u);
                this->computeStressVector(vStress, vStrain, gp, tStep);
            }
        } else if ( nlGeometry == 1 ) {  // First Piola-Kirchhoff stress
            if ( this->domain->giveEngngModel()->giveFormulation() == AL ) { // Cauchy stress
                if ( useUpdatedGpRecord == 1 ) {
                    vStress = matStat->giveCVector();
                } else {
                    this->computeCauchyStressVector(vStress, gp, tStep);
                }

                this->computeBmatrixAt(gp, B);
            } else { // First Piola-Kirchhoff stress
                if ( useUpdatedGpRecord == 1 ) {
                    vStress = matStat->givePVector();
                } else {
                    this->computeFirstPKStressVector(vStress, gp, tStep);
                    ///@todo This is actaully inefficient since it constructs B and twice and collects the nodal unknowns over and over.
                }

                this->computeBHmatrixAt(gp, B);
            }
        }

        if ( vStress.giveSize() == 0 ) { /// @todo is this check really necessary?
            break;
        }

        // Compute nodal internal forces at nodes as f = B^T*Stress dV
        double dV  = this->computeVolumeAround(gp);

        if ( nlGeometry == 1 ) {  // First Piola-Kirchhoff stress
            if ( vStress.giveSize() == 9 ) {
                FloatArray stressTemp;
                StructuralMaterial :: giveReducedVectorForm( stressTemp, vStress, gp->giveMaterialMode() );
                answer.plusProduct(B, stressTemp, dV);
            } else   {
                answer.plusProduct(B, vStress, dV);
            }
        } else   {
            if ( vStress.giveSize() == 6 ) {
                // It may happen that e.g. plane strain is computed
                // using the default 3D implementation. If so,
                // the stress needs to be reduced.
                // (Note that no reduction will take place if
                //  the simulation is actually 3D.)
                FloatArray stressTemp;
                StructuralMaterial :: giveReducedSymVectorForm( stressTemp, vStress, gp->giveMaterialMode() );
                answer.plusProduct(B, stressTemp, dV);
            } else   {
                answer.plusProduct(B, vStress, dV);
            }
        }
    }

    // If inactive: update fields but do not give any contribution to the internal forces
    if ( !this->isActivated(tStep) ) {
        answer.zero();
        return;
    }
}