43 const Foam::vector Foam::boundaryMesh::splitNormal_(3, 2, 1);
46 const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1
E-2;
52 Foam::label Foam::boundaryMesh::nFeatureEdges(label pointI)
const
60 label edgeI = pEdges[pEdgeI];
62 if (edgeToFeature_[edgeI] != -1)
72 Foam::label Foam::boundaryMesh::nextFeatureEdge
82 label nbrEdgeI = pEdges[pEdgeI];
84 if (nbrEdgeI != edgeI)
86 label featI = edgeToFeature_[nbrEdgeI];
105 const label startEdgeI,
111 label edgeI = startEdgeI;
115 label vertI = e.start();
117 while (!isFeaturePoint[vertI])
121 edgeI = nextFeatureEdge(edgeI, vertI);
123 if ((edgeI == -1) || (edgeI == startEdgeI))
132 vertI = e.otherVertex(vertI);
144 labelList featLabels(featureEdges_.size());
146 label featLabelI = 0;
148 label initEdgeI = edgeI;
153 label featI = edgeToFeature_[edgeI];
160 featLabels[featLabelI++] = featI;
162 featVisited[featI] =
true;
168 vertI = e.otherVertex(vertI);
172 edgeI = nextFeatureEdge(edgeI, vertI);
174 if ((edgeI == -1) || (edgeI == initEdgeI))
179 while (!isFeaturePoint[vertI]);
183 featLabels.setSize(featLabelI);
189 void Foam::boundaryMesh::markEdges
191 const label maxDistance,
195 DynamicList<label>& visited
198 if (distance < maxDistance)
202 if (minDistance[edgeI] == -1)
205 visited.append(edgeI);
207 else if (minDistance[edgeI] <= distance)
220 forAll(startEdges, pEdgeI)
250 Foam::label Foam::boundaryMesh::findPatchID
253 const word& patchName
258 if (patches[patchI].
name() == patchName)
274 names[patchI] = patches_[patchI].name();
280 Foam::label Foam::boundaryMesh::whichPatch
288 const polyPatch& pp = patches[patchI];
290 if ((faceI >= pp.start()) && (faceI < (pp.start() + pp.size())))
314 label faceI = changedFaces[i];
320 label edgeI = fEdges[fEdgeI];
322 if (!regionEdge[edgeI] && (edgeRegion[edgeI] == -1))
324 edgeRegion[edgeI] = region;
326 changedEdges[changedI++] = edgeI;
331 changedEdges.setSize(changedI);
350 label edgeI = changedEdges[i];
356 label faceI = eFaces[eFaceI];
358 if (faceRegion[faceI] == -1)
360 faceRegion[faceI] = region;
362 changedFaces[changedI++] = faceI;
367 changedFaces.setSize(changedI);
374 void Foam::boundaryMesh::markZone
382 faceZone[faceI] = currentZone;
394 changedEdges = faceToEdge
404 Pout<<
"From changedFaces:" << changedFaces.size()
405 <<
" to changedEdges:" << changedEdges.size()
409 if (changedEdges.empty())
414 changedFaces = edgeToFace(currentZone, changedEdges, faceZone);
418 Pout<<
"From changedEdges:" << changedEdges.size()
419 <<
" to changedFaces:" << changedFaces.size()
423 if (changedFaces.empty())
473 patches_.setSize(mesh.boundaryMesh().size());
476 label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
480 meshFace_.setSize(nBFaces);
485 forAll(mesh.boundaryMesh(), patchI)
487 const polyPatch& pp = mesh.boundaryMesh()[patchI];
505 meshFace_[bFaceI] = pp.start() + patchFaceI;
507 bFaces[bFaceI] = pp[patchFaceI];
516 Pout<<
"read : patches now:" <<
endl;
520 const boundaryPatch& bp = patches_[patchI];
522 Pout<<
" name : " << bp.
name() << endl
523 <<
" size : " << bp.size() << endl
524 <<
" start : " << bp.start() << endl
525 <<
" type : " << bp.physicalType() << endl
535 PrimitivePatch<face, List, const pointField&> globalPatch
544 meshPtr_ =
new bMesh(globalPatch.localFaces(), globalPatch.localPoints());
549 const bMesh& msh = *meshPtr_;
551 Pout<<
"** Start of Faces **" <<
endl;
555 const face&
f = msh[faceI];
561 ctr += msh.points()[f[fp]];
571 Pout<<
"** End of Faces **" <<
endl;
573 Pout<<
"** Start of Points **" <<
endl;
575 forAll(msh.points(), pointI)
578 <<
" coord:" << msh.points()[pointI]
582 Pout<<
"** End of Points **" <<
endl;
586 featurePoints_.setSize(0);
587 featureEdges_.setSize(0);
589 featureToEdge_.setSize(0);
590 edgeToFeature_.setSize(meshPtr_->nEdges());
593 featureSegments_.setSize(0);
595 extraEdges_.setSize(0);
601 triSurface surf(fName);
609 SortableList<label> regions(surf.size());
613 regions[triI] = surf[triI].region();
618 Map<label> regionToBoundaryPatch;
620 label oldRegion = -1111;
621 label boundPatch = 0;
625 if (regions[i] != oldRegion)
627 regionToBoundaryPatch.insert(regions[i], boundPatch);
629 oldRegion = regions[i];
638 if (surfPatches.size() == regionToBoundaryPatch.size())
643 patches_.setSize(surfPatches.size());
646 forAll(surfPatches, patchI)
648 const geometricSurfacePatch& surfPatch = surfPatches[patchI];
659 surfPatch.geometricType()
668 patches_.setSize(regionToBoundaryPatch.size());
677 "patch" +
name(patchI),
691 const labelList& indices = regions.indices();
695 meshFace_.
setSize(surf.size());
700 label surfRegion = regions[0];
701 label foamRegion = regionToBoundaryPatch[surfRegion];
703 Pout<<
"Surface region " << surfRegion <<
" becomes boundary patch "
704 << foamRegion <<
" with name " << patches_[foamRegion].
name() <<
endl;
708 label startFaceI = 0;
712 label triI = indices[indexI];
714 const labelledTri& tri = surf.localFaces()[triI];
716 if (tri.region() != surfRegion)
719 boundaryPatch& bp = patches_[foamRegion];
721 bp.size() = bFaceI - startFaceI;
722 bp.start() = startFaceI;
724 surfRegion = tri.region();
725 foamRegion = regionToBoundaryPatch[surfRegion];
727 Pout<<
"Surface region " << surfRegion <<
" becomes boundary patch "
728 << foamRegion <<
" with name " << patches_[foamRegion].
name()
734 meshFace_[bFaceI] = triI;
736 bFaces[bFaceI++] = face(tri);
740 boundaryPatch& bp = patches_[foamRegion];
742 bp.size() = bFaceI - startFaceI;
743 bp.start() = startFaceI;
752 meshPtr_ =
new bMesh(bFaces, surf.localPoints());
755 featurePoints_.setSize(0);
756 featureEdges_.setSize(0);
758 featureToEdge_.setSize(0);
759 edgeToFeature_.setSize(meshPtr_->nEdges());
762 featureSegments_.setSize(0);
764 extraEdges_.setSize(0);
774 const boundaryPatch& bp = patches_[patchI];
776 surfPatches[patchI] =
777 geometricSurfacePatch
792 label totalNTris = getNTris(0,
mesh().size(), nTris);
801 startTri[faceI] = triI;
803 triI += nTris[faceI];
809 triangulate(0,
mesh().size(), totalNTris, triVerts);
813 List<labelledTri> tris(totalNTris);
819 const boundaryPatch& bp = patches_[patchI];
823 label faceI = bp.start() + patchFaceI;
825 label triVertI = 3*startTri[faceI];
827 for (label faceTriI = 0; faceTriI < nTris[faceI]; faceTriI++)
829 label v0 = triVerts[triVertI++];
830 label
v1 = triVerts[triVertI++];
831 label
v2 = triVerts[triVertI++];
833 tris[triI++] = labelledTri(v0, v1, v2, patchI);
838 triSurface surf(tris, surfPatches,
mesh().
points());
840 OFstream surfStream(fName);
842 surf.write(surfStream);
859 const primitiveMesh& pMesh,
867 DynamicList<label> leftFaces(
mesh().size()/2);
868 DynamicList<label> rightFaces(
mesh().size()/2);
872 scalar
sign =
mesh().faceNormals()[bFaceI] & splitNormal_;
876 rightFaces.
append(bFaceI);
880 leftFaces.append(bFaceI);
889 Pout<<
"getNearest :"
890 <<
" rightBin:" << rightFaces.size()
891 <<
" leftBin:" << leftFaces.size()
897 UIndirectList<face>(
mesh(), leftFaces),
902 UIndirectList<face>(
mesh(), rightFaces),
908 treeBoundBox overallBb(
mesh().localPoints());
912 scalar tol = 1
E-6 * overallBb.avgDim();
927 treeDataPrimitivePatch<face, UIndirectList, const pointField&>
930 treeDataPrimitivePatch<face, UIndirectList, const pointField&>
942 treeDataPrimitivePatch<face, UIndirectList, const pointField&>
945 treeDataPrimitivePatch<face, UIndirectList, const pointField&>
958 Pout<<
"getNearest : built trees" <<
endl;
969 labelList nearestBFaceI(pMesh.nFaces() - pMesh.nInternalFaces());
971 treeBoundBox tightest;
973 const scalar searchDimSqr =
magSqr(searchSpan);
975 forAll(nearestBFaceI, patchFaceI)
977 label meshFaceI = pMesh.nInternalFaces() + patchFaceI;
979 const point& ctr = pMesh.faceCentres()[meshFaceI];
981 if (debug && (patchFaceI % 1000) == 0)
983 Pout<<
"getNearest : patchFace:" << patchFaceI
984 <<
" meshFaceI:" << meshFaceI <<
" ctr:" << ctr <<
endl;
989 vector n = pMesh.faceAreas()[meshFaceI];
990 scalar area =
mag(n);
993 scalar typDim = -GREAT;
994 const face& f = pMesh.faces()[meshFaceI];
998 typDim =
max(typDim,
mag(pMesh.points()[f[fp]] - ctr));
1002 pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr);
1006 pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr);
1008 if (rightInfo.hit())
1013 label rightFaceI = rightFaces[rightInfo.index()];
1014 label leftFaceI = leftFaces[leftInfo.index()];
1016 label rightDist =
mag(rightInfo.hitPoint()-ctr);
1017 label leftDist =
mag(leftInfo.hitPoint()-ctr);
1019 scalar rightSign = n & ns[rightFaceI];
1020 scalar leftSign = n & ns[leftFaceI];
1024 (rightSign > 0 && leftSign > 0)
1025 || (rightSign < 0 && leftSign < 0)
1029 if (rightDist < leftDist)
1031 nearestBFaceI[patchFaceI] = rightFaceI;
1035 nearestBFaceI[patchFaceI] = leftFaceI;
1047 typDim *= distanceTol_;
1049 if (rightDist < typDim && leftDist < typDim)
1054 nearestBFaceI[patchFaceI] = rightFaceI;
1058 nearestBFaceI[patchFaceI] = leftFaceI;
1064 if (rightDist < leftDist)
1066 nearestBFaceI[patchFaceI] = rightFaceI;
1070 nearestBFaceI[patchFaceI] = leftFaceI;
1079 label rightFaceI = rightFaces[rightInfo.index()];
1080 nearestBFaceI[patchFaceI] = rightFaceI;
1091 nearestBFaceI[patchFaceI] = leftFaces[leftInfo.index()];
1096 nearestBFaceI[patchFaceI] = -1;
1101 return nearestBFaceI;
1108 const polyBoundaryMesh& oldPatches,
1118 HashTable<label> nameToIndex(2*patches_.size());
1120 Map<word> indexToName(2*patches_.size());
1123 label nNewPatches = patches_.size();
1125 forAll(oldPatches, oldPatchI)
1127 const polyPatch& patch = oldPatches[oldPatchI];
1129 label newPatchI = findPatchID(patch.name());
1131 if (newPatchI != -1)
1133 nameToIndex.insert(patch.name(), newPatchI);
1135 indexToName.insert(newPatchI, patch.name());
1141 forAll(patches_, bPatchI)
1143 const boundaryPatch& bp = patches_[bPatchI];
1145 if (!nameToIndex.found(bp.name()))
1147 nameToIndex.insert(bp.name(), bPatchI);
1149 indexToName.insert(bPatchI, bp.name());
1159 List<polyPatch*> newPatchPtrList(nNewPatches);
1161 label meshFaceI = newMesh.nInternalFaces();
1164 label facesToBeDone = newMesh.nFaces() - newMesh.nInternalFaces();
1166 forAll(patches_, bPatchI)
1168 const boundaryPatch& bp = patches_[bPatchI];
1170 label newPatchI = nameToIndex[bp.name()];
1173 label oldPatchI = findPatchID(oldPatches, bp.name());
1175 if (oldPatchI == -1)
1180 Pout<<
"patchify : Creating new polyPatch:" << bp.
name()
1181 <<
" type:" << bp.physicalType() <<
endl;
1191 newMesh.boundaryMesh()
1194 meshFaceI += facesToBeDone;
1202 const polyPatch& oldPatch = oldPatches[oldPatchI];
1206 Pout<<
"patchify : Cloning existing polyPatch:"
1210 newPatchPtrList[newPatchI] = oldPatch.clone
1212 newMesh.boundaryMesh(),
1218 meshFaceI += facesToBeDone;
1228 Pout<<
"Patchify : new polyPatch list:" <<
endl;
1230 forAll(newPatchPtrList, patchI)
1232 const polyPatch& newPatch = *newPatchPtrList[patchI];
1236 Pout<<
"polyPatch:" << newPatch.
name() << endl
1237 <<
" type :" << newPatch.
typeName << endl
1238 <<
" size :" << newPatch.size() << endl
1239 <<
" start:" << newPatch.start() << endl
1240 <<
" index:" << patchI <<
endl;
1246 repatchPolyTopoChanger polyMeshRepatcher(newMesh);
1247 polyMeshRepatcher.changePatches(newPatchPtrList);
1253 if (newPatchPtrList.size())
1255 List<DynamicList<label> > patchFaces(nNewPatches);
1259 (newMesh.nFaces() - newMesh.nInternalFaces())
1262 forAll(patchFaces, newPatchI)
1264 patchFaces[newPatchI].setCapacity(nAvgFaces);
1272 forAll(oldPatches, oldPatchI)
1274 const polyPatch& patch = oldPatches[oldPatchI];
1276 forAll(patch, patchFaceI)
1280 label meshFaceI = patch.start() + patchFaceI;
1282 label bFaceI = meshFaceI - newMesh.nInternalFaces();
1284 patchFaces[whichPatch(nearest[bFaceI])].append(meshFaceI);
1288 forAll(patchFaces, newPatchI)
1290 patchFaces[newPatchI].shrink();
1297 for (label newPatchI = 1; newPatchI < patchFaces.size(); newPatchI++)
1303 polyMeshRepatcher.changePatchID(pFaces[pFaceI], newPatchI);
1307 polyMeshRepatcher.repatch();
1314 edgeToFeature_.setSize(
mesh().nEdges());
1316 edgeToFeature_ = -1;
1321 featureToEdge_.setSize(
mesh().nEdges());
1325 if (minCos >= 0.9999)
1330 edgeToFeature_[edgeI] = featureI;
1331 featureToEdge_[featureI++] = edgeI;
1340 if (eFaces.size() == 2)
1342 label face0I = eFaces[0];
1344 label face1I = eFaces[1];
1355 const vector& n0 =
mesh().faceNormals()[face0I];
1357 const vector& n1 =
mesh().faceNormals()[face1I];
1359 float cosAng = n0 & n1;
1361 if (cosAng < minCos)
1363 edgeToFeature_[edgeI] = featureI;
1364 featureToEdge_[featureI++] = edgeI;
1371 edgeToFeature_[edgeI] = featureI;
1372 featureToEdge_[featureI++] = edgeI;
1378 featureToEdge_.setSize(featureI);
1384 featureEdges_.setSize(featureI);
1391 forAll(featureToEdge_, fEdgeI)
1393 label edgeI = featureToEdge_[fEdgeI];
1397 label start = featToMeshPoint[e.start()];
1401 featToMeshPoint[e.start()] = featPtI;
1403 featurePoints_[featPtI] =
mesh().
points()[e.start()];
1410 label end = featToMeshPoint[e.end()];
1414 featToMeshPoint[e.end()] = featPtI;
1416 featurePoints_[featPtI] =
mesh().
points()[e.end()];
1424 featureEdges_[fEdgeI] = edge(start, end);
1428 featurePoints_.setSize(featPtI);
1440 forAll(featureToEdge_, featI)
1442 label edgeI = featureToEdge_[featI];
1446 if (nFeatureEdges(e.start()) != 2)
1448 isFeaturePoint[e.start()] =
true;
1451 if (nFeatureEdges(e.end()) != 2)
1453 isFeaturePoint[e.end()] =
true;
1463 DynamicList<labelList> segments;
1466 boolList featVisited(featureToEdge_.size(),
false);
1470 label startFeatI = -1;
1472 forAll(featVisited, featI)
1474 if (!featVisited[featI])
1482 if (startFeatI == -1)
1493 featureToEdge_[startFeatI],
1504 featureSegments_.setSize(segments.size());
1506 forAll(featureSegments_, segmentI)
1508 featureSegments_[segmentI] = segments[segmentI];
1518 DynamicList<label> visitedEdges;
1521 markEdges(8, edgeI, 0, minDistance, visitedEdges);
1524 extraEdges_.transfer(visitedEdges);
1528 Foam::label Foam::boundaryMesh::whichPatch(
const label faceI)
const
1532 const boundaryPatch& bp = patches_[patchI];
1534 if ((faceI >= bp.start()) && (faceI < (bp.start() + bp.size())))
1541 <<
"Cannot find face " << faceI <<
" in list of boundaryPatches "
1549 Foam::label Foam::boundaryMesh::findPatchID(
const word& patchName)
const
1553 if (patches_[patchI].
name() == patchName)
1565 patches_.setSize(patches_.size() + 1);
1569 label patchI = patches_.size()-1;
1571 boundaryPatch* bpPtr =
new boundaryPatch
1580 patches_.set(patchI, bpPtr);
1584 Pout<<
"addPatch : patches now:" <<
endl;
1588 const boundaryPatch& bp = patches_[patchI];
1590 Pout<<
" name : " << bp.
name() << endl
1591 <<
" size : " << bp.size() << endl
1592 <<
" start : " << bp.start() << endl
1593 <<
" type : " << bp.physicalType() << endl
1602 label delPatchI = findPatchID(patchName);
1604 if (delPatchI == -1)
1607 <<
"Can't find patch named " << patchName
1611 if (patches_[delPatchI].size())
1614 <<
"Trying to delete non-empty patch " << patchName
1615 << endl <<
"Current size:" << patches_[delPatchI].size()
1619 PtrList<boundaryPatch> newPatches(patches_.size() - 1);
1621 for (label patchI = 0; patchI < delPatchI; patchI++)
1623 newPatches.set(patchI, patches_[patchI].clone());
1628 for (label patchI = delPatchI + 1; patchI < patches_.size(); patchI++)
1630 newPatches.set(patchI - 1, patches_[patchI].clone());
1635 patches_ = newPatches;
1639 Pout<<
"deletePatch : patches now:" <<
endl;
1643 const boundaryPatch& bp = patches_[patchI];
1645 Pout<<
" name : " << bp.
name() << endl
1646 <<
" size : " << bp.size() << endl
1647 <<
" start : " << bp.start() << endl
1648 <<
" type : " << bp.physicalType() << endl
1657 const word& patchName,
1658 const word& patchType
1661 label changeI = findPatchID(patchName);
1665 FatalErrorIn(
"boundaryMesh::changePatchType(const word&, const word&)")
1666 <<
"Can't find patch named " << patchName
1675 PtrList<boundaryPatch> newPatches(patches_.size());
1679 if (patchI == changeI)
1682 const boundaryPatch& oldBp = patches_[patchI];
1684 boundaryPatch* bpPtr =
new boundaryPatch
1693 newPatches.set(patchI, bpPtr);
1698 newPatches.set(patchI, patches_[patchI].clone());
1702 patches_ = newPatches;
1712 if (patchIDs.size() !=
mesh().
size())
1714 FatalErrorIn(
"boundaryMesh::changeFaces(const labelList& patchIDs)")
1715 <<
"List of patchIDs not equal to number of faces." << endl
1716 <<
"PatchIDs size:" << patchIDs.size()
1727 label patchID = patchIDs[faceI];
1729 if (patchID < 0 || patchID >= patches_.size())
1731 FatalErrorIn(
"boundaryMesh::changeFaces(const labelList&)")
1732 <<
"PatchID " << patchID <<
" out of range"
1745 for (label patchI = 1; patchI < patches_.size(); patchI++)
1747 startFace[patchI] = startFace[patchI-1] + nFaces[patchI-1];
1751 PtrList<boundaryPatch> newPatches(patches_.size());
1755 const boundaryPatch& bp = patches_[patchI];
1770 patches_ = newPatches;
1774 Pout<<
"changeFaces : patches now:" <<
endl;
1778 const boundaryPatch& bp = patches_[patchI];
1780 Pout<<
" name : " << bp.
name() << endl
1781 <<
" size : " << bp.size() << endl
1782 <<
" start : " << bp.start() << endl
1783 <<
" type : " << bp.physicalType() << endl
1790 oldToNew.setSize(patchIDs.size());
1794 int patchID = patchIDs[faceI];
1796 oldToNew[faceI] = startFace[patchID]++;
1806 newFaces[oldToNew[faceI]] =
mesh()[faceI];
1807 newMeshFace[oldToNew[faceI]] = meshFace_[faceI];
1814 meshFace_.transfer(newMeshFace);
1821 meshPtr_ = newMeshPtr_;
1827 const face& f =
mesh()[faceI];
1835 const label startFaceI,
1840 label totalNTris = 0;
1842 nTris.setSize(nFaces);
1844 for (label i = 0; i < nFaces; i++)
1846 label faceNTris = getNTris(startFaceI + i);
1848 nTris[i] = faceNTris;
1850 totalNTris += faceNTris;
1860 const label startFaceI,
1862 const label totalNTris,
1867 triVerts.setSize(3*totalNTris);
1871 for (label i = 0; i < nFaces; i++)
1873 label faceI = startFaceI + i;
1875 const face& f =
mesh()[faceI];
1882 f.triangles(
mesh().
points(), nTri, triFaces);
1886 forAll(triFaces, triFaceI)
1888 const face& triF = triFaces[triFaceI];
1890 triVerts[vertI++] = triF[0];
1891 triVerts[vertI++] = triF[1];
1892 triVerts[vertI++] = triF[2];
1901 const label startFaceI,
1905 SubList<face> patchFaces(
mesh(), nFaces, startFaceI);
1909 return patch.nPoints();
1916 const label startFaceI,
1918 const label totalNTris,
1923 SubList<face> patchFaces(
mesh(), nFaces, startFaceI);
1928 localToGlobal = patch.meshPoints();
1931 triVerts.setSize(3*totalNTris);
1935 for (label i = 0; i < nFaces; i++)
1938 const face& f = patch.localFaces()[i];
1941 faceList triFaces(f.nTriangles(patch.localPoints()));
1945 f.triangles(patch.localPoints(), nTri, triFaces);
1949 forAll(triFaces, triFaceI)
1951 const face& triF = triFaces[triFaceI];
1953 triVerts[vertI++] = triF[0];
1954 triVerts[vertI++] = triF[1];
1955 triVerts[vertI++] = triF[2];
1964 const label seedFaceI,
1970 forAll(protectedEdges, i)
1972 protectedEdge[protectedEdges[i]] =
true;
1980 markZone(protectedEdge, seedFaceI, 0, currentZone);
1983 visited.setSize(
mesh().size());
1985 forAll(currentZone, faceI)
1987 if (currentZone[faceI] == 0)
1989 visited[faceI] =
true;
1993 visited[faceI] =
false;