自学教程：C++ FEIElementGeometryWrapper函数代码示例

// Edge load supportvoidDKTPlate :: computeEgdeNMatrixAt(FloatMatrix &answer, int iedge, GaussPoint *gp){    FloatArray n;    this->interp_lin.edgeEvalN( n, iedge, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );    answer.resize(3, 6);    answer.at(1, 1) = n.at(1);    answer.at(1, 4) = n.at(2);    answer.at(2, 2) = answer.at(3, 3) = n.at(1);    answer.at(2, 5) = answer.at(3, 6) = n.at(2);}

void Tet21Stokes :: computeInternalForcesVector(FloatArray &answer, TimeStep *tStep){    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();    FloatArray a_pressure, a_velocity, devStress, epsp, Nh, dN_V(30);    FloatMatrix dN, B(6, 30);    double r_vol, pressure;    this->computeVectorOfVelocities(VM_Total, tStep, a_velocity);    this->computeVectorOfPressures(VM_Total, tStep, a_pressure);    FloatArray momentum, conservation;    B.zero();    for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {        const FloatArray &lcoords = gp->giveNaturalCoordinates();        double detJ = fabs( this->interpolation_quad.evaldNdx( dN, lcoords, FEIElementGeometryWrapper(this) ) );        this->interpolation_lin.evalN( Nh, lcoords, FEIElementGeometryWrapper(this) );        double dV = detJ * gp->giveWeight();        for ( int j = 0, k = 0; j < dN.giveNumberOfRows(); j++, k += 3 ) {            dN_V(k + 0) = B(0, k + 0) = B(3, k + 1) = B(4, k + 2) = dN(j, 0);            dN_V(k + 1) = B(1, k + 1) = B(3, k + 0) = B(5, k + 2) = dN(j, 1);            dN_V(k + 2) = B(2, k + 2) = B(4, k + 0) = B(5, k + 1) = dN(j, 2);        }        epsp.beProductOf(B, a_velocity);        pressure = Nh.dotProduct(a_pressure);        mat->computeDeviatoricStressVector(devStress, r_vol, gp, epsp, pressure, tStep);        momentum.plusProduct(B, devStress, dV);        momentum.add(-pressure * dV, dN_V);        conservation.add(r_vol * dV, Nh);    }    answer.resize(34);    answer.zero();    answer.assemble(momentum, this->momentum_ordering);    answer.assemble(conservation, this->conservation_ordering);}

voidDKTPlate :: computeShearForces(FloatArray &answer, GaussPoint *gp, TimeStep *tStep){    // as shear strains are enforced to be zero (art least on element edges) the shear forces are computed from equlibrium    FloatMatrix m, dndx;    answer.resize(5);    this->computeVertexBendingMoments(m, tStep);    this->interp_lin.evaldNdx( dndx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );    for ( int i = 1; i <= this->numberOfDofMans; i++ ) {        answer.at(4) += m.at(1, i) * dndx.at(i, 1) + m.at(3, i) * dndx.at(i, 2); //dMxdx + dMxydy        answer.at(5) += m.at(2, i) * dndx.at(i, 2) + m.at(3, i) * dndx.at(i, 1); //dMydy + dMxydx    }}

voidStructural3DElement :: computeSurfaceNMatrixAt(FloatMatrix &answer, int iSurf, GaussPoint *sgp){    /* Returns the [ 3 x (nno*3) ] shape function matrix {N} of the receiver,      * evaluated at the given gp.     * {u} = {N}*{a} gives the displacements at the integration point.     */               // Evaluate the shape functions at the position of the gp.     FloatArray N;    static_cast< FEInterpolation3d* > ( this->giveInterpolation() )->        surfaceEvalN( N, iSurf, sgp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );      answer.beNMatrixOf(N, 3);}

voidQuad10_2D_SUPG :: computeDivUMatrix(FloatMatrix &answer, GaussPoint *gp){    FloatMatrix dn;    velocityInterpolation.evaldNdx( dn, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );    answer.resize(1, 8);    answer.zero();    for ( int i = 1; i <= 4; i++ ) {        answer.at(1, 2 * i - 1) = dn.at(i, 1);        answer.at(1, 2 * i) = dn.at(i, 2);    }}

voidQuad10_2D_SUPG :: computeNuMatrix(FloatMatrix &answer, GaussPoint *gp){    FloatArray n;    this->velocityInterpolation.evalN(n, * gp->giveCoordinates(), FEIElementGeometryWrapper(this));    answer.resize(2, 8);    answer.zero();    for ( int i = 1; i <= 4; i++ ) {        answer.at(1, 2 * i - 1) = n.at(i);        answer.at(2, 2 * i - 0) = n.at(i);    }}

void Tr21Stokes :: giveElementFMatrix(FloatMatrix &answer){    IntegrationRule *iRule = integrationRulesArray [ 0 ];    GaussPoint *gp;    double detJ;    FloatArray N, N2, *lcoords;    IntArray col;    FloatMatrix temp;    N2.resize(6);   N2.zero();    col.resize(2);  col.at(1)=1;  col.at(2)=2;    temp.resize(15,2);    temp.zero();    for (int i=0; i<iRule->getNumberOfIntegrationPoints(); i++) {        gp = iRule->getIntegrationPoint(i);        lcoords = gp->giveCoordinates();        this->interpolation_quad.evalN(N, *lcoords, FEIElementGeometryWrapper(this));        detJ = this->interpolation_quad.giveTransformationJacobian(*lcoords, FEIElementGeometryWrapper(this));        N.times(gp->giveWeight()*detJ);        //N.printYourself();        N2.add(N);    }    for (int i=1; i<=6; i++) {        temp.at(i*2-1, 1)=N2.at(i);        temp.at(i*2, 2)=N2.at(i);    }    answer.resize(17,2);    answer.zero();    answer.assemble(temp, this->ordering, col);}

voidQTrPlaneStressGrad :: computeNkappaMatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer)// Returns the displacement interpolation matrix {N} of the receiver, eva-// luated at aGaussPoint.{    FloatArray n(3);    answer.resize(1, 3);    answer.zero();    this->interpolation.evalN(n, * aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this));    for ( int i = 1; i <= 3; i++ ) {        answer.at(1, i) = n.at(i);    }}

voidQ4Axisymm :: computeNmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer)// Returns the displacement interpolation matrix {N} of the receiver,// evaluated at aGaussPoint.{    FloatArray n;    this->interp.evalN(n, *aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this));    answer.resize(2, 16);    answer.zero();    for ( int i = 1; i <= 8; i++ ) {        answer.at(1, 2 * i - 1) = n.at(i);        answer.at(2, 2 * i - 0) = n.at(i);    }}

doubleL4Axisymm :: computeVolumeAround(GaussPoint *aGaussPoint)// Returns the portion of the receiver which is attached to aGaussPoint.{    int i;    double determinant, weight, volume, r, x;    FloatArray n(4);    this->interpolation.evalN( n, * aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this) );    r = 0.;    for ( i = 1; i <= numberOfDofMans; i++ ) {        x  = this->giveNode(i)->giveCoordinate(1);        r += x * n.at(i);    }    determinant = fabs( this->interpolation.giveTransformationJacobian( * aGaussPoint->giveCoordinates(),                                                                       FEIElementGeometryWrapper(this) ) );    weight = aGaussPoint->giveWeight();    volume = determinant * weight * r;    return volume;}

voidStructuralInterfaceElementPhF :: computeBd_matrixAt(GaussPoint *gp, FloatMatrix &answer){    FloatMatrix dNdxi;    FloatMatrix G;    this->computeCovarBaseVectorsAt(gp, G);        this->giveInterpolation( )->evaldNdxi( dNdxi, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper( this ) );    answer.beProductTOf(G,dNdxi);    //answer.beTranspositionOf( dNdx );        }

voidQTrPlaneStressGrad :: computeBkappaMatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer){    FloatMatrix dnx;    this->interpolation.evaldNdx(dnx, * aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this));    answer.resize(2, 3);    answer.zero();    for ( int i = 1; i <= 3; i++ ) {        answer.at(1, i) = dnx.at(i, 1);        answer.at(2, i) = dnx.at(i, 2);    }}

voidIntElLine1PhF :: computeCovarBaseVectorAt(IntegrationPoint *ip, FloatArray &G){    FloatMatrix dNdxi;    FEInterpolation *interp = this->giveInterpolation();    interp->evaldNdxi( dNdxi, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );    G.resize(2);    G.zero();    int numNodes = this->giveNumberOfNodes();    for ( int i = 1; i <= dNdxi.giveNumberOfRows(); i++ ) {        double X1_i = 0.5 * ( this->giveNode(i)->giveCoordinate(1) + this->giveNode(i + numNodes / 2)->giveCoordinate(1) ); // (mean) point on the fictious mid surface        double X2_i = 0.5 * ( this->giveNode(i)->giveCoordinate(2) + this->giveNode(i + numNodes / 2)->giveCoordinate(2) );        G.at(1) += dNdxi.at(i, 1) * X1_i;        G.at(2) += dNdxi.at(i, 1) * X2_i;    }}

doubleQ4Axisymm :: computeVolumeAround(GaussPoint *aGaussPoint)// Returns the portion of the receiver which is attached to aGaussPoint.{    FloatArray n;    double determinant, r;    this->interp.evalN(n, *aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this));    for ( int i = 1; i <= 8; i++ ) {        r += this->giveNode(i)->giveCoordinate(1) * n.at(i);    }    determinant = fabs( this->interp.giveTransformationJacobian(*aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this)) );    return determinant * aGaussPoint->giveWeight() * r;}

double MixedGradientPressureBC :: domainSize(){    int nsd = this->domain->giveNumberOfSpatialDimensions();    double domain_size = 0.0;    // This requires the boundary to be consistent and ordered correctly.    Set *set = this->giveDomain()->giveSet(this->set);    const IntArray &boundaries = set->giveBoundaryList();    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {        Element *e = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );        int boundary = boundaries.at(pos * 2);        FEInterpolation *fei = e->giveInterpolation();        domain_size += fei->evalNXIntegral( boundary, FEIElementGeometryWrapper(e) );    }    return domain_size / nsd;}

double TransportGradientPeriodic :: domainSize(Domain *d, int setNum){    int nsd = d->giveNumberOfSpatialDimensions();    double domain_size = 0.0;    // This requires the boundary to be consistent and ordered correctly.    Set *set = d->giveSet(setNum);    const IntArray &boundaries = set->giveBoundaryList();    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {        Element *e = d->giveElement( boundaries.at(pos * 2 - 1) );        int boundary = boundaries.at(pos * 2);        FEInterpolation *fei = e->giveInterpolation();        domain_size += fei->evalNXIntegral( boundary, FEIElementGeometryWrapper(e) );    }    return fabs(domain_size / nsd);}

void Tet21Stokes :: computeBoundaryLoadVector(FloatArray &answer, BoundaryLoad *load, int iSurf, CharType type, ValueModeType mode, TimeStep *tStep){    if ( type != ExternalForcesVector ) {        answer.clear();        return;    }    answer.resize(34);    answer.zero();    if ( load->giveType() == TransmissionBC ) { // Neumann boundary conditions (traction)        int numberOfSurfaceIPs = ( int ) ceil( ( load->giveApproxOrder() + 1. ) / 2. ) * 2; ///@todo Check this.        GaussIntegrationRule iRule(1, this, 1, 1);        FloatArray N, t, f(18);        f.zero();        iRule.SetUpPointsOnTriangle(numberOfSurfaceIPs, _Unknown);        for ( GaussPoint *gp: iRule ) {            const FloatArray &lcoords = gp->giveNaturalCoordinates();            this->interpolation_quad.surfaceEvalN( N, iSurf, lcoords, FEIElementGeometryWrapper(this) );            double dA = gp->giveWeight() * this->interpolation_quad.surfaceGiveTransformationJacobian( iSurf, lcoords, FEIElementGeometryWrapper(this) );            if ( load->giveFormulationType() == Load :: FT_Entity ) { // load in xi-eta system                load->computeValueAt(t, tStep, lcoords, VM_Total);            } else { // Edge load in x-y system                FloatArray gcoords;                this->interpolation_quad.surfaceLocal2global( gcoords, iSurf, lcoords, FEIElementGeometryWrapper(this) );                load->computeValueAt(t, tStep, gcoords, VM_Total);            }            // Reshape the vector            for ( int j = 0; j < N.giveSize(); j++ ) {                f(3 * j + 0) += N(j) * t(0) * dA;                f(3 * j + 1) += N(j) * t(1) * dA;                f(3 * j + 2) += N(j) * t(2) * dA;            }        }        answer.assemble(f, this->surf_ordering [ iSurf - 1 ]);    } else {        OOFEM_ERROR("Strange boundary condition type");    }}

voidIntElLine1PhF :: computeNmatrixAt(GaussPoint *ip, FloatMatrix &answer){    // Returns the modified N-matrix which multiplied with u give the spatial jump.    FloatArray N;    FEInterpolation *interp = this->giveInterpolation();    interp->evalN( N, ip->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );    answer.resize(2, 8);    answer.zero();    answer.at(1, 1) = answer.at(2, 2) = -N.at(1);    answer.at(1, 3) = answer.at(2, 4) = -N.at(2);    answer.at(1, 5) = answer.at(2, 6) = N.at(1);    answer.at(1, 7) = answer.at(2, 8) = N.at(2);}

intInterfaceElement3dTrLin :: computeGlobalCoordinates(FloatArray &answer, const FloatArray &lcoords){    FloatArray n;    this->interpolation.evalN( n, lcoords, FEIElementGeometryWrapper(this) );    answer.resize(3);    answer.zero();    for ( int i = 1; i <= 3; i++ ) {        answer.at(1) += n.at(i) * this->giveNode(i)->giveCoordinate(1);        answer.at(2) += n.at(i) * this->giveNode(i)->giveCoordinate(2);        answer.at(3) += n.at(i) * this->giveNode(i)->giveCoordinate(3);    }    return 1;}

voidDKTPlate :: computeNmatrixAt(const FloatArray &iLocCoord, FloatMatrix &answer)// Returns the [3x9] displacement interpolation matrix {N} of the receiver,// evaluated at gp.// Note: this interpolation is not available, as the deflection is cubic along the edges,//       but not define in the interior of the element// Note: the interpolation of rotations is quadratic// NOTE: linear interpolation returned instead{    FloatArray N;    answer.resize(3, 9);    answer.zero();    giveInterpolation()->evalN( N, iLocCoord, FEIElementGeometryWrapper(this) );    answer.beNMatrixOf(N, 3);}

void Line2BoundaryElement :: EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(ValueModeType mode,                                                                                   TimeStep *tStep, const FloatArray &lcoords, FloatArray &answer){    FloatArray n;    this->fei.evalN( answer, lcoords, FEIElementGeometryWrapper(this) );    IntArray dofIDs;    this->giveElementDofIDMask(dofIDs);    answer.resize( dofIDs.giveSize() );    answer.zero();    for ( int i = 1; i <= n.giveSize(); ++i ) {        for ( int j = 1; j <= dofIDs.giveSize(); ++j ) {            answer.at(j) += n.at(i) * this->giveNode(i)->giveDofWithID( dofIDs.at(j) )->giveUnknown(mode, tStep);        }    }}

  void Tr21Stokes :: computeBoundarySurfaceLoadVector(FloatArray &answer, BoundaryLoad *load, int boundary, CharType type, ValueModeType mode, TimeStep *tStep, bool global){    if ( type != ExternalForcesVector ) {        answer.clear();        return;    }    if ( load->giveType() == TransmissionBC ) { // Neumann boundary conditions (traction)        int numberOfEdgeIPs = ( int ) ceil( ( load->giveApproxOrder() + 1. ) / 2. ) * 2;        GaussIntegrationRule iRule(1, this, 1, 1);        FloatArray N, t, f(6);        f.zero();        iRule.SetUpPointsOnLine(numberOfEdgeIPs, _Unknown);        for ( GaussPoint *gp: iRule ) {            const FloatArray &lcoords = gp->giveNaturalCoordinates();            this->interpolation_quad.edgeEvalN( N, boundary, lcoords, FEIElementGeometryWrapper(this) );            double detJ = fabs( this->interpolation_quad.boundaryGiveTransformationJacobian( boundary, lcoords, FEIElementGeometryWrapper(this) ) );            double dS = gp->giveWeight() * detJ;            if ( load->giveFormulationType() == Load :: FT_Entity ) { // Edge load in xi-eta system                load->computeValueAt(t, tStep, lcoords, VM_Total);            } else { // Edge load in x-y system                FloatArray gcoords;                this->interpolation_quad.boundaryLocal2Global( gcoords, boundary, lcoords, FEIElementGeometryWrapper(this) );                load->computeValueAt(t, tStep, gcoords, VM_Total);            }            // Reshape the vector            for ( int j = 0; j < 3; j++ ) {                f(2 * j)     += N(j) * t(0) * dS;                f(2 * j + 1) += N(j) * t(1) * dS;            }        }        answer.resize(15);        answer.zero();        answer.assemble(f, this->edge_ordering [ boundary - 1 ]);    } else {        OOFEM_ERROR("Strange boundary condition type");    }}

voidQTRSpace :: computeNmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer)// Returns the displacement interpolation matrix {N} of the receiver, eva-// luated at aGaussPoint.{    FloatArray n(10);    answer.resize(3, 30);    answer.zero();    this->interpolation.evalN(n, * aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this));    for ( int i = 1; i <= 10; i++ ) {        answer.at(1, 3 * i - 2) = n.at(i);        answer.at(2, 3 * i - 1) = n.at(i);        answer.at(3, 3 * i - 0) = n.at(i);    }}

voidTruss3d :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)//// Returns linear part of geometrical equations of the receiver at gp.// Returns the linear part of the B matrix//{    FloatMatrix dN;    this->interp.evaldNdx(dN, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));    answer.resize(1, 6);    answer.at(1, 1) = dN.at(1,1);    answer.at(1, 2) = dN.at(1,2);    answer.at(1, 3) = dN.at(1,3);    answer.at(1, 4) = dN.at(2,1);    answer.at(1, 5) = dN.at(2,2);    answer.at(1, 6) = dN.at(2,3);}

voidQTrPlaneStrain :: computeNmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer)// Returns the displacement interpolation matrix {N} of the receiver, eva-// luated at aGaussPoint.{    int i;    FloatArray n(6);    answer.resize(2, 12);    answer.zero();    this->interpolation.evalN( n, * aGaussPoint->giveCoordinates(), FEIElementGeometryWrapper(this) );    for ( i = 1; i <= 6; i++ ) {        answer.at(1, 2 * i - 1) = n.at(i);        answer.at(2, 2 * i - 0) = n.at(i);    }}

voidLineDistributedSpring :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)// Returns the [3x3] strain-displacement matrix {B} of the receiver,// evaluated at gp.{    FloatArray n;    int ndofs = this->dofs.giveSize();    this->interp_lin.evalN( n, gp->giveNaturalCoordinates(),  FEIElementGeometryWrapper(this) );    answer.resize(ndofs, ndofs*2);    answer.zero();    for (int idof=1; idof<=ndofs; idof++) {      answer.at(idof, idof) = n.at(1);       answer.at(idof, ndofs+idof) = n.at(2);     }}

voidQuad10_2D_SUPG :: computeUDotGradUMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep){    FloatMatrix n, dn;    FloatArray u, un;    this->velocityInterpolation.evaldNdx( dn, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );    this->computeNuMatrix(n, gp);    this->computeVectorOfVelocities(VM_Total, tStep, un);    u.beProductOf(n, un);    answer.resize(2, 8);    answer.zero();    for ( int i = 1; i <= 4; i++ ) {        answer.at(1, 2 * i - 1) = dn.at(i, 1) * u.at(1) + dn.at(i, 2) * u.at(2);        answer.at(2, 2 * i)     = dn.at(i, 1) * u.at(1) + dn.at(i, 2) * u.at(2);    }}