C++ StructuralMaterial类代码示例

void
PerfectlyPlasticMaterial :: giveEffectiveMaterialStiffnessMatrix(FloatMatrix &answer,
                                                                 MatResponseMode mode,
                                                                 GaussPoint *gp,
                                                                 TimeStep *atTime)
//
//
// for case of perfectly plastic material
// computes full elastic constitutive matrix for case of gp stress-strain state.
// if strainIncrement == NULL a loading is assumed
//
// we follow terminology based on paper from R. de Borst:
// "Smeared Cracking, plasticity, creep - Unified Aproach"
//
// if derived material would like to implement failure behaviour
// it must redefine basic Give3dMaterialStiffnessMatrix function
// in order to take possible failure (tension cracking) into account
//
{
    // FloatMatrix *de; // elastic matrix respecting fracture or failure
    StructuralMaterial *lMat = static_cast< StructuralMaterial * >( this->giveLinearElasticMaterial() );

    if ( lMat->hasMaterialModeCapability( gp->giveMaterialMode() ) ) {
        FloatMatrix stiff;
        lMat->giveStiffnessMatrix(stiff, mode, gp, atTime);
        this->giveFullSymMatrixForm(answer, stiff, gp->giveMaterialMode());
    } else {
        OOFEM_ERROR("PerfectlyPlasticMaterial :: giveEffectiveMaterialStiffnessMatrix - unsupported material mode");
    }
}

void
SimpleCrossSection :: give3dShellStiffMtrx(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveMaterial(gp) );

    FloatMatrix mat3d;

    double thickness = this->give(CS_Thickness, gp);
    double thickness3 = thickness * thickness * thickness;

    mat->givePlaneStressStiffMtrx(mat3d, rMode, gp, tStep);

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

    for ( int i = 1; i <= 3; i++ ) {
        for ( int j = 1; j <= 3; j++ ) {
            answer.at(i, j) = mat3d.at(i, j) * thickness;
        }
    }
    for ( int i = 1; i <= 3; i++ ) {
        for ( int j = 1; j <= 3; j++ ) {
            answer.at(i + 3, j + 3) = mat3d.at(i, j) * thickness3 / 12.0;
        }
    }

    answer.at(8, 8) = answer.at(7, 7) = mat3d.at(3, 3) * thickness * ( 5. / 6. );
}

void
SimpleCrossSection :: give2dPlateSubSoilStiffMtrx(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterial *mat;
    mat = dynamic_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    mat->give2dPlateSubSoilStiffMtrx(answer, ElasticStiffness, gp, tStep);
}

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.
  */
    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 :: give3dBeamStiffMtrx(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveMaterial(gp) );

    FloatMatrix mat3d;
    double area, E, G, Iy, Iz, Ik;
    double shearAreay, shearAreaz;

    mat->give1dStressStiffMtrx(mat3d, rMode, gp, tStep);
    E    = mat3d.at(1, 1);
    G    = mat->give('G', gp);
    area = this->give(CS_Area, gp);
    Iy   = this->give(CS_InertiaMomentY, gp);
    Iz   = this->give(CS_InertiaMomentZ, gp);
    Ik   = this->give(CS_TorsionMomentX, gp);

    //shearCoeff = this->give(CS_BeamShearCoeff);
    shearAreay = this->give(CS_ShearAreaY, gp);
    shearAreaz = this->give(CS_ShearAreaZ, gp);

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

    answer.at(1, 1) = E * area;
    ///@todo Do this by using the general 3d tangent matrix instead somehow!!!
    answer.at(2, 2) = shearAreay * G;
    answer.at(3, 3) = shearAreaz * G;
    //answer.at(2, 2) = shearCoeff * G * area;
    //answer.at(3, 3) = shearCoeff * G * area;
    answer.at(4, 4) = G * Ik;
    answer.at(5, 5) = E * Iy;
    answer.at(6, 6) = E * Iz;
}

void
StructuralCrossSection :: giveRealStresses(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
    MaterialMode mode = gp->giveMaterialMode();
    if ( mode == _2dBeam ) {
        this->giveGeneralizedStress_Beam2d(answer, gp, strain, tStep);
    } else if ( mode == _3dBeam ) {
        this->giveGeneralizedStress_Beam3d(answer, gp, strain, tStep);
    } else if ( mode == _2dPlate ) {
        this->giveGeneralizedStress_Plate(answer, gp, strain, tStep);
    } else if ( mode == _3dShell ) {
        this->giveGeneralizedStress_Shell(answer, gp, strain, tStep);
    } else if ( mode == _3dMat ) {
        this->giveRealStress_3d(answer, gp, strain, tStep);
    } else if ( mode == _PlaneStrain ) {
        this->giveRealStress_PlaneStrain(answer, gp, strain, tStep);
    } else if ( mode == _PlaneStress ) {
        this->giveRealStress_PlaneStress(answer, gp, strain, tStep);
    } else if ( mode == _1dMat ) {
        this->giveRealStress_1d(answer, gp, strain, tStep);
    } else if ( mode == _Warping ) {
        this->giveRealStress_Warping(answer, gp, strain, tStep);
    } else {
        // This should never happen ?
        ///@todo this part only works for simple cross section and will be removed soon when new interface elements are done /JB
        StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveMaterial(gp) );
        if ( mat->hasMaterialModeCapability( gp->giveMaterialMode() ) ) {
            mat->giveRealStressVector(answer, gp, strain, tStep);
        } else {
            OOFEM_ERROR("unsupported mode");
        }
    }
}

void
StructuralCrossSection :: giveFullCharacteristicVector(FloatArray &answer,
                                                       GaussPoint *gp,
                                                       const FloatArray &strainVector)
//
// returns full 3d general strain vector from strainVector in reducedMode
// based on StressStrainMode in gp. Included are strains which
// perform nonzero work.
// General strain vector has one of the following forms:
// 1) strainVector3d {eps_x,eps_y,eps_z,gamma_yz,gamma_zx,gamma_xy}
// 2) strainVectorShell {eps_x,eps_y,gamma_xy, kappa_x, kappa_y, kappa_xy, gamma_zx, gamma_zy}
//
// you must assigng your stress strain mode to one of the folloving modes (or add new)
// FullForm of MaterialStiffnessMatrix must have the same form.
//
{
    MaterialMode mode = gp->giveMaterialMode();
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( gp->giveElement()->giveMaterial() );

    if ( ( mode == _3dMat ) || ( mode == _3dMicroplane ) ) {
        answer = strainVector;
        return;
    } else {
        mat->giveFullCharacteristicVector(answer, gp, strainVector);
    }
}

void
FiberedCrossSection :: giveCharMaterialStiffnessMatrix(FloatMatrix &answer,
                                                       MatResponseMode rMode,
                                                       GaussPoint *gp,
                                                       TimeStep *tStep)
{
    MaterialMode mode = gp->giveMaterialMode();
    if ( mode == _2dBeam ) {
        this->give2dBeamStiffMtrx(answer, rMode, gp, tStep);
    } else if ( mode == _3dBeam ) {
        this->give3dBeamStiffMtrx(answer, rMode, gp, tStep);
    } else if ( mode == _2dPlate ) {
        this->give2dPlateStiffMtrx(answer, rMode, gp, tStep);
    } else if ( mode == _3dShell ) {
        this->give3dShellStiffMtrx(answer, rMode, gp, tStep);
    } else {
        OOFEM_ERROR("Not implemented for bulk materials.");
        ///@todo What about the fibers?! Rather give just an error message if the fibers aren't supported than to just silently ignore them.
#if 0
        StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( gp->giveElement()->giveMaterial() );
        if ( mat->hasMaterialModeCapability( gp->giveMaterialMode() ) ) {
            mat->giveStiffnessMatrix(answer, rMode, gp, tStep);
        } else {
            OOFEM_ERROR("unsupported StressStrainMode");
        }
#endif
    }
}

contextIOResultType
FiberedCrossSection :: restoreIPContext(DataStream *stream, ContextMode mode, GaussPoint *masterGp)
//
// restores full material context (saves state variables, that completely describe
// current state)
//
// restores also slaves of master gp
//
{
    contextIOResultType iores;

    if ( ( iores = CrossSection :: restoreIPContext(stream, mode, masterGp) ) != CIO_OK ) {
        THROW_CIOERR(iores);                                                                   // saved masterGp
    }

    for ( int i = 1; i <= numberOfFibers; i++ ) {
        // creates also slaves if they don't exists
        GaussPoint *slaveGP = this->giveSlaveGaussPoint(masterGp, i - 1);
        StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( domain->giveMaterial( fiberMaterials.at(i) ) );
        if ( ( iores = mat->restoreIPContext(stream, mode, slaveGP) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
    }

    return CIO_OK;
}

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
FiberedCrossSection :: giveFiberMaterialStiffnessMatrix(FloatMatrix &fiberMatrix,
                                                        MatResponseMode rMode, GaussPoint *layerGp,
                                                        TimeStep *tStep)
{
    StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( domain->giveMaterial( fiberMaterials.at( layerGp->giveNumber() ) ) );
    mat->giveStiffnessMatrix(fiberMatrix, rMode, layerGp, tStep);
}

void
SimpleCrossSection :: giveReducedCharacteristicVector(FloatArray &answer, GaussPoint *gp,
                                                      const FloatArray &charVector3d)
//
// returns reduced stressVector or strainVector from full 3d vector reduced
// to vector required by gp->giveStressStrainMode()
//
{
    MaterialMode mode = gp->giveMaterialMode();
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( gp->giveElement()->giveMaterial() );
    IntArray indx;
    int size = charVector3d.giveSize();
    int i, j;

    if ( ( mode == _3dShell ) || ( mode == _3dBeam ) || ( mode == _2dPlate ) || ( mode == _2dBeam ) ) {
        if ( size != 12 ) {
            OOFEM_ERROR("SimpleCrossSection :: giveReducedCharacteristicVector - charVector3d size mismatch");
        }

        mat->giveStressStrainMask( indx, ReducedForm, gp->giveMaterialMode() );
        answer.resize( indx.giveSize() );
        answer.zero();

        for ( i = 1; i <= indx.giveSize(); i++ ) {
            if ( ( j = indx.at(i) ) ) {
                answer.at(i) = charVector3d.at(j);
            }
        }

        return;
    } else if ( mode == _3dBeam ) {
        if ( size != 6 ) {
            OOFEM_ERROR("SimpleCrossSection :: giveReducedCharacteristicVector - charVector3d size mismatch");
        }

        answer =  charVector3d;
        return;
    } else if ( mode == _PlaneStressRot ) {
        if ( size != 7 ) {
            OOFEM_ERROR("SimpleCrossSection :: giveReducedCharacteristicVector - charVector3d size mismatch");
        }

        mat->giveStressStrainMask( indx, ReducedForm, gp->giveMaterialMode() );
        answer.resize( indx.giveSize() );
        answer.zero();

        for ( i = 1; i <= indx.giveSize(); i++ ) {
            if ( ( j = indx.at(i) ) ) {
                answer.at(i) = charVector3d.at(j);
            }
        }

        return;
    } else {
        StructuralCrossSection :: giveReducedCharacteristicVector(answer, gp, charVector3d);
        return;
    }
}

void FiberedCrossSection :: createMaterialStatus(GaussPoint &iGP)
{
    for ( int i = 1; i <= numberOfFibers; i++ ) {
        GaussPoint *fiberGp = this->giveSlaveGaussPoint(& iGP, i - 1);
        StructuralMaterial *mat = static_cast< StructuralMaterial * >( domain->giveMaterial( fiberMaterials.at(i) ) );
        MaterialStatus *matStat = mat->CreateStatus(fiberGp);
        iGP.setMaterialStatus(matStat);
    }
}

void
SimpleCrossSection :: giveRealStress_3dDegeneratedShell(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
    StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    IntArray strainControl = {
        1, 2, 4, 5, 6
    };
    mat->giveRealStressVector_ShellStressControl(answer, gp, strain, strainControl, tStep);
}

void
SimpleCrossSection :: giveMembraneRotStiffMtrx(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    StructuralMaterial *mat;
    mat = dynamic_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    mat->givePlaneStressStiffMtrx(answer, ElasticStiffness, gp, tStep);
    answer.resizeWithData(4, 4);
    answer.at(4, 4) = this->give(CS_DrillingStiffness, gp);
}