C++ FloatArray::beDifferenceOf方法代码示例

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
TutorialMaterial :: computeSphDevPartOf(const FloatArray &sigV, FloatArray &sigSph, FloatArray &sigDev)
{
    double volumetricPart = ( sigV.at(1) + sigV.at(2) + sigV.at(3) ) / 3.0;
    sigSph = {volumetricPart, volumetricPart, volumetricPart, 0.0, 0.0, 0.0};
    sigDev.beDifferenceOf(sigV, sigSph);
}

void FE2FluidMaterial :: giveDeviatoricPressureStiffness(FloatArray &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    FE2FluidMaterialStatus *ms = static_cast<FE2FluidMaterialStatus*> (this->giveStatus(gp));
    ms->computeTangents(tStep);
    if ( mode == TangentStiffness ) {
        answer = ms->giveDeviatoricPressureTangent();
#ifdef DEBUG_TANGENT
        // Numerical ATS for debugging
        FloatArray strain(3); strain.zero();
        FloatArray sig, sigh;
        double epspvol, pressure = 0.0;
        double h = 1.00; // Linear problem, size of this doesn't matter.
        computeDeviatoricStressVector (sig, epspvol, gp, strain, pressure, tStep);
        computeDeviatoricStressVector (sigh, epspvol, gp, strain, pressure+h, tStep);

        FloatArray dsigh; dsigh.beDifferenceOf(sigh,sig); dsigh.times(1/h);

        printf("Analytical deviatoric pressure tangent = "); answer.printYourself();
        printf("Numerical deviatoric pressure tangent = "); dsigh.printYourself();
        dsigh.subtract(answer);
        double norm = dsigh.computeNorm();
        if (norm > answer.computeNorm()*DEBUG_ERR && norm > 0.0) {
            OOFEM_ERROR("Error in deviatoric pressure tangent");
        }
#endif
    } else {
        OOFEM_ERROR("Mode not implemented");
    }
}

double Line :: computeTangentialDistanceToEnd(FloatArray *point)
{
    FloatArray projection;
    this->computeProjection(projection);
    FloatArray tmp;
    tmp.beDifferenceOf(* point, mVertices [ 1 ]);
    return tmp.dotProduct(projection) / projection.computeNorm();
}

void PrescribedGradientBCWeak :: giveTractionElNormal(size_t iElInd, FloatArray &oNormal, FloatArray &oTangent) const
{
	FloatArray xS, xE;
	giveTractionElCoord(iElInd, xS, xE);

    oTangent.beDifferenceOf(xE, xS);
    oTangent.normalize();

    oNormal = {
        oTangent [ 1 ], -oTangent [ 0 ]
    };
}

void Line :: transformIntoPolar(FloatArray *point, FloatArray &answer)
{
    FloatArray xp;
    FloatMatrix Qt;
    FloatArray help;
    this->computeTransformationMatrix(Qt);
    help.beDifferenceOf(* point, mVertices [ 1 ]);
    xp.beProductOf(Qt, help);
    answer.resize(2);
    answer.at(1) = xp.computeNorm();
    answer.at(2) = atan2( xp.at(2), xp.at(1) );
}

void
IntElPoint :: computeLocalSlipDir(FloatArray &normal)
{
    normal.resizeWithValues(3);
    if ( this->referenceNode ) {
        // normal
        normal.beDifferenceOf(*domain->giveNode(this->referenceNode)->giveCoordinates(), *this->giveNode(1)->giveCoordinates());
    } else {
        if ( normal.at(1) == 0 && normal.at(2) == 0 && normal.at(3) == 0 ) {
            OOFEM_ERROR("Normal is not defined (referenceNode=0,normal=(0,0,0))");
        }
    }

    normal.normalize();
}

void TransportGradientPeriodic :: computeDofTransformation(ActiveDof *dof, FloatArray &masterContribs)
{
    DofManager *master = this->domain->giveDofManager(this->slavemap[dof->giveDofManager()->giveNumber()]);
    FloatArray *coords = dof->giveDofManager()->giveCoordinates();
    FloatArray *masterCoords = master->giveCoordinates();

    FloatArray dx;
    dx.beDifferenceOf(* coords, * masterCoords );

    masterContribs.resize(dx.giveSize() + 1);

    masterContribs.at(1) = 1.; // Master dof is always weight 1.0
    for ( int i = 1; i <= dx.giveSize(); ++i ) {
        masterContribs.at(i+1) = dx.at(i);
    }
}

void PrescribedGradientBCPeriodic :: findSlaveToMasterMap()
{
    FloatArray coord;
    SpatialLocalizer *sl = this->domain->giveSpatialLocalizer();
    //Set *masterSet = this->domain->giveSet(2);
    const IntArray &nodes = this->domain->giveSet(this->set)->giveNodeList(); // Split into slave set and master set?
    int nsd = jump.giveSize();
    std :: vector< FloatArray > jumps;
    // Construct all the possible jumps;
    jumps.reserve((2 << (nsd-1)) - 1);
    if ( nsd != 3 ) {
        OOFEM_ERROR("Only 3d implemented yet!");
    }
    jumps.emplace_back(jump);
    jumps.emplace_back(FloatArray{jump.at(1), jump.at(2), 0.});
    jumps.emplace_back(FloatArray{jump.at(1), 0., jump.at(3)});
    jumps.emplace_back(FloatArray{0., jump.at(2), jump.at(3)});
    jumps.emplace_back(FloatArray{jump.at(1), 0., 0.});
    jumps.emplace_back(FloatArray{0., jump.at(2), 0.});
    jumps.emplace_back(FloatArray{0., 0., jump.at(3)});

    this->slavemap.clear();
    for ( int inode : nodes ) {
        Node *masterNode = NULL;
        Node *node = this->domain->giveNode(inode);
        const FloatArray &masterCoord = *node->giveCoordinates();
        //printf("node %d\n", node->giveLabel()); masterCoord.printYourself();
        // The difficult part, what offset to subtract to find the master side;
        for ( FloatArray &testJump : jumps ) {
            coord.beDifferenceOf(masterCoord, testJump);
            masterNode = sl->giveNodeClosestToPoint(coord, fabs(jump.at(1))*1e-5);
            if ( masterNode != NULL ) {
                //printf("Found master (%d) to node %d (distance = %e)\n",  masterNode->giveNumber(), node->giveNumber(),
                //       masterNode->giveCoordinates()->distance(coord));
                break;
            }
        }
        if ( masterNode != NULL ) {
            this->slavemap.insert({node->giveNumber(), masterNode->giveNumber()});
        } else {
            OOFEM_ERROR("Couldn't find master node!");
        }
    }
}

double
FEI3dLineLin :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    ///@todo Not clear what this function should return. Just dNds would make sense if the caller defines a local coordinate system.
    FloatArray vec;
    vec.beDifferenceOf(*cellgeo.giveVertexCoordinates(2), *cellgeo.giveVertexCoordinates(1));

    double detJ = vec.computeSquaredNorm() * 0.5;
    double l2_inv = 0.5 / detJ;
    answer.resize(2, 3);

    answer.at(1, 1) = -vec.at(1)*l2_inv;
    answer.at(2, 1) =  vec.at(1)*l2_inv;
    answer.at(1, 2) = -vec.at(2)*l2_inv;
    answer.at(2, 2) =  vec.at(2)*l2_inv;
    answer.at(1, 3) = -vec.at(3)*l2_inv;
    answer.at(2, 3) =  vec.at(3)*l2_inv;

    return detJ;
}

void
StructuralFE2Material :: give3dMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    if ( useNumTangent ) {
        // Numerical tangent
        StructuralFE2MaterialStatus *status = static_cast<StructuralFE2MaterialStatus*>( this->giveStatus( gp ) );
        double h = 1.0e-9;

        const FloatArray &epsRed = status->giveTempStrainVector();
        FloatArray eps;
        StructuralMaterial::giveFullSymVectorForm(eps, epsRed, gp->giveMaterialMode() );


        int dim = eps.giveSize();
        answer.resize(dim, dim);
        answer.zero();

        FloatArray sig, sigPert, epsPert;

        for(int i = 1; i <= dim; i++) {
            // Add a small perturbation to the strain
            epsPert = eps;
            epsPert.at(i) += h;

            giveRealStressVector_3d(sigPert, gp, epsPert, tStep);
            answer.setColumn(sigPert, i);
        }

        giveRealStressVector_3d(sig, gp, eps, tStep);

        for(int i = 1; i <= dim; i++) {
            for(int j = 1; j <= dim; j++) {
                answer.at(j,i) -= sig.at(j);
                answer.at(j,i) /= h;
            }
        }

    } else {

        StructuralFE2MaterialStatus *ms = static_cast< StructuralFE2MaterialStatus * >( this->giveStatus(gp) );
        ms->computeTangent(tStep);
        const FloatMatrix &ans9 = ms->giveTangent();

        StructuralMaterial::giveReducedSymMatrixForm(answer, ans9, _3dMat);

//        const FloatMatrix &ans9 = ms->giveTangent();
//        printf("ans9: "); ans9.printYourself();
//
//        // Compute the (minor) symmetrized tangent:
//        answer.resize(6, 6);
//        for ( int i = 0; i < 6; ++i ) {
//            for ( int j = 0; j < 6; ++j ) {
//                answer(i, j) = ans9(i, j);
//            }
//        }
//        for ( int i = 0; i < 6; ++i ) {
//            for ( int j = 6; j < 9; ++j ) {
//                answer(i, j-3) += ans9(i, j);
//                answer(j-3, i) += ans9(j, i);
//            }
//        }
//        for ( int i = 6; i < 9; ++i ) {
//            for ( int j = 6; j < 9; ++j ) {
//                answer(j-3, i-3) += ans9(j, i);
//            }
//        }
//        for ( int i = 0; i < 6; ++i ) {
//            for ( int j = 3; j < 6; ++j ) {
//                answer(j, i) *= 0.5;
//                answer(i, j) *= 0.5;
//            }
//        }



#if 0
        // Numerical ATS for debugging
        FloatMatrix numericalATS(6, 6);
        FloatArray dsig;
        // Note! We need a copy of the temp strain, since the pertubations might change it.
        FloatArray tempStrain = ms->giveTempStrainVector();

        FloatArray sig, strain, sigPert;
        giveRealStressVector_3d(sig, gp, tempStrain, tStep);
        double hh = 1e-6;
        for ( int k = 1; k <= 6; ++k ) {
            strain = tempStrain;
            strain.at(k) += hh;
            giveRealStressVector_3d(sigPert, gp, strain, tStep);
            dsig.beDifferenceOf(sigPert, sig);
            numericalATS.setColumn(dsig, k);
        }
        numericalATS.times(1. / hh);
        giveRealStressVector_3d(sig, gp, tempStrain, tStep); // Reset

        //answer.printYourself("Analytical deviatoric tangent");
        //numericalATS.printYourself("Numerical deviatoric tangent");

        numericalATS.subtract(answer);
        double norm = numericalATS.computeFrobeniusNorm();
//.........这里部分代码省略.........

void
MisesMat :: giveFirstPKStressVector_3d(FloatArray &answer,
                                       GaussPoint *gp,
                                       const FloatArray &totalDefGradOOFEM,
                                       TimeStep *tStep)
{
    MisesMatStatus *status = static_cast< MisesMatStatus * >( this->giveStatus(gp) );

    double kappa, dKappa, yieldValue, mi;
    FloatMatrix F, oldF, invOldF;
    FloatArray s;
    F.beMatrixForm(totalDefGradOOFEM); //(method assumes full 3D)

    kappa = status->giveCumulativePlasticStrain();
    oldF.beMatrixForm( status->giveFVector() );
    invOldF.beInverseOf(oldF);
    //relative deformation radient
    FloatMatrix f;
    f.beProductOf(F, invOldF);
    //compute elastic predictor
    FloatMatrix trialLeftCauchyGreen, help;

    f.times( 1./cbrt(f.giveDeterminant()) );

    help.beProductOf(f, status->giveTempLeftCauchyGreen());
    trialLeftCauchyGreen.beProductTOf(help, f);
    FloatMatrix E;
    E.beTProductOf(F, F);
    E.at(1, 1) -= 1.0;
    E.at(2, 2) -= 1.0;
    E.at(3, 3) -= 1.0;
    E.times(0.5);

    FloatArray e;
    e.beSymVectorFormOfStrain(E);

    FloatArray leftCauchyGreen;
    FloatArray leftCauchyGreenDev;
    double leftCauchyGreenVol;

    leftCauchyGreen.beSymVectorFormOfStrain(trialLeftCauchyGreen);

    leftCauchyGreenVol = computeDeviatoricVolumetricSplit(leftCauchyGreenDev, leftCauchyGreen);
    FloatArray trialStressDev;
    applyDeviatoricElasticStiffness(trialStressDev, leftCauchyGreenDev, G / 2.);
    s = trialStressDev;

    //check for plastic loading
    double trialS = computeStressNorm(trialStressDev);
    double sigmaY = sig0 + H * kappa;
    //yieldValue = sqrt(3./2.)*trialS-sigmaY;
    yieldValue = trialS - sqrt(2. / 3.) * sigmaY;


    //store deviatoric trial stress(reused by algorithmic stiffness)
    status->letTrialStressDevBe(trialStressDev);
    //the return-mapping algorithm
    double J = F.giveDeterminant();
    mi = leftCauchyGreenVol * G;
    if ( yieldValue > 0 ) {
        //dKappa =sqrt(3./2.)* yieldValue/(H + 3.*mi);
        //kappa = kappa + dKappa;
        //trialStressDev.times(1-sqrt(6.)*mi*dKappa/trialS);
        dKappa = ( yieldValue / ( 2 * mi ) ) / ( 1 + H / ( 3 * mi ) );
        FloatArray n = trialStressDev;
        n.times(2 * mi * dKappa / trialS);
        ////return map
        s.beDifferenceOf(trialStressDev, n);
        kappa += sqrt(2. / 3.) * dKappa;


        //update of intermediate configuration
        trialLeftCauchyGreen.beMatrixForm(s);
        trialLeftCauchyGreen.times(1.0 / G);
        trialLeftCauchyGreen.at(1, 1) += leftCauchyGreenVol;
        trialLeftCauchyGreen.at(2, 2) += leftCauchyGreenVol;
        trialLeftCauchyGreen.at(2, 2) += leftCauchyGreenVol;
        trialLeftCauchyGreen.times(J * J);
    }

    //addition of the elastic mean stress
    FloatMatrix kirchhoffStress;
    kirchhoffStress.beMatrixForm(s);
    kirchhoffStress.at(1, 1) += 1. / 2. *  K * ( J * J - 1 );
    kirchhoffStress.at(2, 2) += 1. / 2. *  K * ( J * J - 1 );
    kirchhoffStress.at(3, 3) += 1. / 2. *  K * ( J * J - 1 );


    FloatMatrix iF, Ep(3, 3), S;
    FloatArray vF, vS, ep;

    //transform Kirchhoff stress into Second Piola - Kirchhoff stress
    iF.beInverseOf(F);
    help.beProductOf(iF, kirchhoffStress);
    S.beProductTOf(help, iF);

    this->computeGLPlasticStrain(F, Ep, trialLeftCauchyGreen, J);

    ep.beSymVectorFormOfStrain(Ep);
    vS.beSymVectorForm(S);
//.........这里部分代码省略.........

void PLHoopStressCirc :: propagateInterfaces(Domain &iDomain, EnrichmentDomain &ioEnrDom)
{
    // Fetch crack tip data
    TipInfo tipInfoStart, tipInfoEnd;
    ioEnrDom.giveTipInfos(tipInfoStart, tipInfoEnd);
    std :: vector< TipInfo >tipInfo = {tipInfoStart, tipInfoEnd};

    SpatialLocalizer *localizer = iDomain.giveSpatialLocalizer();

    for ( size_t tipIndex = 0; tipIndex < tipInfo.size(); tipIndex++ ) {
        // Construct circle points on an arc from -90 to 90 degrees
        double angle = -90.0 + mAngleInc;
        std :: vector< double >angles;
        while ( angle <= ( 90.0 - mAngleInc ) ) {
            angles.push_back(angle * M_PI / 180.0);
            angle += mAngleInc;
        }

        const FloatArray &xT    = tipInfo [ tipIndex ].mGlobalCoord;
        const FloatArray &t             = tipInfo [ tipIndex ].mTangDir;
        const FloatArray &n             = tipInfo [ tipIndex ].mNormalDir;

        // It is meaningless to propagate a tip that is not inside any element
        Element *el = localizer->giveElementContainingPoint(tipInfo [ tipIndex ].mGlobalCoord);
        if ( el != NULL ) {
            std :: vector< FloatArray >circPoints;

            for ( size_t i = 0; i < angles.size(); i++ ) {
                FloatArray tangent(2);
                tangent.zero();
                tangent.add(cos(angles [ i ]), t);
                tangent.add(sin(angles [ i ]), n);
                tangent.normalize();

                FloatArray x(xT);
                x.add(mRadius, tangent);
                circPoints.push_back(x);
            }



            std :: vector< double >sigTTArray, sigRTArray;

            // Loop over circle points
            for ( size_t pointIndex = 0; pointIndex < circPoints.size(); pointIndex++ ) {
                FloatArray stressVec;

                if ( mUseRadialBasisFunc ) {
                    // Interpolate stress with radial basis functions

                    // Choose a cut-off length l:
                    // take the distance between two nodes in the element containing the
                    // crack tip multiplied by a constant factor.
                    // ( This choice implies that we hope that the element has reasonable
                    // aspect ratio.)
                    const FloatArray &x1 = * ( el->giveDofManager(1)->giveCoordinates() );
                    const FloatArray &x2 = * ( el->giveDofManager(2)->giveCoordinates() );
                    const double l = 1.0 * x1.distance(x2);

                    // Use the octree to get all elements that have
                    // at least one Gauss point in a certain region around the tip.
                    const double searchRadius = 3.0 * l;
                    std :: set< int >elIndices;
                    localizer->giveAllElementsWithIpWithinBox(elIndices, circPoints [ pointIndex ], searchRadius);


                    // Loop over the elements and Gauss points obtained.
                    // Evaluate the interpolation.
                    FloatArray sumQiWiVi;
                    double sumWiVi = 0.0;
                    for ( int elIndex: elIndices ) {
                        Element *gpEl = iDomain.giveElement(elIndex);
                        IntegrationRule *iRule = gpEl->giveDefaultIntegrationRulePtr();

                        for ( GaussPoint *gp_i: *iRule ) {

                            ////////////////////////////////////////
                            // Compute global gp coordinates
                            FloatArray N;
                            FEInterpolation *interp = gpEl->giveInterpolation();
                            interp->evalN( N, * ( gp_i->giveCoordinates() ), FEIElementGeometryWrapper(gpEl) );


                            // Compute global coordinates of Gauss point
                            FloatArray globalCoord(2);
                            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;
//.........这里部分代码省略.........

void FE2FluidMaterial :: giveStiffnessMatrices(FloatMatrix &dsdd, FloatArray &dsdp, FloatArray &dedd, double &dedp,
                                               MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    FE2FluidMaterialStatus *ms = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
    ms->computeTangents(tStep);
    if ( mode == TangentStiffness ) {
        dsdd = ms->giveDeviatoricTangent();
        dsdp = ms->giveDeviatoricPressureTangent();
        dedd = ms->giveVolumetricDeviatoricTangent();
        dedp = ms->giveVolumetricPressureTangent();
#if 0
        // Numerical ATS for debugging
        FloatMatrix numericalATS(6, 6);
        FloatArray dsig;
        FloatArray tempStrain(6);

        tempStrain.zero();
        FloatArray sig, strain, sigPert;
        double epspvol;
        computeDeviatoricStressVector(sig, epspvol, gp, tempStrain, 0., tStep);
        double h = 0.001; // Linear problem, size of this doesn't matter.
        for ( int k = 1; k <= 6; ++k ) {
            strain = tempStrain;
            strain.at(k) += h;
            double tmp = strain.at(1) + strain.at(2) + strain.at(3);
            strain.at(1) -= tmp/3.0;
            strain.at(2) -= tmp/3.0;
            strain.at(3) -= tmp/3.0;
            strain.printYourself();
            computeDeviatoricStressVector(sigPert, epspvol, gp, strain, 0., tStep);
            sigPert.printYourself();
            dsig.beDifferenceOf(sigPert, sig);
            numericalATS.setColumn(dsig, k);
        }
        numericalATS.times(1. / h);

        printf("Analytical deviatoric tangent = ");
        dsdd.printYourself();
        printf("Numerical deviatoric tangent = ");
        numericalATS.printYourself();
        numericalATS.subtract(dsdd);
        double norm = numericalATS.computeFrobeniusNorm();
        if ( norm > dsdd.computeFrobeniusNorm() * DEBUG_ERR && norm > 0.0 ) {
            OOFEM_ERROR("Error in deviatoric tangent");
        }
#endif
#if 0
        // Numerical ATS for debugging
        FloatArray strain(3);
        strain.zero();
        FloatArray sig, sigh;
        double epspvol, pressure = 0.0;
        double h = 1.00; // Linear problem, size of this doesn't matter.
        computeDeviatoricStressVector(sig, epspvol, gp, strain, pressure, tStep);
        computeDeviatoricStressVector(sigh, epspvol, gp, strain, pressure + h, tStep);

        FloatArray dsigh;
        dsigh.beDifferenceOf(sigh, sig);
        dsigh.times(1 / h);

        printf("Analytical deviatoric pressure tangent = ");
        dsdp.printYourself();
        printf("Numerical deviatoric pressure tangent = ");
        dsigh.printYourself();
        dsigh.subtract(dsdp);
        double norm = dsigh.computeNorm();
        if ( norm > dsdp.computeNorm() * DEBUG_ERR && norm > 0.0 ) {
            OOFEM_ERROR("Error in deviatoric pressure tangent");
        }
#endif
#if 0
        // Numerical ATS for debugging
        FloatArray tempStrain(3);
        tempStrain.zero();
        FloatArray sig, strain;
        double epspvol, epspvol11, epspvol22, epspvol12, pressure = 0.0;
        double h = 1.0; // Linear problem, size of this doesn't matter.

        computeDeviatoricStressVector(sig, epspvol, gp, tempStrain, pressure, tStep);
        strain = tempStrain;
        strain.at(1) += h;
        computeDeviatoricStressVector(sig, epspvol11, gp, strain, pressure, tStep);
        strain = tempStrain;
        strain.at(2) += h;
        computeDeviatoricStressVector(sig, epspvol22, gp, strain, pressure, tStep);
        strain = tempStrain;
        strain.at(3) += h;
        computeDeviatoricStressVector(sig, epspvol12, gp, strain, pressure, tStep);

        FloatArray dvol(3);
        dvol.at(1) = ( epspvol11 - epspvol ) / h;
        dvol.at(2) = ( epspvol22 - epspvol ) / h;
        dvol.at(3) = ( epspvol12 - epspvol ) / h;
        dvol.at(1) += 1.0;
        dvol.at(2) += 1.0;

        printf("Analytical volumetric deviatoric tangent = ");
        dedd.printYourself();
        printf("Numerical volumetric deviatoric tangent = ");
        dvol.printYourself();
//.........这里部分代码省略.........

void IntMatBilinearCZ :: giveFirstPKTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump,
                                                const FloatMatrix &F, TimeStep *tStep)
{
    IntMatBilinearCZStatus *status = static_cast< IntMatBilinearCZStatus * >( this->giveStatus(gp) );

    status->mJumpNew = jump;

    FloatArray jumpInc;
    jumpInc.beDifferenceOf(status->mJumpNew, status->mJumpOld);

    FloatArray tractionTrial = status->mTractionOld;
    tractionTrial.add(mPenaltyStiffness, jumpInc);

    double TTrNormal    = tractionTrial.at(3);
    double TTrTang              = sqrt( pow(tractionTrial.at(1), 2.0) + pow(tractionTrial.at(2), 2.0) );
    double phiTr = computeYieldFunction(TTrNormal, TTrTang);

    const double damageTol = 1.0e-6;
    if ( status->mDamageOld > ( 1.0 - damageTol ) ) {
        status->mDamageNew = 1.0;
        status->mPlastMultIncNew = 0.0;
        answer.resize(3);
        answer.zero();
        status->mTractionNew = answer;

        status->letTempJumpBe(jump);
        status->letTempFirstPKTractionBe(answer);
        status->letTempTractionBe(answer);

        return;
    }


    answer = tractionTrial;

    if ( phiTr < 0.0 ) {
        status->mDamageNew = status->mDamageOld;
        status->mPlastMultIncNew = 0.0;
        answer.beScaled( ( 1.0 - status->mDamageNew ), answer );

        status->mTractionNew = answer;

        status->letTempJumpBe(jump);
        status->letTempFirstPKTractionBe(answer);
        status->letTempTractionBe(answer);

        return;
    } else {
        // Iterate to find plastic strain increment.
        int maxIter = 50;
        int minIter = 1;
        double absTol = 1.0e-9; // Absolute error tolerance
        double relTol = 1.0e-9; // Relative error tolerance
        double eps = 1.0e-12; // Small value for perturbation when computing numerical Jacobian
        double plastMultInc = 0.0;
        double initialRes = 0.0;

        for ( int iter = 0; iter < maxIter; iter++ ) {
            // Evaluate residual (i.e. yield function)
            computeTraction(answer, tractionTrial, plastMultInc);

            double TNormal      = answer.at(3);
            double TTang                = sqrt( pow(answer.at(1), 2.0) + pow(answer.at(2), 2.0) );
            double phi = computeYieldFunction(TNormal, TTang);

            //          if(iter > 20) {
            //              printf("iter: %d res: %e\n", iter, fabs(phi) );
            //          }

            if ( iter == 0 ) {
                initialRes = fabs(phi);
                initialRes = max(initialRes, 1.0e-12);
            }

            if ( (iter >= minIter && fabs(phi) < absTol) || ( iter >= minIter && ( fabs(phi) / initialRes ) < relTol ) ) {
                // Add damage evolution
                double S = mGIc / mSigmaF;
                status->mPlastMultIncNew = plastMultInc;

                double damageInc = status->mPlastMultIncNew / S;
                status->mDamageNew = status->mDamageOld + damageInc;

                if ( status->mDamageNew > 1.0 ) {
                    status->mDamageNew = 1.0;
                }

                if(mSemiExplicit) {
//                    computeTraction(answer, tractionTrial, status->mPlastMultIncOld);
                    answer.beScaled( ( 1.0 - status->mDamageOld ), answer );
               }
                else {
                    answer.beScaled( ( 1.0 - status->mDamageNew ), answer );
                }

                status->mTractionNew = answer;

                // Jim
                status->letTempJumpBe(jump);
                status->letTempFirstPKTractionBe(answer);
                status->letTempTractionBe(answer);
//.........这里部分代码省略.........