C++ StructuralInterfaceMaterialStatus类代码示例

void
StructuralInterfaceMaterial :: giveStiffnessMatrix_dTdj_Num(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    // Default implementation for computation of the numerical tangent
    // Computes the material stiffness using a central difference method

    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus(gp) );
    if ( status ) {
        FloatMatrix F;
        F = status->giveTempF();
        int dim = F.giveNumberOfRows();
        answer.resize(dim, dim);
        answer.zero();
        const double eps = 1.0e-9;
        FloatArray T, TPlus, TMinus;

        FloatArray jump, jumpPlus, jumpMinus, Kcolumn;
        jump = status->giveTempJump();
        for ( int i = 1; i <= dim; i++ ) {
            jumpPlus = jumpMinus = jump;
            jumpPlus.at(i)  += eps;
            jumpMinus.at(i) -= eps;
            this->giveFirstPKTraction_3d(TPlus, gp, jumpPlus, F, tStep);
            this->giveFirstPKTraction_3d(TMinus, gp, jumpMinus, F, tStep);

            Kcolumn = ( TPlus - TMinus );
            answer.setColumn(Kcolumn, i);
        }
        answer.times( 1.0 / ( 2 * eps ) );
        this->giveFirstPKTraction_3d(T, gp, jump, F, tStep); // reset temp values by recomputing the stress
    }
}

void
StructuralInterfaceMaterial :: give3dStiffnessMatrix_Eng_Num( FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep )
{
    // Default implementation for computation of the numerical tangent d(sig)/d(jump)
    // Computes the material stiffness using a central difference method

    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus( gp ) );
    if(status) {
        const double eps = 1.0e-9;
        FloatArray t, tPlus, tMinus;
        FloatArray jump, jumpPlus, jumpMinus, Kcolumn;
        jump = status->giveTempJump( );
        int dim = jump.giveSize( );
        answer.resize( dim, dim );
        answer.zero( );

        for(int i = 1; i <= dim; i++) {
            jumpPlus = jumpMinus = jump;
            jumpPlus.at( i ) += eps;
            jumpMinus.at( i ) -= eps;
            this->giveEngTraction_3d( tPlus, gp, jumpPlus, tStep );
            this->giveEngTraction_3d( tMinus, gp, jumpMinus, tStep );

            Kcolumn = ( tPlus - tMinus );
            answer.setColumn( Kcolumn, i );
        }
        answer.times( 1.0 / ( 2 * eps ) );
        this->giveEngTraction_3d( t, gp, jump, tStep ); // reset temp values by recomputing the stress
    }
}

void Crack :: AppendCohesiveZoneGaussPoint(GaussPoint *ipGP)
{
    StructuralInterfaceMaterialStatus *matStat = dynamic_cast< StructuralInterfaceMaterialStatus * >( ipGP->giveMaterialStatus() );
    matStat->printYourself();
    if ( matStat != NULL ) {
        // Compute arc length position of the Gauss point
        const FloatArray &coord =  ipGP->giveGlobalCoordinates();
        double tangDist = 0.0, arcPos = 0.0;
        mpEnrichmentDomain->computeTangentialSignDist(tangDist, coord, arcPos);

        // Insert at correct position
        std :: vector< GaussPoint * > :: iterator iteratorGP = mCohesiveZoneGaussPoints.begin();
        std :: vector< double > :: iterator iteratorPos = mCohesiveZoneArcPositions.begin();
        for ( size_t i = 0; i < mCohesiveZoneArcPositions.size(); i++ ) {
            if ( arcPos > mCohesiveZoneArcPositions [ i ] ) {
                iteratorGP++;
                iteratorPos++;
            }
        }

        mCohesiveZoneGaussPoints.insert(iteratorGP, ipGP);
        mCohesiveZoneArcPositions.insert(iteratorPos, arcPos);
    } else   {
        OOFEM_ERROR("matStat == NULL.")
    }
}

void
CohesiveInterfaceMaterial :: give3dStiffnessMatrix_Eng(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    answer.resize(3, 3);
    answer.zero();
    
    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus(gp) );
//     double normalJump = status->giveTempJump().at(1);
    
//     if (normalJump > 0.) {
//         if(normalJump<transitionOpening){ // reduce traction in tension
//             double stiffTmp = kn*stiffCoeffKn + (kn - kn*stiffCoeffKn) * (1. - normalJump/transitionOpening);
//             answer.at(1, 1) = stiffTmp;
//         } else {
//             answer.at(1, 1) = kn * stiffCoeffKn;
//         }
//     } else {
//         // standard part of elastic stress-strain law
//         answer.at(1, 1) = kn;
//     }
//     
//     if ( rMode == SecantStiffness || rMode == TangentStiffness ) {
//         if ( normalJump + transitionOpening <= 0. ) { //local CS
//             answer.at(1, 1) = kn; //compression
//         } else {
//             answer.at(1, 1) = 0*kn*stiffCoeffKn; //tension
//         }
//     } else {
//         answer.at(1, 1) = kn; 
//     }
    
    
    double x = status->giveTempJump().at(1) + transitionOpening;
    
    if (stiffCoeffKn == 1.){//tension stiffness = compression stiffness
        answer.at(1,1) = kn;
    } else {
        //TangentStiffness by derivating traction with regards to x (=relative displacement)
        answer.at(1,1) = (M_PI/2. + atan(smoothMag*x))/M_PI*kn*stiffCoeffKn + (M_PI/2.-atan(smoothMag*x))/M_PI*kn + smoothMag*kn*stiffCoeffKn*x/M_PI/(smoothMag*smoothMag*x*x+1) - smoothMag*kn*x/M_PI/(smoothMag*smoothMag*x*x+1);
    }
    
    //answer.at(1, 1) = kn; 
    answer.at(2, 2) = ks;
    answer.at(3, 3) = ks;
}

int
StructuralInterfaceMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus(gp) );
    if ( type == IST_InterfaceJump ) {
        answer = status->giveJump();
        answer.resizeWithValues(3); // In case some model is not storing all components.
        return 1;
    } else if ( type == IST_InterfaceTraction ) {
        answer = status->giveTraction();
        answer.resizeWithValues(3);
        return 1;
    } else if ( type == IST_InterfaceFirstPKTraction ) {
        answer = status->giveFirstPKTraction();
        answer = status->giveTempFirstPKTraction();
        answer.resizeWithValues(3);
        return 1;
    } else if ( type == IST_DeformationGradientTensor ) {
        answer.beVectorForm( status->giveF() );
        return 1;
    } else if ( type == IST_InterfaceNormal ) {
        answer = status->giveNormal();
        return 1;
    } else {
        return Material :: giveIPValue(answer, gp, type, tStep);
    }
}

void
CohesiveInterfaceMaterial :: giveEngTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
{
    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus(gp) );

    answer.resize(3);

    double x = jump.at(1) + transitionOpening;
    
    if (stiffCoeffKn == 1.){//tension stiffness = compression stiffness
        answer.at(1) = kn*x;
    } else {
        //composed function from two atan's to have smooth intersection between tension and compression
        answer.at(1) = (M_PI/2. + atan(smoothMag*x))/M_PI*kn*stiffCoeffKn*x + (M_PI/2.-atan(smoothMag*x))/M_PI*kn*x;
    }
    
    // shear part of elastic stress-strain law
    answer.at(2) = ks * jump.at(2);
    answer.at(3) = ks * jump.at(3);
    
    // update gp
    status->letTempJumpBe(jump);
    status->letTempTractionBe(answer);
}

//.........这里部分代码省略.........
                        for ( size_t segIndex = 0; segIndex < numSeg; segIndex++ ) {
                            int czRuleNum = 1;
                            mpCZIntegrationRules.emplace_back( new GaussIntegrationRule(czRuleNum, element) );

                            // Add index of current ei
                            mCZEnrItemIndices.push_back(eiIndex);

                            // Add indices of other ei, that cause interaction through
                            // intersection enrichment fronts
                            mCZTouchingEnrItemIndices.push_back(touchingEiIndices);

                            // Compute crack normal
                            FloatArray crackTang;
                            crackTang.beDifferenceOf(crackPolygon [ segIndex + 1 ], crackPolygon [ segIndex ]);

                            if ( crackTang.computeSquaredNorm() > tol2 ) {
                                crackTang.normalize();
                            }

                            FloatArray crackNormal = {
                                -crackTang.at(2), crackTang.at(1)
                            };

                            mpCZIntegrationRules [ segIndex ]->SetUpPointsOn2DEmbeddedLine(mCSNumGaussPoints, matMode,
                                                                                           crackPolygon [ segIndex ], crackPolygon [ segIndex + 1 ]);

                            for ( GaussPoint *gp: *mpCZIntegrationRules [ segIndex ] ) {
                                double gw = gp->giveWeight();
                                double segLength = crackPolygon [ segIndex ].distance(crackPolygon [ segIndex + 1 ]);
                                gw *= 0.5 * segLength;
                                gp->setWeight(gw);

                                // Fetch material status and set normal
                                StructuralInterfaceMaterialStatus *ms = dynamic_cast< StructuralInterfaceMaterialStatus * >( mpCZMat->giveStatus(gp) );
                                if ( ms == NULL ) {
                                    OOFEM_ERROR("Failed to fetch material status.");
                                }

                                ms->letNormalBe(crackNormal);

                                // Give Gauss point reference to the enrichment item
                                // to simplify post processing.
                                crack->AppendCohesiveZoneGaussPoint(gp);
                            }
                        }
                    }



                    partitionSucceeded = true;
                }
            } // if(firstIntersection)
            else {
                // Loop over triangles
                std :: vector< Triangle >allTriCopy;
                for ( size_t triIndex = 0; triIndex < mSubTri.size(); triIndex++ ) {
                    // Call alternative version of XfemElementInterface_prepareNodesForDelaunay
                    std :: vector< std :: vector< FloatArray > >pointPartitionsTri;
                    double startXi, endXi;
                    bool intersection = false;
                    XfemElementInterface_prepareNodesForDelaunay(pointPartitionsTri, startXi, endXi, mSubTri [ triIndex ], eiIndex, intersection);

                    if ( intersection ) {
                        // Use XfemElementInterface_partitionElement to subdivide triangle j
                        for ( int i = 0; i < int ( pointPartitionsTri.size() ); i++ ) {
                            this->XfemElementInterface_partitionElement(allTriCopy, pointPartitionsTri [ i ]);

void GnuplotExportModule::outputXFEM(Crack &iCrack, TimeStep *tStep)
{
    const std::vector<GaussPoint*> &czGaussPoints = iCrack.giveCohesiveZoneGaussPoints();

    std::vector<double> arcLengthPositions, normalJumps, tangJumps, normalTractions;

    const BasicGeometry *bg = iCrack.giveGeometry();

    for( GaussPoint *gp: czGaussPoints ) {

        StructuralInterfaceMaterialStatus *matStat = dynamic_cast<StructuralInterfaceMaterialStatus*> ( gp->giveMaterialStatus() );
        if(matStat != NULL) {

            // Compute arc length position of the Gauss point
            const FloatArray &coord = (gp->giveGlobalCoordinates());
            double tangDist = 0.0, arcPos = 0.0;
            bg->computeTangentialSignDist(tangDist, coord, arcPos);
            arcLengthPositions.push_back(arcPos);

            // Compute displacement jump in normal and tangential direction
            // Local numbering: (tang_z, tang, normal)
            const FloatArray &jumpLoc = matStat->giveJump();

            double normalJump = jumpLoc.at(3);
            normalJumps.push_back(normalJump);


            tangJumps.push_back( jumpLoc.at(2) );


            const FloatArray &trac = matStat->giveFirstPKTraction();
            normalTractions.push_back(trac.at(3));
        }
    }



    Domain *domain = emodel->giveDomain(1);
    XfemManager *xMan = domain->giveXfemManager();
    if ( xMan != NULL ) {
        double time = 0.0;

        TimeStep *ts = emodel->giveCurrentStep();
        if ( ts != NULL ) {
            time = ts->giveTargetTime();
        }

        int eiIndex = iCrack.giveNumber();

        std :: stringstream strNormalJump;
        strNormalJump << "NormalJumpGnuplotEI" << eiIndex << "Time" << time << ".dat";
        std :: string nameNormalJump = strNormalJump.str();
        XFEMDebugTools::WriteArrayToGnuplot(nameNormalJump, arcLengthPositions, normalJumps);

        std :: stringstream strTangJump;
        strTangJump << "TangJumpGnuplotEI" << eiIndex << "Time" << time << ".dat";
        std :: string nameTangJump = strTangJump.str();
        XFEMDebugTools::WriteArrayToGnuplot(nameTangJump, arcLengthPositions, tangJumps);

        std :: stringstream strNormalTrac;
        strNormalTrac << "NormalTracGnuplotEI" << eiIndex << "Time" << time << ".dat";
        std :: string nameNormalTrac = strNormalTrac.str();
        XFEMDebugTools::WriteArrayToGnuplot(nameNormalTrac, arcLengthPositions, normalTractions);


        std::vector<FloatArray> matForcesStart, matForcesEnd;
        std::vector<double> radii;

        // Material forces
        for(double matForceRadius : mMatForceRadii) {

            radii.push_back(matForceRadius);

            EnrichmentFront *efStart = iCrack.giveEnrichmentFrontStart();
            const TipInfo &tipInfoStart = efStart->giveTipInfo();

            FloatArray matForceStart;
            mpMatForceEvaluator->computeMaterialForce(matForceStart, *domain, tipInfoStart, tStep, matForceRadius);

            if(matForceStart.giveSize() > 0) {
                matForcesStart.push_back(matForceStart);
            }
            else {
                matForcesStart.push_back({0.0,0.0});
            }


            EnrichmentFront *efEnd = iCrack.giveEnrichmentFrontEnd();
            const TipInfo &tipInfoEnd = efEnd->giveTipInfo();

            FloatArray matForceEnd;
            mpMatForceEvaluator->computeMaterialForce(matForceEnd, *domain, tipInfoEnd, tStep, matForceRadius);

            if(matForceEnd.giveSize() > 0) {
                matForcesEnd.push_back(matForceEnd);
            }
            else {
                matForcesEnd.push_back({0.0,0.0});
            }

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