C++ polyMesh::points方法代码示例

Foam::labelList Foam::meshToMesh::maskCells
(
    const polyMesh& src,
    const polyMesh& tgt
) const
{
    boundBox intersectBb
    (
        max(src.bounds().min(), tgt.bounds().min()),
        min(src.bounds().max(), tgt.bounds().max())
    );

    intersectBb.inflate(0.01);

    const cellList& srcCells = src.cells();
    const faceList& srcFaces = src.faces();
    const pointField& srcPts = src.points();

    DynamicList<label> cells(src.size());
    forAll(srcCells, srcI)
    {
        boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
        if (intersectBb.overlaps(cellBb))
        {
            cells.append(srcI);
        }
    }

// Write set to VTK readable files
void writeVTK
(
    const polyMesh& mesh,
    const topoSet& currentSet,
    const fileName& vtkName
)
{
    if (isA<faceSet>(currentSet))
    {
        // Faces of set with OpenFOAM faceID as value

        faceList setFaces(currentSet.size());
        labelList faceValues(currentSet.size());
        label setFaceI = 0;

        forAllConstIter(topoSet, currentSet, iter)
        {
            setFaces[setFaceI] = mesh.faces()[iter.key()];
            faceValues[setFaceI] = iter.key();
            setFaceI++;
        }

        primitiveFacePatch fp(setFaces, mesh.points());

        writePatch
        (
            true,
            currentSet.name(),
            fp,
            "faceID",
            faceValues,
            mesh.time().path()/vtkName
        );
    }

void Foam::snappyLayerDriver::averageNeighbours
(
    const polyMesh& mesh,
    const PackedBoolList& isMasterEdge,
    const labelList& meshEdges,
    const labelList& meshPoints,
    const edgeList& edges,
    const scalarField& invSumWeight,
    const Field<Type>& data,
    Field<Type>& average
)
{
    const pointField& pts = mesh.points();

    average = Zero;

    forAll(edges, edgeI)
    {
        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
        {
            const edge& e = edges[edgeI];
            //scalar eWeight = edgeWeights[edgeI];
            //scalar eWeight =  1.0;
            scalar eMag = max
            (
                VSMALL,
                mag
                (
                    pts[meshPoints[e[1]]]
                  - pts[meshPoints[e[0]]]
                )
            );
            scalar eWeight = 1.0/eMag;

            label v0 = e[0];
            label v1 = e[1];

            average[v0] += eWeight*data[v1];
            average[v1] += eWeight*data[v0];
        }
    }

    syncTools::syncPointList
    (
        mesh,
        meshPoints,
        average,
        plusEqOp<Type>(),
        Zero     // null value
    );

    average *= invSumWeight;
}

Foam::Cloud<ParticleType>::Cloud
(
    const polyMesh& pMesh,
    const IDLList<ParticleType>& particles
)
:
    cloud(pMesh),
    IDLList<ParticleType>(particles),
    polyMesh_(pMesh),
    allFaces_(pMesh.faces()),
    points_(pMesh.points()),
    cellFaces_(pMesh.cells()),
    allFaceCentres_(pMesh.faceCentres()),
    owner_(pMesh.faceOwner()),
    neighbour_(pMesh.faceNeighbour()),
    meshInfo_(polyMesh_)
{}

Foam::Cloud<ParticleType>::Cloud
(
    const polyMesh& pMesh,
    const bool checkClass
)
:
    cloud(pMesh),
    polyMesh_(pMesh),
    allFaces_(pMesh.faces()),
    points_(pMesh.points()),
    cellFaces_(pMesh.cells()),
    allFaceCentres_(pMesh.faceCentres()),
    owner_(pMesh.faceOwner()),
    neighbour_(pMesh.faceNeighbour()),
    meshInfo_(polyMesh_)
{
    initCloud(checkClass);
}

void Foam::meshToMeshNew::writeConnectivity
(
    const polyMesh& src,
    const polyMesh& tgt,
    const labelListList& srcToTargetAddr
) const
{
    Pout<< "Source size = " << src.nCells() << endl;
    Pout<< "Target size = " << tgt.nCells() << endl;

    word fName("addressing_" + src.name() + "_to_" + tgt.name());

    if (Pstream::parRun())
    {
        fName = fName +  "_proc" + Foam::name(Pstream::myProcNo());
    }

    OFstream os(src.time().path()/fName + ".obj");

    label vertI = 0;
    forAll(srcToTargetAddr, i)
    {
        const labelList& tgtAddress = srcToTargetAddr[i];
        forAll(tgtAddress, j)
        {
            label tgtI = tgtAddress[j];
            const vector& c0 = src.cellCentres()[i];

            const cell& c = tgt.cells()[tgtI];
            const pointField pts(c.points(tgt.faces(), tgt.points()));
            forAll(pts, j)
            {
                const point& p = pts[j];
                os  << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
                vertI++;
                os  << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
                    << nl;
                vertI++;
                os  << "l " << vertI - 1 << ' ' << vertI << nl;
            }
        }
    }

RBFMeshMotionSolver::RBFMeshMotionSolver(
    const polyMesh & mesh,
    Istream & msData
    )
    :
    motionSolver( mesh ),
    motionCenters( mesh.boundaryMesh().size(), vectorField( 0 ) ),
    staticPatches( lookup( "staticPatches" ) ),
    staticPatchIDs( staticPatches.size() ),
    movingPatches( lookup( "movingPatches" ) ),
    movingPatchIDs( movingPatches.size() ),
    fixedPatches( lookup( "fixedPatches" ) ),
    fixedPatchIDs( fixedPatches.size() ),
    newPoints( mesh.points().size(), vector::zero ),
    rbf( false ),
    nbGlobalFaceCenters( Pstream::nProcs(), 0 ),
    nbGlobalMovingFaceCenters( Pstream::nProcs(), 0 ),
    nbGlobalStaticFaceCenters( Pstream::nProcs(), 0 ),
    nbGlobalFixedFaceCenters( Pstream::nProcs(), 0 ),
    globalMovingPointsLabelList( mesh.boundaryMesh().size(), labelList( 0 ) ),
    twoDCorrector( mesh ),
    nbPoints( 0 ),
    faceCellCenters( true ),
    cpu( false ),
    timeIntegrationScheme( false ),
    corrector( false ),
    k( 0 )
{
    // Find IDs of staticPatches
    forAll( staticPatches, patchI )
    {
        label patchIndex = mesh.boundaryMesh().findPatchID( staticPatches[patchI] );

        assert( patchIndex >= 0 );

        staticPatchIDs[patchI] = patchIndex;
    }

// Naive feature detection. All boundary edges with angle > featureAngle become
// feature edges. All points on feature edges become feature points. All
// boundary faces become feature faces.
void simpleMarkFeatures
(
    const polyMesh& mesh,
    const PackedBoolList& isBoundaryEdge,
    const scalar featureAngle,
    const bool concaveMultiCells,
    const bool doNotPreserveFaceZones,

    labelList& featureFaces,
    labelList& featureEdges,
    labelList& singleCellFeaturePoints,
    labelList& multiCellFeaturePoints
)
{
    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);

    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    // Working sets
    labelHashSet featureEdgeSet;
    labelHashSet singleCellFeaturePointSet;
    labelHashSet multiCellFeaturePointSet;


    // 1. Mark all edges between patches
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
        const labelList& meshEdges = pp.meshEdges();

        // All patch corner edges. These need to be feature points & edges!
        for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
        {
            label meshEdgeI = meshEdges[edgeI];
            featureEdgeSet.insert(meshEdgeI);
            singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][0]);
            singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][1]);
        }
    }



    // 2. Mark all geometric feature edges
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Make distinction between convex features where the boundary point becomes
    // a single cell and concave features where the boundary point becomes
    // multiple 'half' cells.

    // Addressing for all outside faces
    primitivePatch allBoundary
    (
        SubList<face>
        (
            mesh.faces(),
            mesh.nFaces()-mesh.nInternalFaces(),
            mesh.nInternalFaces()
        ),
        mesh.points()
    );

    // Check for non-manifold points (surface pinched at point)
    allBoundary.checkPointManifold(false, &singleCellFeaturePointSet);

    // Check for non-manifold edges (surface pinched at edge)
    const labelListList& edgeFaces = allBoundary.edgeFaces();
    const labelList& meshPoints = allBoundary.meshPoints();

    forAll(edgeFaces, edgeI)
    {
        const labelList& eFaces = edgeFaces[edgeI];

        if (eFaces.size() > 2)
        {
            const edge& e = allBoundary.edges()[edgeI];

            //Info<< "Detected non-manifold boundary edge:" << edgeI
            //    << " coords:"
            //    << allBoundary.points()[meshPoints[e[0]]]
            //    << allBoundary.points()[meshPoints[e[1]]] << endl;

            singleCellFeaturePointSet.insert(meshPoints[e[0]]);
            singleCellFeaturePointSet.insert(meshPoints[e[1]]);
        }
    }

    // Check for features.
    forAll(edgeFaces, edgeI)
    {
        const labelList& eFaces = edgeFaces[edgeI];

        if (eFaces.size() == 2)
        {
            label f0 = eFaces[0];
            label f1 = eFaces[1];

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

// Check the blockMesh topology
void Foam::blockMesh::checkBlockMesh(const polyMesh& bm) const
{
    if (verboseOutput)
    {
        Info<< nl << "Check topology" << endl;
    }

    bool ok = true;

    const pointField& points = bm.points();
    const faceList& faces = bm.faces();
    const cellList& cells = bm.cells();
    const polyPatchList& patches = bm.boundaryMesh();

    label nBoundaryFaces = 0;
    forAll(cells, celli)
    {
        nBoundaryFaces += cells[celli].nFaces();
    }

    nBoundaryFaces -= 2*bm.nInternalFaces();

    label nDefinedBoundaryFaces = 0;
    forAll(patches, patchi)
    {
        nDefinedBoundaryFaces += patches[patchi].size();
    }


    if (verboseOutput)
    {
        Info<< nl << tab << "Basic statistics" << nl
            << tab << tab << "Number of internal faces : "
            << bm.nInternalFaces() << nl
            << tab << tab << "Number of boundary faces : "
            << nBoundaryFaces << nl
            << tab << tab << "Number of defined boundary faces : "
            << nDefinedBoundaryFaces << nl
            << tab << tab << "Number of undefined boundary faces : "
            << nBoundaryFaces - nDefinedBoundaryFaces << nl;

        if ((nBoundaryFaces - nDefinedBoundaryFaces) > 0)
        {
            Info<< tab << tab << tab
                << "(Warning : only leave undefined the front and back planes "
                << "of 2D planar geometries!)" << endl;
        }

        Info<< tab << "Checking patch -> block consistency" << endl;
    }


    forAll(patches, patchi)
    {
        const faceList& Patch = patches[patchi];

        forAll(Patch, patchFacei)
        {
            const face& patchFace = Patch[patchFacei];
            bool patchFaceOK = false;

            forAll(cells, celli)
            {
                const labelList& cellFaces = cells[celli];

                forAll(cellFaces, cellFacei)
                {
                    if (patchFace == faces[cellFaces[cellFacei]])
                    {
                        patchFaceOK = true;

                        if
                        (
                            (
                                patchFace.normal(points)
                              & faces[cellFaces[cellFacei]].normal(points)
                            ) < 0.0
                        )
                        {
                            Info<< tab << tab
                                << "Face " << patchFacei
                                << " of patch " << patchi
                                << " (" << patches[patchi].name() << ")"
                                << " points inwards"
                                << endl;

                            ok = false;
                        }
                    }
                }
            }

            if (!patchFaceOK)
            {
                Info<< tab << tab
                    << "Face " << patchFacei
                    << " of patch " << patchi
                    << " (" << patches[patchi].name() << ")"
                    << " does not match any block faces" << endl;
//.........这里部分代码省略.........

Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
(
    const polyMesh& mesh,
    label fI,
    const point& nCc,
    scalar tol,
    bool report
)
{
    const faceList& pFaces = mesh.faces();
    const pointField& pPts = mesh.points();
    const vectorField& pC = mesh.cellCentres();
    const labelList& pOwner = mesh.faceOwner();

    const face& f = pFaces[fI];

    label oCI = pOwner[fI];

    const point& oCc = pC[oCI];

    List<scalar> tetQualities(2, 0.0);

    forAll(f, faceBasePtI)
    {
        scalar thisBaseMinTetQuality = VGREAT;

        const point& tetBasePt = pPts[f[faceBasePtI]];

        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
        {
            label facePtI = (tetPtI + faceBasePtI) % f.size();
            label otherFacePtI = f.fcIndex(facePtI);

            {
                // owner cell tet
                label ptAI = f[facePtI];
                label ptBI = f[otherFacePtI];

                const point& pA = pPts[ptAI];
                const point& pB = pPts[ptBI];

                tetPointRef tet(oCc, tetBasePt, pA, pB);

                tetQualities[0] = tet.quality();
            }

            {
                // neighbour cell tet
                label ptAI = f[otherFacePtI];
                label ptBI = f[facePtI];

                const point& pA = pPts[ptAI];
                const point& pB = pPts[ptBI];

                tetPointRef tet(nCc, tetBasePt, pA, pB);

                tetQualities[1] = tet.quality();
            }

            if (min(tetQualities) < thisBaseMinTetQuality)
            {
                thisBaseMinTetQuality = min(tetQualities);
            }
        }

        if (thisBaseMinTetQuality > tol)
        {
            return faceBasePtI;
        }
    }

bool Foam::motionSmootherAlgo::checkMesh
(
    const bool report,
    const polyMesh& mesh,
    const dictionary& dict,
    const labelList& checkFaces,
    const List<labelPair>& baffles,
    labelHashSet& wrongFaces
)
{
    const scalar maxNonOrtho
    (
        readScalar(dict.lookup("maxNonOrtho", true))
    );
    const scalar minVol
    (
        readScalar(dict.lookup("minVol", true))
    );
    const scalar minTetQuality
    (
        readScalar(dict.lookup("minTetQuality", true))
    );
    const scalar maxConcave
    (
        readScalar(dict.lookup("maxConcave", true))
    );
    const scalar minArea
    (
        readScalar(dict.lookup("minArea", true))
    );
    const scalar maxIntSkew
    (
        readScalar(dict.lookup("maxInternalSkewness", true))
    );
    const scalar maxBounSkew
    (
        readScalar(dict.lookup("maxBoundarySkewness", true))
    );
    const scalar minWeight
    (
        readScalar(dict.lookup("minFaceWeight", true))
    );
    const scalar minVolRatio
    (
        readScalar(dict.lookup("minVolRatio", true))
    );
    const scalar minTwist
    (
        readScalar(dict.lookup("minTwist", true))
    );
    const scalar minTriangleTwist
    (
        readScalar(dict.lookup("minTriangleTwist", true))
    );
    scalar minFaceFlatness = -1.0;
    dict.readIfPresent("minFaceFlatness", minFaceFlatness, true);
    const scalar minDet
    (
        readScalar(dict.lookup("minDeterminant", true))
    );
    label nWrongFaces = 0;

    Info<< "Checking faces in error :" << endl;
    //Pout.setf(ios_base::left);

    if (maxNonOrtho < 180.0-SMALL)
    {
        polyMeshGeometry::checkFaceDotProduct
        (
            report,
            maxNonOrtho,
            mesh,
            mesh.cellCentres(),
            mesh.faceAreas(),
            checkFaces,
            baffles,
            &wrongFaces
        );

        label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());

        Info<< "    non-orthogonality > "
            << setw(3) << maxNonOrtho
            << " degrees                        : "
            << nNewWrongFaces-nWrongFaces << endl;

        nWrongFaces = nNewWrongFaces;
    }

    if (minVol > -GREAT)
    {
        polyMeshGeometry::checkFacePyramids
        (
            report,
            minVol,
            mesh,
            mesh.cellCentres(),
            mesh.points(),
            checkFaces,
            baffles,
//.........这里部分代码省略.........

//- Returns -1 or cartesian coordinate component (0=x, 1=y, 2=z) of normal
//  in case of 2D mesh
label twoDNess(const polyMesh& mesh)
{
    const pointField& ctrs = mesh.cellCentres();

    if (ctrs.size() < 2)
    {
        return -1;
    }


    //
    // 1. All cell centres on single plane aligned with x, y or z
    //

    // Determine 3 points to base plane on.
    vector vec10 = ctrs[1] - ctrs[0];
    vec10 /= mag(vec10);

    label otherCellI = -1;

    for (label cellI = 2; cellI < ctrs.size(); cellI++)
    {
        vector vec(ctrs[cellI] - ctrs[0]);
        vec /= mag(vec);

        if (mag(vec & vec10) < 0.9)
        {
            // ctrs[cellI] not in line with n
            otherCellI = cellI;

            break;
        }
    }

    if (otherCellI == -1)
    {
        // Cannot find cell to make decent angle with cell0-cell1 vector.
        // Note: what to do here? All cells (almost) in one line. Maybe 1D case?
        return -1;
    }

    plane cellPlane(ctrs[0], ctrs[1], ctrs[otherCellI]);


    forAll (ctrs, cellI)
    {
        const labelList& cEdges = mesh.cellEdges()[cellI];

        scalar minLen = GREAT;

        forAll (cEdges, i)
        {
            minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
        }

        if (cellPlane.distance(ctrs[cellI]) > 1E-6*minLen)
        {
            // Centres not in plane
            return  -1;
        }
    }

Foam::autoPtr<Foam::mapDistribute> Foam::meshToMeshNew::calcProcMap
(
    const polyMesh& src,
    const polyMesh& tgt
) const
{
    // get decomposition of cells on src mesh
    List<boundBox> procBb(Pstream::nProcs());

    if (src.nCells() > 0)
    {
        // bounding box for my mesh - do not parallel reduce
        procBb[Pstream::myProcNo()] = boundBox(src.points(), false);

        // slightly increase size of bounding boxes to allow for cases where
        // bounding boxes are perfectly alligned
        procBb[Pstream::myProcNo()].inflate(0.01);
    }
    else
    {
        procBb[Pstream::myProcNo()] = boundBox();
    }


    Pstream::gatherList(procBb);
    Pstream::scatterList(procBb);


    if (debug)
    {
        Info<< "Determining extent of src mesh per processor:" << nl
            << "\tproc\tbb" << endl;
        forAll(procBb, procI)
        {
            Info<< '\t' << procI << '\t' << procBb[procI] << endl;
        }
    }


    // determine which cells of tgt mesh overlaps src mesh per proc
    const cellList& cells = tgt.cells();
    const faceList& faces = tgt.faces();
    const pointField& points = tgt.points();

    labelListList sendMap;

    {
        // per processor indices into all segments to send
        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
        label iniSize = floor(tgt.nCells()/Pstream::nProcs());

        forAll(dynSendMap, procI)
        {
            dynSendMap[procI].setCapacity(iniSize);
        }

        // work array - whether src processor bb overlaps the tgt cell bounds
        boolList procBbOverlaps(Pstream::nProcs());
        forAll(cells, cellI)
        {
            const cell& c = cells[cellI];

            // determine bounding box of tgt cell
            boundBox cellBb(point::max, point::min);
            forAll(c, faceI)
            {
                const face& f = faces[c[faceI]];
                forAll(f, fp)
                {
                    cellBb.min() = min(cellBb.min(), points[f[fp]]);
                    cellBb.max() = max(cellBb.max(), points[f[fp]]);
                }
            }

            // find the overlapping tgt cells on each src processor
            (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);

            forAll(procBbOverlaps, procI)
            {
                if (procBbOverlaps[procI])
                {
                    dynSendMap[procI].append(cellI);
                }
            }
        }

bool Foam::polyMeshZipUpCells(polyMesh& mesh)
{
    if (polyMesh::debug)
    {
        Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
            << "zipping up topologically open cells" << endl;
    }

    // Algorithm:
    // Take the original mesh and visit all cells.  For every cell
    // calculate the edges of all faces on the cells.  A cell is
    // correctly topologically closed when all the edges are referenced
    // by exactly two faces.  If the edges are referenced only by a
    // single face, additional vertices need to be inserted into some
    // of the faces (topological closedness).  If an edge is
    // referenced by more that two faces, there is an error in
    // topological closedness.
    // Point insertion into the faces is done by attempting to create
    // closed loops and inserting the intermediate points into the
    // defining edge
    // Note:
    // The algorithm is recursive and changes the mesh faces in each
    // pass.  It is therefore essential to discard the addressing
    // after every pass.  The algorithm is completed when the mesh
    // stops changing.

    label nChangedFacesInMesh = 0;
    label nCycles = 0;

    labelHashSet problemCells;

    do
    {
        nChangedFacesInMesh = 0;

        const cellList& Cells = mesh.cells();
        const pointField& Points = mesh.points();

        faceList newFaces = mesh.faces();

        const faceList& oldFaces = mesh.faces();
        const labelListList& pFaces = mesh.pointFaces();

        forAll(Cells, cellI)
        {
            const labelList& curFaces = Cells[cellI];
            const edgeList cellEdges = Cells[cellI].edges(oldFaces);
            const labelList cellPoints = Cells[cellI].labels(oldFaces);

            // Find the edges used only once in the cell

            labelList edgeUsage(cellEdges.size(), 0);

            forAll(curFaces, faceI)
            {
                edgeList curFaceEdges = oldFaces[curFaces[faceI]].edges();

                forAll(curFaceEdges, faceEdgeI)
                {
                    const edge& curEdge = curFaceEdges[faceEdgeI];

                    forAll(cellEdges, cellEdgeI)
                    {
                        if (cellEdges[cellEdgeI] == curEdge)
                        {
                            edgeUsage[cellEdgeI]++;
                            break;
                        }
                    }
                }
            }

            edgeList singleEdges(cellEdges.size());
            label nSingleEdges = 0;

            forAll(edgeUsage, edgeI)
            {
                if (edgeUsage[edgeI] == 1)
                {
                    singleEdges[nSingleEdges] = cellEdges[edgeI];
                    nSingleEdges++;
                }
                else if (edgeUsage[edgeI] != 2)
                {
                    WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
                        << "edge " << cellEdges[edgeI] << " in cell " << cellI
                        << " used " << edgeUsage[edgeI] << " times. " << nl
                        << "Should be 1 or 2 - serious error "
                        << "in mesh structure. " << endl;

#                   ifdef DEBUG_ZIPUP
                    forAll(curFaces, faceI)
                    {
                        Info<< "face: " << oldFaces[curFaces[faceI]]
                            << endl;
                    }

                    Info<< "Cell edges: " << cellEdges << nl
                        << "Edge usage: " << edgeUsage << nl
                        << "Cell points: " << cellPoints << endl;
//.........这里部分代码省略.........

void Foam::cellPointWeight::findTriangle
(
    const polyMesh& mesh,
    const vector& position,
    const label faceIndex
)
{
    if (debug)
    {
        Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
            << "position = " << position << nl
            << "faceIndex = " << faceIndex << endl;
    }

    // Initialise closest triangle variables
    scalar minUVClose = VGREAT;
    label pointIClose = 0;

    // Decompose each face into triangles, making a tet when
    // augmented by the cell centre
    const labelList& facePoints = mesh.faces()[faceIndex];

    const scalar faceArea2 = magSqr(mesh.faceAreas()[faceIndex]);

    label pointI = 1;
    while ((pointI + 1) < facePoints.size())
    {
        // Cartesian co-ordinates of the triangle vertices
        const vector& P1 = mesh.points()[facePoints[0]];
        const vector& P2 = mesh.points()[facePoints[pointI]];
        const vector& P3 = mesh.points()[facePoints[pointI + 1]];

        // Direction vectors
        vector v1 = position - P1;
        const vector v2 = P2 - P1;
        const vector v3 = P3 - P1;

        // Plane normal
        vector n = v2 ^ v3;
        n /= mag(n);

        // Remove any offset to plane
        v1 -= (n & v1)*v1;

        // Helper variables
        const scalar d12 = v1 & v2;
        const scalar d13 = v1 & v3;
        const scalar d22 = v2 & v2;
        const scalar d23 = v2 & v3;
        const scalar d33 = v3 & v3;

        // Determinant of coefficients matrix
        // Note: if det(A) = 0 the triangle is degenerate
        const scalar detA = d22*d33 - d23*d23;

        if (0.25*detA/faceArea2 > tol)
        {
            // Solve using Cramers' rule
            const scalar u = (d12*d33 - d23*d13)/detA;
            const scalar v = (d22*d13 - d12*d23)/detA;

            // Check if point is in triangle
            if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
            {
                // Indices of the cell vertices making up the triangle
                faceVertices_[0] = facePoints[0];
                faceVertices_[1] = facePoints[pointI];
                faceVertices_[2] = facePoints[pointI + 1];

                weights_[0] = u;
                weights_[1] = v;
                weights_[2] = 1.0 - (u + v);
                weights_[3] = 0.0;

                return;
            }
            else
            {
                scalar minU = mag(u);
                scalar minV = mag(v);
                if (minU > 1.0)
                {
                    minU -= 1.0;
                }
                if (minV > 1.0)
                {
                    minV -= 1.0;
                }
                const scalar minUV = mag(minU + minV);

                if (minUV < minUVClose)
                {
                    minUVClose = minUV;
                    pointIClose = pointI;
                }
            }
        }

        pointI++;
    }
//.........这里部分代码省略.........