C++ NLStructuralElement类代码示例

void
GradDpElement :: giveNonlocalInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    double dV, localCumulatedStrain = 0.;
    NLStructuralElement *elem = this->giveNLStructuralElement();
    FloatMatrix stiffKappa, Nk;
    FloatArray fKappa(nlSize), aux(nlSize), dKappa, stress;

    aux.zero();
    int size = nSecVars * nSecNodes;

    //set displacement and nonlocal location array
    this->setDisplacementLocationArray(locU, nPrimNodes, nPrimVars, nSecNodes, nSecVars);
    this->setNonlocalLocationArray(locK, nPrimNodes, nPrimVars, nSecNodes, nSecVars);

    answer.resize(size);
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {
        this->computeNkappaMatrixAt(gp, Nk);
        for ( int j = 1; j <= nlSize; j++ ) {
            fKappa.at(j) = Nk.at(1, j);
        }

        dV  = elem->computeVolumeAround(gp);
        this->computeStressVectorAndLocalCumulatedStrain(stress, localCumulatedStrain, gp, tStep);
        fKappa.times(localCumulatedStrain);
        fKappa.times(-dV);
        aux.add(fKappa);
    }

    this->computeStiffnessMatrix_kk(stiffKappa, TangentStiffness, tStep);
    this->computeNonlocalDegreesOfFreedom(dKappa, tStep);
    answer.beProductOf(stiffKappa, dKappa);
    answer.add(aux);
}

void
GradDpElement :: computeStressVectorAndLocalCumulatedStrain(FloatArray &answer, double localCumulatedStrain, GaussPoint *gp, TimeStep *tStep)
{
    NLStructuralElement *elem = this->giveNLStructuralElement();

    double nlCumulatedStrain;

    int nlGeo = elem->giveGeometryMode();
    StructuralCrossSection *cs = elem->giveStructuralCrossSection();
    GradDpMaterialExtensionInterface *dpmat = static_cast< GradDpMaterialExtensionInterface * >(
        cs->giveMaterialInterface(GradDpMaterialExtensionInterfaceType, gp) );

    if ( !dpmat ) {
        OOFEM_ERROR("Material doesn't implement the required DpGrad interface!");
    }

    this->computeNonlocalCumulatedStrain(nlCumulatedStrain, gp, tStep);
    if ( nlGeo == 0 ) {
        FloatArray Epsilon;
        this->computeLocalStrainVector(Epsilon, gp, tStep);
        dpmat->giveRealStressVectorGrad(answer, localCumulatedStrain, gp, Epsilon, nlCumulatedStrain, tStep);
    } else if ( nlGeo == 1 ) {
        if ( elem->giveDomain()->giveEngngModel()->giveFormulation() == TL ) {
            FloatArray vF;
            this->computeDeformationGradientVector(vF, gp, tStep);
            dpmat->giveFirstPKStressVectorGrad(answer, localCumulatedStrain, gp, vF, nlCumulatedStrain, tStep);
        } else {
            FloatArray vF;
            this->computeDeformationGradientVector(vF, gp, tStep);
            dpmat->giveCauchyStressVectorGrad(answer, localCumulatedStrain, gp, vF, nlCumulatedStrain, tStep);
        }
    }
}

void
PhaseFieldElement :: giveInternalForcesVector_d(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    // Computes int_V ( N^t *Nstress_d  +  B^t * g_c*l*grad(d)  ) dV
    FloatArray NStress, BStress, NS, a_d, grad_d;
    FloatMatrix N, B;
    NLStructuralElement *el = this->giveElement();
    this->computeDamageUnknowns( a_d, VM_Total, tStep );

    for ( auto &gp: el->giveIntegrationRule(0) ) {
        double dV = el->computeVolumeAround(gp);
        
        // compute generalized stress measures
        this->computeNd_matrixAt( *gp->giveNaturalCoordinates(), N);
        computeNStress_d(NStress, gp, tStep, useUpdatedGpRecord);
        NS.beTProductOf(N, NStress);
        answer.add(dV, NS);

        this->computeBd_matrixAt(gp, B);
        grad_d.beProductOf(B, a_d);
        double l = this->giveInternalLength();
        double g_c = this->giveCriticalEnergy();
        BStress = grad_d * l * g_c;
        answer.plusProduct(B, BStress, dV);
    }
}

void
PhaseFieldElement :: computeStiffnessMatrix_ud(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    FloatMatrix B, DB, N_d, DN, S(3,1);
    FloatArray stress;
    NLStructuralElement *el = this->giveElement();
    StructuralCrossSection *cs = dynamic_cast<StructuralCrossSection* > (el->giveCrossSection() );

    answer.clear();

    for ( auto &gp: el->giveIntegrationRule(0) ) {
        StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

        double dV = el->computeVolumeAround(gp);
        // compute int_V ( B^t * D_B * B )dV

        this->computeNd_matrixAt(*gp->giveNaturalCoordinates(), N_d);

        // stress   
        FloatArray reducedStrain, a_u;
        FloatMatrix B_u;
        el->computeBmatrixAt( gp, B_u );
        stress = matStat->giveTempStressVector();
        stress.times( this->computeGPrim( gp, VM_Total, tStep ) );
        
        S.setColumn(stress,1);
        DN.beProductOf(S, N_d);
        answer.plusProductUnsym(B_u, DN, dV);
    }

}

void
GradDpElement :: giveNonlocalInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    double localCumulatedStrain = 0.;
    NLStructuralElement *elem = this->giveNLStructuralElement();
    FloatMatrix stiffKappa;
    FloatArray Nk;
    FloatArray aux, dKappa, stress;

    int size = nSecVars * nSecNodes;

    //set displacement and nonlocal location array
    this->setDisplacementLocationArray();
    this->setNonlocalLocationArray();

    answer.resize(size);
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {
        this->computeNkappaMatrixAt(gp, Nk);

        double dV = elem->computeVolumeAround(gp);
        this->computeStressVectorAndLocalCumulatedStrain(stress, localCumulatedStrain, gp, tStep);
        aux.add(-dV * localCumulatedStrain, Nk);
    }

    this->computeStiffnessMatrix_kk(stiffKappa, TangentStiffness, tStep);
    this->computeNonlocalDegreesOfFreedom(dKappa, tStep);
    answer.beProductOf(stiffKappa, dKappa);
    answer.add(aux);
}

void
PhaseFieldElement :: computeStiffnessMatrix_dd(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    double Delta_t = tStep->giveTimeIncrement();
    double t_star = this->giveRelaxationTime();
    double l = this->giveInternalLength();
    double g_c = this->giveCriticalEnergy();

    FloatMatrix B_d, N_d;
    //StructuralCrossSection *cs = dynamic_cast<StructuralCrossSection* > (this->giveCrossSection() );
    NLStructuralElement *el = this->giveElement();
    answer.clear();

    for ( auto &gp: el->giveIntegrationRule(0) ) {
        double dV = el->computeVolumeAround(gp);

        this->computeNd_matrixAt(*gp->giveNaturalCoordinates(), N_d);
        this->computeBd_matrixAt(gp, B_d, 1, 3);

        double Gprim = this->computeGPrim(gp, VM_Total, tStep);
        double psiBar = this->computeFreeEnergy( gp, tStep );
        //double factorN = t_star/Delta_t + g_c/l + psiBar*Gbis;
        double factorN = t_star / Delta_t + g_c / l + psiBar*(-2.0);
        double factorB = g_c*l;

        answer.plusProductSymmUpper(N_d, N_d, factorN*dV);
        answer.plusProductSymmUpper(B_d, B_d, factorB*dV);   
    }

    answer.symmetrized();
}

void
PhaseFieldElement :: computeStiffnessMatrix_uu(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    // This is the regular stiffness matrix times G
    FloatMatrix B, DB, N, D_B;
    NLStructuralElement *el = this->giveElement();
    StructuralCrossSection *cs = dynamic_cast<StructuralCrossSection* > (el->giveCrossSection() );

    bool matStiffSymmFlag = cs->isCharacteristicMtrxSymmetric(rMode);
    answer.clear();

    for ( auto gp: el->giveIntegrationRule(0) ) {
        double dV = el->computeVolumeAround(gp);
        // compute int_V ( B^t * D_B * B )dV
        el->computeBmatrixAt(gp, B );
        cs->giveCharMaterialStiffnessMatrix(D_B, rMode, gp, tStep);
        D_B.times( computeG(gp, VM_Total, tStep) );
        DB.beProductOf(D_B, B);

        if ( matStiffSymmFlag ) {
            answer.plusProductSymmUpper(B, DB, dV);
        } else {
            answer.plusProductUnsym(B, DB, dV);
        }
    }

    if ( matStiffSymmFlag ) {
        answer.symmetrized();
    }
}

void
GradDpElement :: computeDeformationGradientVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    // Computes the deformation gradient in the Voigt format at the Gauss point gp of
    // the receiver at time step tStep.
    // Order of components: 11, 22, 33, 23, 13, 12, 32, 31, 21 in the 3D.

    // Obtain the current displacement vector of the element and subtract initial displacements (if present)
    FloatArray u;
    FloatMatrix B;
    NLStructuralElement *elem = this->giveNLStructuralElement();

    this->computeDisplacementDegreesOfFreedom(u, tStep);
    // Displacement gradient H = du/dX
    elem->computeBHmatrixAt(gp, B);
    answer.beProductOf(B, u);

    // Deformation gradient F = H + I
    MaterialMode matMode = gp->giveMaterialMode();
    if ( matMode == _3dMat || matMode == _PlaneStrain ) {
        answer.at(1) += 1.0;
        answer.at(2) += 1.0;
        answer.at(3) += 1.0;
    } else if ( matMode == _PlaneStress ) {
        answer.at(1) += 1.0;
        answer.at(2) += 1.0;
    } else if ( matMode == _1dMat ) {
        answer.at(1) += 1.0;
    } else {
        OOFEM_ERROR( "MaterialMode is not supported yet (%s)", __MaterialModeToString(matMode) );
    }
}

void
PhaseFieldElement :: computeNd_matrixAt(const FloatArray &lCoords, FloatMatrix &N)
{
    NLStructuralElement *el = this->giveElement();
    FloatArray Nvec;
    el->giveInterpolation( )->evalN( Nvec, lCoords, FEIElementGeometryWrapper( el ) );
    N.resize(1, Nvec.giveSize());
    N.beNMatrixOf(Nvec,1);

}

void
GradDpElement :: computeLocalStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    FloatArray u;
    FloatMatrix b;
    NLStructuralElement *elem = this->giveNLStructuralElement();

    this->computeDisplacementDegreesOfFreedom(u, tStep);
    elem->computeBmatrixAt(gp, b);
    answer.beProductOf(b, u);
}

double 
PhaseFieldElement :: computeDamageAt(GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
{
    // d = N_d * a_d
    NLStructuralElement *el = this->giveElement();
    FloatArray dVec;
    computeDamageUnknowns(dVec, valueMode, stepN);
    FloatArray Nvec;
    el->giveInterpolation()->evalN(Nvec, *gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(el));
    return Nvec.dotProduct(dVec);
}

void
PhaseFieldElement :: computeBd_matrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
{
    // Returns the [numSpaceDim x nDofs] gradient matrix {B_d} of the receiver,
    // evaluated at gp.

    NLStructuralElement *el = this->giveElement();
    FloatMatrix dNdx;
    el->giveInterpolation( )->evaldNdx( dNdx, *gp->giveNaturalCoordinates( ), FEIElementGeometryWrapper( el ) );
    answer.beTranspositionOf( dNdx );
}

void
GradDpElement :: computeStiffnessMatrix_uk(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    NLStructuralElement *elem = this->giveNLStructuralElement();
    double dV;
    StructuralCrossSection *cs = elem->giveStructuralCrossSection();
    int nlGeo = elem->giveGeometryMode();
    FloatArray Nk;
    FloatMatrix B, SNk, gPSigma;

    answer.clear();

    for ( auto &gp: *elem->giveIntegrationRule(0) ) {

        GradDpMaterialExtensionInterface *dpmat = dynamic_cast< GradDpMaterialExtensionInterface * >(
            cs->giveMaterialInterface(GradDpMaterialExtensionInterfaceType, gp) );
        if ( !dpmat ) {
            OOFEM_ERROR("Material doesn't implement the required DpGrad interface!");
        }
        dpmat->givePDGradMatrix_uk(gPSigma, rMode, gp, tStep);
        this->computeNkappaMatrixAt(gp, Nk);
        elem->computeBmatrixAt(gp, B);
        if ( nlGeo == 1 ) {
            if ( elem->domain->giveEngngModel()->giveFormulation() == AL ) {
                elem->computeBmatrixAt(gp, B);
            } else {
                elem->computeBHmatrixAt(gp, B);
            }
        }
        dV = elem->computeVolumeAround(gp);

        SNk.beProductOf(gPSigma, Nk);
        answer.plusProductUnsym(B, SNk, -dV);
    }
}

void
PhaseFieldElement :: computeBStress_u(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, int useUpdatedGpRecord)
{
    // computes G(d)*sig(u)
    NLStructuralElement *el = this->giveElement();
    FloatArray strain, a_u;
    FloatMatrix B_u;

    el->computeBmatrixAt(gp, B_u);
    this->computeDisplacementUnknowns(a_u, VM_Total, tStep);
    strain.beProductOf(B_u, a_u);
    el->computeStressVector(answer, strain, gp, tStep);
    answer.times( this->computeG(gp, VM_Total, tStep) );

}

void
GradDpElement :: computeStiffnessMatrix_uu(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    NLStructuralElement *elem = this->giveNLStructuralElement();
    StructuralCrossSection *cs = elem->giveStructuralCrossSection();
    FloatMatrix B, D, DB;


    int nlGeo = elem->giveGeometryMode();
    bool matStiffSymmFlag = elem->giveCrossSection()->isCharacteristicMtrxSymmetric(rMode);
    answer.clear();
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {

        GradDpMaterialExtensionInterface *dpmat = dynamic_cast< GradDpMaterialExtensionInterface * >(
                    cs->giveMaterialInterface(GradDpMaterialExtensionInterfaceType, gp) );
        if ( !dpmat ) {
            OOFEM_ERROR("Material doesn't implement the required DpGrad interface!");
        }
        if ( nlGeo == 0 ) {
            elem->computeBmatrixAt(gp, B);
        } else if ( nlGeo == 1 ) {
            if ( elem->domain->giveEngngModel()->giveFormulation() == AL ) {
                elem->computeBmatrixAt(gp, B);
            } else {
                elem->computeBHmatrixAt(gp, B);
            }
        }

        dpmat->givePDGradMatrix_uu(D, rMode, gp, tStep);
        double dV = elem->computeVolumeAround(gp);
        DB.beProductOf(D, B);
        if ( matStiffSymmFlag ) {
            answer.plusProductSymmUpper(B, DB, dV);
        } else {
            answer.plusProductUnsym(B, DB, dV);
        }
    }

    if ( elem->domain->giveEngngModel()->giveFormulation() == AL ) {
        FloatMatrix initialStressMatrix;
        elem->computeInitialStressMatrix(initialStressMatrix, tStep);
        answer.add(initialStressMatrix);
    }

    if ( matStiffSymmFlag ) {
        answer.symmetrized();
    }
}