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 polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
virtual const fileName & name() const
Return the name of the stream.
Class containing data for face removal.
void clearOut()
Clear all geometry and addressing.
void printMeshInfo(const bool, const string &) const
Print some mesh stats.
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.
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Loop across all elements in list.
fileName path() const
Return path.
void checkZoneFaces() const
Debug helper: check faceZones are not on processor patches.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
const word & name() const
Return name.
const labelList & surfaces() const
An STL-conforming const_iterator.
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.
const faceZoneMesh & faceZones() const
Return face zone mesh.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Class describing modification of a face.
A face is a list of labels corresponding to mesh vertices.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
const double e
Elementary charge.
label countHits() const
Count number of intersections (local)
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const labelListList & faceEdges() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
const labelListList & pointEdges() 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.
Wave propagation of information along patch. Every iteration information goes through one layer of fa...
Unit conversion functions.
label nTotalPoints() const
Return total number of points in decomposed mesh. Not.
void size(const label)
Override size to be inconsistent with allocated storage.
const boolList & flipMap() const
Return face flip map.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Given list of cells to remove insert all the topology changes.
Ostream & endl(Ostream &os)
Add newline and flush stream.
autoPtr< mapPolyMesh > createZoneBaffles(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, List< labelPair > &)
Create baffles for faces straddling zoned surfaces. Return.
autoPtr< mapPolyMesh > createBaffles(const labelList &ownPatch, const labelList &neiPatch)
Create baffle for every internal face where ownPatch != -1.
labelList intersectedFaces() const
Get faces with intersection.
static writeType writeLevel()
Get/set write level.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Transport of orientation for use in PatchEdgeFaceWave.
const Field< PointType > & faceCentres() const
Return face centres for patch.
static labelList addCellZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
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.
const cellList & cells() const
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Vector< scalar > vector
A scalar version of the templated Vector.
autoPtr< mapPolyMesh > splitMesh(const label nBufferLayers, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split off (with optional buffer layers) unreachable areas.
const Time & time() const
Return the top-level database.
virtual Ostream & write(const char)
Write character.
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
bool insert(const Key &key)
Insert a new entry.
T & first()
Return the first element of the list.
Class describing modification of a cell.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
const Map< label > & meshPointMap() const
Per point that is to be duplicated the local index.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
label size() const
Return number of elements in table.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
autoPtr< mapPolyMesh > splitMeshRegions(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split mesh. Keep part containing point.
Takes mesh with 'baffles' (= boundary faces sharing points). Determines for selected points on bounda...
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
virtual const pointField & points() const
Return raw points.
static void mapCombineScatter(const List< commsStruct > &comms, Container &Values, const int tag, const label comm)
Scatter data. Reverse of combineGather.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
autoPtr< mapPolyMesh > mergeBaffles(const List< labelPair > &)
Merge baffles. Gets pairs of faces.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
iterator find(const label &)
Find and return an iterator set at the hashedEntry.
const searchableSurfaces & geometry() const
A list of faces which address into the list of points.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
vectorField pointField
pointField is a vectorField.
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
An ordered pair of two objects of type <T> with first() and second() elements.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
dimensionedScalar cos(const dimensionedScalar &ds)
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Simple container to keep together snap specific information.
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
void append(const T &)
Append an element at the end of the list.
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]
Transport of region for use in PatchEdgeFaceWave.
virtual const labelList & faceOwner() const
Return face owner.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
void flip()
Reverse orientation.
Pair< label > labelPair
Label pair.
const globalMeshData & globalData() const
Return parallel info.
static const label labelMax
List< label > labelList
A List of labels.
const word & name() const
Return name.
static void mapCombineGather(const List< commsStruct > &comms, Container &Values, const CombineOp &cop, const int tag, const label comm)
virtual const faceList & faces() const
Return raw faces.
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have 'insidePoint'.
const PtrList< surfaceZonesInfo > & surfZones() const
const vectorField & cellCentres() const
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const Field< PointType > & localPoints() const
Return pointField of points in patch.
label nEdges() const
Return number of edges in patch.
void reverse(UList< T > &, const label n)
OFstream which keeps track of vertices.
void checkData()
Debugging: check that all faces still obey start()>end()
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void updateMesh(const mapPolyMesh &, const labelList &changedFaces)
Update for external change to mesh. changedFaces are in new mesh.
void setInstance(const fileName &)
Set the instance for mesh files.
const Type & second() const
Return second.
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
label size() const
Return the number of elements in the UPtrList.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
const boundBox & bounds() const
Return mesh bounding box.
autoPtr< mapPolyMesh > dupNonManifoldPoints()
Find boundary points that connect to more than one cell.
const wordList & names() const
Names of surfaces.
autoPtr< mapPolyMesh > zonify(const point &keepPoint, const bool allowFreeStandingZoneFaces)
Put faces/cells into zones according to surface specification.
#define WarningInFunction
Report a warning using Foam::Warning.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
A cell is defined as a list of faces with extra functionality.
prefixOSstream Pout(cout, "Pout")
const vectorField & faceAreas() const
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
static labelList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
void setRefinement(const localPointRegion ®ionSide, polyTopoChange &)
Play commands into polyTopoChange to duplicate points. Gets.
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
label nRegions() const
Return total number of regions.
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 labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
const labelListList & pointFaces() const
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
void setSize(const label)
Dummy setSize function.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
virtual Ostream & write(const token &)=0
Write next token to stream.
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
A subset of mesh faces organised as a primitive patch.
List< Key > toc() const
Return the table of contents.
virtual void sync(const polyMesh &mesh)
Sync faceSet across coupled patches.
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.
bool write() const
Write mesh and all data.
A patch is a list of labels that address the faces in the global face list.
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 Type & first() const
Return first.
A List with indirect addressing.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
A HashTable to objects of type <T> with a label key.
void findInside(const labelList &surfacesToTest, const pointField &pt, labelList &insideSurfaces) const
Detect if a point is 'inside' (closed) surfaces.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.