63 x = (
x ==
y) ?
x : value;
71 void Foam::hexRef8::reorder
92 newElems[newI] = elems[i];
96 elems.transfer(newElems);
100 void Foam::hexRef8::getFaceInfo
110 if (!mesh_.isInternalFace(facei))
112 patchID = mesh_.boundaryMesh().whichPatch(facei);
115 zoneID = mesh_.faceZones().whichZone(facei);
121 const faceZone& fZone = mesh_.faceZones()[zoneID];
123 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
129 Foam::label Foam::hexRef8::addFace
131 polyTopoChange& meshMod,
138 label
patchID, zoneID, zoneFlip;
140 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
144 if ((nei == -1) || (own < nei))
147 newFacei = meshMod.setAction
167 newFacei = meshMod.setAction
171 newFace.reverseFace(),
194 Foam::label Foam::hexRef8::addInternalFace
196 polyTopoChange& meshMod,
197 const label meshFacei,
198 const label meshPointi,
204 if (mesh_.isInternalFace(meshFacei))
206 return meshMod.setAction
234 return meshMod.setAction
290 void Foam::hexRef8::modFace
292 polyTopoChange& meshMod,
299 label
patchID, zoneID, zoneFlip;
301 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
305 (own != mesh_.faceOwner()[facei])
307 mesh_.isInternalFace(facei)
308 && (nei != mesh_.faceNeighbour()[facei])
310 || (newFace != mesh_.faces()[facei])
313 if ((nei == -1) || (own < nei))
337 newFace.reverseFace(),
354 Foam::scalar Foam::hexRef8::getLevel0EdgeLength()
const 356 if (cellLevel_.size() != mesh_.nCells())
359 <<
"Number of cells in mesh:" << mesh_.nCells()
360 <<
" does not equal size of cellLevel:" << cellLevel_.size()
362 <<
"This might be because of a restart with inconsistent cellLevel." 369 const scalar GREAT2 =
sqr(GREAT);
371 label nLevels =
gMax(cellLevel_)+1;
386 const label cLevel = cellLevel_[celli];
388 const labelList& cEdges = mesh_.cellEdges(celli);
392 label edgeI = cEdges[i];
394 if (edgeLevel[edgeI] == -1)
396 edgeLevel[edgeI] = cLevel;
398 else if (edgeLevel[edgeI] ==
labelMax)
402 else if (edgeLevel[edgeI] != cLevel)
416 ifEqEqOp<labelMax>(),
424 const label eLevel = edgeLevel[edgeI];
426 if (eLevel >= 0 && eLevel <
labelMax)
428 const edge&
e = mesh_.edges()[edgeI];
430 scalar edgeLenSqr =
magSqr(
e.vec(mesh_.points()));
432 typEdgeLenSqr[eLevel] =
min(typEdgeLenSqr[eLevel], edgeLenSqr);
443 Pout<<
"hexRef8::getLevel0EdgeLength() :" 444 <<
" After phase1: Edgelengths (squared) per refinementlevel:" 445 << typEdgeLenSqr <<
endl;
459 const label cLevel = cellLevel_[celli];
461 const labelList& cEdges = mesh_.cellEdges(celli);
465 const edge&
e = mesh_.edges()[cEdges[i]];
467 scalar edgeLenSqr =
magSqr(
e.vec(mesh_.points()));
469 maxEdgeLenSqr[cLevel] =
max(maxEdgeLenSqr[cLevel], edgeLenSqr);
477 Pout<<
"hexRef8::getLevel0EdgeLength() :" 478 <<
" Poor Edgelengths (squared) per refinementlevel:" 479 << maxEdgeLenSqr <<
endl;
486 forAll(typEdgeLenSqr, levelI)
488 if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
490 typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
496 Pout<<
"hexRef8::getLevel0EdgeLength() :" 497 <<
" Final Edgelengths (squared) per refinementlevel:" 498 << typEdgeLenSqr <<
endl;
502 scalar level0Size = -1;
504 forAll(typEdgeLenSqr, levelI)
506 scalar lenSqr = typEdgeLenSqr[levelI];
514 Pout<<
"hexRef8::getLevel0EdgeLength() :" 515 <<
" For level:" << levelI
517 <<
" with equivalent level0 len:" << level0Size
524 if (level0Size == -1)
536 Foam::label Foam::hexRef8::getAnchorCell
545 if (cellAnchorPoints[celli].size())
547 label index = cellAnchorPoints[celli].find(pointi);
551 return cellAddedCells[celli][index];
558 const face&
f = mesh_.faces()[facei];
562 label index = cellAnchorPoints[celli].find(
f[fp]);
566 return cellAddedCells[celli][index];
572 Perr<<
"cell:" << celli <<
" anchorPoints:" << cellAnchorPoints[celli]
576 <<
"Could not find point " << pointi
577 <<
" in the anchorPoints for cell " << celli <<
endl 578 <<
"Does your original mesh obey the 2:1 constraint and" 579 <<
" did you use consistentRefinement to make your cells to refine" 580 <<
" obey this constraint as well?" 593 void Foam::hexRef8::getFaceNeighbours
609 mesh_.faceOwner()[facei],
614 if (mesh_.isInternalFace(facei))
620 mesh_.faceNeighbour()[facei],
633 Foam::label Foam::hexRef8::findMinLevel(
const labelList&
f)
const 640 label level = pointLevel_[
f[fp]];
642 if (level < minLevel)
654 Foam::label Foam::hexRef8::findMaxLevel(
const labelList&
f)
const 661 label level = pointLevel_[
f[fp]];
663 if (level > maxLevel)
674 Foam::label Foam::hexRef8::countAnchors
677 const label anchorLevel
684 if (pointLevel_[
f[fp]] <= anchorLevel)
693 void Foam::hexRef8::dumpCell(
const label celli)
const 695 OFstream str(mesh_.time().path()/
"cell_" +
Foam::name(celli) +
".obj");
696 Pout<<
"hexRef8 : Dumping cell as obj to " << str.
name() <<
endl;
698 const cell& cFaces = mesh_.cells()[celli];
700 Map<label> pointToObjVert;
705 const face&
f = mesh_.faces()[cFaces[i]];
709 if (pointToObjVert.insert(
f[fp], objVertI))
719 const face&
f = mesh_.faces()[cFaces[i]];
723 label pointi =
f[fp];
726 str <<
"l " << pointToObjVert[pointi]+1
727 <<
' ' << pointToObjVert[nexPointi]+1 <<
nl;
734 Foam::label Foam::hexRef8::findLevel
739 const bool searchForward,
740 const label wantedLevel
747 label pointi =
f[fp];
749 if (pointLevel_[pointi] < wantedLevel)
751 dumpCell(mesh_.faceOwner()[facei]);
752 if (mesh_.isInternalFace(facei))
754 dumpCell(mesh_.faceNeighbour()[facei]);
760 <<
" startFp:" << startFp
761 <<
" wantedLevel:" << wantedLevel
764 else if (pointLevel_[pointi] == wantedLevel)
779 dumpCell(mesh_.faceOwner()[facei]);
780 if (mesh_.isInternalFace(facei))
782 dumpCell(mesh_.faceNeighbour()[facei]);
788 <<
" startFp:" << startFp
789 <<
" wantedLevel:" << wantedLevel
799 const face&
f = mesh_.faces()[facei];
803 return pointLevel_[
f[findMaxLevel(
f)]];
807 label ownLevel = cellLevel_[mesh_.faceOwner()[facei]];
809 if (countAnchors(
f, ownLevel) == 4)
813 else if (countAnchors(
f, ownLevel+1) == 4)
825 void Foam::hexRef8::checkInternalOrientation
827 polyTopoChange& meshMod,
835 face compactFace(
identity(newFace.size()));
836 pointField compactPoints(meshMod.points(), newFace);
838 const vector areaNorm(compactFace.areaNormal(compactPoints));
840 const vector dir(neiPt - ownPt);
842 if ((dir & areaNorm) < 0)
845 <<
"cell:" << celli <<
" old face:" << facei
846 <<
" newFace:" << newFace <<
endl 847 <<
" coords:" << compactPoints
848 <<
" ownPt:" << ownPt
849 <<
" neiPt:" << neiPt
853 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
855 const scalar
s = (fcToOwn & areaNorm) / (dir & areaNorm);
857 if (s < 0.1 || s > 0.9)
860 <<
"cell:" << celli <<
" old face:" << facei
861 <<
" newFace:" << newFace <<
endl 862 <<
" coords:" << compactPoints
863 <<
" ownPt:" << ownPt
864 <<
" neiPt:" << neiPt
871 void Foam::hexRef8::checkBoundaryOrientation
873 polyTopoChange& meshMod,
877 const point& boundaryPt,
881 face compactFace(
identity(newFace.size()));
882 pointField compactPoints(meshMod.points(), newFace);
884 const vector areaNorm(compactFace.areaNormal(compactPoints));
886 const vector dir(boundaryPt - ownPt);
888 if ((dir & areaNorm) < 0)
891 <<
"cell:" << celli <<
" old face:" << facei
892 <<
" newFace:" << newFace
893 <<
" coords:" << compactPoints
894 <<
" ownPt:" << ownPt
895 <<
" boundaryPt:" << boundaryPt
899 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
901 const scalar
s = (fcToOwn & dir) /
magSqr(dir);
903 if (s < 0.7 || s > 1.3)
906 <<
"cell:" << celli <<
" old face:" << facei
907 <<
" newFace:" << newFace
908 <<
" coords:" << compactPoints
909 <<
" ownPt:" << ownPt
910 <<
" boundaryPt:" << boundaryPt
920 void Foam::hexRef8::insertEdgeSplit
925 DynamicList<label>& verts
928 if (
p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
932 if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
934 verts.append(edgeMidPoint[edgeI]);
948 Foam::label Foam::hexRef8::storeMidPointInfo
956 const bool faceOrder,
957 const label edgeMidPointi,
958 const label anchorPointi,
959 const label faceMidPointi,
961 Map<edge>& midPointToAnchors,
962 Map<edge>& midPointToFaceMids,
963 polyTopoChange& meshMod
968 bool changed =
false;
969 bool haveTwoAnchors =
false;
971 auto edgeMidFnd = midPointToAnchors.find(edgeMidPointi);
973 if (!edgeMidFnd.found())
975 midPointToAnchors.insert(edgeMidPointi, edge(anchorPointi, -1));
979 edge&
e = edgeMidFnd.val();
981 if (anchorPointi !=
e[0])
990 if (
e[0] != -1 &&
e[1] != -1)
992 haveTwoAnchors =
true;
996 bool haveTwoFaceMids =
false;
998 auto faceMidFnd = midPointToFaceMids.find(edgeMidPointi);
1000 if (!faceMidFnd.found())
1002 midPointToFaceMids.insert(edgeMidPointi, edge(faceMidPointi, -1));
1006 edge&
e = faceMidFnd.val();
1008 if (faceMidPointi !=
e[0])
1012 e[1] = faceMidPointi;
1017 if (
e[0] != -1 &&
e[1] != -1)
1019 haveTwoFaceMids =
true;
1026 if (changed && haveTwoAnchors && haveTwoFaceMids)
1028 const edge& anchors = midPointToAnchors[edgeMidPointi];
1029 const edge& faceMids = midPointToFaceMids[edgeMidPointi];
1031 label otherFaceMidPointi = faceMids.otherVertex(faceMidPointi);
1038 DynamicList<label> newFaceVerts(4);
1039 if (faceOrder == (mesh_.faceOwner()[facei] == celli))
1041 newFaceVerts.append(faceMidPointi);
1052 newFaceVerts.append(edgeMidPointi);
1062 newFaceVerts.append(otherFaceMidPointi);
1063 newFaceVerts.append(cellMidPoint[celli]);
1067 newFaceVerts.append(otherFaceMidPointi);
1077 newFaceVerts.append(edgeMidPointi);
1087 newFaceVerts.append(faceMidPointi);
1088 newFaceVerts.append(cellMidPoint[celli]);
1092 newFace.transfer(newFaceVerts);
1094 label anchorCell0 = getAnchorCell
1102 label anchorCell1 = getAnchorCell
1108 anchors.otherVertex(anchorPointi)
1115 if (anchorCell0 < anchorCell1)
1120 ownPt = mesh_.points()[anchorPointi];
1121 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1130 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1131 neiPt = mesh_.points()[anchorPointi];
1138 if (anchorCell0 < anchorCell1)
1140 ownPt = mesh_.points()[anchorPointi];
1141 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1145 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1146 neiPt = mesh_.points()[anchorPointi];
1149 checkInternalOrientation
1160 return addInternalFace
1178 void Foam::hexRef8::createInternalFaces
1188 polyTopoChange& meshMod
1194 const cell& cFaces = mesh_.cells()[celli];
1195 const label cLevel = cellLevel_[celli];
1198 Map<edge> midPointToAnchors(24);
1200 Map<edge> midPointToFaceMids(24);
1203 DynamicList<label> storage;
1207 label nFacesAdded = 0;
1211 label facei = cFaces[i];
1213 const face&
f = mesh_.faces()[facei];
1214 const labelList& fEdges = mesh_.faceEdges(facei, storage);
1219 label faceMidPointi = -1;
1221 label nAnchors = countAnchors(
f, cLevel);
1229 label anchorFp = -1;
1233 if (pointLevel_[
f[fp]] <= cLevel)
1241 label edgeMid = findLevel
1249 label faceMid = findLevel
1258 faceMidPointi =
f[faceMid];
1260 else if (nAnchors == 4)
1265 faceMidPointi = faceMidPoint[facei];
1269 dumpCell(mesh_.faceOwner()[facei]);
1270 if (mesh_.isInternalFace(facei))
1272 dumpCell(mesh_.faceNeighbour()[facei]);
1276 <<
"nAnchors:" << nAnchors
1277 <<
" facei:" << facei
1289 label point0 =
f[fp0];
1291 if (pointLevel_[point0] <= cLevel)
1300 label edgeMidPointi = -1;
1304 if (pointLevel_[
f[fp1]] <= cLevel)
1307 label edgeI = fEdges[fp0];
1309 edgeMidPointi = edgeMidPoint[edgeI];
1311 if (edgeMidPointi == -1)
1315 const labelList& cPoints = mesh_.cellPoints(celli);
1318 <<
"cell:" << celli <<
" cLevel:" << cLevel
1319 <<
" cell points:" << cPoints
1322 <<
" face:" << facei
1326 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1327 <<
" faceMidPoint:" << faceMidPoint[facei]
1328 <<
" faceMidPointi:" << faceMidPointi
1336 label edgeMid = findLevel(facei,
f, fp1,
true, cLevel+1);
1338 edgeMidPointi =
f[edgeMid];
1341 label newFacei = storeMidPointInfo
1364 if (nFacesAdded == 12)
1377 if (pointLevel_[
f[fpMin1]] <= cLevel)
1380 label edgeI = fEdges[fpMin1];
1382 edgeMidPointi = edgeMidPoint[edgeI];
1384 if (edgeMidPointi == -1)
1388 const labelList& cPoints = mesh_.cellPoints(celli);
1391 <<
"cell:" << celli <<
" cLevel:" << cLevel
1392 <<
" cell points:" << cPoints
1395 <<
" face:" << facei
1399 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1400 <<
" faceMidPoint:" << faceMidPoint[facei]
1401 <<
" faceMidPointi:" << faceMidPointi
1409 label edgeMid = findLevel
1418 edgeMidPointi =
f[edgeMid];
1421 newFacei = storeMidPointInfo
1444 if (nFacesAdded == 12)
1452 if (nFacesAdded == 12)
1460 void Foam::hexRef8::walkFaceToMid
1465 const label startFp,
1466 DynamicList<label>& faceVerts
1469 const face&
f = mesh_.faces()[facei];
1470 const labelList& fEdges = mesh_.faceEdges(facei);
1479 if (edgeMidPoint[fEdges[fp]] >= 0)
1481 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1486 if (pointLevel_[
f[fp]] <= cLevel)
1492 else if (pointLevel_[
f[fp]] == cLevel+1)
1495 faceVerts.append(
f[fp]);
1499 else if (pointLevel_[
f[fp]] == cLevel+2)
1502 faceVerts.append(
f[fp]);
1509 void Foam::hexRef8::walkFaceFromMid
1514 const label startFp,
1515 DynamicList<label>& faceVerts
1518 const face&
f = mesh_.faces()[facei];
1519 const labelList& fEdges = mesh_.faceEdges(facei);
1525 if (pointLevel_[
f[fp]] <= cLevel)
1530 else if (pointLevel_[
f[fp]] == cLevel+1)
1533 faceVerts.append(
f[fp]);
1536 else if (pointLevel_[
f[fp]] == cLevel+2)
1546 if (edgeMidPoint[fEdges[fp]] >= 0)
1548 faceVerts.append(edgeMidPoint[fEdges[fp]]);
1557 faceVerts.append(
f[fp]);
1564 Foam::label Foam::hexRef8::faceConsistentRefinement
1574 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1576 label own = mesh_.faceOwner()[facei];
1577 label ownLevel = cellLevel[own] + refineCell.get(own);
1579 label nei = mesh_.faceNeighbour()[facei];
1580 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1582 if (ownLevel > (neiLevel+1))
1586 refineCell.set(nei);
1590 refineCell.unset(own);
1594 else if (neiLevel > (ownLevel+1))
1598 refineCell.set(own);
1602 refineCell.unset(nei);
1611 labelList neiLevel(mesh_.nBoundaryFaces());
1615 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1617 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1626 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1627 label ownLevel = cellLevel[own] + refineCell.get(own);
1629 if (ownLevel > (neiLevel[i]+1))
1633 refineCell.unset(own);
1637 else if (neiLevel[i] > (ownLevel+1))
1641 refineCell.set(own);
1652 void Foam::hexRef8::checkWantedRefinementLevels
1658 bitSet refineCell(mesh_.nCells(), cellsToRefine);
1660 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1662 label own = mesh_.faceOwner()[facei];
1663 label ownLevel = cellLevel[own] + refineCell.get(own);
1665 label nei = mesh_.faceNeighbour()[facei];
1666 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1668 if (
mag(ownLevel-neiLevel) > 1)
1674 <<
" current level:" << cellLevel[own]
1675 <<
" level after refinement:" << ownLevel
1677 <<
"neighbour cell:" << nei
1678 <<
" current level:" << cellLevel[nei]
1679 <<
" level after refinement:" << neiLevel
1681 <<
"which does not satisfy 2:1 constraints anymore." 1688 labelList neiLevel(mesh_.nBoundaryFaces());
1692 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1694 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1703 label facei = i + mesh_.nInternalFaces();
1705 label own = mesh_.faceOwner()[facei];
1706 label ownLevel = cellLevel[own] + refineCell.get(own);
1708 if (
mag(ownLevel - neiLevel[i]) > 1)
1710 label patchi = mesh_.boundaryMesh().whichPatch(facei);
1714 <<
"Celllevel does not satisfy 2:1 constraint." 1715 <<
" On coupled face " 1717 <<
" on patch " << patchi <<
" " 1718 << mesh_.boundaryMesh()[patchi].name()
1719 <<
" owner cell " << own
1720 <<
" current level:" << cellLevel[own]
1721 <<
" level after refinement:" << ownLevel
1723 <<
" (coupled) neighbour cell will get refinement " 1736 Pout<<
"hexRef8::setInstance(const fileName& inst) : " 1737 <<
"Resetting file instance to " << inst <<
endl;
1740 cellLevel_.instance() = inst;
1741 pointLevel_.instance() = inst;
1742 level0Edge_.instance() = inst;
1743 history_.instance() = inst;
1747 void Foam::hexRef8::collectLevelPoints
1751 DynamicList<label>&
points 1756 if (pointLevel_[
f[fp]] <= level)
1764 void Foam::hexRef8::collectLevelPoints
1769 DynamicList<label>&
points 1774 label pointi = meshPoints[
f[fp]];
1775 if (pointLevel_[pointi] <= level)
1784 bool Foam::hexRef8::matchHexShape
1787 const label cellLevel,
1788 DynamicList<face>& quads
1791 const cell& cFaces = mesh_.cells()[celli];
1794 DynamicList<label> verts(4);
1802 label facei = cFaces[i];
1803 const face&
f = mesh_.faces()[facei];
1806 collectLevelPoints(
f, cellLevel, verts);
1807 if (verts.size() == 4)
1809 if (mesh_.faceOwner()[facei] != celli)
1813 quads.append(face(0));
1815 quadVerts.transfer(verts);
1820 if (quads.size() < 6)
1822 Map<labelList> pointFaces(2*cFaces.size());
1826 label facei = cFaces[i];
1827 const face&
f = mesh_.faces()[facei];
1835 collectLevelPoints(
f, cellLevel, verts);
1836 if (verts.size() == 1)
1842 label pointi =
f[fp];
1843 if (pointLevel_[pointi] == cellLevel+1)
1845 auto iter = pointFaces.
find(pointi);
1850 pFaces.appendUniq(facei);
1876 label facei =
pFaces[pFacei];
1877 const face&
f = mesh_.faces()[facei];
1878 if (mesh_.faceOwner()[facei] == celli)
1880 fourFaces[pFacei] =
f;
1884 fourFaces[pFacei] =
f.reverseFace();
1890 SubList<face>(fourFaces),
1895 if (edgeLoops.size() == 1)
1901 bigFace.meshPoints(),
1902 bigFace.edgeLoops()[0],
1907 if (verts.size() == 4)
1909 quads.append(face(0));
1911 quadVerts.transfer(verts);
1918 return (quads.size() == 6);
1925 Foam::hexRef8::hexRef8(
const polyMesh&
mesh,
const bool readHistory)
1933 mesh_.facesInstance(),
1946 mesh_.facesInstance(),
1959 mesh_.facesInstance(),
1972 "refinementHistory",
1973 mesh_.facesInstance(),
1980 (readHistory ? mesh_.nCells() : 0)
1982 faceRemover_(mesh_, GREAT),
1983 savedPointLevel_(0),
1998 <<
"History enabled but number of visible cells " 2001 <<
" is not equal to the number of cells in the mesh " 2008 cellLevel_.
size() != mesh_.nCells()
2009 || pointLevel_.
size() != mesh_.nPoints()
2013 <<
"Restarted from inconsistent cellLevel or pointLevel files." 2017 <<
"Number of cells in mesh:" << mesh_.nCells()
2018 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2019 <<
"Number of points in mesh:" << mesh_.nPoints()
2020 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2038 Foam::hexRef8::hexRef8
2040 const polyMesh&
mesh,
2043 const refinementHistory& history,
2044 const scalar level0Edge
2053 mesh_.facesInstance(),
2054 polyMesh::meshSubDir,
2066 mesh_.facesInstance(),
2067 polyMesh::meshSubDir,
2079 mesh_.facesInstance(),
2080 polyMesh::meshSubDir,
2090 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2097 "refinementHistory",
2098 mesh_.facesInstance(),
2099 polyMesh::meshSubDir,
2106 faceRemover_(mesh_, GREAT),
2107 savedPointLevel_(0),
2113 <<
"History enabled but number of visible cells in it " 2115 <<
" is not equal to the number of cells in the mesh " 2121 cellLevel_.
size() != mesh_.nCells()
2122 || pointLevel_.
size() != mesh_.nPoints()
2126 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2127 <<
"Number of cells in mesh:" << mesh_.nCells()
2128 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2129 <<
"Number of points in mesh:" << mesh_.nPoints()
2130 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2147 Foam::hexRef8::hexRef8
2149 const polyMesh&
mesh,
2152 const scalar level0Edge
2161 mesh_.facesInstance(),
2162 polyMesh::meshSubDir,
2174 mesh_.facesInstance(),
2175 polyMesh::meshSubDir,
2187 mesh_.facesInstance(),
2188 polyMesh::meshSubDir,
2198 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2205 "refinementHistory",
2206 mesh_.facesInstance(),
2207 polyMesh::meshSubDir,
2212 List<refinementHistory::splitCell8>(0),
2216 faceRemover_(mesh_, GREAT),
2217 savedPointLevel_(0),
2222 cellLevel_.
size() != mesh_.nCells()
2223 || pointLevel_.
size() != mesh_.nPoints()
2227 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2228 <<
"Number of cells in mesh:" << mesh_.nCells()
2229 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2230 <<
"Number of points in mesh:" << mesh_.nPoints()
2231 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2265 label nChanged = faceConsistentRefinement
2276 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2277 <<
" refinement levels due to 2:1 conflicts." 2292 checkWantedRefinementLevels(cellLevel, newCellsToRefine);
2295 return newCellsToRefine;
2307 const label maxFaceDiff,
2310 const label maxPointDiff,
2314 const labelList& faceOwner = mesh_.faceOwner();
2315 const labelList& faceNeighbour = mesh_.faceNeighbour();
2318 if (maxFaceDiff <= 0)
2321 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2340 forAll(allCellInfo, celli)
2346 maxFaceDiff*(cellLevel_[celli]+1),
2347 maxFaceDiff*cellLevel_[celli]
2354 label celli = cellsToRefine[i];
2356 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2366 int dummyTrackData = 0;
2374 label facei = facesToCheck[i];
2376 if (allFaceInfo[facei].valid(dummyTrackData))
2380 <<
"Argument facesToCheck seems to have duplicate entries!" 2382 <<
"face:" << facei <<
" occurs at positions " 2390 if (mesh_.isInternalFace(facei))
2395 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2398 label faceRefineCount;
2401 faceCount = neiData.
count() + maxFaceDiff;
2402 faceRefineCount = faceCount + maxFaceDiff;
2406 faceCount = ownData.
count() + maxFaceDiff;
2407 faceRefineCount = faceCount + maxFaceDiff;
2410 seedFaces.append(facei);
2411 seedFacesInfo.append
2419 allFaceInfo[facei] = seedFacesInfo.last();
2423 label faceCount = ownData.
count() + maxFaceDiff;
2424 label faceRefineCount = faceCount + maxFaceDiff;
2426 seedFaces.append(facei);
2427 seedFacesInfo.append
2435 allFaceInfo[facei] = seedFacesInfo.last();
2443 forAll(faceNeighbour, facei)
2446 if (!allFaceInfo[facei].valid(dummyTrackData))
2448 label own = faceOwner[facei];
2449 label nei = faceNeighbour[facei];
2453 if (allCellInfo[own].
count() > allCellInfo[nei].
count())
2455 allFaceInfo[facei].updateFace
2465 seedFacesInfo.append(allFaceInfo[facei]);
2467 else if (allCellInfo[own].
count() < allCellInfo[nei].
count())
2469 allFaceInfo[facei].updateFace
2478 seedFaces.append(facei);
2479 seedFacesInfo.append(allFaceInfo[facei]);
2487 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2490 if (!allFaceInfo[facei].valid(dummyTrackData))
2492 label own = faceOwner[facei];
2505 seedFaces.append(facei);
2506 seedFacesInfo.append(faceData);
2524 Pout<<
"hexRef8::consistentSlowRefinement : Seeded " 2525 << seedFaces.size() <<
" faces between cells with different" 2526 <<
" refinement level." <<
endl;
2530 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2532 seedFacesInfo.clear();
2535 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2543 if (maxPointDiff == -1)
2553 forAll(maxPointCount, pointi)
2555 label& pLevel = maxPointCount[pointi];
2557 const labelList& pCells = mesh_.pointCells(pointi);
2561 pLevel =
max(pLevel, allCellInfo[pCells[i]].
count());
2582 label pointi = pointsToCheck[i];
2587 const labelList& pCells = mesh_.pointCells(pointi);
2591 label celli = pCells[pCelli];
2599 maxPointCount[pointi]
2600 >
cellInfo.count() + maxFaceDiff*maxPointDiff
2608 const cell& cFaces = mesh_.cells()[celli];
2612 label facei = cFaces[cFacei];
2625 if (faceData.
count() > allFaceInfo[facei].count())
2627 changedFacesInfo.insert(facei, faceData);
2634 label nChanged = changedFacesInfo.size();
2644 seedFaces.setCapacity(changedFacesInfo.size());
2645 seedFacesInfo.setCapacity(changedFacesInfo.size());
2649 seedFaces.append(iter.key());
2650 seedFacesInfo.append(iter.val());
2657 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2659 label own = mesh_.faceOwner()[facei];
2662 + (allCellInfo[own].isRefined() ? 1 : 0);
2664 label nei = mesh_.faceNeighbour()[facei];
2667 + (allCellInfo[nei].isRefined() ? 1 : 0);
2669 if (
mag(ownLevel-neiLevel) > 1)
2675 <<
" current level:" << cellLevel_[own]
2676 <<
" current refData:" << allCellInfo[own]
2677 <<
" level after refinement:" << ownLevel
2679 <<
"neighbour cell:" << nei
2680 <<
" current level:" << cellLevel_[nei]
2681 <<
" current refData:" << allCellInfo[nei]
2682 <<
" level after refinement:" << neiLevel
2684 <<
"which does not satisfy 2:1 constraints anymore." <<
nl 2685 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2694 labelList neiLevel(mesh_.nBoundaryFaces());
2695 labelList neiCount(mesh_.nBoundaryFaces());
2696 labelList neiRefCount(mesh_.nBoundaryFaces());
2700 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2701 neiLevel[i] = cellLevel_[own];
2702 neiCount[i] = allCellInfo[own].count();
2703 neiRefCount[i] = allCellInfo[own].refinementCount();
2714 label facei = i+mesh_.nInternalFaces();
2716 label own = mesh_.faceOwner()[facei];
2719 + (allCellInfo[own].isRefined() ? 1 : 0);
2723 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2725 if (
mag(ownLevel - nbrLevel) > 1)
2728 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2731 <<
"Celllevel does not satisfy 2:1 constraint." 2732 <<
" On coupled face " 2734 <<
" refData:" << allFaceInfo[facei]
2735 <<
" on patch " << patchi <<
" " 2736 << mesh_.boundaryMesh()[patchi].name() <<
nl 2737 <<
"owner cell " << own
2738 <<
" current level:" << cellLevel_[own]
2739 <<
" current count:" << allCellInfo[own].count()
2740 <<
" current refCount:" 2741 << allCellInfo[own].refinementCount()
2742 <<
" level after refinement:" << ownLevel
2744 <<
"(coupled) neighbour cell" 2745 <<
" has current level:" << neiLevel[i]
2746 <<
" current count:" << neiCount[i]
2747 <<
" current refCount:" << neiRefCount[i]
2748 <<
" level after refinement:" << nbrLevel
2758 forAll(allCellInfo, celli)
2760 if (allCellInfo[celli].isRefined())
2770 forAll(allCellInfo, celli)
2772 if (allCellInfo[celli].isRefined())
2774 newCellsToRefine[nRefined++] = celli;
2780 Pout<<
"hexRef8::consistentSlowRefinement : From " 2781 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
2782 <<
" cells to refine." <<
endl;
2785 return newCellsToRefine;
2791 const label maxFaceDiff,
2796 const labelList& faceOwner = mesh_.faceOwner();
2797 const labelList& faceNeighbour = mesh_.faceNeighbour();
2799 if (maxFaceDiff <= 0)
2802 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2806 const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
2815 List<refinementDistanceData> allCellInfo(mesh_.nCells());
2818 List<refinementDistanceData> allFaceInfo(mesh_.nFaces());
2821 int dummyTrackData = 0;
2827 label celli = cellsToRefine[i];
2829 allCellInfo[celli] = refinementDistanceData
2832 mesh_.cellCentres()[celli],
2837 forAll(allCellInfo, celli)
2839 if (!allCellInfo[celli].valid(dummyTrackData))
2841 allCellInfo[celli] = refinementDistanceData
2844 mesh_.cellCentres()[celli],
2852 DynamicList<label> seedFaces(mesh_.nFaces()/100);
2854 DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100);
2860 label facei = facesToCheck[i];
2862 if (allFaceInfo[facei].valid(dummyTrackData))
2866 <<
"Argument facesToCheck seems to have duplicate entries!" 2868 <<
"face:" << facei <<
" occurs at positions " 2873 label own = faceOwner[facei];
2877 allCellInfo[own].valid(dummyTrackData)
2878 ? allCellInfo[own].originLevel()
2882 if (!mesh_.isInternalFace(facei))
2886 const point& fc = mesh_.faceCentres()[facei];
2888 refinementDistanceData neiData
2895 allFaceInfo[facei].updateFace
2907 label nei = faceNeighbour[facei];
2911 allCellInfo[nei].valid(dummyTrackData)
2912 ? allCellInfo[nei].originLevel()
2916 if (ownLevel == neiLevel)
2919 allFaceInfo[facei].updateFace
2924 refinementDistanceData(level0Size, cc[nei], neiLevel+1),
2928 allFaceInfo[facei].updateFace
2933 refinementDistanceData(level0Size, cc[own], ownLevel+1),
2941 allFaceInfo[facei].updateFace
2946 refinementDistanceData(level0Size, cc[nei], neiLevel),
2950 allFaceInfo[facei].updateFace
2955 refinementDistanceData(level0Size, cc[own], ownLevel),
2962 seedFacesInfo.append(allFaceInfo[facei]);
2969 forAll(faceNeighbour, facei)
2972 if (!allFaceInfo[facei].valid(dummyTrackData))
2974 label own = faceOwner[facei];
2978 allCellInfo[own].valid(dummyTrackData)
2979 ? allCellInfo[own].originLevel()
2983 label nei = faceNeighbour[facei];
2987 allCellInfo[nei].valid(dummyTrackData)
2988 ? allCellInfo[nei].originLevel()
2992 if (ownLevel > neiLevel)
2996 allFaceInfo[facei].updateFace
3001 refinementDistanceData(level0Size, cc[own], ownLevel),
3005 seedFacesInfo.append(allFaceInfo[facei]);
3007 else if (neiLevel > ownLevel)
3009 seedFaces.append(facei);
3010 allFaceInfo[facei].updateFace
3015 refinementDistanceData(level0Size, cc[nei], neiLevel),
3019 seedFacesInfo.append(allFaceInfo[facei]);
3025 seedFacesInfo.shrink();
3028 FaceCellWave<refinementDistanceData, int> levelCalc
3035 mesh_.globalData().nTotalCells()+1,
3076 label celli = cellsToRefine[i];
3078 allCellInfo[celli].originLevel() =
sizeof(label)*8-2;
3079 allCellInfo[celli].origin() = cc[celli];
3085 bitSet refineCell(mesh_.nCells());
3086 forAll(allCellInfo, celli)
3088 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3090 if (wanted > cellLevel_[celli]+1)
3092 refineCell.set(celli);
3095 faceConsistentRefinement(
true, cellLevel_, refineCell);
3099 label nChanged = faceConsistentRefinement(
true, cellLevel_, refineCell);
3101 reduce(nChanged, sumOp<label>());
3105 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3106 <<
" refinement levels due to 2:1 conflicts." 3117 labelList newCellsToRefine(refineCell.toc());
3121 Pout<<
"hexRef8::consistentSlowRefinement2 : From " 3122 << cellsToRefine.
size() <<
" to " << newCellsToRefine.size()
3123 <<
" cells to refine." <<
endl;
3128 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3129 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3130 << cellsIn.size() <<
" to cellSet " 3131 << cellsIn.objectPath() <<
endl;
3135 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3136 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3137 << cellsOut.size() <<
" to cellSet " 3138 << cellsOut.objectPath() <<
endl;
3143 bitSet refineCell(mesh_.nCells(), newCellsToRefine);
3145 const bitSet savedRefineCell(refineCell);
3146 label nChanged = faceConsistentRefinement(
true, cellLevel_, refineCell);
3151 mesh_,
"cellsToRefineOut2", newCellsToRefine.size()
3153 forAll(refineCell, celli)
3155 if (refineCell.test(celli))
3157 cellsOut2.insert(celli);
3160 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3161 << cellsOut2.size() <<
" to cellSet " 3162 << cellsOut2.objectPath() <<
endl;
3168 forAll(refineCell, celli)
3170 if (refineCell.test(celli) && !savedRefineCell.test(celli))
3174 <<
"Cell:" << celli <<
" cc:" 3175 << mesh_.cellCentres()[celli]
3176 <<
" was not marked for refinement but does not obey" 3177 <<
" 2:1 constraints." 3184 return newCellsToRefine;
3197 Pout<<
"hexRef8::setRefinement :" 3198 <<
" Checking initial mesh just to make sure" <<
endl;
3207 savedPointLevel_.clear();
3208 savedCellLevel_.clear();
3212 DynamicList<label> newCellLevel(cellLevel_.size());
3213 forAll(cellLevel_, celli)
3215 newCellLevel.append(cellLevel_[celli]);
3217 DynamicList<label> newPointLevel(pointLevel_.size());
3218 forAll(pointLevel_, pointi)
3220 newPointLevel.append(pointLevel_[pointi]);
3226 Pout<<
"hexRef8::setRefinement :" 3227 <<
" Allocating " << cellLabels.
size() <<
" cell midpoints." 3235 labelList cellMidPoint(mesh_.nCells(), -1);
3239 label celli = cellLabels[i];
3241 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3247 mesh_.cellCentres()[celli],
3254 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3260 cellSet splitCells(mesh_,
"splitCells", cellLabels.
size());
3262 forAll(cellMidPoint, celli)
3264 if (cellMidPoint[celli] >= 0)
3266 splitCells.insert(celli);
3270 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3271 <<
" cells to split to cellSet " << splitCells.objectPath()
3284 Pout<<
"hexRef8::setRefinement :" 3285 <<
" Allocating edge midpoints." 3294 labelList edgeMidPoint(mesh_.nEdges(), -1);
3297 forAll(cellMidPoint, celli)
3299 if (cellMidPoint[celli] >= 0)
3301 const labelList& cEdges = mesh_.cellEdges(celli);
3305 label edgeI = cEdges[i];
3307 const edge&
e = mesh_.edges()[edgeI];
3311 pointLevel_[
e[0]] <= cellLevel_[celli]
3312 && pointLevel_[
e[1]] <= cellLevel_[celli]
3315 edgeMidPoint[edgeI] = 12345;
3342 forAll(edgeMidPoint, edgeI)
3344 if (edgeMidPoint[edgeI] >= 0)
3347 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3355 point(-GREAT, -GREAT, -GREAT)
3360 forAll(edgeMidPoint, edgeI)
3362 if (edgeMidPoint[edgeI] >= 0)
3367 const edge&
e = mesh_.edges()[edgeI];
3380 newPointLevel(edgeMidPoint[edgeI]) =
3393 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3395 forAll(edgeMidPoint, edgeI)
3397 if (edgeMidPoint[edgeI] >= 0)
3399 const edge&
e = mesh_.edges()[edgeI];
3405 Pout<<
"hexRef8::setRefinement :" 3406 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3416 Pout<<
"hexRef8::setRefinement :" 3417 <<
" Allocating face midpoints." 3423 labelList faceAnchorLevel(mesh_.nFaces());
3425 for (label facei = 0; facei < mesh_.nFaces(); facei++)
3427 faceAnchorLevel[facei] = faceLevel(facei);
3432 labelList faceMidPoint(mesh_.nFaces(), -1);
3437 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3439 if (faceAnchorLevel[facei] >= 0)
3441 label own = mesh_.faceOwner()[facei];
3442 label ownLevel = cellLevel_[own];
3443 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3445 label nei = mesh_.faceNeighbour()[facei];
3446 label neiLevel = cellLevel_[nei];
3447 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3451 newOwnLevel > faceAnchorLevel[facei]
3452 || newNeiLevel > faceAnchorLevel[facei]
3455 faceMidPoint[facei] = 12345;
3468 labelList newNeiLevel(mesh_.nBoundaryFaces());
3472 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3473 label ownLevel = cellLevel_[own];
3474 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3476 newNeiLevel[i] = newOwnLevel;
3486 label facei = i+mesh_.nInternalFaces();
3488 if (faceAnchorLevel[facei] >= 0)
3490 label own = mesh_.faceOwner()[facei];
3491 label ownLevel = cellLevel_[own];
3492 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3496 newOwnLevel > faceAnchorLevel[facei]
3497 || newNeiLevel[i] > faceAnchorLevel[facei]
3500 faceMidPoint[facei] = 12345;
3525 mesh_.nBoundaryFaces(),
3526 point(-GREAT, -GREAT, -GREAT)
3531 label facei = i+mesh_.nInternalFaces();
3533 if (faceMidPoint[facei] >= 0)
3535 bFaceMids[i] = mesh_.faceCentres()[facei];
3545 forAll(faceMidPoint, facei)
3547 if (faceMidPoint[facei] >= 0)
3552 const face&
f = mesh_.faces()[facei];
3559 facei < mesh_.nInternalFaces()
3560 ? mesh_.faceCentres()[facei]
3561 : bFaceMids[facei-mesh_.nInternalFaces()]
3571 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3578 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.
size());
3580 forAll(faceMidPoint, facei)
3582 if (faceMidPoint[facei] >= 0)
3584 splitFaces.insert(facei);
3588 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3589 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3610 Pout<<
"hexRef8::setRefinement :" 3611 <<
" Finding cell anchorPoints (8 per cell)" 3624 forAll(cellMidPoint, celli)
3626 if (cellMidPoint[celli] >= 0)
3628 cellAnchorPoints[celli].setSize(8);
3632 forAll(pointLevel_, pointi)
3634 const labelList& pCells = mesh_.pointCells(pointi);
3638 label celli = pCells[pCelli];
3642 cellMidPoint[celli] >= 0
3643 && pointLevel_[pointi] <= cellLevel_[celli]
3646 if (nAnchorPoints[celli] == 8)
3651 <<
" of level " << cellLevel_[celli]
3652 <<
" uses more than 8 points of equal or" 3653 <<
" lower level" <<
nl 3654 <<
"Points so far:" << cellAnchorPoints[celli]
3658 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3663 forAll(cellMidPoint, celli)
3665 if (cellMidPoint[celli] >= 0)
3667 if (nAnchorPoints[celli] != 8)
3671 const labelList& cPoints = mesh_.cellPoints(celli);
3675 <<
" of level " << cellLevel_[celli]
3676 <<
" does not seem to have 8 points of equal or" 3677 <<
" lower level" <<
endl 3678 <<
"cellPoints:" << cPoints <<
endl 3693 Pout<<
"hexRef8::setRefinement :" 3694 <<
" Adding cells (1 per anchorPoint)" 3701 forAll(cellAnchorPoints, celli)
3703 const labelList& cAnchors = cellAnchorPoints[celli];
3705 if (cAnchors.size() == 8)
3707 labelList& cAdded = cellAddedCells[celli];
3714 newCellLevel[celli] = cellLevel_[celli]+1;
3717 for (label i = 1; i < 8; i++)
3727 mesh_.cellZones().whichZone(celli)
3731 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3746 Pout<<
"hexRef8::setRefinement :" 3747 <<
" Marking faces to be handled" 3752 bitSet affectedFace(mesh_.nFaces());
3755 forAll(cellMidPoint, celli)
3757 if (cellMidPoint[celli] >= 0)
3759 const cell& cFaces = mesh_.cells()[celli];
3761 affectedFace.set(cFaces);
3765 forAll(faceMidPoint, facei)
3767 if (faceMidPoint[facei] >= 0)
3769 affectedFace.set(facei);
3773 forAll(edgeMidPoint, edgeI)
3775 if (edgeMidPoint[edgeI] >= 0)
3777 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3779 affectedFace.
set(eFaces);
3790 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3793 forAll(faceMidPoint, facei)
3795 if (faceMidPoint[facei] >= 0 && affectedFace.test(facei))
3801 const face&
f = mesh_.faces()[facei];
3805 bool modifiedFace =
false;
3807 label anchorLevel = faceAnchorLevel[facei];
3813 label pointi =
f[fp];
3815 if (pointLevel_[pointi] <= anchorLevel)
3819 DynamicList<label> faceVerts(4);
3821 faceVerts.append(pointi);
3837 faceVerts.append(faceMidPoint[facei]);
3849 newFace.transfer(faceVerts);
3870 if (mesh_.isInternalFace(facei))
3872 label oldOwn = mesh_.faceOwner()[facei];
3873 label oldNei = mesh_.faceNeighbour()[facei];
3875 checkInternalOrientation
3880 mesh_.cellCentres()[oldOwn],
3881 mesh_.cellCentres()[oldNei],
3887 label oldOwn = mesh_.faceOwner()[facei];
3889 checkBoundaryOrientation
3894 mesh_.cellCentres()[oldOwn],
3895 mesh_.faceCentres()[facei],
3904 modifiedFace =
true;
3906 modFace(meshMod, facei, newFace, own, nei);
3910 addFace(meshMod, facei, newFace, own, nei);
3916 affectedFace.unset(facei);
3926 Pout<<
"hexRef8::setRefinement :" 3927 <<
" Adding edge splits to unsplit faces" 3931 DynamicList<label> eFacesStorage;
3932 DynamicList<label> fEdgesStorage;
3934 forAll(edgeMidPoint, edgeI)
3936 if (edgeMidPoint[edgeI] >= 0)
3940 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3944 label facei = eFaces[i];
3946 if (faceMidPoint[facei] < 0 && affectedFace.test(facei))
3950 const face&
f = mesh_.faces()[facei];
3951 const labelList& fEdges = mesh_.faceEdges
3957 DynamicList<label> newFaceVerts(
f.
size());
3961 newFaceVerts.append(
f[fp]);
3963 label edgeI = fEdges[fp];
3965 if (edgeMidPoint[edgeI] >= 0)
3967 newFaceVerts.append(edgeMidPoint[edgeI]);
3972 newFace.transfer(newFaceVerts);
3976 label anchorFp = findMinLevel(
f);
3993 if (mesh_.isInternalFace(facei))
3995 label oldOwn = mesh_.faceOwner()[facei];
3996 label oldNei = mesh_.faceNeighbour()[facei];
3998 checkInternalOrientation
4003 mesh_.cellCentres()[oldOwn],
4004 mesh_.cellCentres()[oldNei],
4010 label oldOwn = mesh_.faceOwner()[facei];
4012 checkBoundaryOrientation
4017 mesh_.cellCentres()[oldOwn],
4018 mesh_.faceCentres()[facei],
4024 modFace(meshMod, facei, newFace, own, nei);
4027 affectedFace.unset(facei);
4039 Pout<<
"hexRef8::setRefinement :" 4040 <<
" Changing owner/neighbour for otherwise unaffected faces" 4044 forAll(affectedFace, facei)
4046 if (affectedFace.test(facei))
4048 const face&
f = mesh_.faces()[facei];
4052 label anchorFp = findMinLevel(
f);
4066 modFace(meshMod, facei,
f, own, nei);
4069 affectedFace.unset(facei);
4084 Pout<<
"hexRef8::setRefinement :" 4085 <<
" Create new internal faces for split cells" 4089 forAll(cellMidPoint, celli)
4091 if (cellMidPoint[celli] >= 0)
4116 forAll(cellMidPoint, celli)
4118 if (cellMidPoint[celli] >= 0)
4120 minPointi =
min(minPointi, cellMidPoint[celli]);
4121 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4124 forAll(faceMidPoint, facei)
4126 if (faceMidPoint[facei] >= 0)
4128 minPointi =
min(minPointi, faceMidPoint[facei]);
4129 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4132 forAll(edgeMidPoint, edgeI)
4134 if (edgeMidPoint[edgeI] >= 0)
4136 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4137 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4141 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4144 <<
"Added point labels not consecutive to existing mesh points." 4146 <<
"mesh_.nPoints():" << mesh_.nPoints()
4147 <<
" minPointi:" << minPointi
4148 <<
" maxPointi:" << maxPointi
4153 pointLevel_.transfer(newPointLevel);
4154 cellLevel_.transfer(newCellLevel);
4157 setInstance(mesh_.facesInstance());
4164 if (history_.active())
4168 Pout<<
"hexRef8::setRefinement :" 4169 <<
" Updating refinement history to " << cellLevel_.size()
4170 <<
" cells" <<
endl;
4174 history_.resize(cellLevel_.size());
4176 forAll(cellAddedCells, celli)
4178 const labelList& addedCells = cellAddedCells[celli];
4180 if (addedCells.size())
4183 history_.storeSplit(celli, addedCells);
4194 label celli = cellLabels[i];
4196 refinedCells[i].
transfer(cellAddedCells[celli]);
4199 return refinedCells;
4210 savedPointLevel_.clear();
4211 savedPointLevel_.resize(2*pointsToStore.
size());
4214 label pointi = pointsToStore[i];
4215 savedPointLevel_.insert(pointi, pointLevel_[pointi]);
4218 savedCellLevel_.
clear();
4219 savedCellLevel_.resize(2*cellsToStore.
size());
4222 label celli = cellsToStore[i];
4223 savedCellLevel_.insert(celli, cellLevel_[celli]);
4235 updateMesh(map, dummyMap, dummyMap, dummyMap);
4253 Pout<<
"hexRef8::updateMesh :" 4254 <<
" Updating various lists" 4263 Pout<<
"hexRef8::updateMesh :" 4266 <<
" nCells:" << mesh_.nCells()
4268 <<
" cellLevel_:" << cellLevel_.size()
4271 <<
" nPoints:" << mesh_.nPoints()
4273 <<
" pointLevel_:" << pointLevel_.size()
4277 if (reverseCellMap.
size() == cellLevel_.size())
4283 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4291 forAll(cellMap, newCelli)
4293 label oldCelli = cellMap[newCelli];
4297 newCellLevel[newCelli] = -1;
4301 newCellLevel[newCelli] = cellLevel_[oldCelli];
4311 const label newCelli = iter.key();
4312 const label storedCelli = iter.val();
4314 const auto fnd = savedCellLevel_.cfind(storedCelli);
4319 <<
"Problem : trying to restore old value for new cell " 4320 << newCelli <<
" but cannot find old cell " << storedCelli
4321 <<
" in map of stored values " << savedCellLevel_
4324 cellLevel_[newCelli] = fnd.val();
4345 if (reversePointMap.
size() == pointLevel_.size())
4348 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4359 const label oldPointi = pointMap[
newPointi];
4361 if (oldPointi == -1)
4374 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4377 pointLevel_.
transfer(newPointLevel);
4385 const label storedPointi = iter.val();
4387 auto fnd = savedPointLevel_.find(storedPointi);
4392 <<
"Problem : trying to restore old value for new point " 4393 <<
newPointi <<
" but cannot find old point " 4395 <<
" in map of stored values " << savedPointLevel_
4415 if (history_.active())
4417 history_.updateMesh(map);
4421 setInstance(mesh_.facesInstance());
4424 faceRemover_.updateMesh(map);
4427 cellShapesPtr_.clear();
4442 Pout<<
"hexRef8::subset :" 4443 <<
" Updating various lists" 4447 if (history_.active())
4450 <<
"Subsetting will not work in combination with unrefinement." 4452 <<
"Proceed at your own risk." <<
endl;
4460 forAll(cellMap, newCelli)
4462 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4465 cellLevel_.transfer(newCellLevel);
4467 if (cellLevel_.found(-1))
4471 <<
"cellLevel_ contains illegal value -1 after mapping:" 4486 pointLevel_.transfer(newPointLevel);
4488 if (pointLevel_.found(-1))
4492 <<
"pointLevel_ contains illegal value -1 after mapping:" 4499 if (history_.active())
4501 history_.subset(pointMap,
faceMap, cellMap);
4505 setInstance(mesh_.facesInstance());
4511 cellShapesPtr_.clear();
4520 Pout<<
"hexRef8::distribute :" 4521 <<
" Updating various lists" 4531 if (history_.active())
4533 history_.distribute(map);
4537 faceRemover_.distribute(map);
4540 cellShapesPtr_.clear();
4546 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4550 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:" 4551 << smallDim <<
endl;
4564 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4570 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
4574 const polyPatch& pp =
patches[patchi];
4582 label facei = pp.start();
4586 label own = mesh_.faceOwner()[facei];
4587 label bFacei = facei-mesh_.nInternalFaces();
4589 if (!cellToFace.insert(
labelPair(own, nei[bFacei]), facei))
4593 <<
"Faces do not seem to be correct across coupled" 4594 <<
" boundaries" <<
endl 4595 <<
"Coupled face " << facei
4596 <<
" between owner " << own
4597 <<
" on patch " << pp.
name()
4598 <<
" and coupled neighbour " << nei[bFacei]
4599 <<
" has two faces connected to it:" 4601 << cellToFace[
labelPair(own, nei[bFacei])]
4618 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4626 label facei = i+mesh_.nInternalFaces();
4628 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4630 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4632 const face&
f = mesh_.faces()[facei];
4633 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4635 dumpCell(mesh_.faceOwner()[facei]);
4638 <<
"Faces do not seem to be correct across coupled" 4639 <<
" boundaries" <<
endl 4640 <<
"Coupled face " << facei
4641 <<
" on patch " << patchi
4642 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4643 <<
" coords:" << UIndirectList<point>(mesh_.points(),
f)
4644 <<
" has face area:" << magArea
4645 <<
" (coupled) neighbour face area differs:" 4647 <<
" to within tolerance " << smallDim
4657 labelList nVerts(mesh_.nBoundaryFaces());
4661 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4669 label facei = i+mesh_.nInternalFaces();
4671 const face&
f = mesh_.faces()[facei];
4673 if (
f.
size() != nVerts[i])
4675 dumpCell(mesh_.faceOwner()[facei]);
4677 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4680 <<
"Faces do not seem to be correct across coupled" 4681 <<
" boundaries" <<
endl 4682 <<
"Coupled face " << facei
4683 <<
" on patch " << patchi
4684 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4685 <<
" coords:" << UIndirectList<point>(mesh_.points(),
f)
4686 <<
" has size:" <<
f.
size()
4687 <<
" (coupled) neighbour face has size:" 4699 pointField anchorPoints(mesh_.nBoundaryFaces());
4703 label facei = i+mesh_.nInternalFaces();
4704 const point& fc = mesh_.faceCentres()[facei];
4705 const face&
f = mesh_.faces()[facei];
4706 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4708 anchorPoints[i] = anchorVec;
4718 label facei = i+mesh_.nInternalFaces();
4719 const point& fc = mesh_.faceCentres()[facei];
4720 const face&
f = mesh_.faces()[facei];
4721 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4723 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4725 dumpCell(mesh_.faceOwner()[facei]);
4727 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4730 <<
"Faces do not seem to be correct across coupled" 4731 <<
" boundaries" <<
endl 4732 <<
"Coupled face " << facei
4733 <<
" on patch " << patchi
4734 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4735 <<
" coords:" << UIndirectList<point>(mesh_.points(),
f)
4736 <<
" has anchor vector:" << anchorVec
4737 <<
" (coupled) neighbour face anchor vector differs:" 4739 <<
" to within tolerance " << smallDim
4747 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4754 const label maxPointDiff,
4760 Pout<<
"hexRef8::checkRefinementLevels :" 4761 <<
" Checking 2:1 refinement level" <<
endl;
4766 cellLevel_.size() != mesh_.nCells()
4767 || pointLevel_.size() != mesh_.nPoints()
4771 <<
"cellLevel size should be number of cells" 4772 <<
" and pointLevel size should be number of points."<<
nl 4773 <<
"cellLevel:" << cellLevel_.size()
4774 <<
" mesh.nCells():" << mesh_.nCells() <<
nl 4775 <<
"pointLevel:" << pointLevel_.size()
4776 <<
" mesh.nPoints():" << mesh_.nPoints()
4786 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4788 label own = mesh_.faceOwner()[facei];
4789 label nei = mesh_.faceNeighbour()[facei];
4791 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4797 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl 4798 <<
"On face " << facei <<
" owner cell " << own
4799 <<
" has refinement " << cellLevel_[own]
4800 <<
" neighbour cell " << nei <<
" has refinement " 4807 labelList neiLevel(mesh_.nBoundaryFaces());
4811 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4813 neiLevel[i] = cellLevel_[own];
4821 label facei = i+mesh_.nInternalFaces();
4823 label own = mesh_.faceOwner()[facei];
4825 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4829 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4832 <<
"Celllevel does not satisfy 2:1 constraint." 4833 <<
" On coupled face " << facei
4834 <<
" on patch " << patchi <<
" " 4835 << mesh_.boundaryMesh()[patchi].name()
4836 <<
" owner cell " << own <<
" has refinement " 4838 <<
" (coupled) neighbour cell has refinement " 4861 forAll(syncPointLevel, pointi)
4863 if (pointLevel_[pointi] != syncPointLevel[pointi])
4866 <<
"PointLevel is not consistent across coupled patches." 4868 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4869 <<
" has level " << pointLevel_[pointi]
4870 <<
" whereas the coupled point has level " 4871 << syncPointLevel[pointi]
4879 if (maxPointDiff != -1)
4884 forAll(maxPointLevel, pointi)
4886 const labelList& pCells = mesh_.pointCells(pointi);
4888 label& pLevel = maxPointLevel[pointi];
4892 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4908 label pointi = pointsToCheck[i];
4910 const labelList& pCells = mesh_.pointCells(pointi);
4914 label celli = pCells[i];
4918 mag(cellLevel_[celli]-maxPointLevel[pointi])
4925 <<
"Too big a difference between" 4926 <<
" point-connected cells." <<
nl 4928 <<
" cellLevel:" << cellLevel_[celli]
4929 <<
" uses point:" << pointi
4930 <<
" coord:" << mesh_.points()[pointi]
4931 <<
" which is also used by a cell with level:" 4932 << maxPointLevel[pointi]
5007 if (!cellShapesPtr_)
5011 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes." 5012 <<
" cellLevel:" << cellLevel_.size()
5013 <<
" pointLevel:" << pointLevel_.size()
5020 label nSplitHex = 0;
5021 label nUnrecognised = 0;
5023 forAll(cellLevel_, celli)
5025 if (meshShapes[celli].model().index() == 0)
5027 label level = cellLevel_[celli];
5031 bool haveQuads = matchHexShape
5041 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5052 Pout<<
"hexRef8::cellShapes() :" 5053 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl 5054 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5056 <<
" split-hex:" << nSplitHex <<
nl 5057 <<
" poly :" << nUnrecognised <<
nl 5061 return *cellShapesPtr_;
5069 checkRefinementLevels(-1,
labelList(0));
5074 Pout<<
"hexRef8::getSplitPoints :" 5075 <<
" Calculating unrefineable points" <<
endl;
5079 if (!history_.active())
5082 <<
"Only call if constructed with history capability" 5090 labelList splitMaster(mesh_.nPoints(), -1);
5096 for (label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5098 const labelList& pCells = mesh_.pointCells(pointi);
5100 if (pCells.
size() != 8)
5102 splitMaster[pointi] = -2;
5107 const labelList& visibleCells = history_.visibleCells();
5109 forAll(visibleCells, celli)
5111 const labelList& cPoints = mesh_.cellPoints(celli);
5113 if (visibleCells[celli] != -1 && history_.parentIndex(celli) >= 0)
5115 label parentIndex = history_.parentIndex(celli);
5120 label pointi = cPoints[i];
5122 label masterCelli = splitMaster[pointi];
5124 if (masterCelli == -1)
5131 splitMaster[pointi] = parentIndex;
5132 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5134 else if (masterCelli == -2)
5140 (masterCelli != parentIndex)
5141 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5146 splitMaster[pointi] = -2;
5155 label pointi = cPoints[i];
5157 splitMaster[pointi] = -2;
5165 label facei = mesh_.nInternalFaces();
5166 facei < mesh_.nFaces();
5170 const face&
f = mesh_.faces()[facei];
5174 splitMaster[
f[fp]] = -2;
5181 label nSplitPoints = 0;
5183 forAll(splitMaster, pointi)
5185 if (splitMaster[pointi] >= 0)
5194 forAll(splitMaster, pointi)
5196 if (splitMaster[pointi] >= 0)
5198 splitPoints[nSplitPoints++] = pointi;
5276 Pout<<
"hexRef8::consistentUnrefinement :" 5277 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5283 <<
"maxSet not implemented yet." 5293 bitSet unrefinePoint(mesh_.nPoints(), pointsToUnrefine);
5300 bitSet unrefineCell(mesh_.nCells());
5302 forAll(unrefinePoint, pointi)
5304 if (unrefinePoint.test(pointi))
5306 const labelList& pCells = mesh_.pointCells(pointi);
5308 unrefineCell.
set(pCells);
5320 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5322 label own = mesh_.faceOwner()[facei];
5323 label nei = mesh_.faceNeighbour()[facei];
5325 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5326 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5328 if (ownLevel < (neiLevel-1))
5335 unrefineCell.set(nei);
5346 if (!unrefineCell.test(own))
5352 unrefineCell.unset(own);
5356 else if (neiLevel < (ownLevel-1))
5360 unrefineCell.set(own);
5364 if (!unrefineCell.test(nei))
5370 unrefineCell.unset(nei);
5378 labelList neiLevel(mesh_.nBoundaryFaces());
5382 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5384 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5392 label facei = i+mesh_.nInternalFaces();
5393 label own = mesh_.faceOwner()[facei];
5394 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5396 if (ownLevel < (neiLevel[i]-1))
5400 if (!unrefineCell.test(own))
5406 unrefineCell.unset(own);
5410 else if (neiLevel[i] < (ownLevel-1))
5414 if (unrefineCell.test(own))
5420 unrefineCell.set(own);
5430 Pout<<
"hexRef8::consistentUnrefinement :" 5431 <<
" Changed " << nChanged
5432 <<
" refinement levels due to 2:1 conflicts." 5446 forAll(unrefinePoint, pointi)
5448 if (unrefinePoint.test(pointi))
5450 const labelList& pCells = mesh_.pointCells(pointi);
5454 if (!unrefineCell.test(pCells[j]))
5456 unrefinePoint.unset(pointi);
5468 forAll(unrefinePoint, pointi)
5470 if (unrefinePoint.test(pointi))
5479 forAll(unrefinePoint, pointi)
5481 if (unrefinePoint.test(pointi))
5483 newPointsToUnrefine[nSet++] = pointi;
5487 return newPointsToUnrefine;
5497 if (!history_.active())
5500 <<
"Only call if constructed with history capability" 5506 Pout<<
"hexRef8::setUnrefinement :" 5507 <<
" Checking initial mesh just to make sure" <<
endl;
5511 forAll(cellLevel_, celli)
5513 if (cellLevel_[celli] < 0)
5516 <<
"Illegal cell level " << cellLevel_[celli]
5517 <<
" for cell " << celli
5524 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5527 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.
size());
5529 forAll(splitPointLabels, i)
5531 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5533 cSet.insert(pCells);
5537 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5539 << cSet.size() <<
" cells for unrefinement to" <<
nl 5541 <<
" cellSet " << cSet.objectPath()
5553 forAll(splitPointLabels, i)
5557 splitFaces.insert(
pFaces);
5562 faceRemover_.compatibleRemoves
5570 if (facesToRemove.
size() != splitFaces.size())
5574 const_cast<polyMesh&
>(mesh_).setInstance
5576 mesh_.time().timeName()
5579 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5581 faceSet fSet(mesh_,
"splitFaces", splitFaces);
5583 faceSet removeSet(mesh_,
"facesToRemove", facesToRemove);
5588 <<
"Ininitial set of split points to unrefine does not" 5589 <<
" seem to be consistent or not mid points of refined cells" 5590 <<
" splitPoints:" << splitPointLabels.
size()
5591 <<
" splitFaces:" << splitFaces.size()
5592 <<
" facesToRemove:" << facesToRemove.
size()
5601 forAll(splitPointLabels, i)
5603 label pointi = splitPointLabels[i];
5607 const labelList& pCells = mesh_.pointCells(pointi);
5610 if (pCells.
size() != 8)
5613 <<
"splitPoint " << pointi
5614 <<
" should have 8 cells using it. It has " << pCells
5623 label masterCelli =
min(pCells);
5627 label celli = pCells[j];
5629 label region = cellRegion[celli];
5634 <<
"Ininitial set of split points to unrefine does not" 5635 <<
" seem to be consistent or not mid points" 5636 <<
" of refined cells" <<
nl 5637 <<
"cell:" << celli <<
" on splitPoint " << pointi
5638 <<
" has no region to be merged into" 5642 if (masterCelli != cellRegionMaster[region])
5645 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5646 <<
" in region " << region
5647 <<
" has master:" << cellRegionMaster[region]
5648 <<
" which is not the lowest numbered cell" 5649 <<
" among the pointCells:" << pCells
5658 faceRemover_.setRefinement
5669 forAll(splitPointLabels, i)
5671 label pointi = splitPointLabels[i];
5673 const labelList& pCells = mesh_.pointCells(pointi);
5675 label masterCelli =
min(pCells);
5679 cellLevel_[pCells[j]]--;
5682 history_.combineCells(masterCelli, pCells);
5686 setInstance(mesh_.facesInstance());
5696 cellLevel_.write(valid)
5697 && pointLevel_.write(valid)
5698 && level0Edge_.write(valid);
5700 if (history_.active())
5702 writeOk = writeOk && history_.write(valid);
5726 if (
exists(setsDir/
"cellLevel"))
5728 rm(setsDir/
"cellLevel");
5730 if (
exists(setsDir/
"pointLevel"))
5732 rm(setsDir/
"pointLevel");
5734 if (
exists(setsDir/
"level0Edge"))
5736 rm(setsDir/
"level0Edge");
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
A List of labelList.
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
void size(const label n)
Older name for setAddressableSize.
const labelList & reversePointMap() const
Reverse point map.
A class for handling file names.
readOption readOpt() const noexcept
Get the read option.
virtual const fileName & name() const
The name of the stream.
label nOldCells() const
Number of old cells.
const labelList & cellMap() const
Old cell map.
virtual label iterate(const label maxIter)
Iterate until no changes or maxIter reached.
errorManipArg< error, int > exit(error &err, const int errNo=1)
const fileName & facesInstance() const
Return the current instance directory for faces.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
void setUnrefinement(const labelList &splitPointLabels, polyTopoChange &)
Remove some refinement. Needs to be supplied output of.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void append(const T &val)
Append an element at the end of the list.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
bool active() const
Is there unrefinement history?
static void removeFiles(const polyMesh &)
Helper: remove all relevant files from mesh instance.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool updateFace(const polyMesh &, const label thisFacei, const label neighbourCelli, const refinementData &neighbourInfo, const scalar tol, TrackingData &td)
Influence of neighbouring cell.
constexpr char nl
The newline '\n' character (0x0a)
List< face > faceList
A List of faces.
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Wave propagation of information through grid. Every iteration information goes through one layer of c...
void checkMesh() const
Debug: Check coupled mesh for correctness.
void storeData(const labelList &pointsToStore, const labelList &facesToStore, const labelList &cellsToStore)
Signal points/face/cells for which to store data.
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
labelList consistentSlowRefinement2(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck) const
Like consistentSlowRefinement but uses different meshWave.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
label size() const noexcept
The number of elements in table.
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
Holds information regarding type of cell. Used in inside/outside determination in cellClassification...
labelList consistentUnrefinement(const labelList &pointsToUnrefine, const bool maxSet) const
Given proposed.
fileName objectPath() const
The complete path + object name.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
static scalar propagationTol() noexcept
Access to propagation tolerance.
UList< label > labelUList
A UList of labels.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
void subset(const labelList &pointMap, const labelList &faceMap, const labelList &cellMap)
Update local numbering for subsetted mesh.
fileName path() const
The complete path.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
void setInstance(const fileName &inst)
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
void setSize(const label n)
Alias for resize()
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
void clear()
Clear the list, i.e. set size to zero.
Refinement of (split) hexes using polyTopoChange.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
labelList consistentSlowRefinement(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck, const label maxPointDiff, const labelList &pointsToCheck) const
Like consistentRefinement but slower:
bool exists(const fileName &name, const bool checkGzip=true, const bool followLink=true)
Does the name exist (as DIRECTORY or FILE) in the file system?
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
void distributeCellData(List< T > &values) const
Distribute list of cell data.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Container with cells to refine. Refinement given as single direction.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
void updateMesh(const mapPolyMesh &)
Update local numbering for changed mesh.
virtual bool read()
Read object. If global number of visible cells > 0 becomes active.
const labelList & reverseCellMap() const
Reverse cell map.
int debug
Static debugging option.
Pair< label > labelPair
A pair of labels.
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
const labelList & pointMap() const
Old point map.
Reduction class. If x and y are not equal assign value.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::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].reset(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]
T & last()
Access last element of the list, position [size()-1].
const labelList & visibleCells() const
Per cell in the current mesh (i.e. visible) either -1 (unrefined) or an index into splitCells...
Transfers refinement levels such that slow transition between levels is maintained. Used in FaceCellWave.
labelList consistentRefinement(const labelUList &cellLevel, const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
void distributePointData(List< T > &values) const
Distribute list of point data.
UIndirectList< label > labelUIndList
UIndirectList of labels.
static void removeFiles(const polyMesh &)
Helper: remove all sets files from mesh instance.
vector point
Point is a vector.
label faceLevel(const label facei) const
Gets level such that the face has four points <= level.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
bool write(const bool valid=true) const
Force writing refinement+history to polyMesh directory.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
void operator()(label &x, const label y) const
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A cell is defined as a list of faces with extra functionality.
A collection of cell labels.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const polyBoundaryMesh & patches
HashTable< label, labelPair, Foam::Hash< labelPair > > labelPairLookup
This is a Map of a labelPair to a label. Used for e.g. for face1, face2 to shared edge...
labelList getSplitPoints() const
Return the points at the centre of top-level split cells.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
#define DebugVar(var)
Report a variable name and value.
Mesh consisting of general polyhedral cells.
void checkRefinementLevels(const label maxPointDiff, const labelList &pointsToCheck) const
Debug: Check 2:1 consistency across faces.
List< label > labelList
A List of labels.
All refinement history. Used in unrefinement.
virtual bool write(const bool valid=true) const
Write using setting from DB.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
ListType reorder(const labelUList &oldToNew, const ListType &input, const bool prune=false)
Reorder the elements of a list.
Defines the attributes of an object for which implicit objectRegistry management is supported...
const cellShapeList & cellShapes() const
Utility: get hexes as cell shapes.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
void setFaceInfo(const label facei, const Type &faceInfo)
Set single initial changed face.
const volScalarField & p0
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
label nOldPoints() const
Number of old points.
forAllConstIters(mixture.phases(), phase)
A HashTable to objects of type <T> with a label key.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
bool rm(const fileName &file)
Remove a file (or its gz equivalent), returning true if successful.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
static constexpr const zero Zero
Global zero (0)