60 x = (x ==
y) ? x : value;
68 void Foam::hexRef8::reorder
89 newElems[newI] = elems[i];
97 void Foam::hexRef8::getFaceInfo
107 if (!mesh_.isInternalFace(facei))
109 patchID = mesh_.boundaryMesh().whichPatch(facei);
112 zoneID = mesh_.faceZones().whichZone(facei);
118 const faceZone& fZone = mesh_.faceZones()[zoneID];
135 label patchID, zoneID, zoneFlip;
137 getFaceInfo(facei, patchID, zoneID, zoneFlip);
141 if ((nei == -1) || (own < nei))
194 const label meshFacei,
195 const label meshPointi,
201 if (mesh_.isInternalFace(meshFacei))
287 void Foam::hexRef8::modFace
296 label patchID, zoneID, zoneFlip;
298 getFaceInfo(facei, patchID, zoneID, zoneFlip);
302 (own != mesh_.faceOwner()[facei])
304 mesh_.isInternalFace(facei)
305 && (nei != mesh_.faceNeighbour()[facei])
307 || (newFace != mesh_.faces()[facei])
310 if ((nei == -1) || (own < nei))
351 Foam::scalar Foam::hexRef8::getLevel0EdgeLength()
const 353 if (cellLevel_.size() != mesh_.nCells())
356 <<
"Number of cells in mesh:" << mesh_.nCells()
357 <<
" does not equal size of cellLevel:" << cellLevel_.size()
359 <<
"This might be because of a restart with inconsistent cellLevel." 366 const scalar GREAT2 =
sqr(GREAT);
383 const label cLevel = cellLevel_[celli];
385 const labelList& cEdges = mesh_.cellEdges(celli);
389 label edgeI = cEdges[i];
391 if (edgeLevel[edgeI] == -1)
393 edgeLevel[edgeI] = cLevel;
395 else if (edgeLevel[edgeI] ==
labelMax)
399 else if (edgeLevel[edgeI] != cLevel)
421 const label eLevel = edgeLevel[edgeI];
423 if (eLevel >= 0 && eLevel <
labelMax)
425 const edge&
e = mesh_.edges()[edgeI];
427 scalar edgeLenSqr =
magSqr(e.
vec(mesh_.points()));
429 typEdgeLenSqr[eLevel] =
min(typEdgeLenSqr[eLevel], edgeLenSqr);
441 Pout<<
"hexRef8::getLevel0EdgeLength() :" 442 <<
" After phase1: Edgelengths (squared) per refinementlevel:" 443 << typEdgeLenSqr <<
endl;
457 const label cLevel = cellLevel_[celli];
459 const labelList& cEdges = mesh_.cellEdges(celli);
463 const edge&
e = mesh_.edges()[cEdges[i]];
465 scalar edgeLenSqr =
magSqr(e.
vec(mesh_.points()));
467 maxEdgeLenSqr[cLevel] =
max(maxEdgeLenSqr[cLevel], edgeLenSqr);
476 Pout<<
"hexRef8::getLevel0EdgeLength() :" 477 <<
" Crappy Edgelengths (squared) per refinementlevel:" 478 << maxEdgeLenSqr <<
endl;
485 forAll(typEdgeLenSqr, levelI)
487 if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
489 typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
495 Pout<<
"hexRef8::getLevel0EdgeLength() :" 496 <<
" Final Edgelengths (squared) per refinementlevel:" 497 << typEdgeLenSqr <<
endl;
501 scalar level0Size = -1;
503 forAll(typEdgeLenSqr, levelI)
505 scalar lenSqr = typEdgeLenSqr[levelI];
513 Pout<<
"hexRef8::getLevel0EdgeLength() :" 514 <<
" For level:" << levelI
516 <<
" with equivalent level0 len:" << level0Size
523 if (level0Size == -1)
544 if (cellAnchorPoints[celli].size())
550 return cellAddedCells[celli][index];
557 const face& f = mesh_.faces()[facei];
565 return cellAddedCells[celli][index];
571 Perr<<
"cell:" << celli <<
" anchorPoints:" << cellAnchorPoints[celli]
575 <<
"Could not find point " << pointi
576 <<
" in the anchorPoints for cell " << celli << endl
577 <<
"Does your original mesh obey the 2:1 constraint and" 578 <<
" did you use consistentRefinement to make your cells to refine" 579 <<
" obey this constraint as well?" 592 void Foam::hexRef8::getFaceNeighbours
608 mesh_.faceOwner()[facei],
613 if (mesh_.isInternalFace(facei))
619 mesh_.faceNeighbour()[facei],
639 label level = pointLevel_[f[fp]];
641 if (level < minLevel)
660 label level = pointLevel_[f[fp]];
662 if (level > maxLevel)
676 const label anchorLevel
683 if (pointLevel_[f[fp]] <= anchorLevel)
692 void Foam::hexRef8::dumpCell(
const label celli)
const 695 Pout<<
"hexRef8 : Dumping cell as obj to " << str.
name() <<
endl;
697 const cell& cFaces = mesh_.cells()[celli];
704 const face& f = mesh_.faces()[cFaces[i]];
708 if (pointToObjVert.
insert(f[fp], objVertI))
718 const face& f = mesh_.faces()[cFaces[i]];
722 label pointi = f[fp];
725 str <<
"l " << pointToObjVert[pointi]+1
726 <<
' ' << pointToObjVert[nexPointi]+1 <<
nl;
738 const bool searchForward,
739 const label wantedLevel
746 label pointi = f[fp];
748 if (pointLevel_[pointi] < wantedLevel)
750 dumpCell(mesh_.faceOwner()[facei]);
751 if (mesh_.isInternalFace(facei))
753 dumpCell(mesh_.faceNeighbour()[facei]);
759 <<
" startFp:" << startFp
760 <<
" wantedLevel:" << wantedLevel
763 else if (pointLevel_[pointi] == wantedLevel)
778 dumpCell(mesh_.faceOwner()[facei]);
779 if (mesh_.isInternalFace(facei))
781 dumpCell(mesh_.faceNeighbour()[facei]);
787 <<
" startFp:" << startFp
788 <<
" wantedLevel:" << wantedLevel
798 const face& f = mesh_.faces()[facei];
802 return pointLevel_[f[findMaxLevel(f)]];
806 label ownLevel = cellLevel_[mesh_.faceOwner()[facei]];
808 if (countAnchors(f, ownLevel) == 4)
812 else if (countAnchors(f, ownLevel+1) == 4)
824 void Foam::hexRef8::checkInternalOrientation
839 vector dir(neiPt - ownPt);
844 <<
"cell:" << celli <<
" old face:" << facei
845 <<
" newFace:" << newFace <<
endl 846 <<
" coords:" << compactPoints
847 <<
" ownPt:" << ownPt
848 <<
" neiPt:" << neiPt
852 vector fcToOwn(compactFace.
centre(compactPoints) - ownPt);
854 scalar
s = (fcToOwn&
n) / (dir&n);
856 if (s < 0.1 || s > 0.9)
859 <<
"cell:" << celli <<
" old face:" << facei
860 <<
" newFace:" << newFace <<
endl 861 <<
" coords:" << compactPoints
862 <<
" ownPt:" << ownPt
863 <<
" neiPt:" << neiPt
870 void Foam::hexRef8::checkBoundaryOrientation
876 const point& boundaryPt,
885 vector dir(boundaryPt - ownPt);
890 <<
"cell:" << celli <<
" old face:" << facei
891 <<
" newFace:" << newFace
892 <<
" coords:" << compactPoints
893 <<
" ownPt:" << ownPt
894 <<
" boundaryPt:" << boundaryPt
898 vector fcToOwn(compactFace.
centre(compactPoints) - ownPt);
900 scalar
s = (fcToOwn&dir) /
magSqr(dir);
902 if (s < 0.7 || s > 1.3)
905 <<
"cell:" << celli <<
" old face:" << facei
906 <<
" newFace:" << newFace
907 <<
" coords:" << compactPoints
908 <<
" ownPt:" << ownPt
909 <<
" boundaryPt:" << boundaryPt
919 void Foam::hexRef8::insertEdgeSplit
927 if (p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
931 if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
933 verts.
append(edgeMidPoint[edgeI]);
955 const bool faceOrder,
956 const label edgeMidPointi,
957 const label anchorPointi,
958 const label faceMidPointi,
967 bool changed =
false;
968 bool haveTwoAnchors =
false;
972 if (edgeMidFnd == midPointToAnchors.
end())
974 midPointToAnchors.
insert(edgeMidPointi,
edge(anchorPointi, -1));
978 edge&
e = edgeMidFnd();
980 if (anchorPointi != e[0])
989 if (e[0] != -1 && e[1] != -1)
991 haveTwoAnchors =
true;
995 bool haveTwoFaceMids =
false;
999 if (faceMidFnd == midPointToFaceMids.
end())
1001 midPointToFaceMids.
insert(edgeMidPointi,
edge(faceMidPointi, -1));
1005 edge&
e = faceMidFnd();
1007 if (faceMidPointi != e[0])
1011 e[1] = faceMidPointi;
1016 if (e[0] != -1 && e[1] != -1)
1018 haveTwoFaceMids =
true;
1025 if (changed && haveTwoAnchors && haveTwoFaceMids)
1027 const edge& anchors = midPointToAnchors[edgeMidPointi];
1028 const edge& faceMids = midPointToFaceMids[edgeMidPointi];
1038 if (faceOrder == (mesh_.faceOwner()[facei] == celli))
1040 newFaceVerts.
append(faceMidPointi);
1051 newFaceVerts.
append(edgeMidPointi);
1061 newFaceVerts.
append(otherFaceMidPointi);
1062 newFaceVerts.
append(cellMidPoint[celli]);
1066 newFaceVerts.
append(otherFaceMidPointi);
1076 newFaceVerts.
append(edgeMidPointi);
1086 newFaceVerts.
append(faceMidPointi);
1087 newFaceVerts.
append(cellMidPoint[celli]);
1093 label anchorCell0 = getAnchorCell
1101 label anchorCell1 = getAnchorCell
1114 if (anchorCell0 < anchorCell1)
1119 ownPt = mesh_.points()[anchorPointi];
1120 neiPt = mesh_.points()[anchors.
otherVertex(anchorPointi)];
1129 ownPt = mesh_.points()[anchors.
otherVertex(anchorPointi)];
1130 neiPt = mesh_.points()[anchorPointi];
1137 if (anchorCell0 < anchorCell1)
1139 ownPt = mesh_.points()[anchorPointi];
1140 neiPt = mesh_.points()[anchors.
otherVertex(anchorPointi)];
1144 ownPt = mesh_.points()[anchors.
otherVertex(anchorPointi)];
1145 neiPt = mesh_.points()[anchorPointi];
1148 checkInternalOrientation
1159 return addInternalFace
1177 void Foam::hexRef8::createInternalFaces
1193 const cell& cFaces = mesh_.cells()[celli];
1194 const label cLevel = cellLevel_[celli];
1206 label nFacesAdded = 0;
1210 label facei = cFaces[i];
1212 const face& f = mesh_.faces()[facei];
1213 const labelList& fEdges = mesh_.faceEdges(facei, storage);
1218 label faceMidPointi = -1;
1220 label nAnchors = countAnchors(f, cLevel);
1228 label anchorFp = -1;
1232 if (pointLevel_[f[fp]] <= cLevel)
1240 label edgeMid = findLevel
1248 label faceMid = findLevel
1257 faceMidPointi = f[faceMid];
1259 else if (nAnchors == 4)
1264 faceMidPointi = faceMidPoint[facei];
1268 dumpCell(mesh_.faceOwner()[facei]);
1269 if (mesh_.isInternalFace(facei))
1271 dumpCell(mesh_.faceNeighbour()[facei]);
1275 <<
"nAnchors:" << nAnchors
1276 <<
" facei:" << facei
1288 label point0 = f[fp0];
1290 if (pointLevel_[point0] <= cLevel)
1299 label edgeMidPointi = -1;
1303 if (pointLevel_[f[fp1]] <= cLevel)
1306 label edgeI = fEdges[fp0];
1308 edgeMidPointi = edgeMidPoint[edgeI];
1310 if (edgeMidPointi == -1)
1314 const labelList& cPoints = mesh_.cellPoints(celli);
1317 <<
"cell:" << celli <<
" cLevel:" << cLevel
1318 <<
" cell points:" << cPoints
1321 <<
" face:" << facei
1325 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1326 <<
" faceMidPoint:" << faceMidPoint[facei]
1327 <<
" faceMidPointi:" << faceMidPointi
1335 label edgeMid = findLevel(facei, f, fp1,
true, cLevel+1);
1337 edgeMidPointi = f[edgeMid];
1340 label newFacei = storeMidPointInfo
1363 if (nFacesAdded == 12)
1376 if (pointLevel_[f[fpMin1]] <= cLevel)
1379 label edgeI = fEdges[fpMin1];
1381 edgeMidPointi = edgeMidPoint[edgeI];
1383 if (edgeMidPointi == -1)
1387 const labelList& cPoints = mesh_.cellPoints(celli);
1390 <<
"cell:" << celli <<
" cLevel:" << cLevel
1391 <<
" cell points:" << cPoints
1394 <<
" face:" << facei
1398 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1399 <<
" faceMidPoint:" << faceMidPoint[facei]
1400 <<
" faceMidPointi:" << faceMidPointi
1408 label edgeMid = findLevel
1417 edgeMidPointi = f[edgeMid];
1420 newFacei = storeMidPointInfo
1443 if (nFacesAdded == 12)
1451 if (nFacesAdded == 12)
1459 void Foam::hexRef8::walkFaceToMid
1464 const label startFp,
1468 const face& f = mesh_.faces()[facei];
1469 const labelList& fEdges = mesh_.faceEdges(facei);
1478 if (edgeMidPoint[fEdges[fp]] >= 0)
1480 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1485 if (pointLevel_[f[fp]] <= cLevel)
1491 else if (pointLevel_[f[fp]] == cLevel+1)
1498 else if (pointLevel_[f[fp]] == cLevel+2)
1508 void Foam::hexRef8::walkFaceFromMid
1513 const label startFp,
1517 const face& f = mesh_.faces()[facei];
1518 const labelList& fEdges = mesh_.faceEdges(facei);
1524 if (pointLevel_[f[fp]] <= cLevel)
1529 else if (pointLevel_[f[fp]] == cLevel+1)
1535 else if (pointLevel_[f[fp]] == cLevel+2)
1545 if (edgeMidPoint[fEdges[fp]] >= 0)
1547 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1563 Foam::label Foam::hexRef8::faceConsistentRefinement
1572 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1574 label own = mesh_.faceOwner()[facei];
1575 label ownLevel = cellLevel_[own] + refineCell.
get(own);
1577 label nei = mesh_.faceNeighbour()[facei];
1578 label neiLevel = cellLevel_[nei] + refineCell.
get(nei);
1580 if (ownLevel > (neiLevel+1))
1584 refineCell.
set(nei);
1588 refineCell.
unset(own);
1592 else if (neiLevel > (ownLevel+1))
1596 refineCell.
set(own);
1600 refineCell.
unset(nei);
1609 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1613 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1615 neiLevel[i] = cellLevel_[own] + refineCell.
get(own);
1624 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1625 label ownLevel = cellLevel_[own] + refineCell.
get(own);
1627 if (ownLevel > (neiLevel[i]+1))
1631 refineCell.
unset(own);
1635 else if (neiLevel[i] > (ownLevel+1))
1639 refineCell.
set(own);
1650 void Foam::hexRef8::checkWantedRefinementLevels
1658 refineCell.
set(cellsToRefine[i]);
1661 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1663 label own = mesh_.faceOwner()[facei];
1664 label ownLevel = cellLevel_[own] + refineCell.
get(own);
1666 label nei = mesh_.faceNeighbour()[facei];
1667 label neiLevel = cellLevel_[nei] + refineCell.
get(nei);
1669 if (
mag(ownLevel-neiLevel) > 1)
1675 <<
" current level:" << cellLevel_[own]
1676 <<
" level after refinement:" << ownLevel
1678 <<
"neighbour cell:" << nei
1679 <<
" current level:" << cellLevel_[nei]
1680 <<
" level after refinement:" << neiLevel
1682 <<
"which does not satisfy 2:1 constraints anymore." 1689 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1693 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1695 neiLevel[i] = cellLevel_[own] + refineCell.
get(own);
1704 label facei = i + mesh_.nInternalFaces();
1706 label own = mesh_.faceOwner()[facei];
1707 label ownLevel = cellLevel_[own] + refineCell.
get(own);
1709 if (
mag(ownLevel - neiLevel[i]) > 1)
1711 label patchi = mesh_.boundaryMesh().whichPatch(facei);
1715 <<
"Celllevel does not satisfy 2:1 constraint." 1716 <<
" On coupled face " 1718 <<
" on patch " << patchi <<
" " 1719 << mesh_.boundaryMesh()[
patchi].name()
1720 <<
" owner cell " << own
1721 <<
" current level:" << cellLevel_[own]
1722 <<
" level after refinement:" << ownLevel
1724 <<
" (coupled) neighbour cell will get refinement " 1737 Pout<<
"hexRef8::setInstance(const fileName& inst) : " 1738 <<
"Resetting file instance to " << inst <<
endl;
1741 cellLevel_.instance() = inst;
1742 pointLevel_.instance() = inst;
1743 level0Edge_.instance() = inst;
1744 history_.instance() = inst;
1748 void Foam::hexRef8::collectLevelPoints
1757 if (pointLevel_[f[fp]] <= level)
1765 void Foam::hexRef8::collectLevelPoints
1775 label pointi = meshPoints[f[fp]];
1776 if (pointLevel_[pointi] <= level)
1785 bool Foam::hexRef8::matchHexShape
1788 const label cellLevel,
1792 const cell& cFaces = mesh_.cells()[celli];
1803 label facei = cFaces[i];
1804 const face& f = mesh_.faces()[facei];
1807 collectLevelPoints(f, cellLevel, verts);
1808 if (verts.
size() == 4)
1810 if (mesh_.faceOwner()[facei] != celli)
1821 if (quads.
size() < 6)
1827 label facei = cFaces[i];
1828 const face& f = mesh_.faces()[facei];
1836 collectLevelPoints(f, cellLevel, verts);
1837 if (verts.
size() == 1)
1843 label pointi = f[fp];
1844 if (pointLevel_[pointi] == cellLevel+1)
1847 pointFaces.find(pointi);
1848 if (iter != pointFaces.end())
1874 if (pFaces.
size() == 4)
1880 label facei = pFaces[pFacei];
1881 const face& f = mesh_.faces()[facei];
1882 if (mesh_.faceOwner()[facei] == celli)
1884 fourFaces[pFacei] =
f;
1899 if (edgeLoops.size() == 1)
1905 bigFace.meshPoints(),
1906 bigFace.edgeLoops()[0],
1911 if (verts.
size() == 4)
1922 return (quads.
size() == 6);
1929 Foam::hexRef8::hexRef8(
const polyMesh& mesh,
const bool readHistory)
1937 mesh_.facesInstance(),
1950 mesh_.facesInstance(),
1963 mesh_.facesInstance(),
1975 "refinementHistory",
1976 mesh_.facesInstance(),
1982 (readHistory ? mesh_.nCells() : 0)
1984 faceRemover_(mesh_, GREAT),
1985 savedPointLevel_(0),
2006 <<
"History enabled but number of visible cells " 2009 <<
" is not equal to the number of cells in the mesh " 2016 cellLevel_.
size() != mesh_.nCells()
2017 || pointLevel_.
size() != mesh_.nPoints()
2021 <<
"Restarted from inconsistent cellLevel or pointLevel files." 2025 <<
"Number of cells in mesh:" << mesh_.nCells()
2026 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2027 <<
"Number of points in mesh:" << mesh_.nPoints()
2028 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2047 Foam::hexRef8::hexRef8
2053 const scalar level0Edge
2062 mesh_.facesInstance(),
2075 mesh_.facesInstance(),
2088 mesh_.facesInstance(),
2098 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2105 "refinementHistory",
2106 mesh_.facesInstance(),
2114 faceRemover_(mesh_, GREAT),
2115 savedPointLevel_(0),
2121 <<
"History enabled but number of visible cells in it " 2123 <<
" is not equal to the number of cells in the mesh " 2129 cellLevel_.
size() != mesh_.nCells()
2130 || pointLevel_.
size() != mesh_.nPoints()
2134 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2135 <<
"Number of cells in mesh:" << mesh_.nCells()
2136 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2137 <<
"Number of points in mesh:" << mesh_.nPoints()
2138 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2156 Foam::hexRef8::hexRef8
2161 const scalar level0Edge
2170 mesh_.facesInstance(),
2183 mesh_.facesInstance(),
2196 mesh_.facesInstance(),
2206 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2213 "refinementHistory",
2214 mesh_.facesInstance(),
2224 faceRemover_(mesh_, GREAT),
2225 savedPointLevel_(0),
2230 cellLevel_.
size() != mesh_.nCells()
2231 || pointLevel_.
size() != mesh_.nPoints()
2235 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2236 <<
"Number of cells in mesh:" << mesh_.nCells()
2237 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2238 <<
"Number of points in mesh:" << mesh_.nPoints()
2239 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2272 refineCell.
set(cellsToRefine[i]);
2277 label nChanged = faceConsistentRefinement(maxSet, refineCell);
2283 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2284 <<
" refinement levels due to 2:1 conflicts." 2298 forAll(refineCell, celli)
2300 if (refineCell.
get(celli))
2309 forAll(refineCell, celli)
2311 if (refineCell.
get(celli))
2313 newCellsToRefine[nRefined++] = celli;
2319 checkWantedRefinementLevels(newCellsToRefine);
2322 return newCellsToRefine;
2334 const label maxFaceDiff,
2337 const label maxPointDiff,
2341 const labelList& faceOwner = mesh_.faceOwner();
2342 const labelList& faceNeighbour = mesh_.faceNeighbour();
2345 if (maxFaceDiff <= 0)
2348 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2367 forAll(allCellInfo, celli)
2373 maxFaceDiff*(cellLevel_[celli]+1),
2374 maxFaceDiff*cellLevel_[celli]
2381 label celli = cellsToRefine[i];
2383 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2393 int dummyTrackData = 0;
2401 label facei = facesToCheck[i];
2403 if (allFaceInfo[facei].valid(dummyTrackData))
2407 <<
"Argument facesToCheck seems to have duplicate entries!" 2409 <<
"face:" << facei <<
" occurs at positions " 2417 if (mesh_.isInternalFace(facei))
2422 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2425 label faceRefineCount;
2428 faceCount = neiData.
count() + maxFaceDiff;
2429 faceRefineCount = faceCount + maxFaceDiff;
2433 faceCount = ownData.
count() + maxFaceDiff;
2434 faceRefineCount = faceCount + maxFaceDiff;
2437 seedFaces.append(facei);
2438 seedFacesInfo.append
2446 allFaceInfo[facei] = seedFacesInfo.last();
2450 label faceCount = ownData.
count() + maxFaceDiff;
2451 label faceRefineCount = faceCount + maxFaceDiff;
2453 seedFaces.append(facei);
2454 seedFacesInfo.append
2462 allFaceInfo[facei] = seedFacesInfo.last();
2470 forAll(faceNeighbour, facei)
2473 if (!allFaceInfo[facei].valid(dummyTrackData))
2475 label own = faceOwner[facei];
2476 label nei = faceNeighbour[facei];
2480 if (allCellInfo[own].count() > allCellInfo[nei].count())
2482 allFaceInfo[facei].updateFace
2492 seedFacesInfo.append(allFaceInfo[facei]);
2494 else if (allCellInfo[own].count() < allCellInfo[nei].count())
2496 allFaceInfo[facei].updateFace
2505 seedFaces.append(facei);
2506 seedFacesInfo.append(allFaceInfo[facei]);
2514 for (
label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2517 if (!allFaceInfo[facei].valid(dummyTrackData))
2519 label own = faceOwner[facei];
2532 seedFaces.append(facei);
2533 seedFacesInfo.append(faceData);
2551 Pout<<
"hexRef8::consistentSlowRefinement : Seeded " 2552 << seedFaces.size() <<
" faces between cells with different" 2553 <<
" refinement level." <<
endl;
2557 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2559 seedFacesInfo.clear();
2562 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2570 if (maxPointDiff == -1)
2578 labelList maxPointCount(mesh_.nPoints(), 0);
2580 forAll(maxPointCount, pointi)
2582 label& pLevel = maxPointCount[pointi];
2584 const labelList& pCells = mesh_.pointCells(pointi);
2588 pLevel =
max(pLevel, allCellInfo[pCells[i]].count());
2609 label pointi = pointsToCheck[i];
2614 const labelList& pCells = mesh_.pointCells(pointi);
2618 label celli = pCells[pCelli];
2626 maxPointCount[pointi]
2627 > cellInfo.
count() + maxFaceDiff*maxPointDiff
2635 const cell& cFaces = mesh_.cells()[celli];
2639 label facei = cFaces[cFacei];
2652 if (faceData.
count() > allFaceInfo[facei].count())
2654 changedFacesInfo.insert(facei, faceData);
2661 label nChanged = changedFacesInfo.size();
2671 seedFaces.setCapacity(changedFacesInfo.size());
2672 seedFacesInfo.setCapacity(changedFacesInfo.size());
2676 seedFaces.append(iter.key());
2677 seedFacesInfo.append(iter());
2684 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2686 label own = mesh_.faceOwner()[facei];
2689 + (allCellInfo[own].isRefined() ? 1 : 0);
2691 label nei = mesh_.faceNeighbour()[facei];
2694 + (allCellInfo[nei].isRefined() ? 1 : 0);
2696 if (
mag(ownLevel-neiLevel) > 1)
2702 <<
" current level:" << cellLevel_[own]
2703 <<
" current refData:" << allCellInfo[own]
2704 <<
" level after refinement:" << ownLevel
2706 <<
"neighbour cell:" << nei
2707 <<
" current level:" << cellLevel_[nei]
2708 <<
" current refData:" << allCellInfo[nei]
2709 <<
" level after refinement:" << neiLevel
2711 <<
"which does not satisfy 2:1 constraints anymore." <<
nl 2712 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2721 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2722 labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
2723 labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
2727 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2728 neiLevel[i] = cellLevel_[own];
2729 neiCount[i] = allCellInfo[own].count();
2730 neiRefCount[i] = allCellInfo[own].refinementCount();
2741 label facei = i+mesh_.nInternalFaces();
2743 label own = mesh_.faceOwner()[facei];
2746 + (allCellInfo[own].isRefined() ? 1 : 0);
2750 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2752 if (
mag(ownLevel - nbrLevel) > 1)
2755 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2758 <<
"Celllevel does not satisfy 2:1 constraint." 2759 <<
" On coupled face " 2761 <<
" refData:" << allFaceInfo[facei]
2762 <<
" on patch " << patchi <<
" " 2763 << mesh_.boundaryMesh()[
patchi].name() <<
nl 2764 <<
"owner cell " << own
2765 <<
" current level:" << cellLevel_[own]
2766 <<
" current count:" << allCellInfo[own].count()
2767 <<
" current refCount:" 2768 << allCellInfo[own].refinementCount()
2769 <<
" level after refinement:" << ownLevel
2771 <<
"(coupled) neighbour cell" 2772 <<
" has current level:" << neiLevel[i]
2773 <<
" current count:" << neiCount[i]
2774 <<
" current refCount:" << neiRefCount[i]
2775 <<
" level after refinement:" << nbrLevel
2785 forAll(allCellInfo, celli)
2787 if (allCellInfo[celli].isRefined())
2797 forAll(allCellInfo, celli)
2799 if (allCellInfo[celli].isRefined())
2801 newCellsToRefine[nRefined++] = celli;
2807 Pout<<
"hexRef8::consistentSlowRefinement : From " 2808 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
2809 <<
" cells to refine." <<
endl;
2812 return newCellsToRefine;
2818 const label maxFaceDiff,
2823 const labelList& faceOwner = mesh_.faceOwner();
2824 const labelList& faceNeighbour = mesh_.faceNeighbour();
2826 if (maxFaceDiff <= 0)
2829 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2848 int dummyTrackData = 0;
2854 label celli = cellsToRefine[i];
2859 mesh_.cellCentres()[celli],
2864 forAll(allCellInfo, celli)
2866 if (!allCellInfo[celli].valid(dummyTrackData))
2871 mesh_.cellCentres()[celli],
2887 label facei = facesToCheck[i];
2889 if (allFaceInfo[facei].valid(dummyTrackData))
2893 <<
"Argument facesToCheck seems to have duplicate entries!" 2895 <<
"face:" << facei <<
" occurs at positions " 2900 label own = faceOwner[facei];
2904 allCellInfo[own].valid(dummyTrackData)
2905 ? allCellInfo[own].originLevel()
2909 if (!mesh_.isInternalFace(facei))
2913 const point& fc = mesh_.faceCentres()[facei];
2922 allFaceInfo[facei].updateFace
2934 label nei = faceNeighbour[facei];
2938 allCellInfo[nei].valid(dummyTrackData)
2939 ? allCellInfo[nei].originLevel()
2943 if (ownLevel == neiLevel)
2946 allFaceInfo[facei].updateFace
2955 allFaceInfo[facei].updateFace
2968 allFaceInfo[facei].updateFace
2977 allFaceInfo[facei].updateFace
2989 seedFacesInfo.append(allFaceInfo[facei]);
2996 forAll(faceNeighbour, facei)
2999 if (!allFaceInfo[facei].valid(dummyTrackData))
3001 label own = faceOwner[facei];
3005 allCellInfo[own].valid(dummyTrackData)
3006 ? allCellInfo[own].originLevel()
3010 label nei = faceNeighbour[facei];
3014 allCellInfo[nei].valid(dummyTrackData)
3015 ? allCellInfo[nei].originLevel()
3019 if (ownLevel > neiLevel)
3023 allFaceInfo[facei].updateFace
3032 seedFacesInfo.append(allFaceInfo[facei]);
3034 else if (neiLevel > ownLevel)
3036 seedFaces.append(facei);
3037 allFaceInfo[facei].updateFace
3046 seedFacesInfo.append(allFaceInfo[facei]);
3052 seedFacesInfo.shrink();
3062 mesh_.globalData().nTotalCells()+1,
3103 label celli = cellsToRefine[i];
3105 allCellInfo[celli].originLevel() =
sizeof(
label)*8-2;
3106 allCellInfo[celli].origin() = cc[celli];
3113 forAll(allCellInfo, celli)
3115 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3117 if (wanted > cellLevel_[celli]+1)
3119 refineCell.
set(celli);
3122 faceConsistentRefinement(
true, refineCell);
3126 label nChanged = faceConsistentRefinement(
true, refineCell);
3132 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3133 <<
" refinement levels due to 2:1 conflicts." 3146 forAll(refineCell, celli)
3149 if (refineCell[celli])
3158 forAll(refineCell, celli)
3161 if (refineCell[celli])
3163 newCellsToRefine[nRefined++] = celli;
3169 Pout<<
"hexRef8::consistentSlowRefinement2 : From " 3170 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
3171 <<
" cells to refine." <<
endl;
3176 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3177 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3178 << cellsIn.
size() <<
" to cellSet " 3183 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3184 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3185 << cellsOut.
size() <<
" to cellSet " 3192 forAll(newCellsToRefine, i)
3194 refineCell.
set(newCellsToRefine[i]);
3198 label nChanged = faceConsistentRefinement(
true, refineCell);
3203 mesh_,
"cellsToRefineOut2", newCellsToRefine.
size()
3205 forAll(refineCell, celli)
3207 if (refineCell.
get(celli))
3209 cellsOut2.insert(celli);
3212 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3213 << cellsOut2.size() <<
" to cellSet " 3214 << cellsOut2.objectPath() <<
endl;
3220 forAll(refineCell, celli)
3222 if (refineCell.
get(celli) && !savedRefineCell.
get(celli))
3226 <<
"Cell:" << celli <<
" cc:" 3227 << mesh_.cellCentres()[celli]
3228 <<
" was not marked for refinement but does not obey" 3229 <<
" 2:1 constraints." 3236 return newCellsToRefine;
3249 Pout<<
"hexRef8::setRefinement :" 3250 <<
" Checking initial mesh just to make sure" <<
endl;
3259 savedPointLevel_.
clear();
3260 savedCellLevel_.
clear();
3265 forAll(cellLevel_, celli)
3267 newCellLevel.append(cellLevel_[celli]);
3270 forAll(pointLevel_, pointi)
3272 newPointLevel.append(pointLevel_[pointi]);
3278 Pout<<
"hexRef8::setRefinement :" 3279 <<
" Allocating " << cellLabels.
size() <<
" cell midpoints." 3287 labelList cellMidPoint(mesh_.nCells(), -1);
3291 label celli = cellLabels[i];
3293 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3299 mesh_.cellCentres()[celli],
3306 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3312 cellSet splitCells(mesh_,
"splitCells", cellLabels.
size());
3314 forAll(cellMidPoint, celli)
3316 if (cellMidPoint[celli] >= 0)
3318 splitCells.insert(celli);
3322 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3323 <<
" cells to split to cellSet " << splitCells.objectPath()
3336 Pout<<
"hexRef8::setRefinement :" 3337 <<
" Allocating edge midpoints." 3346 labelList edgeMidPoint(mesh_.nEdges(), -1);
3349 forAll(cellMidPoint, celli)
3351 if (cellMidPoint[celli] >= 0)
3353 const labelList& cEdges = mesh_.cellEdges(celli);
3357 label edgeI = cEdges[i];
3359 const edge&
e = mesh_.edges()[edgeI];
3363 pointLevel_[e[0]] <= cellLevel_[celli]
3364 && pointLevel_[e[1]] <= cellLevel_[celli]
3367 edgeMidPoint[edgeI] = 12345;
3394 forAll(edgeMidPoint, edgeI)
3396 if (edgeMidPoint[edgeI] >= 0)
3399 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3407 point(-GREAT, -GREAT, -GREAT)
3412 forAll(edgeMidPoint, edgeI)
3414 if (edgeMidPoint[edgeI] >= 0)
3419 const edge&
e = mesh_.edges()[edgeI];
3432 newPointLevel(edgeMidPoint[edgeI]) =
3445 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3447 forAll(edgeMidPoint, edgeI)
3449 if (edgeMidPoint[edgeI] >= 0)
3451 const edge&
e = mesh_.edges()[edgeI];
3457 Pout<<
"hexRef8::setRefinement :" 3458 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3468 Pout<<
"hexRef8::setRefinement :" 3469 <<
" Allocating face midpoints." 3475 labelList faceAnchorLevel(mesh_.nFaces());
3477 for (
label facei = 0; facei < mesh_.nFaces(); facei++)
3479 faceAnchorLevel[facei] =
faceLevel(facei);
3484 labelList faceMidPoint(mesh_.nFaces(), -1);
3489 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3491 if (faceAnchorLevel[facei] >= 0)
3493 label own = mesh_.faceOwner()[facei];
3494 label ownLevel = cellLevel_[own];
3495 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3497 label nei = mesh_.faceNeighbour()[facei];
3498 label neiLevel = cellLevel_[nei];
3499 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3503 newOwnLevel > faceAnchorLevel[facei]
3504 || newNeiLevel > faceAnchorLevel[facei]
3507 faceMidPoint[facei] = 12345;
3520 labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
3524 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3525 label ownLevel = cellLevel_[own];
3526 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3528 newNeiLevel[i] = newOwnLevel;
3538 label facei = i+mesh_.nInternalFaces();
3540 if (faceAnchorLevel[facei] >= 0)
3542 label own = mesh_.faceOwner()[facei];
3543 label ownLevel = cellLevel_[own];
3544 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3548 newOwnLevel > faceAnchorLevel[facei]
3549 || newNeiLevel[i] > faceAnchorLevel[facei]
3552 faceMidPoint[facei] = 12345;
3577 mesh_.nFaces()-mesh_.nInternalFaces(),
3578 point(-GREAT, -GREAT, -GREAT)
3583 label facei = i+mesh_.nInternalFaces();
3585 if (faceMidPoint[facei] >= 0)
3587 bFaceMids[i] = mesh_.faceCentres()[facei];
3597 forAll(faceMidPoint, facei)
3599 if (faceMidPoint[facei] >= 0)
3604 const face& f = mesh_.faces()[facei];
3611 facei < mesh_.nInternalFaces()
3612 ? mesh_.faceCentres()[facei]
3613 : bFaceMids[facei-mesh_.nInternalFaces()]
3623 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3630 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.
size());
3632 forAll(faceMidPoint, facei)
3634 if (faceMidPoint[facei] >= 0)
3636 splitFaces.insert(facei);
3640 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3641 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3662 Pout<<
"hexRef8::setRefinement :" 3663 <<
" Finding cell anchorPoints (8 per cell)" 3674 labelList nAnchorPoints(mesh_.nCells(), 0);
3676 forAll(cellMidPoint, celli)
3678 if (cellMidPoint[celli] >= 0)
3680 cellAnchorPoints[celli].
setSize(8);
3684 forAll(pointLevel_, pointi)
3686 const labelList& pCells = mesh_.pointCells(pointi);
3690 label celli = pCells[pCelli];
3694 cellMidPoint[celli] >= 0
3695 && pointLevel_[pointi] <= cellLevel_[celli]
3698 if (nAnchorPoints[celli] == 8)
3703 <<
" of level " << cellLevel_[celli]
3704 <<
" uses more than 8 points of equal or" 3705 <<
" lower level" <<
nl 3706 <<
"Points so far:" << cellAnchorPoints[celli]
3710 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3715 forAll(cellMidPoint, celli)
3717 if (cellMidPoint[celli] >= 0)
3719 if (nAnchorPoints[celli] != 8)
3723 const labelList& cPoints = mesh_.cellPoints(celli);
3727 <<
" of level " << cellLevel_[celli]
3728 <<
" does not seem to have 8 points of equal or" 3729 <<
" lower level" <<
endl 3730 <<
"cellPoints:" << cPoints <<
endl 3745 Pout<<
"hexRef8::setRefinement :" 3746 <<
" Adding cells (1 per anchorPoint)" 3753 forAll(cellAnchorPoints, celli)
3755 const labelList& cAnchors = cellAnchorPoints[celli];
3757 if (cAnchors.
size() == 8)
3759 labelList& cAdded = cellAddedCells[celli];
3766 newCellLevel[celli] = cellLevel_[celli]+1;
3769 for (
label i = 1; i < 8; i++)
3779 mesh_.cellZones().whichZone(celli)
3783 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3798 Pout<<
"hexRef8::setRefinement :" 3799 <<
" Marking faces to be handled" 3807 forAll(cellMidPoint, celli)
3809 if (cellMidPoint[celli] >= 0)
3811 const cell& cFaces = mesh_.cells()[celli];
3815 affectedFace.set(cFaces[i]);
3820 forAll(faceMidPoint, facei)
3822 if (faceMidPoint[facei] >= 0)
3824 affectedFace.set(facei);
3828 forAll(edgeMidPoint, edgeI)
3830 if (edgeMidPoint[edgeI] >= 0)
3832 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3836 affectedFace.set(eFaces[i]);
3848 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3851 forAll(faceMidPoint, facei)
3853 if (faceMidPoint[facei] >= 0 && affectedFace.get(facei))
3859 const face& f = mesh_.faces()[facei];
3863 bool modifiedFace =
false;
3865 label anchorLevel = faceAnchorLevel[facei];
3871 label pointi = f[fp];
3873 if (pointLevel_[pointi] <= anchorLevel)
3879 faceVerts.
append(pointi);
3895 faceVerts.
append(faceMidPoint[facei]);
3928 if (mesh_.isInternalFace(facei))
3930 label oldOwn = mesh_.faceOwner()[facei];
3931 label oldNei = mesh_.faceNeighbour()[facei];
3933 checkInternalOrientation
3938 mesh_.cellCentres()[oldOwn],
3939 mesh_.cellCentres()[oldNei],
3945 label oldOwn = mesh_.faceOwner()[facei];
3947 checkBoundaryOrientation
3952 mesh_.cellCentres()[oldOwn],
3953 mesh_.faceCentres()[facei],
3962 modifiedFace =
true;
3964 modFace(meshMod, facei, newFace, own, nei);
3968 addFace(meshMod, facei, newFace, own, nei);
3974 affectedFace.unset(facei);
3984 Pout<<
"hexRef8::setRefinement :" 3985 <<
" Adding edge splits to unsplit faces" 3992 forAll(edgeMidPoint, edgeI)
3994 if (edgeMidPoint[edgeI] >= 0)
3998 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
4002 label facei = eFaces[i];
4004 if (faceMidPoint[facei] < 0 && affectedFace.get(facei))
4008 const face& f = mesh_.faces()[facei];
4009 const labelList& fEdges = mesh_.faceEdges
4019 newFaceVerts.append(f[fp]);
4021 label edgeI = fEdges[fp];
4023 if (edgeMidPoint[edgeI] >= 0)
4025 newFaceVerts.
append(edgeMidPoint[edgeI]);
4034 label anchorFp = findMinLevel(f);
4051 if (mesh_.isInternalFace(facei))
4053 label oldOwn = mesh_.faceOwner()[facei];
4054 label oldNei = mesh_.faceNeighbour()[facei];
4056 checkInternalOrientation
4061 mesh_.cellCentres()[oldOwn],
4062 mesh_.cellCentres()[oldNei],
4068 label oldOwn = mesh_.faceOwner()[facei];
4070 checkBoundaryOrientation
4075 mesh_.cellCentres()[oldOwn],
4076 mesh_.faceCentres()[facei],
4082 modFace(meshMod, facei, newFace, own, nei);
4085 affectedFace.unset(facei);
4097 Pout<<
"hexRef8::setRefinement :" 4098 <<
" Changing owner/neighbour for otherwise unaffected faces" 4102 forAll(affectedFace, facei)
4104 if (affectedFace.get(facei))
4106 const face& f = mesh_.faces()[facei];
4110 label anchorFp = findMinLevel(f);
4124 modFace(meshMod, facei, f, own, nei);
4127 affectedFace.unset(facei);
4142 Pout<<
"hexRef8::setRefinement :" 4143 <<
" Create new internal faces for split cells" 4147 forAll(cellMidPoint, celli)
4149 if (cellMidPoint[celli] >= 0)
4174 forAll(cellMidPoint, celli)
4176 if (cellMidPoint[celli] >= 0)
4178 minPointi =
min(minPointi, cellMidPoint[celli]);
4179 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4182 forAll(faceMidPoint, facei)
4184 if (faceMidPoint[facei] >= 0)
4186 minPointi =
min(minPointi, faceMidPoint[facei]);
4187 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4190 forAll(edgeMidPoint, edgeI)
4192 if (edgeMidPoint[edgeI] >= 0)
4194 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4195 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4199 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4202 <<
"Added point labels not consecutive to existing mesh points." 4204 <<
"mesh_.nPoints():" << mesh_.nPoints()
4205 <<
" minPointi:" << minPointi
4206 <<
" maxPointi:" << maxPointi
4211 pointLevel_.
transfer(newPointLevel);
4226 Pout<<
"hexRef8::setRefinement :" 4227 <<
" Updating refinement history to " << cellLevel_.
size()
4228 <<
" cells" <<
endl;
4234 forAll(cellAddedCells, celli)
4236 const labelList& addedCells = cellAddedCells[celli];
4238 if (addedCells.
size())
4252 label celli = cellLabels[i];
4254 refinedCells[i].
transfer(cellAddedCells[celli]);
4257 return refinedCells;
4268 savedPointLevel_.
resize(2*pointsToStore.
size());
4271 label pointi = pointsToStore[i];
4272 savedPointLevel_.
insert(pointi, pointLevel_[pointi]);
4275 savedCellLevel_.
resize(2*cellsToStore.
size());
4278 label celli = cellsToStore[i];
4279 savedCellLevel_.
insert(celli, cellLevel_[celli]);
4291 updateMesh(map, dummyMap, dummyMap, dummyMap);
4309 Pout<<
"hexRef8::updateMesh :" 4310 <<
" Updating various lists" 4319 Pout<<
"hexRef8::updateMesh :" 4322 <<
" nCells:" << mesh_.nCells()
4324 <<
" cellLevel_:" << cellLevel_.
size()
4327 <<
" nPoints:" << mesh_.nPoints()
4329 <<
" pointLevel_:" << pointLevel_.
size()
4333 if (reverseCellMap.
size() == cellLevel_.
size())
4339 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4347 forAll(cellMap, newCelli)
4349 label oldCelli = cellMap[newCelli];
4353 newCellLevel[newCelli] = -1;
4357 newCellLevel[newCelli] = cellLevel_[oldCelli];
4367 label newCelli = iter.key();
4368 label storedCelli = iter();
4372 if (fnd == savedCellLevel_.
end())
4375 <<
"Problem : trying to restore old value for new cell " 4376 << newCelli <<
" but cannot find old cell " << storedCelli
4377 <<
" in map of stored values " << savedCellLevel_
4380 cellLevel_[newCelli] = fnd();
4401 if (reversePointMap.
size() == pointLevel_.
size())
4404 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4417 if (oldPointi == -1)
4430 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4433 pointLevel_.
transfer(newPointLevel);
4441 label storedPointi = iter();
4445 if (fnd == savedPointLevel_.
end())
4448 <<
"Problem : trying to restore old value for new point " 4449 << newPointi <<
" but cannot find old point " 4451 <<
" in map of stored values " << savedPointLevel_
4483 cellShapesPtr_.clear();
4498 Pout<<
"hexRef8::subset :" 4499 <<
" Updating various lists" 4506 <<
"Subsetting will not work in combination with unrefinement." 4508 <<
"Proceed at your own risk." <<
endl;
4516 forAll(cellMap, newCelli)
4518 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4527 <<
"cellLevel_ contains illegal value -1 after mapping:" 4542 pointLevel_.
transfer(newPointLevel);
4548 <<
"pointLevel_ contains illegal value -1 after mapping:" 4557 history_.
subset(pointMap, faceMap, cellMap);
4567 cellShapesPtr_.clear();
4576 Pout<<
"hexRef8::distribute :" 4577 <<
" Updating various lists" 4596 cellShapesPtr_.clear();
4602 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4606 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:" 4607 << smallDim <<
endl;
4617 labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
4620 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4643 label own = mesh_.faceOwner()[facei];
4644 label bFacei = facei-mesh_.nInternalFaces();
4650 <<
"Faces do not seem to be correct across coupled" 4651 <<
" boundaries" <<
endl 4652 <<
"Coupled face " << facei
4653 <<
" between owner " << own
4654 <<
" on patch " << pp.
name()
4655 <<
" and coupled neighbour " << nei[bFacei]
4656 <<
" has two faces connected to it:" 4672 scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
4675 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4683 label facei = i+mesh_.nInternalFaces();
4685 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4687 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4689 const face& f = mesh_.faces()[facei];
4690 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4692 dumpCell(mesh_.faceOwner()[facei]);
4695 <<
"Faces do not seem to be correct across coupled" 4696 <<
" boundaries" <<
endl 4697 <<
"Coupled face " << facei
4698 <<
" on patch " << patchi
4699 <<
" " << mesh_.boundaryMesh()[
patchi].
name()
4701 <<
" has face area:" << magArea
4702 <<
" (coupled) neighbour face area differs:" 4704 <<
" to within tolerance " << smallDim
4714 labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
4718 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4726 label facei = i+mesh_.nInternalFaces();
4728 const face& f = mesh_.faces()[facei];
4730 if (f.
size() != nVerts[i])
4732 dumpCell(mesh_.faceOwner()[facei]);
4734 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4737 <<
"Faces do not seem to be correct across coupled" 4738 <<
" boundaries" <<
endl 4739 <<
"Coupled face " << facei
4740 <<
" on patch " << patchi
4741 <<
" " << mesh_.boundaryMesh()[
patchi].
name()
4743 <<
" has size:" << f.
size()
4744 <<
" (coupled) neighbour face has size:" 4756 pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
4760 label facei = i+mesh_.nInternalFaces();
4761 const point& fc = mesh_.faceCentres()[facei];
4762 const face& f = mesh_.faces()[facei];
4763 const vector anchorVec(mesh_.points()[f[0]] - fc);
4765 anchorPoints[i] = anchorVec;
4775 label facei = i+mesh_.nInternalFaces();
4776 const point& fc = mesh_.faceCentres()[facei];
4777 const face& f = mesh_.faces()[facei];
4778 const vector anchorVec(mesh_.points()[f[0]] - fc);
4780 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4782 dumpCell(mesh_.faceOwner()[facei]);
4784 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4787 <<
"Faces do not seem to be correct across coupled" 4788 <<
" boundaries" <<
endl 4789 <<
"Coupled face " << facei
4790 <<
" on patch " << patchi
4791 <<
" " << mesh_.boundaryMesh()[
patchi].
name()
4793 <<
" has anchor vector:" << anchorVec
4794 <<
" (coupled) neighbour face anchor vector differs:" 4796 <<
" to within tolerance " << smallDim
4804 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4811 const label maxPointDiff,
4817 Pout<<
"hexRef8::checkRefinementLevels :" 4818 <<
" Checking 2:1 refinement level" <<
endl;
4823 cellLevel_.
size() != mesh_.nCells()
4824 || pointLevel_.
size() != mesh_.nPoints()
4828 <<
"cellLevel size should be number of cells" 4829 <<
" and pointLevel size should be number of points."<<
nl 4830 <<
"cellLevel:" << cellLevel_.
size()
4831 <<
" mesh.nCells():" << mesh_.nCells() <<
nl 4832 <<
"pointLevel:" << pointLevel_.
size()
4833 <<
" mesh.nPoints():" << mesh_.nPoints()
4843 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4845 label own = mesh_.faceOwner()[facei];
4846 label nei = mesh_.faceNeighbour()[facei];
4848 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4854 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl 4855 <<
"On face " << facei <<
" owner cell " << own
4856 <<
" has refinement " << cellLevel_[own]
4857 <<
" neighbour cell " << nei <<
" has refinement " 4864 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
4868 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4870 neiLevel[i] = cellLevel_[own];
4878 label facei = i+mesh_.nInternalFaces();
4880 label own = mesh_.faceOwner()[facei];
4882 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4886 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4889 <<
"Celllevel does not satisfy 2:1 constraint." 4890 <<
" On coupled face " << facei
4891 <<
" on patch " << patchi <<
" " 4892 << mesh_.boundaryMesh()[
patchi].name()
4893 <<
" owner cell " << own <<
" has refinement " 4895 <<
" (coupled) neighbour cell has refinement " 4918 forAll(syncPointLevel, pointi)
4920 if (pointLevel_[pointi] != syncPointLevel[pointi])
4923 <<
"PointLevel is not consistent across coupled patches." 4925 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4926 <<
" has level " << pointLevel_[pointi]
4927 <<
" whereas the coupled point has level " 4928 << syncPointLevel[pointi]
4936 if (maxPointDiff != -1)
4939 labelList maxPointLevel(mesh_.nPoints(), 0);
4941 forAll(maxPointLevel, pointi)
4943 const labelList& pCells = mesh_.pointCells(pointi);
4945 label& pLevel = maxPointLevel[pointi];
4949 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4965 label pointi = pointsToCheck[i];
4967 const labelList& pCells = mesh_.pointCells(pointi);
4971 label celli = pCells[i];
4975 mag(cellLevel_[celli]-maxPointLevel[pointi])
4982 <<
"Too big a difference between" 4983 <<
" point-connected cells." <<
nl 4985 <<
" cellLevel:" << cellLevel_[celli]
4986 <<
" uses point:" << pointi
4987 <<
" coord:" << mesh_.points()[pointi]
4988 <<
" which is also used by a cell with level:" 4989 << maxPointLevel[pointi]
5064 if (cellShapesPtr_.empty())
5068 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes." 5069 <<
" cellLevel:" << cellLevel_.
size()
5070 <<
" pointLevel:" << pointLevel_.
size()
5077 label nSplitHex = 0;
5078 label nUnrecognised = 0;
5080 forAll(cellLevel_, celli)
5082 if (meshShapes[celli].model().index() == 0)
5084 label level = cellLevel_[celli];
5088 bool haveQuads = matchHexShape
5098 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5109 Pout<<
"hexRef8::cellShapes() :" 5110 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl 5111 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5113 <<
" split-hex:" << nSplitHex <<
nl 5114 <<
" poly :" << nUnrecognised <<
nl 5118 return cellShapesPtr_();
5131 Pout<<
"hexRef8::getSplitPoints :" 5132 <<
" Calculating unrefineable points" <<
endl;
5139 <<
"Only call if constructed with history capability" 5147 labelList splitMaster(mesh_.nPoints(), -1);
5148 labelList splitMasterLevel(mesh_.nPoints(), 0);
5153 for (
label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5155 const labelList& pCells = mesh_.pointCells(pointi);
5157 if (pCells.
size() != 8)
5159 splitMaster[pointi] = -2;
5166 forAll(visibleCells, celli)
5168 const labelList& cPoints = mesh_.cellPoints(celli);
5170 if (visibleCells[celli] != -1 && history_.
parentIndex(celli) >= 0)
5177 label pointi = cPoints[i];
5179 label masterCelli = splitMaster[pointi];
5181 if (masterCelli == -1)
5188 splitMaster[pointi] = parentIndex;
5189 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5191 else if (masterCelli == -2)
5197 (masterCelli != parentIndex)
5198 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5203 splitMaster[pointi] = -2;
5212 label pointi = cPoints[i];
5214 splitMaster[pointi] = -2;
5222 label facei = mesh_.nInternalFaces();
5223 facei < mesh_.nFaces();
5227 const face& f = mesh_.faces()[facei];
5231 splitMaster[f[fp]] = -2;
5238 label nSplitPoints = 0;
5240 forAll(splitMaster, pointi)
5242 if (splitMaster[pointi] >= 0)
5251 forAll(splitMaster, pointi)
5253 if (splitMaster[pointi] >= 0)
5255 splitPoints[nSplitPoints++] = pointi;
5333 Pout<<
"hexRef8::consistentUnrefinement :" 5334 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5340 <<
"maxSet not implemented yet." 5352 forAll(pointsToUnrefine, i)
5354 label pointi = pointsToUnrefine[i];
5356 unrefinePoint.set(pointi);
5367 forAll(unrefinePoint, pointi)
5369 if (unrefinePoint.get(pointi))
5371 const labelList& pCells = mesh_.pointCells(pointi);
5375 unrefineCell.set(pCells[j]);
5388 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5390 label own = mesh_.faceOwner()[facei];
5391 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5393 label nei = mesh_.faceNeighbour()[facei];
5394 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5396 if (ownLevel < (neiLevel-1))
5403 unrefineCell.set(nei);
5414 if (unrefineCell.get(own) == 0)
5420 unrefineCell.unset(own);
5424 else if (neiLevel < (ownLevel-1))
5428 unrefineCell.set(own);
5432 if (unrefineCell.get(nei) == 0)
5438 unrefineCell.unset(nei);
5446 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
5450 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5452 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5460 label facei = i+mesh_.nInternalFaces();
5461 label own = mesh_.faceOwner()[facei];
5462 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5464 if (ownLevel < (neiLevel[i]-1))
5468 if (unrefineCell.get(own) == 0)
5474 unrefineCell.unset(own);
5478 else if (neiLevel[i] < (ownLevel-1))
5482 if (unrefineCell.get(own) == 1)
5488 unrefineCell.set(own);
5498 Pout<<
"hexRef8::consistentUnrefinement :" 5499 <<
" Changed " << nChanged
5500 <<
" refinement levels due to 2:1 conflicts." 5514 forAll(unrefinePoint, pointi)
5516 if (unrefinePoint.get(pointi))
5518 const labelList& pCells = mesh_.pointCells(pointi);
5522 if (!unrefineCell.get(pCells[j]))
5524 unrefinePoint.unset(pointi);
5536 forAll(unrefinePoint, pointi)
5538 if (unrefinePoint.get(pointi))
5547 forAll(unrefinePoint, pointi)
5549 if (unrefinePoint.get(pointi))
5551 newPointsToUnrefine[nSet++] = pointi;
5555 return newPointsToUnrefine;
5568 <<
"Only call if constructed with history capability" 5574 Pout<<
"hexRef8::setUnrefinement :" 5575 <<
" Checking initial mesh just to make sure" <<
endl;
5579 forAll(cellLevel_, celli)
5581 if (cellLevel_[celli] < 0)
5584 <<
"Illegal cell level " << cellLevel_[celli]
5585 <<
" for cell " << celli
5592 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5595 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.
size());
5597 forAll(splitPointLabels, i)
5599 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5603 cSet.insert(pCells[j]);
5608 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5610 << cSet.size() <<
" cells for unrefinement to" <<
nl 5612 <<
" cellSet " << cSet.objectPath()
5624 forAll(splitPointLabels, i)
5630 splitFaces.insert(pFaces[j]);
5644 if (facesToRemove.
size() != splitFaces.size())
5647 <<
"Ininitial set of split points to unrefine does not" 5648 <<
" seem to be consistent or not mid points of refined cells" 5657 forAll(splitPointLabels, i)
5659 label pointi = splitPointLabels[i];
5663 const labelList& pCells = mesh_.pointCells(pointi);
5666 if (pCells.
size() != 8)
5669 <<
"splitPoint " << pointi
5670 <<
" should have 8 cells using it. It has " << pCells
5683 label celli = pCells[j];
5685 label region = cellRegion[celli];
5690 <<
"Ininitial set of split points to unrefine does not" 5691 <<
" seem to be consistent or not mid points" 5692 <<
" of refined cells" <<
nl 5693 <<
"cell:" << celli <<
" on splitPoint " << pointi
5694 <<
" has no region to be merged into" 5698 if (masterCelli != cellRegionMaster[region])
5701 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5702 <<
" in region " << region
5703 <<
" has master:" << cellRegionMaster[region]
5704 <<
" which is not the lowest numbered cell" 5705 <<
" among the pointCells:" << pCells
5725 forAll(splitPointLabels, i)
5727 label pointi = splitPointLabels[i];
5729 const labelList& pCells = mesh_.pointCells(pointi);
5735 cellLevel_[pCells[j]]--;
5753 && pointLevel_.
write()
5754 && level0Edge_.
write();
5758 writeOk = writeOk && history_.
write();
const DynamicList< point > & points() const
Points. Shrunk after constructing mesh (or calling of compact())
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
const labelList & visibleCells() const
Per cell in the current mesh (i.e. visible) either -1 (unrefined)
point centre(const pointField &) const
Centre point of face.
fileName objectPath() const
Return complete path + object name.
scalar level0EdgeLength() const
Typical edge length between unrefined points.
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.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A class for handling file names.
label iterate(const label maxIter)
Iterate until no changes or maxIter reached. Returns actual.
errorManipArg< error, int > exit(error &err, const int errNo=1)
const boolList & flipMap() const
Return face flip map.
Class describing modification of a face.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
A face is a list of labels corresponding to mesh vertices.
void setFaceInfo(const labelList &changedFaces, const List< Type > &changedFacesInfo)
Set initial changed faces.
const double e
Elementary charge.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void setUnrefinement(const labelList &splitPointLabels, polyTopoChange &)
Remove some refinement. Needs to be supplied output of.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
point centre(const pointField &) const
Return centre (centroid)
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.
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
void size(const label)
Override size to be inconsistent with allocated storage.
fileCheckTypes
Types of communications.
label nOldPoints() const
Number of old points.
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 updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
const labelList & pointMap() const
Old point map.
void storeData(const labelList &pointsToStore, const labelList &facesToStore, const labelList &cellsToStore)
Signal points/face/cells for which to store data.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
Class containing data for point addition.
label otherVertex(const label a) const
Given one vertex, return the other.
label size() const
Return number of elements in table.
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
void distributeCellData(List< T > &lst) const
Distribute list of cell data.
Holds information regarding type of cell. Used in inside/outside determination in cellClassification...
vector normal(const pointField &) const
Vector normal; magnitude is equal to area of face.
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurences of given element. Linear search.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
labelList getSplitPoints() const
Return the points at the centre of top-level split cells.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
A face addition data class. A face can be inflated either from a point or from another face and can e...
void subset(const labelList &pointMap, const labelList &faceMap, const labelList &cellMap)
Update local numbering for subsetted mesh.
bool write() const
Force writing refinement+history to polyMesh directory.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
iterator find(const label &)
Find and return an iterator set at the hashedEntry.
void setInstance(const fileName &inst)
static fileCheckTypes fileModificationChecking
void distribute(const mapDistributePolyMesh &)
Force recalculation of locally stored data for mesh distribution.
A topoSetSource to select a faceSet from cells.
An STL-conforming iterator.
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Reduction class. If x and y are not equal assign value.
Class containing data for cell addition.
A list of faces which address into the list of points.
A List obtained as a section of another List.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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))
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
An ordered pair of two objects of type <T> with first() and second() elements.
void combineCells(const label masterCelli, const labelList &combinedCells)
Store combining 8 cells into master.
void set(const PackedList< 1 > &)
Set specified bits.
label start() const
Return start label of this patch in the polyMesh face list.
void operator()(label &x, const label y) const
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
label compatibleRemoves(const labelList &inPiercedFaces, labelList &cellRegion, labelList &cellRegionMaster, labelList &outPiercedFaces) const
Find faces including those with cells which have the same mastercell.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
bool active() const
Is there unrefinement history?
labelList consistentUnrefinement(const labelList &pointsToUnrefine, const bool maxSet) const
Given proposed.
label parentIndex(const label celli) const
Get parent of cell.
void clear()
Clear all entries from table.
virtual const fileName & name() const
Return the name of the stream.
label faceLevel(const label facei) const
Gets level such that the face has four points <= level.
void append(const T &)
Append an element at the end of the list.
void subset(const labelList &pointMap, const labelList &faceMap, const labelList &cellMap)
Update numbering for subsetting.
Xfer< List< T > > xfer()
Transfer contents to the Xfer container as a plain 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]
const cellShapeList & cellShapes() const
Utility: get hexes as cell shapes.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Pair< label > labelPair
Label pair.
labelList consistentSlowRefinement2(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck) const
Like consistentSlowRefinement but uses different meshWave.
static const label labelMax
List< label > labelList
A List of labels.
void distributePointData(List< T > &lst) const
Distribute list of point data.
labelList consistentRefinement(const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
const word & name() const
Return name.
An STL-conforming hash table.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
void updateMesh(const mapPolyMesh &)
Update local numbering for changed mesh.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
void reverse(UList< T > &, const label n)
virtual const fileName & name() const
Return the name of the stream.
virtual bool read()
Read object. If global number of visible cells > 0 becomes active.
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
label refinementCount() const
labelList consistentSlowRefinement(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck, const label maxPointDiff, const labelList &pointsToCheck) const
Like consistentRefinement but slower:
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
void storeSplit(const label celli, const labelList &addedCells)
Store splitting of cell into 8.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
prefixOSstream Perr(cerr,"Perr")
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
label nOldCells() const
Number of old cells.
void updateMesh(const mapPolyMesh &)
Update numbering for mesh changes.
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
Transfers refinement levels such that slow transition between levels is maintained. Used in FaceCellWave.
void setSize(const label)
Reset size of List.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
void checkMesh() const
Debug: Check coupled mesh for correctness.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
const labelList & cellMap() const
Old cell map.
A cell is defined as a list of faces with extra functionality.
A collection of cell labels.
A List with indirect addressing.
Direct mesh changes based on v1.3 polyTopoChange syntax.
void flip()
Flip the face in-place.
readOption readOpt() const
const labelList & reverseCellMap() const
Reverse cell map.
vector vec(const pointField &) const
Return the vector (end - start)
virtual bool write() const
Write using setting from DB.
void unset(const PackedList< 1 > &)
Unset specified bits.
void setRefinement(const labelList &piercedFaces, const labelList &cellRegion, const labelList &cellRegionMaster, polyTopoChange &) const
Play commands into polyTopoChange to remove faces.
dimensioned< scalar > mag(const dimensioned< Type > &)
bool headerOk()
Read and check header info.
face reverseFace() const
Return face with reverse direction.
void resize(const label newSize)
Resize the hash table for efficiency.
virtual Ostream & write(const token &)=0
Write next token to stream.
Mesh consisting of general polyhedral cells.
A subset of mesh faces organised as a primitive patch.
All refinement history. Used in unrefinement.
void checkRefinementLevels(const label maxPointDiff, const labelList &pointsToCheck) const
Debug: Check 2:1 consistency across faces.
A patch is a list of labels that address the faces in the global face list.
T & last()
Return the last element of the list.
void clear()
Clear the addressed list, i.e. set the size to zero.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
void resize(const label nCells)
Extend/shrink storage. additional visibleCells_ elements get.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
static const label labelMin
Transfers refinement levels such that slow transition between levels is maintained. Used in FaceCellWave.
const labelList & reversePointMap() const
Reverse point map.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
unsigned int get(const label) const
Get value at index I.