C++ GaussPoint::giveNaturalCoordinates方法代码示例

bool Inclusion :: isMaterialModified(GaussPoint &iGP, Element &iEl, CrossSection * &opCS) const
{
    // Check if the point is located inside the inclusion

    FloatArray N;
    FEInterpolation *interp = iEl.giveInterpolation();
    interp->evalN( N, * iGP.giveNaturalCoordinates(), FEIElementGeometryWrapper(& iEl) );

    const IntArray &elNodes = iEl.giveDofManArray();

    double levelSetGP = 0.0;
    interpLevelSet(levelSetGP, N, elNodes);

    if ( levelSetGP < 0.0 ) {
        opCS = mpCrossSection;
        return true;
    }

    return false;
}

void
tet21ghostsolid :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{

    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix Kf, G, Kx, D, B, Ed, EdB, dNx;
    FloatArray Nlin, dNv;

    for (int j = 0; j<iRule->giveNumberOfIntegrationPoints(); j++) {
        GaussPoint *gp = iRule->getIntegrationPoint(j);

        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30); // dNv = [dN1/dx dN1/dy dN1/dz dN2/dx dN2/dy dN2/dz ... dN10/dz]

        for (int k = 0; k<dNx.giveNumberOfRows(); k++) {
            dNv.at(k*3+1) = dNx.at(k+1,1);
            dNv.at(k*3+2) = dNx.at(k+1,2);
            dNv.at(k*3+3) = dNx.at(k+1,3);
        }

        if (nlGeometry == 0) {

            this->computeBmatrixAt(gp, B);

            // Fluid part
            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->giveDeviatoricStiffnessMatrix(Ed, TangentStiffness, gp, tStep);
#else
            OOFEM_ERROR("Fluid module missing\n");
#endif
            gp->setMaterialMode(_3dMat);

            EdB.beProductOf(Ed, B);
            Kf.plusProductSymmUpper(B, EdB, detJ*weight);

            // Ghost solid part
            EdB.beProductOf(Dghost, B);
            Kx.plusProductSymmUpper(B, EdB, detJ*weight);

            // Incompressibility part
            G.plusDyadUnsym(dNv, Nlin, -detJ*weight);

        } else {
            OOFEM_ERROR ("No support for large deformations yet!");
        }

    }

    FloatMatrix GT;

    GT.beTranspositionOf(G);
    //GTdeltat.beTranspositionOf(G);
    //GTdeltat.times(deltat);
    Kf.symmetrized();
    Kx.symmetrized();
    //    Kf.printYourself();
    //    G.printYourself();
    //    GT.printYourself();
    //    Kx.printYourself();

    answer.resize(64, 64);
    answer.zero();
#define USEUNCOUPLED 0

#if USEUNCOUPLED == 1
    // Totaly uncoupled
    answer.assemble(Kf, momentum_ordering, momentum_ordering);
    answer.assemble(G, momentum_ordering, conservation_ordering);
    answer.assemble(GT, conservation_ordering, momentum_ordering);
    answer.assemble(Kx, ghostdisplacement_ordering, ghostdisplacement_ordering);
#else
    answer.assemble(Kf, ghostdisplacement_ordering, ghostdisplacement_ordering);
    answer.assemble(Kf, ghostdisplacement_ordering, momentum_ordering);
    answer.assemble(G, ghostdisplacement_ordering, conservation_ordering);
    answer.assemble(GT, conservation_ordering, ghostdisplacement_ordering);
    answer.assemble(GT, conservation_ordering, momentum_ordering);
    answer.assemble(Kx, momentum_ordering, ghostdisplacement_ordering);
#endif

    //answer.printYourself();

}

void
tet21ghostsolid :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{

    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix Kf, G, Kx, B, Ed, dNx;
    FloatArray Strain, Stress, Nlin, dNv, a, aVelocity, aPressure, aGhostDisplacement, fluidStress, epsf;
    FloatArray momentum, conservation, auxstress;
    double pressure, epsvol;

    this->computeVectorOf( VM_Total, tStep, a);
    if (!tStep->isTheFirstStep()) {
        // a.printYourself();
    }

    aVelocity.beSubArrayOf(a, momentum_ordering);
    aPressure.beSubArrayOf(a, conservation_ordering);
    aGhostDisplacement.beSubArrayOf(a, ghostdisplacement_ordering);

    for (int j = 0; j<iRule->giveNumberOfIntegrationPoints(); j++) {
        GaussPoint *gp = iRule->getIntegrationPoint(j);

        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for (int k = 0; k<dNx.giveNumberOfColumns(); k++) {
            dNv.at(k*3+1) = dNx.at(1,k+1);
            dNv.at(k*3+2) = dNx.at(2,k+1);
            dNv.at(k*3+3) = dNx.at(3,k+1);
        }

        if (nlGeometry == 0) {

            this->computeBmatrixAt(gp, B);
            epsf.beProductOf(B, aVelocity);
            pressure = Nlin.dotProduct(aPressure);

            // Fluid part
            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->computeDeviatoricStressVector(fluidStress, epsvol, gp, epsf, pressure, tStep);
#else
            OOFEM_ERROR("Missing FM module");
#endif
            gp->setMaterialMode(_3dMat);

            momentum.plusProduct(B, fluidStress, detJ*weight);
            momentum.add(-pressure * detJ * weight, dNv);
            conservation.add(epsvol * detJ * weight, Nlin);

            // Ghost solid part
            Strain.beProductOf(B, aGhostDisplacement);
            Stress.beProductOf(Dghost, Strain);
            auxstress.plusProduct(B, Stress, detJ * weight);

        } else {
            OOFEM_ERROR("No support for large deformations yet!");
        }

    }
    answer.resize(64);
    answer.zero();

#if USEUNCOUPLED == 1
    // Totaly uncoupled
    answer.assemble(momentum, momentum_ordering);
    answer.assemble(conservation, conservation_ordering);
    answer.assemble(auxstress, ghostdisplacement_ordering);
#else
    answer.assemble(momentum, ghostdisplacement_ordering);
    answer.assemble(conservation, conservation_ordering);
    answer.assemble(auxstress, momentum_ordering);
#endif

    // Test linear

    /*    if (this->giveNumber() == 364) {
        FloatMatrix K;
        FloatArray ans;
        this->computeStiffnessMatrix(K, TangentStiffness, tStep);
        ans.beProductOf(K, a);
        ans.printYourself();
        answer.printYourself();
    } */


}

void
tet21ghostsolid :: computeLoadVector(FloatArray &answer, Load *load, CharType type, ValueModeType mode, TimeStep *tStep)
{

    answer.resize(64);
    answer.zero();

    if ( type != ExternalForcesVector ) {
        answer.resize(0);
        return;
    }

#ifdef __FM_MODULE
    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    FloatArray N, gVector, temparray(30), dNv, u, inc, u_prev, vload;
    FloatMatrix dNx, G;

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();

    vload.resize(4);
    vload.zero();

    this->giveDisplacementsIncrementData(u_prev, u, inc, tStep);

    for ( int k = 0; k < iRule->giveNumberOfIntegrationPoints(); k++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(k);
        FloatArray *lcoords = gp->giveNaturalCoordinates();

        double detJ = fabs( this->interpolation.giveTransformationJacobian( * lcoords, FEIElementGeometryWrapper(this) ) );
        double dA = detJ * gp->giveWeight();

        // Body load
        if ( gVector.giveSize() ) {
#ifdef __FM_MODULE
            double rho = mat->give('d', gp);
#else
            OOFEM_ERROR("Missing FM module");
            double rho = 1.0;
#endif
            this->interpolation.evalN( N, * lcoords, FEIElementGeometryWrapper(this) );

            for ( int j = 0; j < N.giveSize(); j++ ) {
                temparray(3 * j + 0) += N(j) * rho * gVector(0) * dA;
                temparray(3 * j + 1) += N(j) * rho * gVector(1) * dA;
                temparray(3 * j + 2) += N(j) * rho * gVector(2) * dA;
            }
        }

        // "load" from previous step
        this->interpolation.evaldNdx( dNx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( N, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for (int k = 0; k<dNx.giveNumberOfColumns(); k++) {
            dNv.at(k*3+1) = dNx.at(1,k+1);
            dNv.at(k*3+2) = dNx.at(2,k+1);
            dNv.at(k*3+3) = dNx.at(3,k+1);
        }

        G.plusDyadUnsym(N, dNv, dA);
        FloatMatrix GT;
        GT.beTranspositionOf(G);

        vload.plusProduct(GT, u_prev, -0.0 );
        //vload.printYourself();

    }

#if USEUNCOUPLED == 1
    // Totaly uncoupled
    answer.assemble(temparray, this->momentum_ordering);
#else
    answer.assemble(temparray, this->ghostdisplacement_ordering);
    answer.assemble(vload, this->conservation_ordering);
#endif



    // answer.printYourself();
}

int
PatchIntegrationRule :: SetUpPointsOnWedge(int nPointsTri, int nPointsDepth, MaterialMode mode)
{
    //int pointsPassed = 0;

    // TODO: set properly
    firstLocalStrainIndx = 1;
    lastLocalStrainIndx = 3;

    double totArea = 0.0;
    for ( size_t i = 0; i < mTriangles.size(); i++ ) {
        totArea += mTriangles [ i ].getArea();
    }

    std :: vector< int >triToKeep;
    const double triTol = ( 1.0e-6 ) * totArea;

    for ( size_t i = 0; i < mTriangles.size(); i++ ) {
        if ( mTriangles [ i ].getArea() > triTol ) {
            triToKeep.push_back(i);
        }
    }

    int nPointsTot = nPointsTri * nPointsDepth * triToKeep.size();
    FloatArray coords_xi1, coords_xi2, coords_xi3, weightsTri, weightsDepth;
    this->giveTriCoordsAndWeights(nPointsTri, coords_xi1, coords_xi2, weightsTri);
    this->giveLineCoordsAndWeights(nPointsDepth, coords_xi3, weightsDepth);
    this->gaussPoints.resize(nPointsTot);

    std :: vector< FloatArray >newGPCoord;

    double parentArea = this->elem->computeArea();
    int count = 0;

    // Loop over triangles
    for ( int i = 0; i < int( triToKeep.size() ); i++ ) {

        Triangle triangle = mTriangles [ triToKeep [ i ] ];
        
        // global coords of the the triangle verticies
        std::vector< FloatArray > gCoords( triangle.giveNrVertices() );
        for ( int j = 0; j < triangle.giveNrVertices(); j++ ) {
            gCoords[j] = (triangle.giveVertex(j + 1));
        }
        

        for ( int k = 1; k <= nPointsTri; k++ ) {
            for ( int m = 1; m <= nPointsDepth; m++ ) {
                // local coords in the parent triangle
                FloatArray *lCoords = new FloatArray(3);
                lCoords->at(1) = coords_xi1.at(k);
                lCoords->at(2) = coords_xi2.at(k);
                lCoords->at(3) = coords_xi3.at(m);

                double refElArea = 0.5;
                double oldWeight = weightsTri.at(k) * weightsDepth.at(m);
                double newWeight = 2.0 * refElArea * oldWeight * triangle.getArea() / parentArea; 
                
                GaussPoint *gp = new GaussPoint(this, count + 1, lCoords, newWeight, mode);
                this->gaussPoints[count] = gp;
                count++;
                
                
                // Compute global gp coordinate in the element from local gp coord in the sub triangle
                FloatArray global;
                mTriInterp.local2global( global, * gp->giveNaturalCoordinates(),
                                    FEIVertexListGeometryWrapper(gCoords) );
                 
                
                // Compute local gp coordinate in the element from global gp coord in the element
                FloatArray local;
                this->elem->computeLocalCoordinates(local, global);
                local.at(3) = coords_xi3.at(m); // manually set third coordinate
                // compute global coords again, since interpolator dosn't give the z-coord 
                this->elem->computeGlobalCoordinates(global, local);
                
                gp->setGlobalCoordinates(global);
                gp->setNaturalCoordinates(local);
                gp->setSubPatchCoordinates(local);

                // Store new global gp coord for vtk output
                newGPCoord.push_back(global);
            }
        }


        //for ( int k = 0; k < mTriangles [ triToKeep [ i ] ].giveNrVertices(); k++ ) {
        //    delete gCoords [ k ];
        //}

        //delete [] gCoords;
    }

    XfemManager *xMan = elem->giveDomain()->giveXfemManager();
    if ( xMan != NULL ) {
        if ( xMan->giveVtkDebug() ) {
            double time = 0.0;

            Element *el = this->elem;
            if ( el != NULL ) {
//.........这里部分代码省略.........