C++ FEInterpolation::boundaryEvalNormal方法代码示例

void MixedGradientPressureWeakPeriodic :: computeStress(FloatArray &sigmaDev, FloatArray &tractions, double rve_size)
{
    FloatMatrix mMatrix;
    FloatArray normal, coords, t;

    int nsd = domain->giveNumberOfSpatialDimensions();
    Set *set = this->giveDomain()->giveSet(this->set);
    const IntArray &boundaries = set->giveBoundaryList();
    // Reminder: sigma = int t * n dA, where t = sum( C_i * n t_i ).
    // This loop will construct sigma in matrix form.

    FloatMatrix sigma;

    for ( int pos = 1; pos <= boundaries.giveSize() / 2; ++pos ) {
        Element *el = this->giveDomain()->giveElement( boundaries.at(pos * 2 - 1) );
        int boundary = boundaries.at(pos * 2);

        FEInterpolation *interp = el->giveInterpolation(); // Geometry interpolation. The displacements or velocities must have the same interpolation scheme (on the boundary at least).

        int maxorder = this->order + interp->giveInterpolationOrder() * 3;
        std :: unique_ptr< IntegrationRule >ir( interp->giveBoundaryIntegrationRule(maxorder, boundary) );

        for ( GaussPoint *gp: *ir ) {
            const FloatArray &lcoords = gp->giveNaturalCoordinates();
            FEIElementGeometryWrapper cellgeo(el);

            double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
            // Compute v_m = d_dev . x
            interp->boundaryLocal2Global(coords, boundary, lcoords, cellgeo);

            this->constructMMatrix(mMatrix, coords, normal);
            t.beProductOf(mMatrix, tractions);
            sigma.plusDyadUnsym(t, coords, detJ * gp->giveWeight());
        }
    }
    sigma.times(1. / rve_size);

    double pressure = 0.;
    for ( int i = 1; i <= nsd; i++ ) {
        pressure += sigma.at(i, i);
    }
    pressure /= 3; // Not 100% sure about this for 2D cases.
    if ( nsd == 3 ) {
        sigmaDev.resize(6);
        sigmaDev.at(1) = sigma.at(1, 1) - pressure;
        sigmaDev.at(2) = sigma.at(2, 2) - pressure;
        sigmaDev.at(3) = sigma.at(3, 3) - pressure;
        sigmaDev.at(4) = 0.5 * ( sigma.at(2, 3) + sigma.at(3, 2) );
        sigmaDev.at(5) = 0.5 * ( sigma.at(1, 3) + sigma.at(3, 1) );
        sigmaDev.at(6) = 0.5 * ( sigma.at(1, 2) + sigma.at(2, 1) );
    } else if ( nsd == 2 ) {
        sigmaDev.resize(3);
        sigmaDev.at(1) = sigma.at(1, 1) - pressure;
        sigmaDev.at(2) = sigma.at(2, 2) - pressure;
        sigmaDev.at(3) = 0.5 * ( sigma.at(1, 2) + sigma.at(2, 1) );
    } else {
        sigmaDev.resize(1);
        sigmaDev.at(1) = sigma.at(1, 1) - pressure;
    }
}

void SurfaceTensionBoundaryCondition :: computeLoadVectorFromElement(FloatArray &answer, Element *e, int side, TimeStep *tStep)
{
    FEInterpolation *fei = e->giveInterpolation();
    if ( !fei ) {
        OOFEM_ERROR("No interpolation or default integration available for element.");
    }
    std :: unique_ptr< IntegrationRule > iRule( fei->giveBoundaryIntegrationRule(fei->giveInterpolationOrder()-1, side) );

    int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();

    if ( side == -1 ) {
        side = 1;
    }

    answer.clear();

    if ( nsd == 2 ) {

        FEInterpolation2d *fei2d = static_cast< FEInterpolation2d * >(fei);

        ///@todo More of this grunt work should be moved to the interpolation classes
        IntArray bnodes;
        fei2d->boundaryGiveNodes(bnodes, side);
        int nodes = bnodes.giveSize();
        FloatMatrix xy(2, nodes);
        for ( int i = 1; i <= nodes; i++ ) {
            Node *node = e->giveNode(bnodes.at(i));
            ///@todo This should rely on the xindex and yindex in the interpolator..
            xy.at(1, i) = node->giveCoordinate(1);
            xy.at(2, i) = node->giveCoordinate(2);
        }

        FloatArray tmp; // Integrand
        FloatArray es; // Tangent vector to curve
        FloatArray dNds;

        if ( e->giveDomain()->isAxisymmetric() ) {
            FloatArray N;
            FloatArray gcoords;
            for ( GaussPoint *gp: *iRule ) {
                fei2d->edgeEvaldNds( dNds, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                fei->boundaryEvalN( N, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                double J = fei->boundaryGiveTransformationJacobian( side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                fei->boundaryLocal2Global( gcoords, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                double r = gcoords(0); // First coordinate is the radial coord.

                es.beProductOf(xy, dNds);

                tmp.resize( 2 * nodes);
                for ( int i = 0; i < nodes; i++ ) {
                    tmp(2 * i)   = dNds(i) * es(0) * r + N(i);
                    tmp(2 * i + 1) = dNds(i) * es(1) * r;
                }

                answer.add(- 2 * M_PI * gamma * J * gp->giveWeight(), tmp);
            }
        } else {
            for ( GaussPoint *gp: *iRule ) {
                double t = e->giveCrossSection()->give(CS_Thickness, gp);
                fei2d->edgeEvaldNds( dNds, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                double J = fei->boundaryGiveTransformationJacobian( side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
                es.beProductOf(xy, dNds);

                tmp.resize( 2 * nodes);
                for ( int i = 0; i < nodes; i++ ) {
                    tmp(2 * i)   = dNds(i) * es(0);
                    tmp(2 * i + 1) = dNds(i) * es(1);
                    //B.at(1, 1+i*2) = B.at(2, 2+i*2) = dNds(i);
                }
                //tmp.beTProductOf(B, es);

                answer.add(- t * gamma * J * gp->giveWeight(), tmp);
            }
        }
    } else if ( nsd ==  3 ) {

        FEInterpolation3d *fei3d = static_cast< FEInterpolation3d * >(fei);
        FloatArray n, surfProj;
        FloatMatrix dNdx, B;
        for ( GaussPoint *gp: *iRule ) {
            fei3d->surfaceEvaldNdx( dNdx, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );
            double J = fei->boundaryEvalNormal( n, side, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e) );

            // [I - n(x)n]  in voigt form:
            surfProj = {1. - n(0)*n(0), 1. - n(1)*n(1), 1. - n(2)*n(2),
                        - n(1)*n(2), - n(0)*n(2), - n(0)*n(1),
            };

            // Construct B matrix of the surface nodes
            B.resize(6, 3 * dNdx.giveNumberOfRows());
            for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
                B.at(1, 3 * i - 2) = dNdx.at(i, 1);
                B.at(2, 3 * i - 1) = dNdx.at(i, 2);
                B.at(3, 3 * i - 0) = dNdx.at(i, 3);

                B.at(5, 3 * i - 2) = B.at(4, 3 * i - 1) = dNdx.at(i, 3);
                B.at(6, 3 * i - 2) = B.at(4, 3 * i - 0) = dNdx.at(i, 2);
                B.at(6, 3 * i - 1) = B.at(5, 3 * i - 0) = dNdx.at(i, 1);
            }

//.........这里部分代码省略.........

void PrescribedGradientBCNeumann :: integrateTangent(FloatMatrix &oTangent, Element *e, int iBndIndex)
{
    FloatArray normal, n;
    FloatMatrix nMatrix, E_n;
    FloatMatrix contrib;

    Domain *domain = e->giveDomain();
    XfemElementInterface *xfemElInt = dynamic_cast< XfemElementInterface * >( e );

    FEInterpolation *interp = e->giveInterpolation(); // Geometry interpolation

    int nsd = e->giveDomain()->giveNumberOfSpatialDimensions();

    // Interpolation order
    int order = interp->giveInterpolationOrder();
    IntegrationRule *ir = NULL;

    IntArray edgeNodes;
    FEInterpolation2d *interp2d = dynamic_cast< FEInterpolation2d * >( interp );
    if ( interp2d == NULL ) {
        OOFEM_ERROR("failed to cast to FEInterpolation2d.")
    }
    interp2d->computeLocalEdgeMapping(edgeNodes, iBndIndex);

    const FloatArray &xS = * ( e->giveDofManager( edgeNodes.at(1) )->giveCoordinates() );
    const FloatArray &xE = * ( e->giveDofManager( edgeNodes.at( edgeNodes.giveSize() ) )->giveCoordinates() );

    if ( xfemElInt != NULL && domain->hasXfemManager() ) {
        std :: vector< Line >segments;
        std :: vector< FloatArray >intersecPoints;
        xfemElInt->partitionEdgeSegment(iBndIndex, segments, intersecPoints);
        MaterialMode matMode = e->giveMaterialMode();
        ir = new DiscontinuousSegmentIntegrationRule(1, e, segments, xS, xE);
        int numPointsPerSeg = 1;
        ir->SetUpPointsOnLine(numPointsPerSeg, matMode);
    } else   {
        ir = interp->giveBoundaryIntegrationRule(order, iBndIndex);
    }

    oTangent.clear();

    for ( GaussPoint *gp: *ir ) {
        FloatArray &lcoords = * gp->giveNaturalCoordinates();
        FEIElementGeometryWrapper cellgeo(e);

        // Evaluate the normal;
        double detJ = interp->boundaryEvalNormal(normal, iBndIndex, lcoords, cellgeo);

        interp->boundaryEvalN(n, iBndIndex, lcoords, cellgeo);
        // If cracks cross the edge, special treatment is necessary.
        // Exploit the XfemElementInterface to minimize duplication of code.
        if ( xfemElInt != NULL && domain->hasXfemManager() ) {
            // Compute global coordinates of Gauss point
            FloatArray globalCoord;

            interp->boundaryLocal2Global(globalCoord, iBndIndex, lcoords, cellgeo);

            // Compute local coordinates on the element
            FloatArray locCoord;
            e->computeLocalCoordinates(locCoord, globalCoord);

            xfemElInt->XfemElementInterface_createEnrNmatrixAt(nMatrix, locCoord, * e, false);
        } else {
            // Evaluate the velocity/displacement coefficients
            nMatrix.beNMatrixOf(n, nsd);
        }

        if ( nsd == 3 ) {
            OOFEM_ERROR("not implemented for nsd == 3.")
        } else if ( nsd == 2 ) {
            E_n.resize(4, 2);
            E_n.at(1, 1) = normal.at(1);
            E_n.at(1, 2) = 0.0;

            E_n.at(2, 1) = 0.0;
            E_n.at(2, 2) = normal.at(2);

            E_n.at(3, 1) = normal.at(2);
            E_n.at(3, 2) = 0.0;

            E_n.at(4, 1) = 0.0;
            E_n.at(4, 2) = normal.at(1);
        } else {
            E_n.clear();
        }

        contrib.beProductOf(E_n, nMatrix);

        oTangent.add(detJ * gp->giveWeight(), contrib);
    }
    delete ir;
}

void
SolutionbasedShapeFunction :: initializeSurfaceData(modeStruct *mode)
{
    EngngModel *m = mode->myEngngModel;
    double TOL2 = 1e-5;
    IntArray pNodes, mNodes, zNodes;

    Set *mySet = this->domain->giveSet( this->giveSetNumber() );
    IntArray BoundaryList = mySet->giveBoundaryList();

    // First add all nodes to pNodes or nNodes respectively depending on coordinate and normal.
    for ( int i = 0; i < BoundaryList.giveSize() / 2; i++ ) {
        int ElementID = BoundaryList(2 * i);
        int Boundary = BoundaryList(2 * i + 1);

        Element *e = m->giveDomain(1)->giveElement(ElementID);
        FEInterpolation *geoInterpolation = e->giveInterpolation();

        // Check all sides of element
        IntArray bnodes;

#define usePoints 1
#if usePoints == 1
        // Check if all nodes are on the boundary
        geoInterpolation->boundaryGiveNodes(bnodes, Boundary);
        for ( int k = 1; k <= bnodes.giveSize(); k++ ) {
            DofManager *dman = e->giveDofManager( bnodes.at(k) );
            for ( int l = 1; l <= dman->giveCoordinates()->giveSize(); l++ ) {
                if ( fabs( dman->giveCoordinates()->at(l) - maxCoord.at(l) ) < TOL2 ) {
                    pNodes.insertOnce( dman->giveNumber() );
                }
                if ( fabs( dman->giveCoordinates()->at(l) - minCoord.at(l) ) < TOL2 ) {
                    mNodes.insertOnce( dman->giveNumber() );
                }
            }
        }
#else
        // Check normal
        FloatArray lcoords;
        lcoords.resize(2);
        lcoords.at(1) = 0.33333;
        lcoords.at(2) = 0.33333;

        FloatArray normal;
        geoInterpolation->boundaryEvalNormal( normal, j, lcoords, FEIElementGeometryWrapper(e) );
        geoInterpolation->boundaryGiveNodes(bnodes, j);

        printf( "i=%u\tj=%u\t(%f\t%f\t%f)\n", i, j, normal.at(1), normal.at(2), normal.at(3) );
        for ( int k = 1; k <= normal.giveSize(); k++ ) {
            if ( fabs( ( fabs( normal.at(k) ) - 1 ) ) < 1e-4 ) { // Points in x, y or z direction
                addTo = NULL;
                if ( normal.at(k) > 0.5 ) {
                    addTo = & pNodes;
                }
                if ( normal.at(k) < -0.5 ) {
                    addTo = & mNodes;
                }
                if ( addTo != NULL ) {
                    for ( int l = 1; l <= bnodes.giveSize(); l++ ) {
                        bool isSurface = false;
                        DofManager *dman = e->giveDofManager( bnodes.at(l) );
                        dman->giveCoordinates()->printYourself();
                        for ( int m = 1; m <= dman->giveCoordinates()->giveSize(); m++ ) {
                            if ( ( fabs( dman->giveCoordinates()->at(m) - maxCoord.at(m) ) < TOL2 ) || ( fabs( dman->giveCoordinates()->at(m) - minCoord.at(m) ) < TOL2 ) ) {
                                isSurface = true;
                            }
                        }

                        if ( isSurface ) {
                            addTo->insertOnce( e->giveDofManagerNumber( bnodes.at(l) ) );
                        }
                    }
                }
            }
        }
#endif
    }

#if 0
    printf("p=[");
    for ( int i = 1; i < pNodes.giveSize(); i++ ) {
        printf( "%u, ", pNodes.at(i) );
    }
    printf("];\n");
    printf("m=[");
    for ( int i = 1; i < mNodes.giveSize(); i++ ) {
        printf( "%u, ", mNodes.at(i) );
    }
    printf("];\n");
#endif
    //The intersection of pNodes and mNodes constitutes zNodes
    {
        int i = 1, j = 1;
        while ( i <= pNodes.giveSize() ) {
            j = 1;
            while ( j <= mNodes.giveSize() && ( i <= pNodes.giveSize() ) ) {
                //printf("%u == %u?\n", pNodes.at(i), mNodes.at(j));
                if ( pNodes.at(i) == mNodes.at(j) ) {
                    zNodes.insertOnce( pNodes.at(i) );
                    pNodes.erase(i);
//.........这里部分代码省略.........

void
SolutionbasedShapeFunction :: computeCorrectionFactors(modeStruct &myMode, IntArray *Dofs, double *am, double *ap)
{
    /*
     * *Compute c0, cp, cm, Bp, Bm, Ap and Am
     */

    double A0p = 0.0, App = 0.0, A0m = 0.0, Amm = 0.0, Bp = 0.0, Bm = 0.0, c0 = 0.0, cp = 0.0, cm = 0.0;

    EngngModel *m = myMode.myEngngModel;
    Set *mySet = m->giveDomain(1)->giveSet(externalSet);

    IntArray BoundaryList = mySet->giveBoundaryList();

    for ( int i = 0; i < BoundaryList.giveSize() / 2; i++ ) {
        int ElementID = BoundaryList(2 * i);
        int Boundary = BoundaryList(2 * i + 1);

        Element *thisElement = m->giveDomain(1)->giveElement(ElementID);
        FEInterpolation *geoInterpolation = thisElement->giveInterpolation();
        IntArray bnodes, zNodes, pNodes, mNodes;
        FloatMatrix nodeValues;

        geoInterpolation->boundaryGiveNodes(bnodes, Boundary);

        nodeValues.resize( this->dofs.giveSize(), bnodes.giveSize() );
        nodeValues.zero();

        // Change to global ID for bnodes and identify the intersection of bnodes and the zero boundary
        splitBoundaryNodeIDs(myMode, * thisElement, bnodes, pNodes, mNodes, zNodes, nodeValues);

        std :: unique_ptr< IntegrationRule >iRule(geoInterpolation->giveBoundaryIntegrationRule(order, Boundary));

        for ( GaussPoint *gp: *iRule ) {
            FloatArray *lcoords = gp->giveCoordinates();
            FloatArray gcoords, normal, N;
            FloatArray Phi;

            double detJ = fabs( geoInterpolation->boundaryGiveTransformationJacobian( Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) ) ) * gp->giveWeight();

            geoInterpolation->boundaryEvalNormal( normal, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );
            geoInterpolation->boundaryEvalN( N, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );
            geoInterpolation->boundaryLocal2Global( gcoords, Boundary, * lcoords, FEIElementGeometryWrapper(thisElement) );

            FloatArray pPhi, mPhi, zPhi;
            pPhi.resize( Dofs->giveSize() );
            pPhi.zero();
            mPhi.resize( Dofs->giveSize() );
            mPhi.zero();
            zPhi.resize( Dofs->giveSize() );
            zPhi.zero();

            // Build phi (analytical averaging, not projected onto the mesh)
            computeBaseFunctionValueAt(Phi, gcoords, * Dofs, * myMode.myEngngModel);

            // Build zPhi for this DofID
            for ( int l = 1; l <= zNodes.giveSize(); l++ ) {
                int nodeID = zNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    zPhi.at(m) = zPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }


            // Build pPhi for this DofID
            for ( int l = 1; l <= pNodes.giveSize(); l++ ) {
                int nodeID = pNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    pPhi.at(m) = pPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }

            // Build mPhi for this DofID
            for ( int l = 1; l <= mNodes.giveSize(); l++ ) {
                int nodeID = mNodes.at(l);
                for ( int m = 1; m <= this->dofs.giveSize(); m++ ) {
                    mPhi.at(m) = mPhi.at(m) + N.at(nodeID) * nodeValues.at(m, nodeID);
                }
            }

            c0 = c0 + zPhi.dotProduct(normal, 3) * detJ;
            cp = cp + pPhi.dotProduct(normal, 3) * detJ;
            cm = cm + mPhi.dotProduct(normal, 3) * detJ;

            App = App + pPhi.dotProduct(pPhi, 3) * detJ;
            Amm = Amm + mPhi.dotProduct(mPhi, 3) * detJ;
            A0p = A0p + zPhi.dotProduct(pPhi, 3) * detJ;
            A0m = A0m + zPhi.dotProduct(mPhi, 3) * detJ;

            Bp = Bp + Phi.dotProduct(pPhi, 3) * detJ;
            Bm = Bm + Phi.dotProduct(mPhi, 3) * detJ;
        }
    }

    * am = -( A0m * cp * cp - Bm * cp * cp - A0p * cm * cp + App * c0 * cm + Bp * cm * cp ) / ( App * cm * cm + Amm * cp * cp );
    * ap = -( A0p * cm * cm - Bp * cm * cm - A0m * cm * cp + Amm * c0 * cp + Bm * cm * cp ) / ( App * cm * cm + Amm * cp * cp );
}