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

void
DKTPlate :: computeVertexBendingMoments(FloatMatrix &answer, TimeStep *tStep)
{
#ifdef DKT_EnableVertexMomentsCache
    if ( stateCounter == tStep->giveSolutionStateCounter() ) {
        answer = vertexMoments;
        return;
    }
#endif

    // the results should be cached somehow, as computing on the fly is highly inefficient
    // due to multiple requests
    FloatMatrix dndx;
    answer.resize(5, 3);

    FloatMatrix b;
    FloatArray eps, m;
    FloatArray coords [ 3 ]; // vertex local coordinates
    coords [ 0 ] = {
        1.0, 0.0
    };
    coords [ 1 ] = {
        0.0, 1.0
    };
    coords [ 2 ] = {
        0.0, 0.0
    };

    GaussIntegrationRule iRule = GaussIntegrationRule(1, this, 1, 1); // dummy rule used to evaluate B at vertices
    iRule.SetUpPointsOnTriangle(1, _Unknown);
    GaussPoint *vgp = iRule.getIntegrationPoint(0);

    for ( int i = 1; i <= this->numberOfDofMans; i++ ) {
        vgp->setNaturalCoordinates(coords [ i - 1 ]);
        this->computeStrainVector(eps, vgp, tStep);
        this->giveStructuralCrossSection()->giveGeneralizedStress_Plate(m, vgp, eps, tStep);
        answer.setColumn(m, i);
    }

#ifdef DKT_EnableVertexMomentsCache
    this->vertexMoments = answer;
    this->stateCounter = tStep->giveSolutionStateCounter();
#endif
}

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 ) {
//.........这里部分代码省略.........