C++ StructuralElement类代码示例

void
FreeWarping :: computeResultAtCenterOfGravity(TimeStep *tStep)
{
    int noCS = this->giveDomain(1)->giveNumberOfCrossSectionModels(); //number of warping Crosssections
    SolutionAtCG.resize(noCS);
    Element *closestElement;
    FloatArray lcoords,  closest, lcg;
    SpatialLocalizer *sp = this->giveDomain(1)->giveSpatialLocalizer();
    sp->init();
    lcoords.resize(2);
    closest.resize(2);
    lcg.resize(2);

    for ( int j = 1; j <= noCS; ++j ) {
        lcg.at(1) = CG.at(j, 1);
        lcg.at(2) = CG.at(j, 2);
        closestElement = sp->giveElementClosestToPoint(lcoords, closest, lcg, 0);

        StructuralElement *sE = dynamic_cast< StructuralElement * >(closestElement);
        FloatArray u, r, b;
        FloatMatrix N;
        sE->computeNmatrixAt(lcoords, N);
        sE->computeVectorOf(VM_Total, tStep, u);
        u.resizeWithValues(3);
        r.beProductOf(N, u);

        SolutionAtCG.at(j) = r.at(1);
    }
}

void
GradDpElement :: computeStiffnessMatrix_kk(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    StructuralElement *elem = this->giveStructuralElement();
    double dV;
    FloatMatrix lStiff;
    FloatArray Nk;
    FloatMatrix Bk, LBk;
    StructuralCrossSection *cs = elem->giveStructuralCrossSection();

    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!");
        }

        this->computeNkappaMatrixAt(gp, Nk);
        this->computeBkappaMatrixAt(gp, Bk);
        dV = elem->computeVolumeAround(gp);

        dpmat->givePDGradMatrix_kk(lStiff, rMode, gp, tStep);
        answer.plusProductUnsym(Nk, Nk, dV);
        if ( dpmat->giveAveragingType() == 0 || dpmat->giveAveragingType() == 1 ) {
            double l = lStiff.at(1, 1);
            answer.plusProductUnsym(Bk, Bk, l * l * dV);
        } else if ( dpmat->giveAveragingType() == 2 ) {
            LBk.beProductOf(lStiff, Bk);
            answer.plusProductUnsym(Bk, LBk, dV);
        }
    }
}

// computes Btransposed*eta for the given Gauss point
// (where eta is the derivative of cum. plastic strain wrt final strain)
void
RankineMatNl :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                        FloatArray &rcontrib, TimeStep *atTime)
{
    RankineMatNlStatus *status = ( RankineMatNlStatus * ) this->giveStatus(gp);
    StructuralElement *elem = ( StructuralElement * )( gp->giveElement() );
    elem->giveLocationArray(rloc, EID_MomentumBalance, s);
    FloatMatrix b;
    elem->computeBmatrixAt(gp, b);
    int ncols = b.giveNumberOfColumns();
    rcontrib.resize(ncols);

    double kappa = status->giveCumulativePlasticStrain();
    double tempKappa = status->giveTempCumulativePlasticStrain();
    if ( tempKappa <= kappa ) {
        rcontrib.zero();
        return;
    }

    int i, j;
    double sum;
    int nsize = 3;
    FloatArray eta(3);
    computeEta(eta, status);
    for ( i = 1; i <= ncols; i++ ) {
        sum = 0.;
        for ( j = 1; j <= nsize; j++ ) {
            sum += eta.at(j) * b.at(j, i);
        }

        rcontrib.at(i) = sum;
    }
}

void
TrabBoneNL3D :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                        FloatArray &rcontrib, TimeStep *tStep)
{
    TrabBoneNL3DStatus *nlStatus = static_cast< TrabBoneNL3DStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );

    FloatMatrix b;

    elem->giveLocationArray(rloc, s);
    elem->computeBmatrixAt(gp, b);

    double kappa = nlStatus->giveKappa();
    double tempKappa = nlStatus->giveTempKappa();
    double dKappa = tempKappa - kappa;
    if ( dKappa < 10.e-9 ) {
        dKappa = 0;
    }

    if ( dKappa > 0.0 ) {
        FloatArray remoteNu, prodTensor;
        const FloatArray &plasFlowDirec = nlStatus->givePlasFlowDirec();
        const FloatMatrix &SSaTensor = nlStatus->giveSSaTensor();
        double beta = nlStatus->giveBeta();

        prodTensor.beTProductOf(SSaTensor, plasFlowDirec);
        remoteNu = 1 / beta * prodTensor;
        rcontrib.beTProductOf(b, remoteNu);
    } else {
        rcontrib.resize(b.giveNumberOfColumns());
        rcontrib.zero();
    }
}

int
MisesMatNl :: giveLocalNonlocalStiffnessContribution(GaussPoint *gp, IntArray &loc, const UnknownNumberingScheme &s,
        FloatArray &lcontrib, TimeStep *atTime)
{
    int nrows, nsize, i, j;
    double sum, nlKappa, damage, tempDamage, dDamF;
    MisesMatNlStatus *status = ( MisesMatNlStatus * ) this->giveStatus(gp);
    StructuralElement *elem = ( StructuralElement * )( gp->giveElement() );
    FloatMatrix b;
    FloatArray stress;

    this->computeCumPlasticStrain(nlKappa, gp, atTime);
    damage = status->giveDamage();
    tempDamage = status->giveTempDamage();
    if ( ( tempDamage - damage ) > 0 ) {
        elem->giveLocationArray(loc, EID_MomentumBalance, s);
        status->giveTempEffectiveStress(stress);
        elem->computeBmatrixAt(gp, b);
        dDamF = computeDamageParamPrime(nlKappa);
        nrows = b.giveNumberOfColumns();
        nsize = stress.giveSize();
        lcontrib.resize(nrows);

        for ( i = 1; i <= nrows; i++ ) {
            sum = 0.0;
            for ( j = 1; j <= nsize; j++ ) {
                sum += b.at(j, i) * stress.at(j);
            }

            lcontrib.at(i) = sum * mm * dDamF;
        }
    }

    return 1;
}

void
MisesMatNl :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                      FloatArray &rcontrib, TimeStep *tStep)
{
    double kappa, tempKappa;
    MisesMatNlStatus *status = static_cast< MisesMatNlStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );
    FloatMatrix b;
    LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
    double E = lmat->give('E', gp);

    elem->giveLocationArray(rloc, s);
    elem->computeBmatrixAt(gp, b);
    kappa = status->giveCumulativePlasticStrain();
    tempKappa = status->giveTempCumulativePlasticStrain();

    rcontrib.clear();
    if ( ( tempKappa - kappa ) > 0 ) {
        const FloatArray &stress = status->giveTempEffectiveStress();
        if ( gp->giveMaterialMode() == _1dMat ) {
            double coeff = sgn( stress.at(1) ) * E / ( E + H );
            rcontrib.plusProduct(b, stress, coeff);
            return;
        }
    }
    rcontrib.resize(b.giveNumberOfColumns());
    rcontrib.zero();
}

void
SimpleCrossSection :: giveTemperatureVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    Element *elem = gp->giveElement();
    answer.clear();
    //sum up all prescribed temperatures over an element
    StructuralElement *selem = dynamic_cast< StructuralElement * >(elem);
    selem->computeResultingIPTemperatureAt(answer, tStep, gp, VM_Total);

    /* add external source, if provided */
    FieldManager *fm = this->domain->giveEngngModel()->giveContext()->giveFieldManager();
    FM_FieldPtr tf;

    if ( ( tf = fm->giveField(FT_Temperature) ) ) {
        // temperature field registered
        FloatArray gcoords, et2;
        int err;
        elem->computeGlobalCoordinates( gcoords, gp->giveNaturalCoordinates() );
        if ( ( err = tf->evaluateAt(et2, gcoords, VM_Total, tStep) ) ) {
            OOFEM_ERROR("tf->evaluateAt failed, element %d, error code %d", elem->giveNumber(), err);
        }

        if ( et2.isNotEmpty() ) {
            if ( answer.isEmpty() ) {
                answer = et2;
            } else {
                answer.at(1) += et2.at(1);
            }
        }
    }
}

void
FiberedCrossSection :: giveGeneralizedStress_Beam3d(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
    double fiberThick, fiberWidth, fiberZCoord, fiberYCoord;
    FloatArray fiberStrain, reducedFiberStress;
    StructuralElement *element = static_cast< StructuralElement * >( gp->giveElement() );
    FiberedCrossSectionInterface *interface;

    if ( ( interface = static_cast< FiberedCrossSectionInterface * >( element->giveInterface(FiberedCrossSectionInterfaceType) ) ) == NULL ) {
        OOFEM_ERROR("element with no fiber support encountered");
    }

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

    for ( int i = 1; i <= numberOfFibers; i++ ) {
        GaussPoint *fiberGp = this->giveSlaveGaussPoint(gp, i - 1);
        StructuralMaterial *fiberMat = static_cast< StructuralMaterial * >( domain->giveMaterial( fiberMaterials.at(i) ) );
        // the question is whether this function should exist ?
        // if yes the element details will be hidden.
        // good idea also should be existence of element::GiveBmatrixOfLayer
        // and computing strains here - but first idea looks better
        // but treating of geometric non-linearities may become more complicated
        // another approach - use several functions with assumed kinematic constraints

        // resolve current layer z-coordinate
        fiberThick  = this->fiberThicks.at(i);
        fiberWidth  = this->fiberWidths.at(i);
        fiberYCoord = fiberGp->giveNaturalCoordinate(1);
        fiberZCoord = fiberGp->giveNaturalCoordinate(2);

        interface->FiberedCrossSectionInterface_computeStrainVectorInFiber(fiberStrain, strain, fiberGp, tStep);

        fiberMat->giveRealStressVector_Fiber(reducedFiberStress, fiberGp, fiberStrain, tStep);

        // perform integration
        // 1) membrane terms N, Qz, Qy
        answer.at(1) += reducedFiberStress.at(1) * fiberWidth * fiberThick;
        answer.at(2) += reducedFiberStress.at(2) * fiberWidth * fiberThick;
        answer.at(3) += reducedFiberStress.at(3) * fiberWidth * fiberThick;
        // 2) bending terms mx, my, mxy
        answer.at(4) += ( reducedFiberStress.at(2) * fiberWidth * fiberThick * fiberYCoord -
                          reducedFiberStress.at(3) * fiberWidth * fiberThick * fiberZCoord );
        answer.at(5) += reducedFiberStress.at(1) * fiberWidth * fiberThick * fiberZCoord;
        answer.at(6) -= reducedFiberStress.at(1) * fiberWidth * fiberThick * fiberYCoord;
    }

    // now we must update master gp ///@ todo simply chosen the first fiber material as master material /JB
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >
                                       ( domain->giveMaterial( fiberMaterials.at(1) )->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}

void GradDpElement :: computeNonlocalDegreesOfFreedom(FloatArray &answer, TimeStep *tStep)
{
    ///@todo Just use this instead! (should work, but I don't have any tests right now)
#if 0
    this->giveStructuralElement()->computeVectorOf({G_0}, VM_Total, tStep, answer);
#else
    StructuralElement *elem = this->giveStructuralElement();
    FloatArray u;
    answer.resize(nlSize);

    elem->computeVectorOf(VM_Total, tStep, u);
    u.resizeWithValues(totalSize);
    for  ( int i = 1; i <= nlSize; i++ ) {
        answer.at(i) = u.at( locK.at(i) );
    }
#endif
}

int
TrabBoneNL3D :: giveLocalNonlocalStiffnessContribution(GaussPoint *gp, IntArray &loc, const UnknownNumberingScheme &s,
                                                       FloatArray &lcontrib, TimeStep *tStep)
{
    TrabBoneNL3DStatus *nlStatus = static_cast< TrabBoneNL3DStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );

    int nrows, nsize;
    double sum, nlKappa, dDamFunc, dam, tempDam;
    FloatArray localNu;
    FloatMatrix b;

    this->computeCumPlastStrain(nlKappa, gp, tStep);
    dam = nlStatus->giveDam();
    tempDam = nlStatus->giveTempDam();

    if ( ( tempDam - dam ) > 0.0 ) {
        elem->giveLocationArray(loc, s);
        localNu = nlStatus->giveTempEffectiveStress();

        elem->giveLocationArray(loc, EModelDefaultEquationNumbering() );
        elem->computeBmatrixAt(gp, b);
        dDamFunc = expDam * critDam * exp(-expDam * nlKappa);

        nrows = b.giveNumberOfColumns();
        nsize = localNu.giveSize();
        lcontrib.resize(nrows);
        for ( int i = 1; i <= nrows; i++ ) {
            sum = 0.0;
            for ( int j = 1; j <= nsize; j++ ) {
                sum += b.at(j, i) * localNu.at(j);
            }

            lcontrib.at(i) = mParam * dDamFunc * sum;
        }

        return 1;
    } else {
        loc.clear();
        return 0;
    }
}

void
MisesMatNl :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
        FloatArray &rcontrib, TimeStep *atTime)
{
    int ncols, nsize, i, j;
    double sum, kappa, tempKappa;
    MisesMatNlStatus *status = ( MisesMatNlStatus * ) this->giveStatus(gp);
    StructuralElement *elem = ( StructuralElement * )( gp->giveElement() );
    FloatMatrix b;
    FloatArray stress;
    LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
    double E = lmat->give('E', gp);

    elem->giveLocationArray(rloc, EID_MomentumBalance, s);
    elem->computeBmatrixAt(gp, b);
    ncols = b.giveNumberOfColumns();
    rcontrib.resize(ncols);
    kappa = status->giveCumulativePlasticStrain();
    tempKappa = status->giveTempCumulativePlasticStrain();

    if ( ( tempKappa - kappa ) > 0 ) {
        status->giveTempEffectiveStress(stress);
        if ( gp->giveMaterialMode() == _1dMat ) {
            nsize = stress.giveSize();
            double coeff = sgn( stress.at(1) ) * E / ( E + H );
            for ( i = 1; i <= ncols; i++ ) {
                sum = 0.;
                for ( j = 1; j <= nsize; j++ ) {
                    sum += stress.at(j) * coeff * b.at(j, i);
                }

                rcontrib.at(i) = sum;
            }
        }
    } else {
        rcontrib.zero();
    }
}

// computes m*gprime*Btransposed*sigmaeff for the given Gauss point
// and returns 0 of damage is not growing, 1 if it is growing
// (if damage is not growing, the contribution is not considered at all)
int
RankineMatNl :: giveLocalNonlocalStiffnessContribution(GaussPoint *gp, IntArray &loc, const UnknownNumberingScheme &s,
                                                       FloatArray &lcontrib, TimeStep *atTime)
{
    int nrows, nsize, i, j;
    double sum, nlKappa, damage, tempDamage;
    RankineMatNlStatus *status = ( RankineMatNlStatus * ) this->giveStatus(gp);
    StructuralElement *elem = ( StructuralElement * )( gp->giveElement() );
    FloatMatrix b;
    FloatArray stress;

    damage = status->giveDamage();
    tempDamage = status->giveTempDamage();
    if ( tempDamage <= damage ) {
        return 0; // no contribution if damage is not growing
    }

    elem->giveLocationArray(loc, EID_MomentumBalance, s);
    status->giveTempEffectiveStress(stress);
    elem->computeBmatrixAt(gp, b);
    this->computeCumPlasticStrain(nlKappa, gp, atTime);
    double factor = computeDamageParamPrime(nlKappa);
    factor *= mm; // this factor is m*gprime
    nrows = b.giveNumberOfColumns();
    nsize = stress.giveSize();
    lcontrib.resize(nrows);
    // compute the product Btransposed*stress and multiply by factor
    for ( i = 1; i <= nrows; i++ ) {
        sum = 0.0;
        for ( j = 1; j <= nsize; j++ ) {
            sum += b.at(j, i) * stress.at(j);
        }

        lcontrib.at(i) = sum * factor;
    }

    return 1; // contribution will be considered
}

int
MisesMatNl :: giveLocalNonlocalStiffnessContribution(GaussPoint *gp, IntArray &loc, const UnknownNumberingScheme &s,
                                                     FloatArray &lcontrib, TimeStep *tStep)
{
    double nlKappa, damage, tempDamage, dDamF;
    MisesMatNlStatus *status = static_cast< MisesMatNlStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );
    FloatMatrix b;

    this->computeCumPlasticStrain(nlKappa, gp, tStep);
    damage = status->giveDamage();
    tempDamage = status->giveTempDamage();
    if ( ( tempDamage - damage ) > 0 ) {
        const FloatArray &stress = status->giveTempEffectiveStress();
        elem->giveLocationArray(loc, s);
        elem->computeBmatrixAt(gp, b);
        dDamF = computeDamageParamPrime(nlKappa);
        lcontrib.clear();
        lcontrib.plusProduct(b, stress, mm * dDamF);
    }

    return 1;
}

 MorphoFilter<T>::MorphoFilter(const StructuralElement& element)
   : NonLinearFilter<T>(&(element.getKernel()[0]), element.getWidth(), element.getHeight(), element.getDepth())
 {}

void LSPrimaryVariableMapper :: mapPrimaryVariables(FloatArray &oU, Domain &iOldDom, Domain &iNewDom, ValueModeType iMode, TimeStep &iTStep)
{
    EngngModel *engngMod = iNewDom.giveEngngModel();
    EModelDefaultEquationNumbering num;


    const int dim = iNewDom.giveNumberOfSpatialDimensions();

    int numElNew = iNewDom.giveNumberOfElements();

    // Count dofs
    int numDofsNew = engngMod->giveNumberOfDomainEquations( 1, num );


    oU.resize(numDofsNew);
    oU.zero();

    FloatArray du(numDofsNew);
    du.zero();

    FloatArray res(numDofsNew);

    std :: unique_ptr< SparseMtrx > K;
    std :: unique_ptr< SparseLinearSystemNM > solver;

    solver.reset( classFactory.createSparseLinSolver(ST_Petsc, & iOldDom, engngMod) );
    if (!solver) {
        solver.reset( classFactory.createSparseLinSolver(ST_Direct, & iOldDom, engngMod) );
    }
    K.reset( classFactory.createSparseMtrx(solver->giveRecommendedMatrix(true)) );

    K->buildInternalStructure( engngMod, iNewDom.giveNumber(), num );

    int maxIter = 1;

    for ( int iter = 0; iter < maxIter; iter++ ) {
        K->zero();
        res.zero();


        // Contribution from elements
        for ( int elIndex = 1; elIndex <= numElNew; elIndex++ ) {
            StructuralElement *elNew = dynamic_cast< StructuralElement * >( iNewDom.giveElement(elIndex) );
            if ( elNew == NULL ) {
                OOFEM_ERROR("Failed to cast Element new to StructuralElement.");
            }

            ///////////////////////////////////
            // Compute residual

            // Count element dofs
            int numElNodes = elNew->giveNumberOfDofManagers();
            int numElDofs = 0;
            for ( int i = 1; i <= numElNodes; i++ ) {
                numElDofs += elNew->giveDofManager(i)->giveNumberOfDofs();
            }

            FloatArray elRes(numElDofs);
            elRes.zero();

            IntArray elDofsGlob;
            elNew->giveLocationArray( elDofsGlob, num );


            // Loop over Gauss points
            for ( int intRuleInd = 0; intRuleInd < elNew->giveNumberOfIntegrationRules(); intRuleInd++ ) {

                for ( GaussPoint *gp: *elNew->giveIntegrationRule(intRuleInd) ) {

                    // New N-matrix
                    FloatMatrix NNew;
                    elNew->computeNmatrixAt(gp->giveNaturalCoordinates(), NNew);


                    //////////////
                    // Global coordinates of GP
                    const FloatArray &localCoord = gp->giveNaturalCoordinates();
                    FloatArray globalCoord;
                    elNew->computeGlobalCoordinates(globalCoord, localCoord);
                    //////////////


                    // Localize element and point in the old domain
                    FloatArray localCoordOld(dim), pointCoordOld(dim);
                    StructuralElement *elOld = dynamic_cast< StructuralElement * >( iOldDom.giveSpatialLocalizer()->giveElementClosestToPoint(localCoordOld, pointCoordOld, globalCoord, 0) );
                    if ( elOld == NULL ) {
                        OOFEM_ERROR("Failed to cast Element old to StructuralElement.");
                    }


                    // Compute N-Matrix for the old element
                    FloatMatrix NOld;
                    elOld->computeNmatrixAt(localCoordOld, NOld);

                    // Fetch nodal displacements for the new element
                    FloatArray nodeDispNew( elDofsGlob.giveSize() );


                    int dofsPassed = 1;
                    for ( int elNode: elNew->giveDofManArray() ) {
//.........这里部分代码省略.........