// initialize the nodes of the outer face
// and the corresponding faces
void ComputeTricOrder::initOuterNodes(node v1, node v2){
	// set nodes v1 and v2 on the "base"
	m_v1 = v1;
	m_v2 = v2;

	adjEntry firstAdj = m_outerFace->firstAdj();
	// set firstAdj, so outerface is on the right
	if (m_pEmbedding->rightFace(firstAdj) == m_outerFace)
		firstAdj = firstAdj->cyclicSucc();

	 adjEntry adjRun = firstAdj;
	 // traverse all nodes of the outer face
	 do {
	 	 	node v = adjRun->theNode();
			// now traverse the faces f of v
			// and increase outv[f] and add v to outerNodes[f]
			for(adjEntry adjV : v->adjEntries){
				face f = m_pEmbedding->rightFace(adjV);
				if (f != m_outerFace){
					m_outv[f]++;
					m_outerNodes[f].pushBack(v);
				}
			}
			adjRun = adjRun->twin()->cyclicSucc();
	} while (adjRun != firstAdj);
}

    forAllConstIter(Map<label>, meshFaceMap, iter)
    {
        label faceI = iter.key();
        label localI = iter();

        // Replace points with duplicated ones.
        const face& fRegion = faceRegions[localI];
        const face& f = mesh_.faces()[faceI];

        newFace.setSize(f.size());
        forAll(f, fp)
        {
            label pointI = f[fp];

            Map<label>::const_iterator iter = meshPointMap.find(pointI);

            if (iter != meshPointMap.end())
            {
                // Point has been duplicated. Find correct one for my
                // region.

                // Get the regions and added points for this point
                const labelList& regions = pointRegions[iter()];
                const labelList& dupPoints = duplicates_[iter()];

                // Look up index of my region in the regions for this point
                label index = findIndex(regions, fRegion[fp]);
                // Get the corresponding added point
                newFace[fp] = dupPoints[index];
            }
            else
            {
                newFace[fp] = pointI;
            }
        }

bool Foam::extrudedMesh::sameOrder(const face& f, const edge& e)
{
    label i = findIndex(f, e[0]);

    label nextI = (i == f.size()-1 ? 0 : i+1);

    return f[nextI] == e[1];
}

// Calculate (normalized) edge vectors.
// edges[i] gives edge between point i+1 and i.
Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
(
    const face& f,
    const pointField& points
)
{
    tmp<vectorField> tedges(new vectorField(f.size()));
    vectorField& edges = tedges();

    forAll(f, i)
    {
        point thisPt = points[f[i]];
        point nextPt = points[f[f.fcIndex(i)]];

        vector vec(nextPt - thisPt);
        vec /= mag(vec) + VSMALL;

        edges[i] = vec;
    }

// Test face for (almost) convexeness. Allows certain convexity before
// complaining.
// For every two consecutive edges calculate the normal. If it differs too
// much from the average face normal complain.
bool Foam::combineFaces::convexFace
(
    const scalar minConcaveCos,
    const pointField& points,
    const face& f
)
{
    // Get outwards pointing normal of f.
    vector n = f.normal(points);
    n /= mag(n);

    // Get edge from f[0] to f[size-1];
    vector ePrev(points[f.first()] - points[f.last()]);
    scalar magEPrev = mag(ePrev);
    ePrev /= magEPrev + VSMALL;

    forAll(f, fp0)
    {
        // Get vertex after fp
        label fp1 = f.fcIndex(fp0);

        // Normalized vector between two consecutive points
        vector e10(points[f[fp1]] - points[f[fp0]]);
        scalar magE10 = mag(e10);
        e10 /= magE10 + VSMALL;

        if (magEPrev > SMALL && magE10 > SMALL)
        {
            vector edgeNormal = ePrev ^ e10;

            if ((edgeNormal & n) < 0)
            {
                // Concave. Check angle.
                if ((ePrev & e10) < minConcaveCos)
                {
                    return false;
                }
            }
        }

        ePrev = e10;
        magEPrev = magE10;
    }

void Foam::pointMVCWeight::calcWeights
(
    const Map<label>& toLocal,
    const face& f,
    const DynamicList<point>& u,
    const scalarField& dist,
    scalarField& weights
) const
{
    weights.setSize(toLocal.size());
    weights = 0.0;

    scalarField theta(f.size());

    // recompute theta, the theta computed previously are not robust
    forAll(f, j)
    {
        label jPlus1 = f.fcIndex(j);
        scalar l = mag(u[j] - u[jPlus1]);
        theta[j] = 2.0*Foam::asin(l/2.0);
    }

        face operator()
        (
            const face& x,
            const label offset
        ) const
        {
            face result(x.size());

            forAll(x, xI)
            {
                result[xI] = x[xI] + offset;
            }

void Foam::sweptFaceAreaWeightAMI<SourcePatch, TargetPatch>::writeFaceOBJ
(
    const face& f,
    const pointField& ps,
    const string& name
) const
{
    OFstream obj(typeName + "_" + name + ".obj");

    for (label i = 0; i < f.size(); ++ i)
    {
        const vector& x = ps[f[i]];
        obj << "v " << x.x() << ' ' << x.y() << ' ' << x.z() << endl;
    }

    obj << 'f';
    for (label i = 0; i < f.size(); ++ i)
    {
        obj << ' ' << i + 1;
    }
    obj << endl;
}

Foam::treeBoundBox Foam::treeDataPrimitivePatch<PatchType>::calcBb
(
    const pointField& points,
    const face& f
)
{
    treeBoundBox bb(points[f[0]], points[f[0]]);

    for (label fp = 1; fp < f.size(); fp++)
    {
        const point& p = points[f[fp]];

        bb.min() = min(bb.min(), p);
        bb.max() = max(bb.max(), p);
    }
    return bb;
}

// initialize the edges of the external face
// and the corresponding faces
void ComputeTricOrder::initOuterEdges()
{
	adjEntry firstAdj = m_outerFace->firstAdj();
	// set firstAdj, so outerface is on the right
	if (m_pEmbedding->rightFace(firstAdj) == m_outerFace)
		firstAdj = firstAdj->cyclicSucc();
	adjEntry adjRun = firstAdj;
	// traverse all edges of the outer face
	do {
		edge e = adjRun->theEdge();
		face f = m_pEmbedding->rightFace(adjRun);
		// verify that actual edge is not edge (v1,v2)
		if (!((e->source() == m_v1 && e->target() == m_v2) || ((e->source() == m_v2 && e->target() == m_v1)))){
			m_oute[f]++;
			m_outerEdges[f].pushBack(e);
		}
		adjRun = adjRun->twin()->cyclicSucc();
	} while (adjRun != firstAdj);
}

	void FixEdgeInserterCore::insertEdgesIntoDualAfterRemove(const CombinatorialEmbedding &E, face f)
	{
		node vRight = m_nodeOf[f];

		adjEntry adj1 = f->firstAdj(), adj = adj1;
		do {
			if(m_pForbidden && (*m_pForbidden)[m_pr.original(adj->theEdge())] == true)
				continue;

			node vLeft = m_nodeOf[E.leftFace(adj)];

			edge eLR = m_dual.newEdge(vLeft,vRight);
			m_primalAdj[eLR] = adj;

			edge eRL = m_dual.newEdge(vRight,vLeft);
			m_primalAdj[eRL] = adj->twin();
		}
		while((adj = adj->faceCycleSucc()) != adj1);
	}

void Foam::faceCracker::detachFaceCracker
(
    polyTopoChange& ref
) const
{
    if (debug)
    {
        Pout<< "void faceCracker::detachFaceCracker("
            << "polyTopoChange& ref) const "
            << " for object " << name() << " : "
            << "Detaching faces" << endl;
    }

    const polyMesh& mesh = topoChanger().mesh();
    const faceZoneMesh& zoneMesh = mesh.faceZones();

    const primitiveFacePatch& masterFaceLayer =
        zoneMesh[crackZoneID_.index()]();
    const pointField& points = mesh.points();
    const labelListList& meshEdgeFaces = mesh.edgeFaces();

    const labelList& mp = masterFaceLayer.meshPoints();
    const edgeList& zoneLocalEdges = masterFaceLayer.edges();

    const labelList& meshEdges = zoneMesh[crackZoneID_.index()].meshEdges();

    // Create the points

    labelList addedPoints(mp.size(), -1);

    // Go through boundary edges of the master patch.  If all the faces from
    // this patch are internal, mark the points in the addedPoints lookup
    // with their original labels to stop duplication
    label nIntEdges = masterFaceLayer.nInternalEdges();

    for
    (
        label curEdgeID = nIntEdges;
        curEdgeID < meshEdges.size();
        curEdgeID++
    )
    {
        const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];

        bool edgeIsInternal = true;

        forAll (curFaces, faceI)
        {
            if (!mesh.isInternalFace(curFaces[faceI]))
            {
                // The edge belongs to a boundary face
                edgeIsInternal = false;
                break;
            }
        }

        if (edgeIsInternal)
        {
            // Reset the point creation
            addedPoints[zoneLocalEdges[curEdgeID].start()] =
                mp[zoneLocalEdges[curEdgeID].start()];

            addedPoints[zoneLocalEdges[curEdgeID].end()] =
                mp[zoneLocalEdges[curEdgeID].end()];
        }
    }

    // Create new points for face zone
    forAll (addedPoints, pointI)
    {
        if (addedPoints[pointI] < 0)
        {
            addedPoints[pointI] =
                ref.setAction
                (
                    polyAddPoint
                    (
                        points[mp[pointI]],        // point
                        mp[pointI],                // master point
                        -1,                        // zone ID
                        true                       // supports a cell
                    )
                );
        }
    }

    // Modify faces in the master zone and duplicate for the slave zone

    const labelList& mf = zoneMesh[crackZoneID_.index()];
    const boolList& mfFlip = zoneMesh[crackZoneID_.index()].flipMap();
    const faceList& zoneFaces = masterFaceLayer.localFaces();

    const faceList& faces = mesh.faces();
    const labelList& own = mesh.faceOwner();
    const labelList& nei = mesh.faceNeighbour();

    forAll (mf, faceI)
    {
        const label curFaceID = mf[faceI];

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

 
 points[pointI] = p;
 }
 
 //- read the number of faces
 file >> counter;
 
 faceListPMG& faces = meshModifier.facesAccess();
 
 //- read faces from file
 faces.setSize(counter);
 forAll(faces, faceI)
 {
     file >> counter;
 
 face f;
 f.setSize(counter);
 
 forAll(f, pI)
   file >> f[pI];
 
 faces[faceI] = f.reverseFace();
 }
 
 //- read the number of cells
 file >> counter;
 
 //- read cells from file
 cellListPMG& cells = meshModifier.cellsAccess();
 cells.setSize(counter);
 

 inline unsigned Hash<face>::operator()(const face& t, unsigned seed) const
 {
     return Hasher(t.cdata(),t.size()*sizeof(label), seed);
 }

void Foam::attachDetach::detachInterface
(
    polyTopoChange& ref
) const
{
    // Algorithm:
    // 1. Create new points for points of the master face zone
    // 2. Modify all faces of the master zone, by putting them into the master
    //    patch (look for orientation) and their renumbered mirror images
    //    into the slave patch
    // 3. Create a point renumbering list, giving a new point index for original
    //    points in the face patch
    // 4. Grab all faces in cells on the master side and renumber them 
    //    using the point renumbering list.  Exclude the ones that belong to
    //    the master face zone
    //
    // Note on point creation:
    // In order to take into account the issues related to partial
    // blocking in an attach/detach mesh modifier, special treatment
    // is required for the duplication of points on the edge of the
    // face zone.  Points which are shared only by internal edges need
    // not to be duplicated, as this would propagate the discontinuity
    // in the mesh beyond the face zone.  Therefore, before creating
    // the new points, check the external edge loop.  For each edge
    // check if the edge is internal (i.e. does not belong to a
    // patch); if so, exclude both of its points from duplication.

    if (debug)
    {
        Pout<< "void attachDetach::detachInterface("
            << "polyTopoChange& ref) const "
            << " for object " << name() << " : "
            << "Detaching interface" << endl;
    }

    const polyMesh& mesh = topoChanger().mesh();
    const faceZoneMesh& zoneMesh = mesh.faceZones();

    const primitiveFacePatch& masterFaceLayer = zoneMesh[faceZoneID_.index()]();
    const pointField& points = mesh.points();
    const labelListList& meshEdgeFaces = mesh.edgeFaces();

    const labelList& mp = masterFaceLayer.meshPoints();
    const edgeList& zoneLocalEdges = masterFaceLayer.edges();

    const labelList& meshEdges = zoneMesh[faceZoneID_.index()].meshEdges();

    // Create the points

    labelList addedPoints(mp.size(), -1);

    // Go through boundary edges of the master patch.  If all the faces from
    // this patch are internal, mark the points in the addedPoints lookup
    // with their original labels to stop duplication
    label nIntEdges = masterFaceLayer.nInternalEdges();

    for (label curEdgeID = nIntEdges; curEdgeID < meshEdges.size(); curEdgeID++)
    {
        const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];

        bool edgeIsInternal = true;

        forAll (curFaces, faceI)
        {
            if (!mesh.isInternalFace(curFaces[faceI]))
            {
                // The edge belongs to a boundary face
                edgeIsInternal = false;
                break;
            }
        }

        if (edgeIsInternal)
        {
// Pout<< "Internal edge found: (" << mp[zoneLocalEdges[curEdgeID].start()] << " " << mp[zoneLocalEdges[curEdgeID].end()] << ")" << endl;

            // Reset the point creation
            addedPoints[zoneLocalEdges[curEdgeID].start()] =
                mp[zoneLocalEdges[curEdgeID].start()];

            addedPoints[zoneLocalEdges[curEdgeID].end()] =
                mp[zoneLocalEdges[curEdgeID].end()];
        }
    }
// Pout << "addedPoints before point creation: " << addedPoints << endl;

    // Create new points for face zone
    forAll (addedPoints, pointI)
    {
        if (addedPoints[pointI] < 0)
        {
            addedPoints[pointI] =
                ref.setAction
                (
                    polyAddPoint
                    (
                        points[mp[pointI]],        // point
                        mp[pointI],                // master point
                        -1,                        // zone ID
                        true                       // supports a cell
//.........这里部分代码省略.........