C++ FEIElementGeometryWrapper函数代码示例

// Edge load support
void
DKTPlate :: 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);
}

void
DKTPlate :: 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
    }
}

void
Structural3DElement :: 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);
}

void
Quad10_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);
    }
}

void
Quad10_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);

}

void
QTrPlaneStressGrad :: 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);
    }
}

void
Q4Axisymm :: 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);
    }
}

double
L4Axisymm :: 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;
}

void
StructuralInterfaceElementPhF :: 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 );
    
    
}

void
QTrPlaneStressGrad :: 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);
    }
}

void
IntElLine1PhF :: 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;
    }
}

double
Q4Axisymm :: 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;
}