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(),
1983 (readHistory ? mesh_.nCells() : 0)
1985 faceRemover_(mesh_, GREAT),
1986 savedPointLevel_(0),
2007 <<
"History enabled but number of visible cells " 2010 <<
" is not equal to the number of cells in the mesh " 2017 cellLevel_.
size() != mesh_.nCells()
2018 || pointLevel_.
size() != mesh_.nPoints()
2022 <<
"Restarted from inconsistent cellLevel or pointLevel files." 2026 <<
"Number of cells in mesh:" << mesh_.nCells()
2027 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2028 <<
"Number of points in mesh:" << mesh_.nPoints()
2029 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2048 Foam::hexRef8::hexRef8
2054 const scalar level0Edge
2063 mesh_.facesInstance(),
2076 mesh_.facesInstance(),
2089 mesh_.facesInstance(),
2099 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2106 "refinementHistory",
2107 mesh_.facesInstance(),
2115 faceRemover_(mesh_, GREAT),
2116 savedPointLevel_(0),
2122 <<
"History enabled but number of visible cells in it " 2124 <<
" is not equal to the number of cells in the mesh " 2130 cellLevel_.
size() != mesh_.nCells()
2131 || pointLevel_.
size() != mesh_.nPoints()
2135 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2136 <<
"Number of cells in mesh:" << mesh_.nCells()
2137 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2138 <<
"Number of points in mesh:" << mesh_.nPoints()
2139 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2157 Foam::hexRef8::hexRef8
2162 const scalar level0Edge
2171 mesh_.facesInstance(),
2184 mesh_.facesInstance(),
2197 mesh_.facesInstance(),
2207 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2214 "refinementHistory",
2215 mesh_.facesInstance(),
2225 faceRemover_(mesh_, GREAT),
2226 savedPointLevel_(0),
2231 cellLevel_.
size() != mesh_.nCells()
2232 || pointLevel_.
size() != mesh_.nPoints()
2236 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2237 <<
"Number of cells in mesh:" << mesh_.nCells()
2238 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2239 <<
"Number of points in mesh:" << mesh_.nPoints()
2240 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2273 refineCell.
set(cellsToRefine[i]);
2278 label nChanged = faceConsistentRefinement(maxSet, refineCell);
2284 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2285 <<
" refinement levels due to 2:1 conflicts." 2299 forAll(refineCell, celli)
2301 if (refineCell.
get(celli))
2310 forAll(refineCell, celli)
2312 if (refineCell.
get(celli))
2314 newCellsToRefine[nRefined++] = celli;
2320 checkWantedRefinementLevels(newCellsToRefine);
2323 return newCellsToRefine;
2335 const label maxFaceDiff,
2338 const label maxPointDiff,
2342 const labelList& faceOwner = mesh_.faceOwner();
2343 const labelList& faceNeighbour = mesh_.faceNeighbour();
2346 if (maxFaceDiff <= 0)
2349 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2368 forAll(allCellInfo, celli)
2374 maxFaceDiff*(cellLevel_[celli]+1),
2375 maxFaceDiff*cellLevel_[celli]
2382 label celli = cellsToRefine[i];
2384 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2394 int dummyTrackData = 0;
2402 label facei = facesToCheck[i];
2404 if (allFaceInfo[facei].valid(dummyTrackData))
2408 <<
"Argument facesToCheck seems to have duplicate entries!" 2410 <<
"face:" << facei <<
" occurs at positions " 2418 if (mesh_.isInternalFace(facei))
2423 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2426 label faceRefineCount;
2429 faceCount = neiData.
count() + maxFaceDiff;
2430 faceRefineCount = faceCount + maxFaceDiff;
2434 faceCount = ownData.
count() + maxFaceDiff;
2435 faceRefineCount = faceCount + maxFaceDiff;
2438 seedFaces.append(facei);
2439 seedFacesInfo.append
2447 allFaceInfo[facei] = seedFacesInfo.last();
2451 label faceCount = ownData.
count() + maxFaceDiff;
2452 label faceRefineCount = faceCount + maxFaceDiff;
2454 seedFaces.append(facei);
2455 seedFacesInfo.append
2463 allFaceInfo[facei] = seedFacesInfo.last();
2471 forAll(faceNeighbour, facei)
2474 if (!allFaceInfo[facei].valid(dummyTrackData))
2476 label own = faceOwner[facei];
2477 label nei = faceNeighbour[facei];
2481 if (allCellInfo[own].count() > allCellInfo[nei].count())
2483 allFaceInfo[facei].updateFace
2493 seedFacesInfo.append(allFaceInfo[facei]);
2495 else if (allCellInfo[own].count() < allCellInfo[nei].count())
2497 allFaceInfo[facei].updateFace
2506 seedFaces.append(facei);
2507 seedFacesInfo.append(allFaceInfo[facei]);
2515 for (
label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2518 if (!allFaceInfo[facei].valid(dummyTrackData))
2520 label own = faceOwner[facei];
2533 seedFaces.append(facei);
2534 seedFacesInfo.append(faceData);
2552 Pout<<
"hexRef8::consistentSlowRefinement : Seeded " 2553 << seedFaces.size() <<
" faces between cells with different" 2554 <<
" refinement level." <<
endl;
2558 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2560 seedFacesInfo.clear();
2563 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2571 if (maxPointDiff == -1)
2579 labelList maxPointCount(mesh_.nPoints(), 0);
2581 forAll(maxPointCount, pointi)
2583 label& pLevel = maxPointCount[pointi];
2585 const labelList& pCells = mesh_.pointCells(pointi);
2589 pLevel =
max(pLevel, allCellInfo[pCells[i]].count());
2610 label pointi = pointsToCheck[i];
2615 const labelList& pCells = mesh_.pointCells(pointi);
2619 label celli = pCells[pCelli];
2627 maxPointCount[pointi]
2628 > cellInfo.
count() + maxFaceDiff*maxPointDiff
2636 const cell& cFaces = mesh_.cells()[celli];
2640 label facei = cFaces[cFacei];
2653 if (faceData.
count() > allFaceInfo[facei].count())
2655 changedFacesInfo.insert(facei, faceData);
2662 label nChanged = changedFacesInfo.size();
2672 seedFaces.setCapacity(changedFacesInfo.size());
2673 seedFacesInfo.setCapacity(changedFacesInfo.size());
2677 seedFaces.append(iter.key());
2678 seedFacesInfo.append(iter());
2685 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2687 label own = mesh_.faceOwner()[facei];
2690 + (allCellInfo[own].isRefined() ? 1 : 0);
2692 label nei = mesh_.faceNeighbour()[facei];
2695 + (allCellInfo[nei].isRefined() ? 1 : 0);
2697 if (
mag(ownLevel-neiLevel) > 1)
2703 <<
" current level:" << cellLevel_[own]
2704 <<
" current refData:" << allCellInfo[own]
2705 <<
" level after refinement:" << ownLevel
2707 <<
"neighbour cell:" << nei
2708 <<
" current level:" << cellLevel_[nei]
2709 <<
" current refData:" << allCellInfo[nei]
2710 <<
" level after refinement:" << neiLevel
2712 <<
"which does not satisfy 2:1 constraints anymore." <<
nl 2713 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2722 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2723 labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
2724 labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
2728 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2729 neiLevel[i] = cellLevel_[own];
2730 neiCount[i] = allCellInfo[own].count();
2731 neiRefCount[i] = allCellInfo[own].refinementCount();
2742 label facei = i+mesh_.nInternalFaces();
2744 label own = mesh_.faceOwner()[facei];
2747 + (allCellInfo[own].isRefined() ? 1 : 0);
2751 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2753 if (
mag(ownLevel - nbrLevel) > 1)
2756 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2759 <<
"Celllevel does not satisfy 2:1 constraint." 2760 <<
" On coupled face " 2762 <<
" refData:" << allFaceInfo[facei]
2763 <<
" on patch " << patchi <<
" " 2764 << mesh_.boundaryMesh()[
patchi].name() <<
nl 2765 <<
"owner cell " << own
2766 <<
" current level:" << cellLevel_[own]
2767 <<
" current count:" << allCellInfo[own].count()
2768 <<
" current refCount:" 2769 << allCellInfo[own].refinementCount()
2770 <<
" level after refinement:" << ownLevel
2772 <<
"(coupled) neighbour cell" 2773 <<
" has current level:" << neiLevel[i]
2774 <<
" current count:" << neiCount[i]
2775 <<
" current refCount:" << neiRefCount[i]
2776 <<
" level after refinement:" << nbrLevel
2786 forAll(allCellInfo, celli)
2788 if (allCellInfo[celli].isRefined())
2798 forAll(allCellInfo, celli)
2800 if (allCellInfo[celli].isRefined())
2802 newCellsToRefine[nRefined++] = celli;
2808 Pout<<
"hexRef8::consistentSlowRefinement : From " 2809 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
2810 <<
" cells to refine." <<
endl;
2813 return newCellsToRefine;
2819 const label maxFaceDiff,
2824 const labelList& faceOwner = mesh_.faceOwner();
2825 const labelList& faceNeighbour = mesh_.faceNeighbour();
2827 if (maxFaceDiff <= 0)
2830 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2849 int dummyTrackData = 0;
2855 label celli = cellsToRefine[i];
2860 mesh_.cellCentres()[celli],
2865 forAll(allCellInfo, celli)
2867 if (!allCellInfo[celli].valid(dummyTrackData))
2872 mesh_.cellCentres()[celli],
2888 label facei = facesToCheck[i];
2890 if (allFaceInfo[facei].valid(dummyTrackData))
2894 <<
"Argument facesToCheck seems to have duplicate entries!" 2896 <<
"face:" << facei <<
" occurs at positions " 2901 label own = faceOwner[facei];
2905 allCellInfo[own].valid(dummyTrackData)
2906 ? allCellInfo[own].originLevel()
2910 if (!mesh_.isInternalFace(facei))
2914 const point& fc = mesh_.faceCentres()[facei];
2923 allFaceInfo[facei].updateFace
2935 label nei = faceNeighbour[facei];
2939 allCellInfo[nei].valid(dummyTrackData)
2940 ? allCellInfo[nei].originLevel()
2944 if (ownLevel == neiLevel)
2947 allFaceInfo[facei].updateFace
2956 allFaceInfo[facei].updateFace
2969 allFaceInfo[facei].updateFace
2978 allFaceInfo[facei].updateFace
2990 seedFacesInfo.append(allFaceInfo[facei]);
2997 forAll(faceNeighbour, facei)
3000 if (!allFaceInfo[facei].valid(dummyTrackData))
3002 label own = faceOwner[facei];
3006 allCellInfo[own].valid(dummyTrackData)
3007 ? allCellInfo[own].originLevel()
3011 label nei = faceNeighbour[facei];
3015 allCellInfo[nei].valid(dummyTrackData)
3016 ? allCellInfo[nei].originLevel()
3020 if (ownLevel > neiLevel)
3024 allFaceInfo[facei].updateFace
3033 seedFacesInfo.append(allFaceInfo[facei]);
3035 else if (neiLevel > ownLevel)
3037 seedFaces.append(facei);
3038 allFaceInfo[facei].updateFace
3047 seedFacesInfo.append(allFaceInfo[facei]);
3053 seedFacesInfo.shrink();
3063 mesh_.globalData().nTotalCells()+1,
3104 label celli = cellsToRefine[i];
3106 allCellInfo[celli].originLevel() =
sizeof(
label)*8-2;
3107 allCellInfo[celli].origin() = cc[celli];
3114 forAll(allCellInfo, celli)
3116 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3118 if (wanted > cellLevel_[celli]+1)
3120 refineCell.
set(celli);
3123 faceConsistentRefinement(
true, refineCell);
3127 label nChanged = faceConsistentRefinement(
true, refineCell);
3133 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3134 <<
" refinement levels due to 2:1 conflicts." 3147 forAll(refineCell, celli)
3150 if (refineCell[celli])
3159 forAll(refineCell, celli)
3162 if (refineCell[celli])
3164 newCellsToRefine[nRefined++] = celli;
3170 Pout<<
"hexRef8::consistentSlowRefinement2 : From " 3171 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
3172 <<
" cells to refine." <<
endl;
3177 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3178 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3179 << cellsIn.
size() <<
" to cellSet " 3184 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3185 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3186 << cellsOut.
size() <<
" to cellSet " 3193 forAll(newCellsToRefine, i)
3195 refineCell.
set(newCellsToRefine[i]);
3199 label nChanged = faceConsistentRefinement(
true, refineCell);
3204 mesh_,
"cellsToRefineOut2", newCellsToRefine.
size()
3206 forAll(refineCell, celli)
3208 if (refineCell.
get(celli))
3210 cellsOut2.insert(celli);
3213 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3214 << cellsOut2.size() <<
" to cellSet " 3215 << cellsOut2.objectPath() <<
endl;
3221 forAll(refineCell, celli)
3223 if (refineCell.
get(celli) && !savedRefineCell.
get(celli))
3227 <<
"Cell:" << celli <<
" cc:" 3228 << mesh_.cellCentres()[celli]
3229 <<
" was not marked for refinement but does not obey" 3230 <<
" 2:1 constraints." 3237 return newCellsToRefine;
3250 Pout<<
"hexRef8::setRefinement :" 3251 <<
" Checking initial mesh just to make sure" <<
endl;
3260 savedPointLevel_.
clear();
3261 savedCellLevel_.
clear();
3266 forAll(cellLevel_, celli)
3268 newCellLevel.append(cellLevel_[celli]);
3271 forAll(pointLevel_, pointi)
3273 newPointLevel.append(pointLevel_[pointi]);
3279 Pout<<
"hexRef8::setRefinement :" 3280 <<
" Allocating " << cellLabels.
size() <<
" cell midpoints." 3288 labelList cellMidPoint(mesh_.nCells(), -1);
3292 label celli = cellLabels[i];
3294 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3300 mesh_.cellCentres()[celli],
3307 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3313 cellSet splitCells(mesh_,
"splitCells", cellLabels.
size());
3315 forAll(cellMidPoint, celli)
3317 if (cellMidPoint[celli] >= 0)
3319 splitCells.insert(celli);
3323 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3324 <<
" cells to split to cellSet " << splitCells.objectPath()
3337 Pout<<
"hexRef8::setRefinement :" 3338 <<
" Allocating edge midpoints." 3347 labelList edgeMidPoint(mesh_.nEdges(), -1);
3350 forAll(cellMidPoint, celli)
3352 if (cellMidPoint[celli] >= 0)
3354 const labelList& cEdges = mesh_.cellEdges(celli);
3358 label edgeI = cEdges[i];
3360 const edge&
e = mesh_.edges()[edgeI];
3364 pointLevel_[e[0]] <= cellLevel_[celli]
3365 && pointLevel_[e[1]] <= cellLevel_[celli]
3368 edgeMidPoint[edgeI] = 12345;
3395 forAll(edgeMidPoint, edgeI)
3397 if (edgeMidPoint[edgeI] >= 0)
3400 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3408 point(-GREAT, -GREAT, -GREAT)
3413 forAll(edgeMidPoint, edgeI)
3415 if (edgeMidPoint[edgeI] >= 0)
3420 const edge&
e = mesh_.edges()[edgeI];
3433 newPointLevel(edgeMidPoint[edgeI]) =
3446 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3448 forAll(edgeMidPoint, edgeI)
3450 if (edgeMidPoint[edgeI] >= 0)
3452 const edge&
e = mesh_.edges()[edgeI];
3458 Pout<<
"hexRef8::setRefinement :" 3459 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3469 Pout<<
"hexRef8::setRefinement :" 3470 <<
" Allocating face midpoints." 3476 labelList faceAnchorLevel(mesh_.nFaces());
3478 for (
label facei = 0; facei < mesh_.nFaces(); facei++)
3480 faceAnchorLevel[facei] =
faceLevel(facei);
3485 labelList faceMidPoint(mesh_.nFaces(), -1);
3490 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3492 if (faceAnchorLevel[facei] >= 0)
3494 label own = mesh_.faceOwner()[facei];
3495 label ownLevel = cellLevel_[own];
3496 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3498 label nei = mesh_.faceNeighbour()[facei];
3499 label neiLevel = cellLevel_[nei];
3500 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3504 newOwnLevel > faceAnchorLevel[facei]
3505 || newNeiLevel > faceAnchorLevel[facei]
3508 faceMidPoint[facei] = 12345;
3521 labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
3525 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3526 label ownLevel = cellLevel_[own];
3527 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3529 newNeiLevel[i] = newOwnLevel;
3539 label facei = i+mesh_.nInternalFaces();
3541 if (faceAnchorLevel[facei] >= 0)
3543 label own = mesh_.faceOwner()[facei];
3544 label ownLevel = cellLevel_[own];
3545 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3549 newOwnLevel > faceAnchorLevel[facei]
3550 || newNeiLevel[i] > faceAnchorLevel[facei]
3553 faceMidPoint[facei] = 12345;
3578 mesh_.nFaces()-mesh_.nInternalFaces(),
3579 point(-GREAT, -GREAT, -GREAT)
3584 label facei = i+mesh_.nInternalFaces();
3586 if (faceMidPoint[facei] >= 0)
3588 bFaceMids[i] = mesh_.faceCentres()[facei];
3598 forAll(faceMidPoint, facei)
3600 if (faceMidPoint[facei] >= 0)
3605 const face& f = mesh_.faces()[facei];
3612 facei < mesh_.nInternalFaces()
3613 ? mesh_.faceCentres()[facei]
3614 : bFaceMids[facei-mesh_.nInternalFaces()]
3624 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3631 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.
size());
3633 forAll(faceMidPoint, facei)
3635 if (faceMidPoint[facei] >= 0)
3637 splitFaces.insert(facei);
3641 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3642 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3663 Pout<<
"hexRef8::setRefinement :" 3664 <<
" Finding cell anchorPoints (8 per cell)" 3675 labelList nAnchorPoints(mesh_.nCells(), 0);
3677 forAll(cellMidPoint, celli)
3679 if (cellMidPoint[celli] >= 0)
3681 cellAnchorPoints[celli].
setSize(8);
3685 forAll(pointLevel_, pointi)
3687 const labelList& pCells = mesh_.pointCells(pointi);
3691 label celli = pCells[pCelli];
3695 cellMidPoint[celli] >= 0
3696 && pointLevel_[pointi] <= cellLevel_[celli]
3699 if (nAnchorPoints[celli] == 8)
3704 <<
" of level " << cellLevel_[celli]
3705 <<
" uses more than 8 points of equal or" 3706 <<
" lower level" <<
nl 3707 <<
"Points so far:" << cellAnchorPoints[celli]
3711 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3716 forAll(cellMidPoint, celli)
3718 if (cellMidPoint[celli] >= 0)
3720 if (nAnchorPoints[celli] != 8)
3724 const labelList& cPoints = mesh_.cellPoints(celli);
3728 <<
" of level " << cellLevel_[celli]
3729 <<
" does not seem to have 8 points of equal or" 3730 <<
" lower level" <<
endl 3731 <<
"cellPoints:" << cPoints <<
endl 3746 Pout<<
"hexRef8::setRefinement :" 3747 <<
" Adding cells (1 per anchorPoint)" 3754 forAll(cellAnchorPoints, celli)
3756 const labelList& cAnchors = cellAnchorPoints[celli];
3758 if (cAnchors.
size() == 8)
3760 labelList& cAdded = cellAddedCells[celli];
3767 newCellLevel[celli] = cellLevel_[celli]+1;
3770 for (
label i = 1; i < 8; i++)
3780 mesh_.cellZones().whichZone(celli)
3784 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3799 Pout<<
"hexRef8::setRefinement :" 3800 <<
" Marking faces to be handled" 3808 forAll(cellMidPoint, celli)
3810 if (cellMidPoint[celli] >= 0)
3812 const cell& cFaces = mesh_.cells()[celli];
3816 affectedFace.set(cFaces[i]);
3821 forAll(faceMidPoint, facei)
3823 if (faceMidPoint[facei] >= 0)
3825 affectedFace.set(facei);
3829 forAll(edgeMidPoint, edgeI)
3831 if (edgeMidPoint[edgeI] >= 0)
3833 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3837 affectedFace.set(eFaces[i]);
3849 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3852 forAll(faceMidPoint, facei)
3854 if (faceMidPoint[facei] >= 0 && affectedFace.get(facei))
3860 const face& f = mesh_.faces()[facei];
3864 bool modifiedFace =
false;
3866 label anchorLevel = faceAnchorLevel[facei];
3872 label pointi = f[fp];
3874 if (pointLevel_[pointi] <= anchorLevel)
3880 faceVerts.
append(pointi);
3896 faceVerts.
append(faceMidPoint[facei]);
3929 if (mesh_.isInternalFace(facei))
3931 label oldOwn = mesh_.faceOwner()[facei];
3932 label oldNei = mesh_.faceNeighbour()[facei];
3934 checkInternalOrientation
3939 mesh_.cellCentres()[oldOwn],
3940 mesh_.cellCentres()[oldNei],
3946 label oldOwn = mesh_.faceOwner()[facei];
3948 checkBoundaryOrientation
3953 mesh_.cellCentres()[oldOwn],
3954 mesh_.faceCentres()[facei],
3963 modifiedFace =
true;
3965 modFace(meshMod, facei, newFace, own, nei);
3969 addFace(meshMod, facei, newFace, own, nei);
3975 affectedFace.unset(facei);
3985 Pout<<
"hexRef8::setRefinement :" 3986 <<
" Adding edge splits to unsplit faces" 3993 forAll(edgeMidPoint, edgeI)
3995 if (edgeMidPoint[edgeI] >= 0)
3999 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
4003 label facei = eFaces[i];
4005 if (faceMidPoint[facei] < 0 && affectedFace.get(facei))
4009 const face& f = mesh_.faces()[facei];
4010 const labelList& fEdges = mesh_.faceEdges
4020 newFaceVerts.append(f[fp]);
4022 label edgeI = fEdges[fp];
4024 if (edgeMidPoint[edgeI] >= 0)
4026 newFaceVerts.
append(edgeMidPoint[edgeI]);
4035 label anchorFp = findMinLevel(f);
4052 if (mesh_.isInternalFace(facei))
4054 label oldOwn = mesh_.faceOwner()[facei];
4055 label oldNei = mesh_.faceNeighbour()[facei];
4057 checkInternalOrientation
4062 mesh_.cellCentres()[oldOwn],
4063 mesh_.cellCentres()[oldNei],
4069 label oldOwn = mesh_.faceOwner()[facei];
4071 checkBoundaryOrientation
4076 mesh_.cellCentres()[oldOwn],
4077 mesh_.faceCentres()[facei],
4083 modFace(meshMod, facei, newFace, own, nei);
4086 affectedFace.unset(facei);
4098 Pout<<
"hexRef8::setRefinement :" 4099 <<
" Changing owner/neighbour for otherwise unaffected faces" 4103 forAll(affectedFace, facei)
4105 if (affectedFace.get(facei))
4107 const face& f = mesh_.faces()[facei];
4111 label anchorFp = findMinLevel(f);
4125 modFace(meshMod, facei, f, own, nei);
4128 affectedFace.unset(facei);
4143 Pout<<
"hexRef8::setRefinement :" 4144 <<
" Create new internal faces for split cells" 4148 forAll(cellMidPoint, celli)
4150 if (cellMidPoint[celli] >= 0)
4175 forAll(cellMidPoint, celli)
4177 if (cellMidPoint[celli] >= 0)
4179 minPointi =
min(minPointi, cellMidPoint[celli]);
4180 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4183 forAll(faceMidPoint, facei)
4185 if (faceMidPoint[facei] >= 0)
4187 minPointi =
min(minPointi, faceMidPoint[facei]);
4188 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4191 forAll(edgeMidPoint, edgeI)
4193 if (edgeMidPoint[edgeI] >= 0)
4195 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4196 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4200 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4203 <<
"Added point labels not consecutive to existing mesh points." 4205 <<
"mesh_.nPoints():" << mesh_.nPoints()
4206 <<
" minPointi:" << minPointi
4207 <<
" maxPointi:" << maxPointi
4212 pointLevel_.
transfer(newPointLevel);
4227 Pout<<
"hexRef8::setRefinement :" 4228 <<
" Updating refinement history to " << cellLevel_.
size()
4229 <<
" cells" <<
endl;
4235 forAll(cellAddedCells, celli)
4237 const labelList& addedCells = cellAddedCells[celli];
4239 if (addedCells.
size())
4253 label celli = cellLabels[i];
4255 refinedCells[i].
transfer(cellAddedCells[celli]);
4258 return refinedCells;
4269 savedPointLevel_.
resize(2*pointsToStore.
size());
4272 label pointi = pointsToStore[i];
4273 savedPointLevel_.
insert(pointi, pointLevel_[pointi]);
4276 savedCellLevel_.
resize(2*cellsToStore.
size());
4279 label celli = cellsToStore[i];
4280 savedCellLevel_.
insert(celli, cellLevel_[celli]);
4292 updateMesh(map, dummyMap, dummyMap, dummyMap);
4310 Pout<<
"hexRef8::updateMesh :" 4311 <<
" Updating various lists" 4320 Pout<<
"hexRef8::updateMesh :" 4323 <<
" nCells:" << mesh_.nCells()
4325 <<
" cellLevel_:" << cellLevel_.
size()
4328 <<
" nPoints:" << mesh_.nPoints()
4330 <<
" pointLevel_:" << pointLevel_.
size()
4334 if (reverseCellMap.
size() == cellLevel_.
size())
4340 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4348 forAll(cellMap, newCelli)
4350 label oldCelli = cellMap[newCelli];
4354 newCellLevel[newCelli] = -1;
4358 newCellLevel[newCelli] = cellLevel_[oldCelli];
4368 label newCelli = iter.key();
4369 label storedCelli = iter();
4373 if (fnd == savedCellLevel_.
end())
4376 <<
"Problem : trying to restore old value for new cell " 4377 << newCelli <<
" but cannot find old cell " << storedCelli
4378 <<
" in map of stored values " << savedCellLevel_
4381 cellLevel_[newCelli] = fnd();
4402 if (reversePointMap.
size() == pointLevel_.
size())
4405 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4418 if (oldPointi == -1)
4431 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4434 pointLevel_.
transfer(newPointLevel);
4442 label storedPointi = iter();
4446 if (fnd == savedPointLevel_.
end())
4449 <<
"Problem : trying to restore old value for new point " 4450 << newPointi <<
" but cannot find old point " 4452 <<
" in map of stored values " << savedPointLevel_
4484 cellShapesPtr_.clear();
4499 Pout<<
"hexRef8::subset :" 4500 <<
" Updating various lists" 4507 <<
"Subsetting will not work in combination with unrefinement." 4509 <<
"Proceed at your own risk." <<
endl;
4517 forAll(cellMap, newCelli)
4519 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4528 <<
"cellLevel_ contains illegal value -1 after mapping:" 4543 pointLevel_.
transfer(newPointLevel);
4549 <<
"pointLevel_ contains illegal value -1 after mapping:" 4558 history_.
subset(pointMap, faceMap, cellMap);
4568 cellShapesPtr_.clear();
4577 Pout<<
"hexRef8::distribute :" 4578 <<
" Updating various lists" 4597 cellShapesPtr_.clear();
4603 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4607 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:" 4608 << smallDim <<
endl;
4618 labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
4621 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4644 label own = mesh_.faceOwner()[facei];
4645 label bFacei = facei-mesh_.nInternalFaces();
4651 <<
"Faces do not seem to be correct across coupled" 4652 <<
" boundaries" <<
endl 4653 <<
"Coupled face " << facei
4654 <<
" between owner " << own
4655 <<
" on patch " << pp.
name()
4656 <<
" and coupled neighbour " << nei[bFacei]
4657 <<
" has two faces connected to it:" 4673 scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
4676 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4684 label facei = i+mesh_.nInternalFaces();
4686 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4688 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4690 const face& f = mesh_.faces()[facei];
4691 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4693 dumpCell(mesh_.faceOwner()[facei]);
4696 <<
"Faces do not seem to be correct across coupled" 4697 <<
" boundaries" <<
endl 4698 <<
"Coupled face " << facei
4699 <<
" on patch " << patchi
4700 <<
" " << mesh_.boundaryMesh()[
patchi].
name()
4702 <<
" has face area:" << magArea
4703 <<
" (coupled) neighbour face area differs:" 4705 <<
" to within tolerance " << smallDim
4715 labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
4719 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4727 label facei = i+mesh_.nInternalFaces();
4729 const face& f = mesh_.faces()[facei];
4731 if (f.
size() != nVerts[i])
4733 dumpCell(mesh_.faceOwner()[facei]);
4735 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4738 <<
"Faces do not seem to be correct across coupled" 4739 <<
" boundaries" <<
endl 4740 <<
"Coupled face " << facei
4741 <<
" on patch " << patchi
4742 <<
" " << mesh_.boundaryMesh()[
patchi].
name()
4744 <<
" has size:" << f.
size()
4745 <<
" (coupled) neighbour face has size:" 4757 pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
4761 label facei = i+mesh_.nInternalFaces();
4762 const point& fc = mesh_.faceCentres()[facei];
4763 const face& f = mesh_.faces()[facei];
4764 const vector anchorVec(mesh_.points()[f[0]] - fc);
4766 anchorPoints[i] = anchorVec;
4776 label facei = i+mesh_.nInternalFaces();
4777 const point& fc = mesh_.faceCentres()[facei];
4778 const face& f = mesh_.faces()[facei];
4779 const vector anchorVec(mesh_.points()[f[0]] - fc);
4781 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4783 dumpCell(mesh_.faceOwner()[facei]);
4785 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4788 <<
"Faces do not seem to be correct across coupled" 4789 <<
" boundaries" <<
endl 4790 <<
"Coupled face " << facei
4791 <<
" on patch " << patchi
4792 <<
" " << mesh_.boundaryMesh()[
patchi].
name()
4794 <<
" has anchor vector:" << anchorVec
4795 <<
" (coupled) neighbour face anchor vector differs:" 4797 <<
" to within tolerance " << smallDim
4805 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4812 const label maxPointDiff,
4818 Pout<<
"hexRef8::checkRefinementLevels :" 4819 <<
" Checking 2:1 refinement level" <<
endl;
4824 cellLevel_.
size() != mesh_.nCells()
4825 || pointLevel_.
size() != mesh_.nPoints()
4829 <<
"cellLevel size should be number of cells" 4830 <<
" and pointLevel size should be number of points."<<
nl 4831 <<
"cellLevel:" << cellLevel_.
size()
4832 <<
" mesh.nCells():" << mesh_.nCells() <<
nl 4833 <<
"pointLevel:" << pointLevel_.
size()
4834 <<
" mesh.nPoints():" << mesh_.nPoints()
4844 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4846 label own = mesh_.faceOwner()[facei];
4847 label nei = mesh_.faceNeighbour()[facei];
4849 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4855 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl 4856 <<
"On face " << facei <<
" owner cell " << own
4857 <<
" has refinement " << cellLevel_[own]
4858 <<
" neighbour cell " << nei <<
" has refinement " 4865 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
4869 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4871 neiLevel[i] = cellLevel_[own];
4879 label facei = i+mesh_.nInternalFaces();
4881 label own = mesh_.faceOwner()[facei];
4883 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4887 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4890 <<
"Celllevel does not satisfy 2:1 constraint." 4891 <<
" On coupled face " << facei
4892 <<
" on patch " << patchi <<
" " 4893 << mesh_.boundaryMesh()[
patchi].name()
4894 <<
" owner cell " << own <<
" has refinement " 4896 <<
" (coupled) neighbour cell has refinement " 4919 forAll(syncPointLevel, pointi)
4921 if (pointLevel_[pointi] != syncPointLevel[pointi])
4924 <<
"PointLevel is not consistent across coupled patches." 4926 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4927 <<
" has level " << pointLevel_[pointi]
4928 <<
" whereas the coupled point has level " 4929 << syncPointLevel[pointi]
4937 if (maxPointDiff != -1)
4940 labelList maxPointLevel(mesh_.nPoints(), 0);
4942 forAll(maxPointLevel, pointi)
4944 const labelList& pCells = mesh_.pointCells(pointi);
4946 label& pLevel = maxPointLevel[pointi];
4950 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4966 label pointi = pointsToCheck[i];
4968 const labelList& pCells = mesh_.pointCells(pointi);
4972 label celli = pCells[i];
4976 mag(cellLevel_[celli]-maxPointLevel[pointi])
4983 <<
"Too big a difference between" 4984 <<
" point-connected cells." <<
nl 4986 <<
" cellLevel:" << cellLevel_[celli]
4987 <<
" uses point:" << pointi
4988 <<
" coord:" << mesh_.points()[pointi]
4989 <<
" which is also used by a cell with level:" 4990 << maxPointLevel[pointi]
5065 if (cellShapesPtr_.empty())
5069 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes." 5070 <<
" cellLevel:" << cellLevel_.
size()
5071 <<
" pointLevel:" << pointLevel_.
size()
5078 label nSplitHex = 0;
5079 label nUnrecognised = 0;
5081 forAll(cellLevel_, celli)
5083 if (meshShapes[celli].model().index() == 0)
5085 label level = cellLevel_[celli];
5089 bool haveQuads = matchHexShape
5099 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5110 Pout<<
"hexRef8::cellShapes() :" 5111 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl 5112 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5114 <<
" split-hex:" << nSplitHex <<
nl 5115 <<
" poly :" << nUnrecognised <<
nl 5119 return cellShapesPtr_();
5132 Pout<<
"hexRef8::getSplitPoints :" 5133 <<
" Calculating unrefineable points" <<
endl;
5140 <<
"Only call if constructed with history capability" 5148 labelList splitMaster(mesh_.nPoints(), -1);
5149 labelList splitMasterLevel(mesh_.nPoints(), 0);
5154 for (
label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5156 const labelList& pCells = mesh_.pointCells(pointi);
5158 if (pCells.
size() != 8)
5160 splitMaster[pointi] = -2;
5167 forAll(visibleCells, celli)
5169 const labelList& cPoints = mesh_.cellPoints(celli);
5171 if (visibleCells[celli] != -1 && history_.
parentIndex(celli) >= 0)
5178 label pointi = cPoints[i];
5180 label masterCelli = splitMaster[pointi];
5182 if (masterCelli == -1)
5189 splitMaster[pointi] = parentIndex;
5190 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5192 else if (masterCelli == -2)
5198 (masterCelli != parentIndex)
5199 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5204 splitMaster[pointi] = -2;
5213 label pointi = cPoints[i];
5215 splitMaster[pointi] = -2;
5223 label facei = mesh_.nInternalFaces();
5224 facei < mesh_.nFaces();
5228 const face& f = mesh_.faces()[facei];
5232 splitMaster[f[fp]] = -2;
5239 label nSplitPoints = 0;
5241 forAll(splitMaster, pointi)
5243 if (splitMaster[pointi] >= 0)
5252 forAll(splitMaster, pointi)
5254 if (splitMaster[pointi] >= 0)
5256 splitPoints[nSplitPoints++] = pointi;
5334 Pout<<
"hexRef8::consistentUnrefinement :" 5335 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5341 <<
"maxSet not implemented yet." 5353 forAll(pointsToUnrefine, i)
5355 label pointi = pointsToUnrefine[i];
5357 unrefinePoint.set(pointi);
5368 forAll(unrefinePoint, pointi)
5370 if (unrefinePoint.get(pointi))
5372 const labelList& pCells = mesh_.pointCells(pointi);
5376 unrefineCell.set(pCells[j]);
5389 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5391 label own = mesh_.faceOwner()[facei];
5392 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5394 label nei = mesh_.faceNeighbour()[facei];
5395 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5397 if (ownLevel < (neiLevel-1))
5404 unrefineCell.set(nei);
5415 if (unrefineCell.get(own) == 0)
5421 unrefineCell.unset(own);
5425 else if (neiLevel < (ownLevel-1))
5429 unrefineCell.set(own);
5433 if (unrefineCell.get(nei) == 0)
5439 unrefineCell.unset(nei);
5447 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
5451 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5453 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5461 label facei = i+mesh_.nInternalFaces();
5462 label own = mesh_.faceOwner()[facei];
5463 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5465 if (ownLevel < (neiLevel[i]-1))
5469 if (unrefineCell.get(own) == 0)
5475 unrefineCell.unset(own);
5479 else if (neiLevel[i] < (ownLevel-1))
5483 if (unrefineCell.get(own) == 1)
5489 unrefineCell.set(own);
5499 Pout<<
"hexRef8::consistentUnrefinement :" 5500 <<
" Changed " << nChanged
5501 <<
" refinement levels due to 2:1 conflicts." 5515 forAll(unrefinePoint, pointi)
5517 if (unrefinePoint.get(pointi))
5519 const labelList& pCells = mesh_.pointCells(pointi);
5523 if (!unrefineCell.get(pCells[j]))
5525 unrefinePoint.unset(pointi);
5537 forAll(unrefinePoint, pointi)
5539 if (unrefinePoint.get(pointi))
5548 forAll(unrefinePoint, pointi)
5550 if (unrefinePoint.get(pointi))
5552 newPointsToUnrefine[nSet++] = pointi;
5556 return newPointsToUnrefine;
5569 <<
"Only call if constructed with history capability" 5575 Pout<<
"hexRef8::setUnrefinement :" 5576 <<
" Checking initial mesh just to make sure" <<
endl;
5580 forAll(cellLevel_, celli)
5582 if (cellLevel_[celli] < 0)
5585 <<
"Illegal cell level " << cellLevel_[celli]
5586 <<
" for cell " << celli
5593 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5596 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.
size());
5598 forAll(splitPointLabels, i)
5600 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5604 cSet.insert(pCells[j]);
5609 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5611 << cSet.size() <<
" cells for unrefinement to" <<
nl 5613 <<
" cellSet " << cSet.objectPath()
5625 forAll(splitPointLabels, i)
5631 splitFaces.insert(pFaces[j]);
5645 if (facesToRemove.
size() != splitFaces.size())
5648 <<
"Ininitial set of split points to unrefine does not" 5649 <<
" seem to be consistent or not mid points of refined cells" 5658 forAll(splitPointLabels, i)
5660 label pointi = splitPointLabels[i];
5664 const labelList& pCells = mesh_.pointCells(pointi);
5667 if (pCells.
size() != 8)
5670 <<
"splitPoint " << pointi
5671 <<
" should have 8 cells using it. It has " << pCells
5684 label celli = pCells[j];
5686 label region = cellRegion[celli];
5691 <<
"Ininitial set of split points to unrefine does not" 5692 <<
" seem to be consistent or not mid points" 5693 <<
" of refined cells" <<
nl 5694 <<
"cell:" << celli <<
" on splitPoint " << pointi
5695 <<
" has no region to be merged into" 5699 if (masterCelli != cellRegionMaster[region])
5702 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5703 <<
" in region " << region
5704 <<
" has master:" << cellRegionMaster[region]
5705 <<
" which is not the lowest numbered cell" 5706 <<
" among the pointCells:" << pCells
5726 forAll(splitPointLabels, i)
5728 label pointi = splitPointLabels[i];
5730 const labelList& pCells = mesh_.pointCells(pointi);
5736 cellLevel_[pCells[j]]--;
5753 cellLevel_.
write(valid)
5754 && pointLevel_.
write(valid)
5755 && level0Edge_.
write(valid);
5759 writeOk = writeOk && history_.
write(valid);
fileCheckTypes
Enumeration defining the file checking options.
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
virtual const fileName & name() const
Return the name of the stream.
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.
const labelList & reversePointMap() const
Reverse point map.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
A class for handling file names.
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. Returns actual.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Class describing modification of a face.
unsigned int get(const label) const
Get value at index I.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
A face is a list of labels corresponding to mesh vertices.
void distributePointData(List< T > &lst) const
Distribute list of point data.
void setFaceInfo(const labelList &changedFaces, const List< Type > &changedFacesInfo)
Set initial changed faces.
const double e
Elementary charge.
void operator()(label &x, const label y) const
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.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
bool active() const
Is there unrefinement history?
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.
const boolList & flipMap() const
Return face flip map.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
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...
bool headerOk()
Read and check header info.
void checkMesh() const
Debug: Check coupled mesh for correctness.
point centre(const pointField &) const
Centre point of face.
void updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
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...
labelList consistentSlowRefinement2(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck) const
Like consistentSlowRefinement but uses different meshWave.
Class containing data for point addition.
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
label otherVertex(const label a) const
Given one vertex, return the other.
Holds information regarding type of cell. Used in inside/outside determination in cellClassification...
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 consistentUnrefinement(const labelList &pointsToUnrefine, const bool maxSet) const
Given proposed.
void distributeCellData(List< T > &lst) const
Distribute list of cell data.
label size() const
Return number of elements in table.
scalar level0EdgeLength() const
Typical edge length between unrefined points.
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...
point centre(const pointField &) const
Return centre (centroid)
void subset(const labelList &pointMap, const labelList &faceMap, const labelList &cellMap)
Update local numbering for subsetted mesh.
vector normal(const pointField &) const
Vector normal; magnitude is equal to area of face.
label refinementCount() const
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)
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.
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.
face reverseFace() const
Return face with reverse direction.
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.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
void combineCells(const label masterCelli, const labelList &combinedCells)
Store combining 8 cells into master.
void set(const PackedList< 1 > &)
Set specified bits.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
labelList consistentRefinement(const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
virtual const fileName & name() const
Return the name of the stream.
const DynamicList< point > & points() const
Points. Shrunk after constructing mesh (or calling of compact())
void clear()
Clear all entries from table.
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]
labelList consistentSlowRefinement(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck, const label maxPointDiff, const labelList &pointsToCheck) const
Like consistentRefinement but slower:
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Pair< label > labelPair
Label pair.
static const label labelMax
List< label > labelList
A List of labels.
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 > &)
vector vec(const pointField &) const
Return the vector (end - start)
void reverse(UList< T > &, const label n)
virtual bool read()
Read object. If global number of visible cells > 0 becomes active.
const labelList & reverseCellMap() const
Reverse cell map.
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
static fileCheckTypes fileModificationChecking
Type of file modification checking.
prefixOSstream Perr(cerr, "Perr")
const labelList & pointMap() const
Old point map.
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.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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.
void updateMesh(const mapPolyMesh &)
Update numbering for mesh changes.
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)
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.
label faceLevel(const label facei) const
Gets level such that the face has four points <= level.
#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)
A cell is defined as a list of faces with extra functionality.
A collection of cell labels.
prefixOSstream Pout(cout, "Pout")
label parentIndex(const label celli) const
Get parent of cell.
A List with indirect addressing.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
void flip()
Flip the face in-place.
labelList getSplitPoints() const
Return the points at the centre of top-level split cells.
void setRefinement(const labelList &piercedFaces, const labelList &cellRegion, const labelList &cellRegionMaster, polyTopoChange &) const
Play commands into polyTopoChange to remove faces.
void unset(const PackedList< 1 > &)
Unset specified bits.
dimensioned< scalar > mag(const dimensioned< Type > &)
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.
void checkRefinementLevels(const label maxPointDiff, const labelList &pointsToCheck) const
Debug: Check 2:1 consistency across faces.
A subset of mesh faces organised as a primitive patch.
All refinement history. Used in unrefinement.
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.
virtual bool write(const bool valid=true) const
Write using setting from DB.
void clear()
Clear the addressed list, i.e. set the size to zero.
readOption readOpt() const
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
fileName objectPath() const
Return complete path + object name.
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.
const cellShapeList & cellShapes() const
Utility: get hexes as cell shapes.
static const label labelMin
label nOldPoints() const
Number of old points.
Transfers refinement levels such that slow transition between levels is maintained. Used in FaceCellWave.
label compatibleRemoves(const labelList &inPiercedFaces, labelList &cellRegion, labelList &cellRegionMaster, labelList &outPiercedFaces) const
Find faces including those with cells which have the same mastercell.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.