59 polyTopoChange& meshMod
62 const face& f = mesh_.
faces()[faceI];
64 bool zoneFlip =
false;
68 const faceZone& fZone = mesh_.
faceZones()[zoneID];
69 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
96 <<
"No neighbour patch for internal face " << faceI
101 bool reverseFlip =
false;
104 reverseFlip = !zoneFlip;
107 dupFaceI = meshMod.setAction
128 void Foam::meshRefinement::getBafflePatches
137 autoPtr<OFstream> str;
139 if (debug&OBJINTERSECTIONS)
150 Pout<<
"getBafflePatches : Writing surface intersections to file " 162 ownPatch.setSize(mesh_.
nFaces());
164 neiPatch.setSize(mesh_.
nFaces());
181 label faceI = testFaces[i];
187 start[i] = cellCentres[own];
192 start[i] = cellCentres[own];
208 List<pointIndexHit> hit1;
211 List<pointIndexHit> hit2;
230 label faceI = testFaces[i];
232 if (hit1[i].hit() && hit2[i].hit())
244 str()<<
"l " << vertI-3 <<
' ' << vertI-2 <<
nl;
245 str()<<
"l " << vertI-2 <<
' ' << vertI-1 <<
nl;
246 str()<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
250 ownPatch[faceI] = globalToMasterPatch
254 neiPatch[faceI] = globalToMasterPatch
259 if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
282 const bool allowBoundary,
287 Map<labelPair> bafflePatch(mesh_.
nFaces()/1000);
289 const PtrList<surfaceZonesInfo>& surfZones = surfaces_.
surfZones();
294 const word& faceZoneName = surfZones[surfI].faceZoneName();
296 if (faceZoneName.size())
299 label zoneI = fZones.findZoneID(faceZoneName);
301 const faceZone& fZone = fZones[zoneI];
307 globalToMasterPatch[globalRegionI],
308 globalToSlavePatch[globalRegionI]
311 Info<<
"For zone " << fZone.name() <<
" found patches " 318 label faceI = fZone[i];
323 if (fZone.flipMap()[i])
328 if (!bafflePatch.insert(faceI, patches))
333 <<
" in zone " << fZone.name()
334 <<
" is in multiple zones!" 359 <<
" ownPatch:" << ownPatch.
size()
360 <<
" neiPatch:" << neiPatch.
size()
361 <<
". Should be number of faces:" << mesh_.
nFaces()
372 forAll(syncedOwnPatch, faceI)
376 (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
377 || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
381 <<
"Non synchronised at face:" << faceI
384 <<
"ownPatch:" << ownPatch[faceI]
385 <<
" syncedOwnPatch:" << syncedOwnPatch[faceI]
386 <<
" neiPatch:" << neiPatch[faceI]
387 <<
" syncedNeiPatch:" << syncedNeiPatch[faceI]
400 if (ownPatch[faceI] != -1)
422 if (map().hasMotionPoints())
438 faceSet baffledFacesSet(mesh_,
"baffledFacesSet", 2*nBaffles);
440 const labelList& reverseFaceMap = map().reverseFaceMap();
444 forAll(ownPatch, oldFaceI)
446 label faceI = reverseFaceMap[oldFaceI];
448 if (ownPatch[oldFaceI] != -1 && faceI >= 0)
454 baffledFacesSet.
insert(ownFaces[i]);
461 label oldFaceI = faceMap[faceI];
463 if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
469 baffledFacesSet.
insert(ownFaces[i]);
473 baffledFacesSet.
sync(mesh_);
491 if (isA<processorPolyPatch>(pp))
501 <<
"face:" << faceI <<
" on patch " << pp.
name()
502 <<
" is in zone " << fZones[zoneI].
name()
526 if (zonedSurfaces.
size())
529 Info<<
"Converting zoned faces into baffles ..." <<
endl;
551 ownPatch[iter.key()] = iter().
first();
552 neiPatch[iter.key()] = iter().second();
566 const labelList& reverseFaceMap = map().reverseFaceMap();
570 label oldFaceI = faceMap[faceI];
578 if (iter != faceToPatch.
end())
580 label masterFaceI = reverseFaceMap[oldFaceI];
581 if (faceI != masterFaceI)
583 baffles[baffleI++] =
labelPair(masterFaceI, faceI);
588 if (baffleI != faceToPatch.
size())
591 <<
"Had " << faceToPatch.
size() <<
" patches to create " 592 <<
" but encountered " << baffleI
593 <<
" slave faces originating from patcheable faces." 600 Pout<<
"Writing zone-baffled mesh to time " <<
timeName()
610 Info<<
"Created " << nZoneFaces <<
" baffles in = " 620 const scalar planarAngle
647 const label baffleValue = 1000000;
670 nBafflesPerEdge[fEdges[fEdgeI]]++;
692 nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
697 label f1 = couples[i].second();
701 nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
737 label edgeI = fEdges[fEdgeI];
739 if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
741 filteredCouples[filterI++] = couple;
747 filteredCouples.
setSize(filterI);
752 Info<<
"freeStandingBaffles : detected " 754 <<
" free-standing baffles out of " 769 forAll(filteredCouples, i)
771 const labelPair& couple = filteredCouples[i];
822 forAll(filteredCouples, i)
824 const labelPair& couple = filteredCouples[i];
831 surface1[i] != surface2[i]
832 || hit1[i].index() != hit2[i].index()
844 if ((normal1[i]&normal2[i]) > planarAngleCos)
848 scalar magN =
mag(n);
851 filteredCouples[filterI++] = couple;
855 else if (hit1[i].hit() || hit2[i].hit())
861 filteredCouples.
setSize(filterI);
863 Info<<
"freeStandingBaffles : detected " 865 <<
" planar (within " << planarAngle
866 <<
" degrees) free-standing baffles out of " 871 return filteredCouples;
890 label face1 = couples[i].second();
894 label own0 = faceOwner[face0];
895 label own1 = faceOwner[face1];
897 if (face1 < 0 || own0 < own1)
901 bool zoneFlip =
false;
905 const faceZone& fZone = faceZones[zoneID];
909 label nei = (face1 < 0 ? -1 : own1);
932 bool zoneFlip =
false;
936 const faceZone& fZone = faceZones[zoneID];
966 if (map().hasMotionPoints())
987 label newFace0 = map().reverseFaceMap()[couples[i].
first()];
990 newExposedFaces[newI++] = newFace0;
993 label newFace1 = map().reverseFaceMap()[couples[i].second()];
996 newExposedFaces[newI++] = newFace1;
999 newExposedFaces.setSize(newI);
1007 void Foam::meshRefinement::findCellZoneGeometric
1025 closedNamedSurfaces,
1030 forAll(insideSurfaces, cellI)
1032 if (cellToZone[cellI] == -2)
1034 label surfI = insideSurfaces[cellI];
1038 cellToZone[cellI] = surfaceToCellZone[surfI];
1051 label nCandidates = 0;
1052 forAll(namedSurfaceIndex, faceI)
1054 label surfI = namedSurfaceIndex[faceI];
1072 forAll(namedSurfaceIndex, faceI)
1074 label surfI = namedSurfaceIndex[faceI];
1078 label own = faceOwner[faceI];
1079 const point& ownCc = cellCentres[own];
1083 label nei = faceNeighbour[faceI];
1084 const point& neiCc = cellCentres[nei];
1086 const vector d = 1
e-4*(neiCc - ownCc);
1087 candidatePoints[nCandidates++] = ownCc-d;
1088 candidatePoints[nCandidates++] = neiCc+d;
1096 const vector d = 1
e-4*(neiFc - ownCc);
1097 candidatePoints[nCandidates++] = ownCc-d;
1107 closedNamedSurfaces,
1116 forAll(namedSurfaceIndex, faceI)
1118 label surfI = namedSurfaceIndex[faceI];
1122 label own = faceOwner[faceI];
1126 label ownSurfI = insideSurfaces[nCandidates++];
1129 cellToZone[own] = surfaceToCellZone[ownSurfI];
1132 label neiSurfI = insideSurfaces[nCandidates++];
1135 label nei = faceNeighbour[faceI];
1137 cellToZone[nei] = surfaceToCellZone[neiSurfI];
1142 label ownSurfI = insideSurfaces[nCandidates++];
1145 cellToZone[own] = surfaceToCellZone[ownSurfI];
1161 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1167 max(ownZone, neiZone)
1203 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1209 max(ownZone, neiZone)
1221 void Foam::meshRefinement::findCellZoneInsideWalk
1234 forAll(namedSurfaceIndex, faceI)
1236 if (namedSurfaceIndex[faceI] == -1)
1238 blockedFace[faceI] =
false;
1242 blockedFace[faceI] =
true;
1249 blockedFace.clear();
1258 forAll(locationSurfaces, i)
1260 label surfI = locationSurfaces[i];
1262 const point& insidePoint = surfZones[surfI].zoneInsidePoint();
1264 Info<<
"For surface " << surfaces_.
names()[surfI]
1265 <<
" finding inside point " << insidePoint
1273 mergeDistance_*
vector(1,1,1),
1277 Info<<
"For surface " << surfaces_.
names()[surfI]
1278 <<
" found point " << insidePoint
1279 <<
" in global region " << keepRegionI
1280 <<
" out of " << cellRegion.
nRegions() <<
" regions." <<
endl;
1282 if (keepRegionI == -1)
1285 <<
"Point " << insidePoint
1286 <<
" is not inside the mesh." <<
nl 1287 <<
"Bounding box of the mesh:" << mesh_.
bounds()
1292 forAll(cellRegion, cellI)
1294 if (cellRegion[cellI] == keepRegionI)
1296 if (cellToZone[cellI] == -2)
1298 cellToZone[cellI] = surfaceToCellZone[surfI];
1300 else if (cellToZone[cellI] != surfaceToCellZone[surfI])
1305 <<
" is inside surface " << surfaces_.
names()[surfI]
1306 <<
" but already marked as being in zone " 1307 << cellToZone[cellI] << endl
1308 <<
"This can happen if your surfaces are not" 1309 <<
" (sufficiently) closed." 1318 bool Foam::meshRefinement::calcRegionToZone
1320 const label surfZoneI,
1321 const label ownRegion,
1322 const label neiRegion,
1327 bool changed =
false;
1330 if (ownRegion != neiRegion)
1337 if (regionToCellZone[ownRegion] == -2)
1339 if (regionToCellZone[neiRegion] == surfZoneI)
1343 regionToCellZone[ownRegion] = -1;
1346 else if (regionToCellZone[neiRegion] != -2)
1350 regionToCellZone[ownRegion] = surfZoneI;
1354 else if (regionToCellZone[neiRegion] == -2)
1356 if (regionToCellZone[ownRegion] == surfZoneI)
1360 regionToCellZone[neiRegion] = -1;
1363 else if (regionToCellZone[ownRegion] != -2)
1367 regionToCellZone[neiRegion] = surfZoneI;
1376 void Foam::meshRefinement::findCellZoneTopo
1378 const point& keepPoint,
1392 forAll(namedSurfaceIndex, faceI)
1394 if (namedSurfaceIndex[faceI] == -1)
1396 blockedFace[faceI] =
false;
1400 blockedFace[faceI] =
true;
1407 blockedFace.clear();
1418 forAll(cellToZone, cellI)
1420 if (cellToZone[cellI] != -2)
1422 regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
1432 mergeDistance_*
vector(1,1,1),
1436 Info<<
"Found point " << keepPoint
1437 <<
" in global region " << keepRegionI
1438 <<
" out of " << cellRegion.
nRegions() <<
" regions." <<
endl;
1440 if (keepRegionI == -1)
1443 <<
"Point " << keepPoint
1444 <<
" is not inside the mesh." <<
nl 1445 <<
"Bounding box of the mesh:" << mesh_.
bounds()
1450 if (regionToCellZone[keepRegionI] == -2)
1452 regionToCellZone[keepRegionI] = -1;
1471 bool changed =
false;
1477 label surfI = namedSurfaceIndex[faceI];
1484 bool changedCell = calcRegionToZone
1486 surfaceToCellZone[surfI],
1492 changed = changed | changedCell;
1530 label surfI = namedSurfaceIndex[faceI];
1535 bool changedCell = calcRegionToZone
1537 surfaceToCellZone[surfI],
1543 changed = changed | changedCell;
1556 forAll(regionToCellZone, regionI)
1558 label zoneI = regionToCellZone[regionI];
1563 <<
"For region " << regionI <<
" haven't set cell zone." 1570 forAll(regionToCellZone, regionI)
1572 Pout<<
"Region " << regionI
1573 <<
" becomes cellZone:" << regionToCellZone[regionI]
1579 forAll(cellToZone, cellI)
1581 cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
1586 void Foam::meshRefinement::makeConsistentFaceIndex
1597 label ownZone = cellToZone[faceOwner[faceI]];
1598 label neiZone = cellToZone[faceNeighbour[faceI]];
1600 if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1602 namedSurfaceIndex[faceI] = -1;
1604 else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1607 <<
"Different cell zones on either side of face " << faceI
1609 <<
" but face not marked with a surface." 1645 label ownZone = cellToZone[faceOwner[faceI]];
1648 if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1650 namedSurfaceIndex[faceI] = -1;
1652 else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1655 <<
"Different cell zones on either side of face " 1657 <<
" but face not marked with a surface." 1668 namedSurfaceIndex[faceI] = -1;
1675 void Foam::meshRefinement::handleSnapProblems
1678 const bool useTopologicalSnapDetection,
1679 const bool removeEdgeConnectedCells,
1688 <<
"Introducing baffles to block off problem cells" <<
nl 1689 <<
"----------------------------------------------" <<
nl 1693 if (useTopologicalSnapDetection)
1695 facePatch = markFacesOnProblemCells
1698 removeEdgeConnectedCells,
1705 facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
1707 Info<<
"Analyzed problem cells in = " 1712 faceSet problemFaces(mesh_,
"problemFaces", mesh_.
nFaces()/100);
1716 if (facePatch[faceI] != -1)
1718 problemFaces.insert(faceI);
1721 problemFaces.instance() =
timeName();
1722 Pout<<
"Dumping " << problemFaces.size()
1723 <<
" problem faces to " << problemFaces.objectPath() <<
endl;
1724 problemFaces.
write();
1727 Info<<
"Introducing baffles to delete problem cells." <<
nl <<
endl;
1742 Info<<
"Created baffles in = " 1749 Pout<<
"Writing extra baffled mesh to time " 1755 runTime.
path()/
"extraBaffles" 1757 Pout<<
"Dumped debug data in = " 1793 if (faceToZone[faceI] != -1)
1796 label ownZone = cellToZone[faceOwner[faceI]];
1797 label neiZone = cellToZone[faceNeighbour[faceI]];
1798 if (ownZone == neiZone)
1800 faceLabels.
append(faceI);
1811 if (faceToZone[faceI] != -1)
1814 label ownZone = cellToZone[faceOwner[faceI]];
1816 if (ownZone == neiZone)
1818 faceLabels.
append(faceI);
1823 return faceLabels.shrink();
1827 void Foam::meshRefinement::calcPatchNumMasterFaces
1836 nMasterFacesPerEdge = 0;
1838 forAll(patch.addressing(), faceI)
1840 const label meshFaceI = patch.addressing()[faceI];
1842 if (isMasterFace[meshFaceI])
1847 nMasterFacesPerEdge[fEdges[fEdgeI]]++;
1856 nMasterFacesPerEdge,
1876 label nProtected = 0;
1878 forAll(nMasterFacesPerEdge, edgeI)
1880 if (nMasterFacesPerEdge[edgeI] > 2)
1882 allEdgeInfo[edgeI] = -2;
1901 >::propagationTol();
1909 label currentZoneI = 0;
1915 for (; faceI < allFaceInfo.size(); faceI++)
1917 if (!allFaceInfo[faceI].valid(dummyTrackData))
1919 globalSeed = globalFaces.toGlobal(faceI);
1931 label procI = globalFaces.whichProcID(globalSeed);
1932 label seedFaceI = globalFaces.toLocal(procI, globalSeed);
1944 faceInfo = currentZoneI;
1950 label edgeI = fEdges[fEdgeI];
1968 changedEdges.
append(edgeI);
1969 changedInfo.
append(edgeInfo);
2001 faceToZone.
setSize(patch.size());
2002 forAll(allFaceInfo, faceI)
2004 if (!allFaceInfo[faceI].valid(dummyTrackData))
2007 <<
"Problem: unvisited face " << faceI
2011 faceToZone[faceI] = allFaceInfo[faceI].region();
2014 return currentZoneI;
2018 void Foam::meshRefinement::consistentOrientation
2038 label nProtected = 0;
2040 forAll(patch.addressing(), faceI)
2042 const label meshFaceI = patch.addressing()[faceI];
2048 && bm[patchI].coupled()
2049 && !isMasterFace[meshFaceI]
2062 label nProtected = 0;
2064 forAll(nMasterFacesPerEdge, edgeI)
2066 if (nMasterFacesPerEdge[edgeI] > 2)
2073 Info<<
"Protected from visiting " 2075 <<
" non-manifold edges" <<
nl <<
endl;
2087 >::propagationTol();
2097 forAll(allFaceInfo, faceI)
2101 globalSeed = globalFaces.toGlobal(faceI);
2113 label procI = globalFaces.whichProcID(globalSeed);
2114 label seedFaceI = globalFaces.toLocal(procI, globalSeed);
2128 if (zoneToOrientation[faceToZone[seedFaceI]] < 0)
2137 label edgeI = fEdges[fEdgeI];
2155 changedEdges.
append(edgeI);
2156 changedInfo.
append(edgeInfo);
2195 forAll(patch.addressing(), i)
2197 const label meshFaceI = patch.addressing()[i];
2201 allFaceInfo[i].flipStatus();
2206 forAll(patch.addressing(), i)
2208 const label meshFaceI = patch.addressing()[i];
2214 && bm[patchI].coupled()
2215 && !isMasterFace[meshFaceI]
2232 <<
"Incorrect status for face " << meshFaceI
2243 meshFlipMap =
false;
2245 forAll(allFaceInfo, faceI)
2247 label meshFaceI = patch.addressing()[faceI];
2251 meshFlipMap[meshFaceI] =
false;
2255 meshFlipMap[meshFaceI] =
true;
2260 <<
"Problem : unvisited face " << faceI
2272 const bool doHandleSnapProblems,
2274 const bool useTopologicalSnapDetection,
2275 const bool removeEdgeConnectedCells,
2277 const bool mergeFreeStanding,
2278 const scalar planarAngle,
2283 const point& keepPoint
2292 Info<<
"Introducing baffles for " 2294 <<
" faces that are intersected by the surface." <<
nl <<
endl;
2299 calcNeighbourData(neiLevel, neiCc);
2304 globalToMasterPatch,
2320 Info<<
"Created baffles in = " 2333 runTime.
path()/
"baffles" 2335 Pout<<
"Dumped debug data in = " 2345 if (doHandleSnapProblems)
2350 useTopologicalSnapDetection,
2351 removeEdgeConnectedCells,
2355 globalToMasterPatch,
2365 <<
"Remove unreachable sections of mesh" <<
nl 2366 <<
"-----------------------------------" <<
nl 2381 Info<<
"Split mesh in = " 2396 Pout<<
"Dumped debug data in = " 2404 if (mergeFreeStanding)
2407 <<
"Merge free-standing baffles" <<
nl 2408 <<
"---------------------------" <<
nl 2425 Info<<
"Detected free-standing baffles : " << nCouples <<
endl;
2439 useTopologicalSnapDetection,
2440 removeEdgeConnectedCells,
2444 globalToMasterPatch,
2454 Info<<
"Merged free-standing baffles in = " 2462 const label nBufferLayers,
2465 const point& keepPoint
2474 calcNeighbourData(neiLevel, neiCc);
2479 globalToMasterPatch,
2492 if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
2494 blockedFace[faceI] =
true;
2501 blockedFace.clear();
2508 mergeDistance_*
vector(1,1,1),
2512 Info<<
"Found point " << keepPoint
2513 <<
" in global region " << keepRegionI
2514 <<
" out of " << cellRegion.
nRegions() <<
" regions." <<
endl;
2516 if (keepRegionI == -1)
2519 <<
"Point " << keepPoint
2520 <<
" is not inside the mesh." <<
nl 2521 <<
"Bounding box of the mesh:" << mesh_.
bounds()
2536 label defaultPatch = 0;
2537 if (globalToMasterPatch.
size())
2539 defaultPatch = globalToMasterPatch[0];
2542 for (
label i = 0; i < nBufferLayers; i++)
2548 forAll(faceNeighbour, faceI)
2552 label ownRegion = cellRegion[faceOwner[faceI]];
2553 label neiRegion = cellRegion[faceNeighbour[faceI]];
2555 if (ownRegion == keepRegionI && neiRegion != keepRegionI)
2562 pointBaffle[f[fp]] =
max(defaultPatch, ownPatch[faceI]);
2565 else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
2567 label newPatchI = neiPatch[faceI];
2568 if (newPatchI == -1)
2570 newPatchI =
max(defaultPatch, ownPatch[faceI]);
2574 pointBaffle[f[fp]] = newPatchI;
2587 label ownRegion = cellRegion[faceOwner[faceI]];
2589 if (ownRegion == keepRegionI)
2593 pointBaffle[f[fp]] =
max(defaultPatch, ownPatch[faceI]);
2612 forAll(pointFaces, pointI)
2614 if (pointBaffle[pointI] != -1)
2620 label faceI = pFaces[pFaceI];
2622 if (ownPatch[faceI] == -1)
2624 ownPatch[faceI] = pointBaffle[pointI];
2638 if (ownPatch[faceI] != -1)
2640 label own = faceOwner[faceI];
2642 if (cellRegion[own] != keepRegionI)
2644 cellRegion[own] = keepRegionI;
2646 const cell& ownFaces = mesh_.
cells()[own];
2649 if (ownPatch[ownFaces[j]] == -1)
2651 newOwnPatch[ownFaces[j]] = ownPatch[faceI];
2657 label nei = faceNeighbour[faceI];
2659 if (cellRegion[nei] != keepRegionI)
2661 cellRegion[nei] = keepRegionI;
2663 const cell& neiFaces = mesh_.
cells()[nei];
2666 if (ownPatch[neiFaces[j]] == -1)
2668 newOwnPatch[neiFaces[j]] = ownPatch[faceI];
2688 forAll(cellRegion, cellI)
2690 if (cellRegion[cellI] != keepRegionI)
2692 cellsToRemove.append(cellI);
2695 cellsToRemove.shrink();
2697 label nCellsToKeep = mesh_.
nCells() - cellsToRemove.size();
2700 Info<<
"Keeping all cells in region " << keepRegionI
2701 <<
" containing point " << keepPoint <<
endl 2702 <<
"Selected for keeping : " << nCellsToKeep
2703 <<
" cells." <<
endl;
2711 labelList exposedPatches(exposedFaces.size());
2715 label faceI = exposedFaces[i];
2717 if (ownPatch[faceI] != -1)
2719 exposedPatches[i] = ownPatch[faceI];
2724 <<
"For exposed face " << faceI
2726 <<
" found no patch." << endl
2727 <<
" Taking patch " << defaultPatch
2728 <<
" instead." <<
endl;
2729 exposedPatches[i] = defaultPatch;
2733 return doRemoveCells
2757 Info<<
"dupNonManifoldPoints : Found : " << nNonManifPoints
2758 <<
" non-manifold points (out of " 2775 if (map().hasMotionPoints())
2807 const point& keepPoint,
2808 const bool allowFreeStandingZoneFaces
2817 label surfI = namedSurfaces[i];
2819 Info<<
"Surface : " << surfaces_.
names()[surfI] <<
nl 2820 <<
" faceZone : " << surfZones[surfI].faceZoneName() <<
nl 2821 <<
" cellZone : " << surfZones[surfI].cellZoneName() <<
endl;
2852 calcNeighbourData(neiLevel, neiCc);
2885 label faceI = testFaces[i];
2889 start[i] = cellCentres[faceOwner[faceI]];
2890 end[i] = cellCentres[faceNeighbour[faceI]];
2894 start[i] = cellCentres[faceOwner[faceI]];
2901 const vectorField smallVec(ROOTSMALL*(end-start));
2941 label faceI = testFaces[i];
2944 if (surface1[i] != -1)
2951 magSqr(hit2[i].hitPoint())
2952 <
magSqr(hit1[i].hitPoint())
2956 namedSurfaceIndex[faceI] = surface2[i];
2957 posOrientation[faceI] = ((area&normal2[i]) > 0);
2958 nSurfFaces[surface2[i]]++;
2962 namedSurfaceIndex[faceI] = surface1[i];
2963 posOrientation[faceI] = ((area&normal1[i]) > 0);
2964 nSurfFaces[surface1[i]]++;
2967 else if (surface2[i] != -1)
2969 namedSurfaceIndex[faceI] = surface2[i];
2970 posOrientation[faceI] = ((area&normal2[i]) > 0);
2971 nSurfFaces[surface2[i]]++;
2990 forAll(nSurfFaces, surfI)
2993 << surfaces_.
names()[surfI]
2994 <<
" nZoneFaces:" << nSurfFaces[surfI] <<
nl;
3025 if (closedNamedSurfaces.size())
3027 Info<<
"Found " << closedNamedSurfaces.size()
3028 <<
" closed, named surfaces. Assigning cells in/outside" 3029 <<
" these surfaces to the corresponding cellZone." 3032 findCellZoneGeometric
3035 closedNamedSurfaces,
3052 if (locationSurfaces.
size())
3054 Info<<
"Found " << locationSurfaces.
size()
3055 <<
" named surfaces with a provided inside point." 3056 <<
" Assigning cells inside these surfaces" 3057 <<
" to the corresponding cellZone." 3060 findCellZoneInsideWalk
3075 Info<<
"Walking from location-in-mesh " << keepPoint
3076 <<
" to assign cellZones " 3077 <<
"- crossing a faceZone face changes cellZone" <<
nl <<
endl;
3092 if (!allowFreeStandingZoneFaces)
3094 Info<<
"Only keeping zone faces inbetween different cellZones." 3097 makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
3107 forAll(namedSurfaceIndex, faceI)
3109 label surfI = namedSurfaceIndex[faceI];
3112 faceToZone[faceI] = surfaceToFaceZone[surfI];
3134 neiCellZone[bFaceI++] = -1;
3164 freeStandingBaffleFaces
3175 if (nFreeStanding > 0)
3177 Info<<
"Detected " << nFreeStanding <<
" free-standing zone faces" 3189 calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);
3194 const label nZones = markPatchZones
3197 nMasterFacesPerEdge,
3202 for (
label zoneI = 0; zoneI < nZones; zoneI++)
3204 nPosOrientation.
insert(zoneI, 0);
3210 consistentOrientation
3214 nMasterFacesPerEdge,
3215 faceToConnectedZone,
3223 forAll(patch.addressing(), faceI)
3225 label meshFaceI = patch.addressing()[faceI];
3227 if (isMasterFace[meshFaceI])
3232 bool(posOrientation[meshFaceI])
3233 == meshFlipMap[meshFaceI]
3239 nPosOrientation.
find(faceToConnectedZone[faceI])() += n;
3246 Info<<
"Split " << nFreeStanding <<
" free-standing zone faces" 3247 <<
" into " << nZones <<
" disconnected regions with size" 3248 <<
" (negative denotes wrong orientation) :" 3251 for (
label zoneI = 0; zoneI < nZones; zoneI++)
3253 Info<<
" " << zoneI <<
"\t" << nPosOrientation[zoneI]
3261 consistentOrientation
3265 nMasterFacesPerEdge,
3266 faceToConnectedZone,
3280 label faceZoneI = faceToZone[faceI];
3282 if (faceZoneI != -1)
3288 label ownZone = cellToZone[faceOwner[faceI]];
3289 label neiZone = cellToZone[faceNeighbour[faceI]];
3293 if (ownZone == neiZone)
3296 flip = meshFlipMap[faceI];
3303 || (neiZone != -1 && ownZone > neiZone)
3311 mesh_.
faces()[faceI],
3314 faceNeighbour[faceI],
3335 label faceZoneI = faceToZone[faceI];
3337 if (faceZoneI != -1)
3339 label ownZone = cellToZone[faceOwner[faceI]];
3344 if (ownZone == neiZone)
3347 flip = meshFlipMap[faceI];
3354 || (neiZone != -1 && ownZone > neiZone)
3362 mesh_.
faces()[faceI],
3382 forAll(cellToZone, cellI)
3384 label zoneI = cellToZone[cellI];
3407 if (map().hasMotionPoints())
const vectorField & faceAreas() const
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
const cellZoneMesh & cellZones() const
Return cell zone mesh.
A List with indirect addressing.
void calc(const argList &args, const Time &runTime, const fvMesh &mesh)
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have 'insidePoint'.
Transport of region for use in PatchEdgeFaceWave.
OFstream which keeps track of vertices.
label countHits() const
Count number of intersections (local)
void findNearestIntersection(const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList ®ion1, labelList &surface2, List< pointIndexHit > &hit2, labelList ®ion2) const
Find intersection nearest to the endpoints. surface1,2 are.
const boolList & flipMap() const
Return face flip map.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
const faceZoneMesh & faceZones() const
Return face zone mesh.
autoPtr< mapPolyMesh > createBaffles(const labelList &ownPatch, const labelList &neiPatch)
Create baffle for every internal face where ownPatch != -1.
static PackedBoolList getMasterFaces(const polyMesh &)
Get per face whether it is uncoupled or a master of a.
bool write() const
Write mesh and all data.
label size() const
Return the number of elements in the PtrList.
const searchableSurfaces & geometry() const
static label findRegion(const polyMesh &, const labelList &cellRegion, const vector &perturbVec, const point &p)
Find region point is in. Uses optional perturbation to re-test.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
void baffleAndSplitMesh(const bool handleSnapProblems, const snapParameters &snapParams, const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField &perpendicularAngle, const bool mergeFreeStanding, const scalar freeStandingAngle, const dictionary &motionDict, Time &runTime, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split off unreachable areas of mesh.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
An ordered pair of two objects of type <T> with first() and second() elements.
void setInstance(const fileName &)
Set the instance for mesh files.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
fileName path() const
Return path.
const word & name() const
Return name.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
void checkZoneFaces() const
Debug helper: check faceZones are not on processor patches.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
dimensioned< scalar > mag(const dimensioned< Type > &)
An STL-conforming const_iterator.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
A cell is defined as a list of faces with extra functionality.
const labelListList & pointEdges() const
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
void flip()
Reverse orientation.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
label nRegions() const
Return total number of regions.
A HashTable to objects of type <T> with a label key.
A subset of mesh faces organised as a primitive patch.
autoPtr< mapPolyMesh > splitMesh(const label nBufferLayers, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split off (with optional buffer layers) unreachable areas.
const Type & second() const
Return second.
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void updateMesh(const mapPolyMesh &, const labelList &changedFaces)
Update for external change to mesh. changedFaces are in new mesh.
void size(const label)
Override size to be inconsistent with allocated storage.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Given list of cells to remove insert all the topology changes.
void setSize(const label)
Dummy setSize function.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
const cellList & cells() const
autoPtr< mapPolyMesh > splitMeshRegions(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split mesh. Keep part containing point.
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
const vectorField & cellCentres() const
void printMeshInfo(const bool, const string &) const
Print some mesh stats.
A list of keyword definitions, which are a keyword followed by any number of values (e...
bool updateEdge(const polyMesh &mesh, const indirectPrimitivePatch &patch, const label edgeI, const label faceI, const patchFaceOrientation &faceInfo, const scalar tol, TrackingData &td)
Influence of face on edge.
Simple container to keep together snap specific information.
static labelList getClosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are closed and.
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
T & first()
Return the first element of the list.
const labelListList & faceEdges() const
vectorField pointField
pointField is a vectorField.
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
Class containing data for face removal.
A patch is a list of labels that address the faces in the global face list.
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
const wordList & names() const
Names of surfaces.
const Time & time() const
Return the top-level database.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
A face is a list of labels corresponding to mesh vertices.
void findInside(const labelList &surfacesToTest, const pointField &pt, labelList &insideSurfaces) const
Detect if a point is 'inside' (closed) surfaces.
static labelList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
static void mapCombineScatter(const List< commsStruct > &comms, Container &Values, const int tag, const label comm)
Scatter data. Reverse of combineGather.
void setRefinement(const localPointRegion ®ionSide, polyTopoChange &)
Play commands into polyTopoChange to duplicate points. Gets.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
Wave propagation of information along patch. Every iteration information goes through one layer of fa...
const double e
Elementary charge.
#define WarningInFunction
Report a warning using Foam::Warning.
const Field< PointType > & faceCentres() const
Return face centres for patch.
void setSize(const label)
Reset size of List.
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual void sync(const polyMesh &mesh)
Sync faceSet across coupled patches.
virtual const pointField & points() const
Return raw points.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
virtual Ostream & write(const char)
Write character.
autoPtr< mapPolyMesh > createZoneBaffles(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, List< labelPair > &)
Create baffles for faces straddling zoned surfaces. Return.
virtual Ostream & write(const token &)=0
Write next token to stream.
dimensionedScalar cos(const dimensionedScalar &ds)
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Unit conversion functions.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Class describing modification of a face.
static labelList addCellZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
label size() const
Return number of elements in table.
List< Key > toc() const
Return the table of contents.
static void mapCombineGather(const List< commsStruct > &comms, Container &Values, const CombineOp &cop, const int tag, const label comm)
const labelListList & pointFaces() const
errorManip< error > abort(error &err)
const word & name() const
Return name.
virtual const labelList & faceOwner() const
Return face owner.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
const word & name() const
Return name.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
virtual const fileName & name() const
Return the name of the stream.
label nTotalPoints() const
Return total number of points in decomposed mesh. Not.
Transport of orientation for use in PatchEdgeFaceWave.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
A list of faces which address into the list of points.
static writeType writeLevel()
Get/set write level.
label start() const
Return start label of this patch in the polyMesh face list.
label nInternalFaces() const
labelList intersectedFaces() const
Get faces with intersection.
const Map< label > & meshPointMap() const
Per point that is to be duplicated the local index.
const PtrList< surfaceZonesInfo > & surfZones() const
const boundBox & bounds() const
Return mesh bounding box.
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
List< label > labelList
A List of labels.
void clearOut()
Clear all geometry and addressing.
autoPtr< mapPolyMesh > dupNonManifoldPoints()
Find boundary points that connect to more than one cell.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const labelListList & faceEdges() const
Return face-edge addressing.
Pair< label > labelPair
Label pair.
const Type & first() const
Return first.
Vector< scalar > vector
A scalar version of the templated Vector.
void reverse(UList< T > &, const label n)
Direct mesh changes based on v1.3 polyTopoChange syntax.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
Class describing modification of a cell.
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
void append(const T &)
Append an element at the end of the list.
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
void checkData()
Debugging: check that all faces still obey start()>end()
Takes mesh with 'baffles' (= boundary faces sharing points). Determines for selected points on bounda...
virtual const faceList & faces() const
Return raw faces.
iterator find(const label &)
Find and return an iterator set at the hashedEntry.
bool updateEdge(const polyMesh &mesh, const indirectPrimitivePatch &patch, const label edgeI, const label faceI, const patchEdgeFaceRegion &faceInfo, const scalar tol, TrackingData &td)
Influence of face on edge.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
const labelList & surfaces() const
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
autoPtr< mapPolyMesh > mergeBaffles(const List< labelPair > &)
Merge baffles. Gets pairs of faces.
label nEdges() const
Return number of edges in patch.
virtual const labelList & faceNeighbour() const
Return face neighbour.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
const globalMeshData & globalData() const
Return parallel info.
autoPtr< mapPolyMesh > zonify(const point &keepPoint, const bool allowFreeStandingZoneFaces)
Put faces/cells into zones according to surface specification.
prefixOSstream Pout(cout,"Pout")
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const vectorField & faceCentres() const
static const label labelMax
bool insert(const Key &key)
Insert a new entry.