61bool Foam::medialAxisMeshMover::isMaxEdge
63 const List<PointData<vector>>& pointWallDist,
66 const bool disableWallEdges
76 scalar magV0(
mag(v0));
83 scalar magV1(
mag(v1));
117 if ((pointWallDist[
e[0]].data() & pointWallDist[
e[1]].data()) < minCos)
129 <<
" : Calculating distance to Medial Axis ..." <<
endl;
134 const labelList& meshPoints = pp.meshPoints();
141 const label nSmoothSurfaceNormals
143 meshRefinement::get<label>
146 "nSmoothSurfaceNormals",
154 scalar minMedialAxisAngle(
Zero);
157 !coeffDict.readCompat
159 "minMedialAxisAngle",
160 {{
"minMedianAxisAngle", 1712 }},
168 <<
"Entry '" <<
"minMedialAxisAngle"
169 <<
"' not found in dictionary " << coeffDict.name() <<
endl;
175 const scalar featureAngle
177 meshRefinement::get<scalar>(coeffDict,
"featureAngle", dryRun_)
181 const scalar slipFeatureAngle
183 coeffDict.getOrDefault<scalar>(
"slipFeatureAngle", (0.5*featureAngle))
187 const label nSmoothNormals
189 meshRefinement::get<label>(coeffDict,
"nSmoothNormals", dryRun_)
193 const label nMedialAxisIter = coeffDict.getOrDefault<label>
199 const bool disableWallEdges = coeffDict.getOrDefault<
bool>
224 const bitSet isPatchMasterPoint
232 const bitSet isPatchMasterEdge
247 fieldSmoother_.smoothPatchNormals
249 nSmoothSurfaceNormals,
261 List<PointData<vector>> pointWallDist(
mesh().
nPoints());
264 int dummyTrackData = 0;
270 List<PointData<vector>> wallInfo(meshPoints.size());
272 forAll(meshPoints, patchPointI)
274 label pointI = meshPoints[patchPointI];
275 wallInfo[patchPointI] = PointData<vector>
279 pointNormals[patchPointI]
284 List<PointData<vector>> edgeWallDist(
mesh().
nEdges());
285 PointEdgeWave<PointData<vector>> wallDistCalc
295 wallDistCalc.iterate(nMedialAxisIter);
299 wallDistCalc.nUnvisitedPoints(),
305 if (nMedialAxisIter > 0)
308 <<
" : Limited walk to " << nMedialAxisIter
309 <<
" steps. Not visited " << nUnvisit
311 <<
" points" <<
endl;
316 <<
"Walking did not visit all points." <<
nl
317 <<
" Did not visit " << nUnvisit
319 <<
" points. This is not necessarily a problem" <<
nl
320 <<
" and might be due to faceZones splitting of part"
321 <<
" of the domain." <<
nl <<
endl;
330 List<pointEdgePoint> pointMedialDist(
mesh().
nPoints());
331 List<pointEdgePoint> edgeMedialDist(
mesh().
nEdges());
334 DynamicList<pointEdgePoint> maxInfo(meshPoints.size());
335 DynamicList<label> maxPoints(meshPoints.size());
343 const edge&
e = edges[edgeI];
347 !pointWallDist[
e[0]].valid(dummyTrackData)
348 || !pointWallDist[
e[1]].valid(dummyTrackData)
354 label pointI =
e[ep];
356 if (!pointMedialDist[pointI].valid(dummyTrackData))
367 pointMedialDist[pointI] = maxInfo.last();
378 minMedialAxisAngleCos,
389 scalar eMag =
mag(eVec);
396 const point& origin0 = pointWallDist[
e[0]].origin();
398 const point& origin1 = pointWallDist[
e[1]].origin();
399 scalar dist0 = (
p0-origin0) & eVec;
400 scalar dist1 = (origin1-p1) & eVec;
401 scalar
s = 0.5*(dist1+eMag+dist0);
413 else if (
s >= dist0+eMag)
424 medialAxisPt =
p0+(
s-dist0)*eVec;
431 label pointI =
e[ep];
433 if (!pointMedialDist[pointI].valid(dummyTrackData))
444 pointMedialDist[pointI] = maxInfo.last();
461 const polyPatch& pp =
patches[patchI];
463 pointDisplacement().boundaryField()[patchI];
468 && !isA<emptyPolyPatch>(pp)
469 && !adaptPatches.found(patchI)
472 const labelList& meshPoints = pp.meshPoints();
479 if (pvf.fixesValue())
483 <<
" : Inserting all points on patch " << pp.name()
488 label pointI = meshPoints[i];
489 if (!pointMedialDist[pointI].valid(dummyTrackData))
491 maxPoints.append(pointI);
500 pointMedialDist[pointI] = maxInfo.last();
512 <<
" : Inserting points on patch " << pp.name()
513 <<
" if angle to nearest layer patch > "
514 << slipFeatureAngle <<
" degrees." <<
endl;
527 label pointI = meshPoints[i];
531 pointWallDist[pointI].valid(dummyTrackData)
532 && !pointMedialDist[pointI].valid(dummyTrackData)
538 -pointWallDist[pointI].data()
541 if (cosAngle > slipFeatureAngleCos)
546 maxPoints.append(pointI);
555 pointMedialDist[pointI] = maxInfo.last();
575 OBJstream str(
mesh().time().timePath()/
"medialSurfacePoints.obj");
580 "medialSurfacePoints",
585 <<
" : Writing estimated medial surface:" <<
nl <<
incrIndent
590 for (
const auto& info : maxInfo)
592 str.
write(info.origin());
599 PointEdgeWave<pointEdgePoint> medialDistCalc
610 medialDistCalc.iterate(2*nMedialAxisIter);
614 forAll(pointMedialDist, pointI)
616 if (pointMedialDist[pointI].valid(dummyTrackData))
620 pointMedialDist[pointI].distSqr()
622 medialVec_[pointI] = pointMedialDist[pointI].origin();
627 medialDist_[pointI] = 0.0;
628 medialVec_[pointI] =
point(1, 0, 0);
636 if (!pointWallDist[i].valid(dummyTrackData))
638 dispVec_[i] =
vector(1, 0, 0);
642 dispVec_[i] = pointWallDist[i].data();
647 fieldSmoother_.smoothNormals
657 forAll(medialRatio_, pointI)
659 if (!pointWallDist[pointI].valid(dummyTrackData))
661 medialRatio_[pointI] = 0.0;
665 scalar wDist2 = pointWallDist[pointI].distSqr();
666 scalar mDist = medialDist_[pointI];
668 if (wDist2 <
sqr(SMALL) && mDist < SMALL)
675 medialRatio_[pointI] = 0.0;
679 medialRatio_[pointI] = mDist / (
Foam::sqrt(wDist2) + mDist);
689 <<
indent <<
"ratio of medial distance to wall distance : "
690 << medialRatio_.
name() <<
nl
691 <<
indent <<
"distance to nearest medial axis : "
692 << medialDist_.
name() <<
nl
693 <<
indent <<
"nearest medial axis location : "
694 << medialVec_.
name() <<
nl
695 <<
indent <<
"normal at nearest wall : "
700 medialRatio_.write();
707bool Foam::medialAxisMeshMover::unmarkExtrusion
709 const label patchPointI,
711 List<snappyLayerDriver::extrudeMode>& extrudeStatus
717 patchDisp[patchPointI] =
Zero;
723 patchDisp[patchPointI] =
Zero;
731void Foam::medialAxisMeshMover::syncPatchDisplacement
735 List<snappyLayerDriver::extrudeMode>& extrudeStatus
739 const labelList& meshPoints = pp.meshPoints();
741 label nChangedTotal = 0;
753 minMagSqrEqOp<vector>(),
760 if (
mag(patchDisp[i]) < minThickness[i])
762 if (unmarkExtrusion(i, patchDisp, extrudeStatus))
833 nChangedTotal += nChanged;
849void Foam::medialAxisMeshMover::
850handleFeatureAngleLayerTerminations
853 const bitSet& isPatchMasterPoint,
855 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
865 boolList extrudedFaces(pp.size(),
true);
867 forAll(pp.localFaces(), faceI)
869 const face&
f = pp.localFaces()[faceI];
875 extrudedFaces[faceI] =
false;
890 List<List<point>> edgeFaceNormals(pp.nEdges());
891 List<List<bool>> edgeFaceExtrude(pp.nEdges());
898 const labelList& eFaces = edgeFaces[edgeI];
900 edgeFaceNormals[edgeI].
setSize(eFaces.size());
901 edgeFaceExtrude[edgeI].setSize(eFaces.size());
904 label faceI = eFaces[i];
906 edgeFaceExtrude[edgeI][i] = extrudedFaces[faceI];
915 ListOps::appendEqOp<point>(),
924 ListOps::appendEqOp<bool>(),
929 forAll(edgeFaceNormals, edgeI)
931 const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
932 const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
934 if (eFaceNormals.size() == 2)
936 const edge&
e = pp.edges()[edgeI];
946 if (!eFaceExtrude[0] || !eFaceExtrude[1])
948 const vector& n0 = eFaceNormals[0];
949 const vector& n1 = eFaceNormals[1];
951 if ((n0 & n1) < minCos)
953 if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
955 if (isPatchMasterPoint[v0])
960 if (unmarkExtrusion(v1, patchDisp, extrudeStatus))
962 if (isPatchMasterPoint[v1])
982void Foam::medialAxisMeshMover::findIsolatedRegions
984 const scalar minCosLayerTermination,
985 const bool detectExtrusionIsland,
986 const bitSet& isPatchMasterPoint,
987 const bitSet& isPatchMasterEdge,
990 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
996 const labelList& meshPoints = pp.meshPoints();
999 Info<< typeName <<
" : Removing isolated regions ..." <<
nl
1000 <<
indent <<
"- if partially extruded faces make angle < "
1002 if (detectExtrusionIsland)
1004 Info<<
indent <<
"- if exclusively surrounded by non-extruded points"
1009 Info<<
indent <<
"- if exclusively surrounded by non-extruded faces"
1014 label nPointCounter = 0;
1020 if (minCosLayerTermination > -1)
1022 handleFeatureAngleLayerTerminations
1024 minCosLayerTermination,
1033 syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
1044 boolList keptPoints(pp.nPoints(),
false);
1046 if (detectExtrusionIsland)
1055 const face&
f = pp.localFaces()[faceI];
1061 if (islandPoint[faceI] == -1)
1064 islandPoint[faceI] =
f[fp];
1066 else if (islandPoint[faceI] != -2)
1069 islandPoint[faceI] = -2;
1081 forAll(pointFaces, patchPointI)
1090 if (islandPoint[faceI] != patchPointI)
1092 keptPoints[patchPointI] =
true;
1105 boolList extrudedFaces(pp.size(),
true);
1106 forAll(pp.localFaces(), faceI)
1108 const face&
f = pp.localFaces()[faceI];
1113 extrudedFaces[faceI] =
false;
1121 forAll(keptPoints, patchPointI)
1128 if (extrudedFaces[faceI])
1130 keptPoints[patchPointI] =
true;
1149 forAll(keptPoints, patchPointI)
1151 if (!keptPoints[patchPointI])
1153 if (unmarkExtrusion(patchPointI, patchDisp, extrudeStatus))
1168 const edgeList& edges = pp.edges();
1178 if (isPatchMasterEdge[edgeI])
1180 const edge&
e = edges[edgeI];
1187 isolatedPoint[v0] += 1;
1191 isolatedPoint[v1] += 1;
1209 const face&
f = pp.localFaces()[faceI];
1210 bool failed =
false;
1213 if (isolatedPoint[
f[fp]] > 2)
1219 bool allPointsExtruded =
true;
1226 allPointsExtruded =
false;
1231 if (allPointsExtruded)
1252 reduce(nPointCounter, sumOp<label>());
1254 <<
" : Number of isolated points extrusion stopped : "<< nPointCounter
1270 adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
1271 adaptPatchPtr_(getPatch(
mesh(), adaptPatchIDs_)),
1277 pointDisplacement.time().
timeName(),
1278 pointDisplacement.db(),
1298 fieldSmoother_(
mesh()),
1304 pointDisplacement.time().
timeName(),
1305 pointDisplacement.db(),
1318 pointDisplacement.time().
timeName(),
1319 pointDisplacement.db(),
1332 pointDisplacement.time().
timeName(),
1333 pointDisplacement.db(),
1346 pointDisplacement.time().
timeName(),
1347 pointDisplacement.db(),
1368void Foam::medialAxisMeshMover::calculateDisplacement
1376 Info<< typeName <<
" : Smoothing using Medial Axis ..." <<
endl;
1386 const label nSmoothDisplacement =
1390 const scalar maxThicknessToMedialRatio =
1391 coeffDict.
get<scalar>(
"maxThicknessToMedialRatio");
1394 const scalar featureAngle = coeffDict.
get<scalar>(
"featureAngle");
1397 scalar layerTerminationAngle = coeffDict.
getOrDefault<scalar>
1399 "layerTerminationAngle",
1405 const label nSmoothPatchThickness =
1406 coeffDict.
get<label>(
"nSmoothThickness");
1409 const label nMedialAxisIter = coeffDict.
getOrDefault<label>
1416 const bool detectExtrusionIsland = coeffDict.
getOrDefault
1418 "detectExtrusionIsland",
1437 const bitSet isPatchMasterPoint
1445 const bitSet isPatchMasterEdge
1457 forAll(thickness, patchPointI)
1461 thickness[patchPointI] = 0.0;
1465 label numThicknessRatioExclude = 0;
1470 autoPtr<OBJstream> str;
1478 /
"thicknessRatioExcludePoints_"
1484 <<
" : Writing points with too large an extrusion distance to "
1485 << str().name() <<
endl;
1488 autoPtr<OBJstream> medialVecStr;
1496 /
"thicknessRatioExcludeMedialVec_"
1502 <<
" : Writing medial axis vectors on points with too large"
1503 <<
" an extrusion distance to " << medialVecStr().name() <<
endl;
1506 forAll(meshPoints, patchPointI)
1510 label pointI = meshPoints[patchPointI];
1514 scalar mDist = medialDist_[pointI];
1515 scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1526 thicknessRatio *= (
n & mVec);
1528 if (thicknessRatio > maxThicknessToMedialRatio)
1533 Pout<<
"truncating displacement at "
1535 <<
" from " << thickness[patchPointI]
1539 minThickness[patchPointI]
1540 +thickness[patchPointI]
1542 <<
" medial direction:" << mVec
1543 <<
" extrusion direction:" <<
n
1544 <<
" with thicknessRatio:" << thicknessRatio
1548 thickness[patchPointI] =
1549 0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1551 patchDisp[patchPointI] = thickness[patchPointI]*
n;
1553 if (isPatchMasterPoint[patchPointI])
1555 numThicknessRatioExclude++;
1561 str().write(
linePointRef(pt, pt+patchDisp[patchPointI]));
1566 medialVecStr().write
1579 reduce(numThicknessRatioExclude, sumOp<label>());
1580 Info<< typeName <<
" : Reducing layer thickness at "
1581 << numThicknessRatioExclude
1582 <<
" nodes where thickness to medial axis distance is large " <<
endl;
1589 minCosLayerTermination,
1590 detectExtrusionIsland,
1602 forAll(thickness, patchPointI)
1606 thickness[patchPointI] = 0.0;
1615 fieldSmoother_.minSmoothField
1617 nSmoothPatchThickness,
1627 int dummyTrackData = 0;
1629 List<PointData<scalar>> pointWallDist(
mesh().
nPoints());
1636 List<PointData<scalar>> edgeWallDist(
mesh().
nEdges());
1640 List<PointData<scalar>> wallInfo(meshPoints.
size());
1642 forAll(meshPoints, patchPointI)
1644 label pointI = meshPoints[patchPointI];
1645 wallPoints[patchPointI] = pointI;
1646 wallInfo[patchPointI] = PointData<scalar>
1650 thickness[patchPointI]
1655 PointEdgeWave<PointData<scalar>> wallDistCalc
1665 wallDistCalc.iterate(nMedialAxisIter);
1670 pointField& displacement = pointDisplacement_;
1672 forAll(displacement, pointI)
1674 if (!pointWallDist[pointI].valid(dummyTrackData))
1676 displacement[pointI] =
Zero;
1685 displacement[pointI] =
1686 -medialRatio_[pointI]
1687 *pointWallDist[pointI].data()
1694 if (nSmoothDisplacement > 0)
1696 bitSet isToBeSmoothed(displacement.size(),
false);
1700 if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
1702 isToBeSmoothed.set(i);
1706 fieldSmoother_.smoothLambdaMuDisplacement
1708 nSmoothDisplacement,
1718bool Foam::medialAxisMeshMover::shrinkMesh
1720 const dictionary& meshQualityDict,
1721 const label nAllowableErrors,
1726 const label nSnap = meshQualityDict.get<label>(
"nRelaxIter");
1731 meshMover_.setDisplacementPatchFields();
1733 Info<< typeName <<
" : Moving mesh ..." <<
endl;
1734 scalar oldErrorReduction = -1;
1736 bool meshOk =
false;
1738 for (label iter = 0; iter < 2*nSnap ; iter++)
1741 <<
" : Iteration " << iter <<
endl;
1745 <<
" : Displacement scaling for error reduction set to 0."
1747 oldErrorReduction = meshMover_.setErrorReduction(0.0);
1752 meshMover_.scaleMesh
1756 meshMover_.paramDict(),
1763 Info<< typeName <<
" : Successfully moved mesh" <<
endl;
1769 if (oldErrorReduction >= 0)
1771 meshMover_.setErrorReduction(oldErrorReduction);
1774 Info<< typeName <<
" : Finished moving mesh ..." <<
endl;
1783 const label nAllowableErrors,
1791 const word minThicknessName = moveDict.
get<
word>(
"minThicknessName");
1798 if (minThicknessName ==
"none")
1804 (minThicknessName ==
"none")
1812 pointDisplacement_.internalField(),
1821 forAll(extrudeStatus, pointI)
1823 if (
mag(patchDisp[pointI]) <= minThickness[pointI]+SMALL)
1831 calculateDisplacement(moveDict, minThickness, extrudeStatus, patchDisp);
1848 meshMover_.movePoints();
1851 meshMover_.correct();
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
virtual const fileName & name() const
Return the name of the stream.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void setSize(const label n)
Alias for resize()
void append(const T &val)
Append an element at the end of the list.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
A list of faces which address into the list of points.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
const labelList & meshPoints() const
Return labelList of mesh points in patch.
void size(const label n)
Older name for setAddressableSize.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
virtual const vectorField & pointNormals() const
Return point unit normals.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Virtual base class for mesh movers with externally provided displacement field giving the boundary co...
const polyMesh & mesh() const
label nTotalPoints() const noexcept
Return total number of points in decomposed mesh. Not.
virtual bool update()
Update the mesh for both mesh motion and topology change.
@ REGEX
Regular expression.
static bitSet getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
static bitSet getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
void movePoints()
Update for new mesh geometry.
const Type & lookupObject(const word &name, const bool recursive=false) const
static const complex rootMax
complex (ROOTVGREAT, ROOTVGREAT)
static const complex max
complex (VGREAT,VGREAT)
Mesh representing a set of points created from polyMesh.
Mesh consisting of general polyhedral cells.
const globalMeshData & globalData() const
Return parallel info.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const pointField & points() const
Return raw points.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
@ NOEXTRUDE
Do not extrude. No layers added.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
const volScalarField & p0
const polyBoundaryMesh & patches
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimless
Dimensionless.
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
line< point, const point & > linePointRef
A line using referred points.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ostream & indent(Ostream &os)
Indent stream.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
List< labelList > labelListList
A List of labelList.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Field< vector > vectorField
Specialisation of Field<T> for vector.
pointPatchField< vector > pointPatchVectorField
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< edge > edgeList
A List of edges.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
#define forAll(list, i)
Loop across all elements in list.
Unit conversion functions.