static triSurface pack
(
    const triSurface& surf,
    const pointField& localPoints,
    const labelList& pointMap
)
{
    List<labelledTri> newTriangles(surf.size());
    label newTriangleI = 0;

    forAll(surf, faceI)
    {
        const labelledTri& f = surf.localFaces()[faceI];

        label newA = pointMap[f[0]];
        label newB = pointMap[f[1]];
        label newC = pointMap[f[2]];

        if ((newA != newB) && (newA != newC) && (newB != newC))
        {
            newTriangles[newTriangleI++] =
                labelledTri(newA, newB, newC, f.region());
        }
    }
    newTriangles.setSize(newTriangleI);

    return triSurface(newTriangles, surf.patches(), localPoints);
}

Foam::vector Foam::surfaceLocation::normal(const triSurface& s) const
{
    const vectorField& n = s.faceNormals();

    if (elementType_ == triPointRef::NONE)
    {
        return n[index()];
    }
    else if (elementType_ == triPointRef::EDGE)
    {
        const labelList& eFaces = s.edgeFaces()[index()];

        if (eFaces.size() == 1)
        {
            return n[eFaces[0]];
        }
        else
        {
            vector edgeNormal(vector::zero);

            forAll(eFaces, i)
            {
                edgeNormal += n[eFaces[i]];
            }
            return edgeNormal/(mag(edgeNormal) + VSMALL);
        }
    }
    else
    {
        return s.pointNormals()[index()];

// Update intersections for selected edges.
void Foam::edgeIntersections::intersectEdges
(
    const triSurface& surf1,
    const pointField& points1,          // surf1 meshPoints (not localPoints!)
    const triSurfaceSearch& querySurf2,
    const scalarField& surf1PointTol,   // surf1 tolerance per point
    const labelList& edgeLabels
)
{
    const triSurface& surf2 = querySurf2.surface();
    const vectorField& normals2 = surf2.faceNormals();

    const labelList& meshPoints = surf1.meshPoints();

    if (debug)
    {
        Pout<< "Calculating intersection of " << edgeLabels.size() << " edges"
            << " out of " << surf1.nEdges() << " with " << surf2.size()
            << " triangles ..." << endl;
    }

    pointField start(edgeLabels.size());
    pointField end(edgeLabels.size());
    vectorField edgeDirs(edgeLabels.size());

    // Go through all edges, calculate intersections
    forAll(edgeLabels, i)
    {
        label edgeI = edgeLabels[i];

        if (debug)// && (i % 1000 == 0))
        {
            Pout<< "Intersecting edge " << edgeI << " with surface" << endl;
        }

        const edge& e = surf1.edges()[edgeI];

        const point& pStart = points1[meshPoints[e.start()]];
        const point& pEnd = points1[meshPoints[e.end()]];

        const vector eVec(pEnd - pStart);
        const vector n(eVec/(mag(eVec) + VSMALL));

        // Start tracking somewhat before pStart and up to somewhat after p1.
        // Note that tolerances here are smaller than those used to classify
        // hit below.
        // This will cause this hit to be marked as degenerate and resolved
        // later on.
        start[i] = pStart - 0.5*surf1PointTol[e[0]]*n;
        end[i] = pEnd + 0.5*surf1PointTol[e[1]]*n;

        edgeDirs[i] = n;
    }

// Dump collapse region to .obj file
static void writeRegionOBJ
(
    const triSurface& surf,
    const label regionI,
    const labelList& collapseRegion,
    const labelList& outsideVerts
)
{
    fileName dir("regions");

    mkDir(dir);
    fileName regionName(dir / "region_" + name(regionI) + ".obj");

    Pout<< "Dumping region " << regionI << " to file " << regionName << endl;

    boolList include(surf.size(), false);

    forAll(collapseRegion, faceI)
    {
        if (collapseRegion[faceI] == regionI)
        {
            include[faceI] = true;
        }
    }

    labelList pointMap, faceMap;

    triSurface regionSurf(surf.subsetMesh(include, pointMap, faceMap));

    Pout<< "Region " << regionI << " surface:" << nl;
    regionSurf.writeStats(Pout);

    regionSurf.write(regionName);


    // Dump corresponding outside vertices.
    fileName pointsName(dir / "regionPoints_" + name(regionI) + ".obj");

    Pout<< "Dumping region " << regionI << " points to file " << pointsName
        << endl;

    OFstream str(pointsName);

    forAll(outsideVerts, i)
    {
        meshTools::writeOBJ(str, surf.localPoints()[outsideVerts[i]]);
    }

void Foam::edgeIntersections::checkEdges(const triSurface& surf)
{
    const pointField& localPoints = surf.localPoints();
    const edgeList& edges = surf.edges();
    const labelListList& edgeFaces = surf.edgeFaces();

    treeBoundBox bb(localPoints);

    scalar minSize = SMALL * bb.minDim();

    forAll(edges, edgeI)
    {
        const edge& e = edges[edgeI];

        scalar eMag = e.mag(localPoints);

        if (eMag < minSize)
        {
            WarningIn
            (
                "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
            )   << "Edge " << edgeI << " vertices " << e
                << " coords:" << localPoints[e[0]] << ' '
                << localPoints[e[1]] << " is very small compared to bounding"
                << " box dimensions " << bb << endl
                << "This might lead to problems in intersection"
                << endl;
        }

        if (edgeFaces[edgeI].size() == 1)
        {
            WarningIn
            (
                "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
            )   << "Edge " << edgeI << " vertices " << e
                << " coords:" << localPoints[e[0]] << ' '
                << localPoints[e[1]] << " has only one face connected to it:"
                << edgeFaces[edgeI] << endl
                << "This might lead to problems in intersection"
                << endl;
        }
    }
}

//- Sets point neighbours of face to val
static void markPointNbrs
(
    const triSurface& surf,
    const label faceI,
    const bool val,
    boolList& okToCollapse
)
{
    const triSurface::FaceType& f = surf.localFaces()[faceI];

    forAll(f, fp)
    {
        const labelList& pFaces = surf.pointFaces()[f[fp]];

        forAll(pFaces, i)
        {
            okToCollapse[pFaces[i]] = false;
        }
    }
}

// Print on master all the per-processor surface stats.
void writeProcStats
(
    const triSurface& s,
    const List<List<treeBoundBox> >& meshBb
)
{
    // Determine surface bounding boxes, faces, points
    List<treeBoundBox> surfBb(Pstream::nProcs());
    {
        surfBb[Pstream::myProcNo()] = treeBoundBox(s.points());
        Pstream::gatherList(surfBb);
        Pstream::scatterList(surfBb);
    }

    labelList nPoints(Pstream::nProcs());
    nPoints[Pstream::myProcNo()] = s.points().size();
    Pstream::gatherList(nPoints);
    Pstream::scatterList(nPoints);

    labelList nFaces(Pstream::nProcs());
    nFaces[Pstream::myProcNo()] = s.size();
    Pstream::gatherList(nFaces);
    Pstream::scatterList(nFaces);

    forAll(surfBb, procI)
    {
        const List<treeBoundBox>& bbs = meshBb[procI];

        Info<< "processor" << procI << nl
            << "\tMesh bounds          : " << bbs[0] << nl;
        for (label i = 1; i < bbs.size(); i++)
        {
            Info<< "\t                       " << bbs[i]<< nl;
        }
        Info<< "\tSurface bounding box : " << surfBb[procI] << nl
            << "\tTriangles            : " << nFaces[procI] << nl
            << "\tVertices             : " << nPoints[procI]
            << endl;
    }
    Info<< endl;
}

void Foam::momentOfInertia::massPropertiesSolid
(
    const triSurface& surf,
    scalar density,
    scalar& mass,
    vector& cM,
    tensor& J
)
{
    triFaceList faces(surf.size());

    forAll(surf, i)
    {
        faces[i] = triFace(surf[i]);
    }

scalarField calcCurvature(const triSurface& surf)
{
    scalarField k(surf.points().size(), 0);

    Polyhedron P;

    buildCGALPolyhedron convert(surf);
    P.delegate(convert);

    // Info<< "Created CGAL Polyhedron with " << label(P.size_of_vertices())
    //     << " vertices and " << label(P.size_of_facets())
    //     << " facets. " << endl;

    // The rest of this function adapted from
    //     CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp

     //Vertex property map, with std::map
    typedef std::map<Vertex*, int> Vertex2int_map_type;
    typedef boost::associative_property_map< Vertex2int_map_type >
        Vertex_PM_type;
    typedef T_PolyhedralSurf_rings<Polyhedron, Vertex_PM_type > Poly_rings;

    typedef CGAL::Monge_via_jet_fitting<Kernel>         Monge_via_jet_fitting;
    typedef Monge_via_jet_fitting::Monge_form           Monge_form;

    std::vector<Point_3> in_points;  //container for data points

    // default parameter values and global variables
    unsigned int d_fitting = 2;
    unsigned int d_monge = 2;
    unsigned int min_nb_points = (d_fitting + 1)*(d_fitting + 2)/2;

    //initialize the tag of all vertices to -1
    Vertex_iterator vitb = P.vertices_begin();
    Vertex_iterator vite = P.vertices_end();

    Vertex2int_map_type vertex2props;
    Vertex_PM_type vpm(vertex2props);

    CGAL_For_all(vitb, vite)
    {
        put(vpm, &(*vitb), -1);
    }

Foam::labelList Foam::orientedSurface::faceToEdge
(
    const triSurface& s,
    const labelList& changedFaces
)
{
    labelList changedEdges(3*changedFaces.size());
    label changedI = 0;

    forAll(changedFaces, i)
    {
        const labelList& fEdges = s.faceEdges()[changedFaces[i]];

        forAll(fEdges, j)
        {
            changedEdges[changedI++] = fEdges[j];
        }
    }
    changedEdges.setSize(changedI);

    return changedEdges;
}

// Does face use valid vertices?
bool validTri
(
    const bool verbose,
    const triSurface& surf,
    const label faceI
)
{
    // Simple check on indices ok.

    const labelledTri& f = surf[faceI];

    forAll(f, fp)
    {
        if (f[fp] < 0 || f[fp] >= surf.points().size())
        {
            WarningIn("validTri(const triSurface&, const label)")
                << "triangle " << faceI << " vertices " << f
                << " uses point indices outside point range 0.."
                << surf.points().size()-1 << endl;
            return false;
        }
    }

    if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
    {
        WarningIn("validTri(const triSurface&, const label)")
            << "triangle " << faceI
            << " uses non-unique vertices " << f
            << " coords:" << f.points(surf.points())
            << endl;
        return false;
    }

    // duplicate triangle check

    const labelList& fFaces = surf.faceFaces()[faceI];

    // Check if faceNeighbours use same points as this face.
    // Note: discards normal information - sides of baffle are merged.
    forAll(fFaces, i)
    {
        label nbrFaceI = fFaces[i];

        if (nbrFaceI <= faceI)
        {
            // lower numbered faces already checked
            continue;
        }

        const labelledTri& nbrF = surf[nbrFaceI];

        if
        (
            ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
         && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
         && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
        )
        {
            WarningIn("validTri(const triSurface&, const label)")
                << "triangle " << faceI << " vertices " << f
                << " has the same vertices as triangle " << nbrFaceI
                << " vertices " << nbrF
                << " coords:" << f.points(surf.points())
                << endl;

            return false;
        }
    }

// Update intersections for selected edges.
void Foam::edgeIntersections::intersectEdges
(
    const triSurface& surf1,
    const pointField& points1,          // surf1 meshPoints (not localPoints!)
    const triSurfaceSearch& querySurf2,
    const scalarField& surf1PointTol,   // surf1 tolerance per point
    const labelList& edgeLabels
)
{
    const triSurface& surf2 = querySurf2.surface();
    const vectorField& normals2 = surf2.faceNormals();

    const labelList& meshPoints = surf1.meshPoints();

    if (debug)
    {
        Pout<< "Calculating intersection of " << edgeLabels.size() << " edges"
            << " out of " << surf1.nEdges() << " with " << surf2.size()
            << " triangles ..." << endl;
    }

    // Construct octree.
    const indexedOctree<treeDataTriSurface>& tree = querySurf2.tree();


    label nHits = 0;


    // Go through all edges, calculate intersections
    forAll(edgeLabels, i)
    {
        label edgeI = edgeLabels[i];

        if (debug && (i % 1000 == 0))
        {
            Pout<< "Intersecting edge " << edgeI << " with surface" << endl;
        }

        const edge& e = surf1.edges()[edgeI];

        const point& pStart = points1[meshPoints[e.start()]];
        const point& pEnd = points1[meshPoints[e.end()]];

        const vector eVec(pEnd - pStart);
        const scalar eMag = mag(eVec);
        const vector n(eVec/(eMag + VSMALL));

        // Smallish length for intersection calculation errors.
        const point tolVec = 1e-6*eVec;

        // Start tracking somewhat before pStart and upto somewhat after p1.
        // Note that tolerances here are smaller than those used to classify
        // hit below.
        // This will cause this hit to be marked as degenerate and resolved
        // later on.
        point p0 = pStart - 0.5*surf1PointTol[e[0]]*n;
        const point p1 = pEnd + 0.5*surf1PointTol[e[1]]*n;
        const scalar maxS = mag(p1 - pStart);

        // Get all intersections of the edge with the surface

        DynamicList<pointIndexHit> currentIntersections(100);
        DynamicList<label> currentIntersectionTypes(100);

        while (true)
        {
            pointIndexHit pHit = tree.findLine(p0, p1);

            if (pHit.hit())
            {
                nHits++;

                currentIntersections.append(pHit);

                // Classify point on surface1 edge.
                label edgeEnd = -1;

                if (mag(pHit.hitPoint() - pStart) < surf1PointTol[e[0]])
                {
                    edgeEnd = 0;
                }
                else if (mag(pHit.hitPoint() - pEnd) < surf1PointTol[e[1]])
                {
                    edgeEnd = 1;
                }
                else if (mag(n & normals2[pHit.index()]) < alignedCos_)
                {
                    Pout<< "Flat angle edge:" << edgeI
                        << " face:" << pHit.index()
                        << " cos:" << mag(n & normals2[pHit.index()])
                        << endl;
                    edgeEnd = 2;
                }

                currentIntersectionTypes.append(edgeEnd);

                if (edgeEnd == 1)
                {
                    // Close to end
//.........这里部分代码省略.........

// Checks if there exists a special topological situation that causes
// edge and the face it hit not to be recognized.
//
// For now if the face shares a point with the edge
bool Foam::surfaceIntersection::excludeEdgeHit
(
    const triSurface& surf,
    const label edgeI,
    const label faceI,
    const scalar
)
{
    const labelledTri& f = surf.localFaces()[faceI];

    const edge& e = surf.edges()[edgeI];

    if
    (
        (f[0] == e.start())
     || (f[0] == e.end())
     || (f[1] == e.start())
     || (f[1] == e.end())
     || (f[2] == e.start())
     || (f[2] == e.end())
    )
    {
        return true;

//        // Get edge vector
//        vector eVec = e.vec(surf.localPoints());
//        eVec /= mag(eVec) + VSMALL;
//
//        const labelList& eLabels = surf.faceEdges()[faceI];
//
//        // Get edge vector of 0th edge of face
//        vector e0Vec = surf.edges()[eLabels[0]].vec(surf.localPoints());
//        e0Vec /= mag(e0Vec) + VSMALL;
//
//        vector n = e0Vec ^ eVec;
//
//        if (mag(n) < SMALL)
//        {
//            // e0 is aligned with e. Choose next edge of face.
//            vector e1Vec = surf.edges()[eLabels[1]].vec(surf.localPoints());
//            e1Vec /= mag(e1Vec) + VSMALL;
//
//            n = e1Vec ^ eVec;
//
//            if (mag(n) < SMALL)
//            {
//                // Problematic triangle. Two edges aligned with edgeI. Give
//                // up.
//                return true;
//            }
//        }
//
//        // Check if same as faceNormal
//        if (mag(n & surf.faceNormals()[faceI]) > 1-tol)
//        {
//
//            Pout<< "edge:" << e << "  face:" << faceI
//                << "  e0Vec:" << e0Vec << "  n:" << n
//                << "  normalComponent:" << (n & surf.faceNormals()[faceI])
//                << "  tol:" << tol << endl;
//
//            return true;
//        }
//        else
//        {
//            return false;
//        }
    }
    else
    {
        return false;
    }
}

// Collapses small edge to point, thus removing triangle.
label collapseEdge(triSurface& surf, const scalar minLen)
{
    label nTotalCollapsed = 0;

    while (true)
    {
        const pointField& localPoints = surf.localPoints();
        const List<labelledTri>& localFaces = surf.localFaces();


        // Mapping from old to new points
        labelList pointMap(surf.nPoints());
        forAll(pointMap, i)
        {
            pointMap[i] = i;
        }

        // Storage for new points.
        pointField newPoints(localPoints);

        // To protect neighbours of collapsed faces.
        boolList okToCollapse(surf.size(), true);
        label nCollapsed = 0;

        forAll(localFaces, faceI)
        {
            if (okToCollapse[faceI])
            {
                // Check edge lengths.
                const triSurface::FaceType& f = localFaces[faceI];

                forAll(f, fp)
                {
                    label v = f[fp];
                    label v1 = f[f.fcIndex(fp)];

                    if (mag(localPoints[v1] - localPoints[v]) < minLen)
                    {
                        // Collapse f[fp1] onto f[fp].
                        pointMap[v1] = v;
                        newPoints[v] = 0.5*(localPoints[v1] + localPoints[v]);

                        Pout<< "Collapsing triange " << faceI << " to edge mid "
                            << newPoints[v] << endl;

                        nCollapsed++;
                        okToCollapse[faceI] = false;

                        // Protect point neighbours from collapsing.
                        markPointNbrs(surf, faceI, false, okToCollapse);

                        break;
                    }
                }
            }
        }

        Pout<< "collapseEdge : collapsing " << nCollapsed << " triangles"
            << endl;

        nTotalCollapsed += nCollapsed;

        if (nCollapsed == 0)
        {
            break;
        }

        // Pack the triangles
        surf = pack(surf, newPoints, pointMap);
    }

Foam::labelList Foam::orientedSurface::edgeToFace
(
    const triSurface& s,
    const labelList& changedEdges,
    labelList& flip
)
{
    labelList changedFaces(2*changedEdges.size());
    label changedI = 0;

    // 1.6.x merge: using local faces.  Reconsider
    // Rewrite uses cached local faces for efficiency
    // HJ, 24/Aug/2010
    const List<labelledTri> lf =  s.localFaces();

    forAll(changedEdges, i)
    {
        label edgeI = changedEdges[i];

        const labelList& eFaces = s.edgeFaces()[edgeI];

        if (eFaces.size() < 2)
        {
            // Do nothing, faces was already visited.
        }
        else if (eFaces.size() == 2)
        {
            label face0 = eFaces[0];
            label face1 = eFaces[1];

            const labelledTri& f0 = lf[face0];
            const labelledTri& f1 = lf[face1];

            // Old.  HJ, 24/Aug/2010
//            const labelledTri& f0 = s[face0];
//            const labelledTri& f1 = s[face1];

            if (flip[face0] == UNVISITED)
            {
                if (flip[face1] == UNVISITED)
                {
                    FatalErrorIn("orientedSurface::edgeToFace") << "Problem"
                            << abort(FatalError);
                }
                else
                {
                    // Face1 has a flip state, face0 hasn't
                    if (consistentEdge(s.edges()[edgeI], f0, f1))
                    {
                        // Take over flip status
                        flip[face0] = (flip[face1] == FLIP ? FLIP : NOFLIP);
                    }
                    else
                    {
                        // Invert
                        flip[face0] = (flip[face1] == FLIP ? NOFLIP : FLIP);
                    }
                    changedFaces[changedI++] = face0;
                }
            }
            else
            {
                if (flip[face1] == UNVISITED)
                {
                    // Face0 has a flip state, face1 hasn't
                    if (consistentEdge(s.edges()[edgeI], f0, f1))
                    {
                        flip[face1] = (flip[face0] == FLIP ? FLIP : NOFLIP);
                    }
                    else
                    {
                        flip[face1] = (flip[face0] == FLIP ? NOFLIP : FLIP);
                    }
                    changedFaces[changedI++] = face1;
                }
            }
        }
        else
        {
            // Multiply connected. Do what?
        }
    }