42 const Foam::vector Foam::boundaryMesh::splitNormal_(3, 2, 1);
45 const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1
E-2;
51 Foam::label Foam::boundaryMesh::nFeatureEdges(label pointI)
const
59 label edgeI = pEdges[pEdgeI];
61 if (edgeToFeature_[edgeI] != -1)
71 Foam::label Foam::boundaryMesh::nextFeatureEdge
81 label nbrEdgeI = pEdges[pEdgeI];
83 if (nbrEdgeI != edgeI)
85 label featI = edgeToFeature_[nbrEdgeI];
104 const label startEdgeI,
110 label edgeI = startEdgeI;
114 label vertI = e.start();
116 while (!isFeaturePoint[vertI])
120 edgeI = nextFeatureEdge(edgeI, vertI);
122 if ((edgeI == -1) || (edgeI == startEdgeI))
131 vertI = e.otherVertex(vertI);
143 labelList featLabels(featureEdges_.size());
145 label featLabelI = 0;
147 label initEdgeI = edgeI;
152 label featI = edgeToFeature_[edgeI];
159 featLabels[featLabelI++] = featI;
161 featVisited[featI] =
true;
167 vertI = e.otherVertex(vertI);
171 edgeI = nextFeatureEdge(edgeI, vertI);
173 if ((edgeI == -1) || (edgeI == initEdgeI))
178 while (!isFeaturePoint[vertI]);
182 featLabels.setSize(featLabelI);
188 void Foam::boundaryMesh::markEdges
190 const label maxDistance,
194 DynamicList<label>& visited
197 if (distance < maxDistance)
201 if (minDistance[edgeI] == -1)
204 visited.append(edgeI);
206 else if (minDistance[edgeI] <= distance)
219 forAll(startEdges, pEdgeI)
249 Foam::label Foam::boundaryMesh::findPatchID
252 const word& patchName
257 if (patches[patchI].
name() == patchName)
273 names[patchI] = patches_[patchI].name();
279 Foam::label Foam::boundaryMesh::whichPatch
289 if ((faceI >= pp.
start()) && (faceI < (pp.
start() + pp.size())))
313 label faceI = changedFaces[i];
319 label edgeI = fEdges[fEdgeI];
321 if (!regionEdge[edgeI] && (edgeRegion[edgeI] == -1))
323 edgeRegion[edgeI] = region;
325 changedEdges[changedI++] = edgeI;
330 changedEdges.setSize(changedI);
349 label edgeI = changedEdges[i];
355 label faceI = eFaces[eFaceI];
357 if (faceRegion[faceI] == -1)
359 faceRegion[faceI] = region;
361 changedFaces[changedI++] = faceI;
366 changedFaces.setSize(changedI);
373 void Foam::boundaryMesh::markZone
381 faceZone[faceI] = currentZone;
393 changedEdges = faceToEdge
403 Pout<<
"From changedFaces:" << changedFaces.size()
404 <<
" to changedEdges:" << changedEdges.size()
408 if (changedEdges.empty())
413 changedFaces = edgeToFace(currentZone, changedEdges, faceZone);
417 Pout<<
"From changedEdges:" << changedEdges.size()
418 <<
" to changedFaces:" << changedFaces.size()
422 if (changedFaces.empty())
479 meshFace_.setSize(nBFaces);
504 meshFace_[bFaceI] = pp.
start() + patchFaceI;
506 bFaces[bFaceI] = pp[patchFaceI];
515 Pout<<
"read : patches now:" <<
endl;
521 Pout<<
" name : " << bp.
name() << endl
522 <<
" size : " << bp.
size() << endl
523 <<
" start : " << bp.
start() << endl
543 meshPtr_ =
new bMesh(globalPatch.localFaces(), globalPatch.localPoints());
548 const bMesh& msh = *meshPtr_;
550 Pout<<
"** Start of Faces **" <<
endl;
554 const face&
f = msh[faceI];
560 ctr += msh.
points()[f[fp]];
570 Pout<<
"** End of Faces **" <<
endl;
572 Pout<<
"** Start of Points **" <<
endl;
577 <<
" coord:" << msh.
points()[pointI]
581 Pout<<
"** End of Points **" <<
endl;
585 featurePoints_.setSize(0);
586 featureEdges_.setSize(0);
588 featureToEdge_.setSize(0);
589 edgeToFeature_.setSize(meshPtr_->nEdges());
592 featureSegments_.setSize(0);
594 extraEdges_.setSize(0);
612 regions[triI] = surf[triI].region();
619 label oldRegion = -1111;
620 label boundPatch = 0;
624 if (regions[i] != oldRegion)
626 regionToBoundaryPatch.
insert(regions[i], boundPatch);
628 oldRegion = regions[i];
637 if (surfPatches.
size() == regionToBoundaryPatch.
size())
642 patches_.setSize(surfPatches.
size());
645 forAll(surfPatches, patchI)
676 "patch" +
name(patchI),
690 const labelList& indices = regions.indices();
699 label surfRegion = regions[0];
700 label foamRegion = regionToBoundaryPatch[surfRegion];
702 Pout<<
"Surface region " << surfRegion <<
" becomes boundary patch "
703 << foamRegion <<
" with name " << patches_[foamRegion].
name() <<
endl;
707 label startFaceI = 0;
711 label triI = indices[indexI];
715 if (tri.
region() != surfRegion)
720 bp.
size() = bFaceI - startFaceI;
721 bp.
start() = startFaceI;
723 surfRegion = tri.
region();
724 foamRegion = regionToBoundaryPatch[surfRegion];
726 Pout<<
"Surface region " << surfRegion <<
" becomes boundary patch "
727 << foamRegion <<
" with name " << patches_[foamRegion].
name()
733 meshFace_[bFaceI] = triI;
735 bFaces[bFaceI++] =
face(tri);
741 bp.
size() = bFaceI - startFaceI;
742 bp.
start() = startFaceI;
754 featurePoints_.setSize(0);
755 featureEdges_.setSize(0);
757 featureToEdge_.setSize(0);
758 edgeToFeature_.setSize(meshPtr_->nEdges());
761 featureSegments_.setSize(0);
763 extraEdges_.setSize(0);
775 surfPatches[patchI] =
791 label totalNTris = getNTris(0,
mesh().size(), nTris);
800 startTri[faceI] = triI;
802 triI += nTris[faceI];
808 triangulate(0,
mesh().size(), totalNTris, triVerts);
822 label faceI = bp.
start() + patchFaceI;
824 label triVertI = 3*startTri[faceI];
826 for (label faceTriI = 0; faceTriI < nTris[faceI]; faceTriI++)
828 label v0 = triVerts[triVertI++];
829 label
v1 = triVerts[triVertI++];
830 label
v2 = triVerts[triVertI++];
841 surf.write(surfStream);
871 scalar
sign =
mesh().faceNormals()[bFaceI] & splitNormal_;
875 rightFaces.
append(bFaceI);
888 Pout<<
"getNearest :"
889 <<
" rightBin:" << rightFaces.
size()
890 <<
" leftBin:" << leftFaces.
size()
900 scalar tol = 1
E-6 * overallBb.
avgDim();
940 Pout<<
"getNearest : built trees" <<
endl;
955 const scalar searchDim =
mag(searchSpan);
957 forAll(nearestBFaceI, patchFaceI)
963 if (debug && (patchFaceI % 1000) == 0)
965 Pout<<
"getNearest : patchFace:" << patchFaceI
966 <<
" meshFaceI:" << meshFaceI <<
" ctr:" << ctr <<
endl;
972 scalar area =
mag(n);
975 scalar typDim = -GREAT;
984 tightest.min() = ctr - searchSpan;
985 tightest.
max() = ctr + searchSpan;
986 scalar rightDist = searchDim;
987 label rightI = rightTree.
findNearest(ctr, tightest, rightDist);
992 tightest.min() = ctr - searchSpan;
993 tightest.
max() = ctr + searchSpan;
994 scalar leftDist = searchDim;
995 label leftI = leftTree.
findNearest(ctr, tightest, leftDist);
1005 nearestBFaceI[patchFaceI] = -1;
1011 nearestBFaceI[patchFaceI] = leftFaces[leftI];
1020 nearestBFaceI[patchFaceI] = rightFaces[rightI];
1026 scalar rightSign = n & ns[rightFaces[rightI]];
1027 scalar leftSign = n & ns[leftFaces[leftI]];
1031 (rightSign > 0 && leftSign > 0)
1032 || (rightSign < 0 && leftSign < 0)
1036 if (rightDist < leftDist)
1038 nearestBFaceI[patchFaceI] = rightFaces[rightI];
1042 nearestBFaceI[patchFaceI] = leftFaces[leftI];
1054 typDim *= distanceTol_;
1056 if (rightDist < typDim && leftDist < typDim)
1061 nearestBFaceI[patchFaceI] = rightFaces[rightI];
1065 nearestBFaceI[patchFaceI] = leftFaces[leftI];
1071 if (rightDist < leftDist)
1073 nearestBFaceI[patchFaceI] = rightFaces[rightI];
1077 nearestBFaceI[patchFaceI] = leftFaces[leftI];
1085 return nearestBFaceI;
1104 Map<word> indexToName(2*patches_.size());
1107 label nNewPatches = patches_.
size();
1109 forAll(oldPatches, oldPatchI)
1111 const polyPatch& patch = oldPatches[oldPatchI];
1113 label newPatchI = findPatchID(patch.
name());
1115 if (newPatchI != -1)
1117 nameToIndex.insert(patch.
name(), newPatchI);
1119 indexToName.insert(newPatchI, patch.
name());
1125 forAll(patches_, bPatchI)
1129 if (!nameToIndex.found(bp.
name()))
1131 nameToIndex.insert(bp.
name(), bPatchI);
1133 indexToName.insert(bPatchI, bp.
name());
1150 forAll(patches_, bPatchI)
1154 label newPatchI = nameToIndex[bp.
name()];
1157 label oldPatchI = findPatchID(oldPatches, bp.
name());
1159 if (oldPatchI == -1)
1164 Pout<<
"patchify : Creating new polyPatch:" << bp.
name()
1178 meshFaceI += facesToBeDone;
1186 const polyPatch& oldPatch = oldPatches[oldPatchI];
1190 Pout<<
"patchify : Cloning existing polyPatch:"
1194 newPatchPtrList[newPatchI] = oldPatch.
clone
1202 meshFaceI += facesToBeDone;
1212 Pout<<
"Patchify : new polyPatch list:" <<
endl;
1214 forAll(newPatchPtrList, patchI)
1216 const polyPatch& newPatch = *newPatchPtrList[patchI];
1220 Pout<<
"polyPatch:" << newPatch.
name() << endl
1221 <<
" type :" << newPatch.typeName << endl
1222 <<
" size :" << newPatch.size() << endl
1223 <<
" start:" << newPatch.
start() << endl
1224 <<
" index:" << patchI <<
endl;
1237 if (newPatchPtrList.
size())
1246 forAll(patchFaces, newPatchI)
1248 patchFaces[newPatchI].setCapacity(nAvgFaces);
1256 forAll(oldPatches, oldPatchI)
1258 const polyPatch& patch = oldPatches[oldPatchI];
1260 forAll(patch, patchFaceI)
1264 label meshFaceI = patch.
start() + patchFaceI;
1268 patchFaces[whichPatch(nearest[bFaceI])].
append(meshFaceI);
1272 forAll(patchFaces, newPatchI)
1274 patchFaces[newPatchI].shrink();
1281 for (label newPatchI = 1; newPatchI < patchFaces.
size(); newPatchI++)
1298 edgeToFeature_.setSize(
mesh().nEdges());
1300 edgeToFeature_ = -1;
1305 featureToEdge_.setSize(
mesh().nEdges());
1309 if (minCos >= 0.9999)
1314 edgeToFeature_[edgeI] = featureI;
1315 featureToEdge_[featureI++] = edgeI;
1324 if (eFaces.
size() == 2)
1326 label face0I = eFaces[0];
1328 label face1I = eFaces[1];
1339 const vector& n0 =
mesh().faceNormals()[face0I];
1341 const vector& n1 =
mesh().faceNormals()[face1I];
1343 float cosAng = n0 & n1;
1345 if (cosAng < minCos)
1347 edgeToFeature_[edgeI] = featureI;
1348 featureToEdge_[featureI++] = edgeI;
1355 edgeToFeature_[edgeI] = featureI;
1356 featureToEdge_[featureI++] = edgeI;
1362 featureToEdge_.setSize(featureI);
1368 featureEdges_.setSize(featureI);
1375 forAll(featureToEdge_, fEdgeI)
1377 label edgeI = featureToEdge_[fEdgeI];
1381 label start = featToMeshPoint[e.
start()];
1385 featToMeshPoint[e.
start()] = featPtI;
1394 label end = featToMeshPoint[e.
end()];
1398 featToMeshPoint[e.
end()] = featPtI;
1408 featureEdges_[fEdgeI] =
edge(start, end);
1412 featurePoints_.setSize(featPtI);
1424 forAll(featureToEdge_, featI)
1426 label edgeI = featureToEdge_[featI];
1430 if (nFeatureEdges(e.
start()) != 2)
1432 isFeaturePoint[e.
start()] =
true;
1435 if (nFeatureEdges(e.
end()) != 2)
1437 isFeaturePoint[e.
end()] =
true;
1450 boolList featVisited(featureToEdge_.size(),
false);
1454 label startFeatI = -1;
1456 forAll(featVisited, featI)
1458 if (!featVisited[featI])
1466 if (startFeatI == -1)
1477 featureToEdge_[startFeatI],
1488 featureSegments_.setSize(segments.
size());
1490 forAll(featureSegments_, segmentI)
1492 featureSegments_[segmentI] = segments[segmentI];
1505 markEdges(8, edgeI, 0, minDistance, visitedEdges);
1508 extraEdges_.transfer(visitedEdges);
1512 Foam::label Foam::boundaryMesh::whichPatch(
const label faceI)
const
1525 <<
"Cannot find face " << faceI <<
" in list of boundaryPatches "
1533 Foam::label Foam::boundaryMesh::findPatchID(
const word& patchName)
const
1537 if (patches_[patchI].
name() == patchName)
1549 patches_.setSize(patches_.size() + 1);
1553 label patchI = patches_.size()-1;
1564 patches_.set(patchI, bpPtr);
1568 Pout<<
"addPatch : patches now:" <<
endl;
1574 Pout<<
" name : " << bp.
name() << endl
1575 <<
" size : " << bp.
size() << endl
1576 <<
" start : " << bp.
start() << endl
1586 label delPatchI = findPatchID(patchName);
1588 if (delPatchI == -1)
1591 <<
"Can't find patch named " << patchName
1595 if (patches_[delPatchI].size())
1598 <<
"Trying to delete non-empty patch " << patchName
1599 <<
endl <<
"Current size:" << patches_[delPatchI].size()
1605 for (label patchI = 0; patchI < delPatchI; patchI++)
1607 newPatches.
set(patchI, patches_[patchI].clone());
1612 for (label patchI = delPatchI + 1; patchI < patches_.size(); patchI++)
1614 newPatches.set(patchI - 1, patches_[patchI].clone());
1619 patches_ = newPatches;
1623 Pout<<
"deletePatch : patches now:" <<
endl;
1629 Pout<<
" name : " << bp.
name() << endl
1630 <<
" size : " << bp.
size() << endl
1631 <<
" start : " << bp.
start() << endl
1641 const word& patchName,
1642 const word& patchType
1645 label changeI = findPatchID(patchName);
1649 FatalErrorIn(
"boundaryMesh::changePatchType(const word&, const word&)")
1650 <<
"Can't find patch named " << patchName
1663 if (patchI == changeI)
1677 newPatches.set(patchI, bpPtr);
1682 newPatches.set(patchI, patches_[patchI].clone());
1686 patches_ = newPatches;
1698 FatalErrorIn(
"boundaryMesh::changeFaces(const labelList& patchIDs)")
1699 <<
"List of patchIDs not equal to number of faces." <<
endl
1700 <<
"PatchIDs size:" << patchIDs.
size()
1711 label patchID = patchIDs[faceI];
1713 if (patchID < 0 || patchID >= patches_.size())
1715 FatalErrorIn(
"boundaryMesh::changeFaces(const labelList&)")
1716 <<
"PatchID " << patchID <<
" out of range"
1729 for (label patchI = 1; patchI < patches_.size(); patchI++)
1731 startFace[patchI] = startFace[patchI-1] + nFaces[patchI-1];
1754 patches_ = newPatches;
1758 Pout<<
"changeFaces : patches now:" <<
endl;
1764 Pout<<
" name : " << bp.
name() << endl
1765 <<
" size : " << bp.
size() << endl
1766 <<
" start : " << bp.
start() << endl
1778 int patchID = patchIDs[faceI];
1780 oldToNew[faceI] = startFace[patchID]++;
1790 newFaces[oldToNew[faceI]] =
mesh()[faceI];
1791 newMeshFace[oldToNew[faceI]] = meshFace_[faceI];
1798 meshFace_.transfer(newMeshFace);
1805 meshPtr_ = newMeshPtr_;
1819 const label startFaceI,
1824 label totalNTris = 0;
1828 for (label i = 0; i < nFaces; i++)
1830 label faceNTris = getNTris(startFaceI + i);
1832 nTris[i] = faceNTris;
1834 totalNTris += faceNTris;
1844 const label startFaceI,
1846 const label totalNTris,
1851 triVerts.
setSize(3*totalNTris);
1855 for (label i = 0; i < nFaces; i++)
1857 label faceI = startFaceI + i;
1870 forAll(triFaces, triFaceI)
1872 const face& triF = triFaces[triFaceI];
1874 triVerts[vertI++] = triF[0];
1875 triVerts[vertI++] = triF[1];
1876 triVerts[vertI++] = triF[2];
1885 const label startFaceI,
1900 const label startFaceI,
1902 const label totalNTris,
1915 triVerts.
setSize(3*totalNTris);
1919 for (label i = 0; i < nFaces; i++)
1933 forAll(triFaces, triFaceI)
1935 const face& triF = triFaces[triFaceI];
1937 triVerts[vertI++] = triF[0];
1938 triVerts[vertI++] = triF[1];
1939 triVerts[vertI++] = triF[2];
1948 const label seedFaceI,
1954 forAll(protectedEdges, i)
1956 protectedEdge[protectedEdges[i]] =
true;
1964 markZone(protectedEdge, seedFaceI, 0, currentZone);
1969 forAll(currentZone, faceI)
1971 if (currentZone[faceI] == 0)
1973 visited[faceI] =
true;
1977 visited[faceI] =
false;