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

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 PrescribedGradientBCWeak :: computeField(FloatArray &sigma, TimeStep *tStep)
{
	double Lx = mUC[0] - mLC[0];
	double Ly = mUC[1] - mLC[1];
	double dSize = Lx*Ly;
//    printf("dSize: %e\n", dSize);

    const int dim = domain->giveNumberOfSpatialDimensions();
    FloatMatrix stressMatrix(dim, dim);

    for( auto *el: mpTracElNew ) {

        for ( GaussPoint *gp: *(el->mIntRule.get()) ) {

            // For now, assume piecewise constant approx
            FloatArray Ntrac = FloatArray { 1.0 };

            // N-matrix
            FloatMatrix Nmat;
            Nmat.beNMatrixOf(Ntrac, dim);

        	// Fetch global coordinate x
        	const FloatArray &x = gp->giveGlobalCoordinates();

            // Compute vector of traction unknowns
            FloatArray tracUnknowns;
            el->mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), VM_Total, tStep);


            FloatArray traction;
            traction.beProductOf(Nmat, tracUnknowns);

            FloatArray tmp;
            giveBoundaryCoordVector(tmp, x);

            FloatMatrix contrib;
            contrib.beDyadicProductOf(traction, tmp);

            double detJ = 0.5 * el->giveLength();
            contrib.times( detJ * gp->giveWeight() );

            for ( int m = 0; m < dim; m++ ) {
                for ( int n = 0; n < dim; n++ ) {
                    stressMatrix(m, n) += contrib(m, n);
                }
            }

        }

    }

    if ( dim == 2 ) {
        sigma = {
            stressMatrix(0, 0), stressMatrix(1, 1), 0.0, 0.0, 0.0, 0.5*(stressMatrix(0, 1) + stressMatrix(1, 0))
        };
    } else   {
        sigma.beVectorForm(stressMatrix);
    }

    sigma.times(1.0 / dSize);

}

//.........这里部分代码省略.........
    stress.changeComponentOrder();
    strain.changeComponentOrder();
    strainIncrement.changeComponentOrder();

    OOFEM_LOG_DEBUG("AbaqusUserMaterial :: giveRealStressVector - Calling subroutine");
    this->umat(stress.givePointer(), // STRESS
               state.givePointer(), // STATEV
               jacobian.givePointer(), // DDSDDE
               & sse, // SSE
               & spd, // SPD
               & scd, // SCD
               & rpl, // RPL
               ddsddt.givePointer(), // DDSDDT
               drplde.givePointer(), // DRPLDE
               & drpldt, // DRPLDT
               strain.givePointer(), // STRAN
               strainIncrement.givePointer(), // DSTRAN
               time, // TIME
               & dtime, // DTIME
               & temp, // TEMP
               & dtemp, // DTEMP
               & predef, // PREDEF
               & dpred, // DPRED
               this->cmname, // CMNAME
               & ndi, // NDI
               & nshr, // NSHR
               & ntens, // NTENS
               & numState, // NSTATV
               properties.givePointer(), // PROPS
               & numProperties, // NPROPS
               coords.givePointer(), // COORDS
               drot.givePointer(), // DROT
               & pnewdt, // PNEWDT
               & celent, // CELENT
               dfgrd0.givePointer(), // DFGRD0
               dfgrd1.givePointer(), // DFGRD1
               & noel, // NOEL
               & npt, // NPT
               & layer, // LAYER
               & kspt, // KSPT
               & kstep, // KSTEP
               & kinc); // KINC


    // Change to OOFEM's component order
    stress.changeComponentOrder();
    strain.changeComponentOrder();
    strainIncrement.changeComponentOrder();
    jacobian.changeComponentOrder();


    if ( mStressInterpretation == 0 ) {
        FloatMatrix P, cauchyStress;
        P.beMatrixForm(stress);

        FloatArray vP;
        vP.beVectorForm(P);

        cauchyStress.beProductTOf(P, F);
        cauchyStress.times( 1.0 / F.giveDeterminant() );

        FloatArray vCauchyStress;
        vCauchyStress.beSymVectorForm(cauchyStress);

        ms->letTempStrainVectorBe(strain);
        ms->letTempStressVectorBe(vCauchyStress);
        ms->letTempStateVectorBe(state);
        ms->letTempTangentBe(jacobian);
        ms->letTempPVectorBe(vP);
        ms->letTempFVectorBe(vF);

        answer = vP;
    } else   {
        FloatArray vP;
        FloatMatrix P, sigma, F_inv;
        F_inv.beInverseOf(F);

        sigma.beMatrixForm(stress);
        P.beProductOf(F, sigma);
        vP.beVectorForm(P);
        answer = vP;

        // Convert from sigma to S
        FloatMatrix S;
        S.beProductOf(F_inv, P);
        vS.beSymVectorForm(S);

        // @todo Should convert from dsigmadE to dSdE here
        // L2=detF*matmul( matmul ( inv(op_a_V9(F,F), cm-op_a_V9(ident,Tau)-op_b_V9(Tau,ident) ), inv(op_a_V9(Ft,Ft)))

        ms->letTempStrainVectorBe(strain);
        ms->letTempStressVectorBe(vS);
        ms->letTempStateVectorBe(state);
        ms->letTempTangentBe(jacobian);
        ms->letTempPVectorBe(vP);
        ms->letTempFVectorBe(vF);
    }

    OOFEM_LOG_DEBUG("AbaqusUserMaterial :: giveRealStressVector_3d - Calling subroutine was successful");
}

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

    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);
    vF.beVectorForm(F);
    answer.beVectorForm(kirchhoffStress);

    status->setTrialStressVol(mi);
    status->letTempLeftCauchyGreenBe(trialLeftCauchyGreen);
    status->setTempCumulativePlasticStrain(kappa);
    status->letTempStressVectorBe(answer);
    status->letTempStrainVectorBe(e);
    status->letTempPlasticStrainBe(ep);
    status->letTempPVectorBe(answer);
    status->letTempFVectorBe(vF);
}