60 DynamicList<label> adaptPatchIDs;
63 const pointPatchField<vector>& patchFld =
66 if (
isA<valuePointPatchField<vector>>(patchFld))
68 if (
isA<zeroFixedValuePointPatchField<vector>>(patchFld))
75 adaptPatchIDs.append(
patchi);
85 Foam::medialAxisMeshMover::getPatch
91 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
98 const polyPatch& pp =
patches[patchIDs[i]];
109 const polyPatch& pp =
patches[patchIDs[i]];
111 label meshFacei = pp.start();
115 addressing[nFaces++] = meshFacei++;
119 return autoPtr<indirectPrimitivePatch>
123 IndirectList<face>(mesh.faces(), addressing),
130 void Foam::medialAxisMeshMover::smoothPatchNormals
132 const label nSmoothDisp,
133 const PackedBoolList& isPatchMasterPoint,
134 const PackedBoolList& isPatchMasterEdge,
140 const labelList& meshPoints = pp.meshPoints();
143 Info<< typeName <<
" : Smoothing normals ..." <<
endl;
159 for (
label iter = 0; iter < nSmoothDisp; iter++)
174 if ((iter % 10) == 0)
181 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
189 normals[pointi] =
average[pointi];
190 normals[pointi] /=
mag(normals[pointi]) + vSmall;
197 void Foam::medialAxisMeshMover::smoothNormals
199 const label nSmoothDisp,
200 const PackedBoolList& isMeshMasterPoint,
201 const PackedBoolList& isMeshMasterEdge,
208 <<
" : Smoothing normals in interior ..." <<
endl;
210 const edgeList& edges = mesh().edges();
213 PackedBoolList isFixedPoint(mesh().
nPoints());
218 label meshPointi = fixedPoints[i];
219 isFixedPoint.set(meshPointi, 1);
245 for (
label iter = 0; iter < nSmoothDisp; iter++)
260 if ((iter % 10) == 0)
267 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
274 if (isFixedPoint.get(pointi) == 0)
278 normals[pointi] =
average[pointi];
279 normals[pointi] /=
mag(normals[pointi]) + vSmall;
288 bool Foam::medialAxisMeshMover::isMaxEdge
290 const List<pointData>& pointWallDist,
299 const edge&
e = mesh().edges()[edgei];
302 scalar magV0(
mag(v0));
310 scalar magV1(
mag(v1));
335 if ((pointWallDist[
e[0]].v() & pointWallDist[
e[1]].v()) < minCos)
346 void Foam::medialAxisMeshMover::update(
const dictionary& coeffDict)
349 <<
" : Calculating distance to Medial Axis ..." <<
endl;
354 const labelList& meshPoints = pp.meshPoints();
361 const label nSmoothSurfaceNormals =
362 coeffDict.lookup<
label>(
"nSmoothSurfaceNormals");
365 word angleKey =
"minMedialAxisAngle";
366 if (!coeffDict.found(angleKey))
369 angleKey =
"minMedianAxisAngle";
373 degToRad(coeffDict.lookup<scalar>(angleKey))
377 const scalar featureAngle = coeffDict.lookup<scalar>(
"featureAngle");
380 const scalar slipFeatureAngle =
382 coeffDict.found(
"slipFeatureAngle")
383 ? coeffDict.lookup<scalar>(
"slipFeatureAngle")
388 const label nSmoothNormals =
389 coeffDict.lookup<
label>(
"nSmoothNormals");
392 const label nMedialAxisIter = coeffDict.lookupOrDefault<
label>
395 mesh().globalData().nTotalPoints()
417 const PackedBoolList isPatchMasterPoint
425 const PackedBoolList isPatchMasterEdge
442 nSmoothSurfaceNormals,
453 List<pointData> pointWallDist(mesh().
nPoints());
456 int dummyTrackData = 0;
462 List<pointData> wallInfo(meshPoints.size());
464 forAll(meshPoints, patchPointi)
466 label pointi = meshPoints[patchPointi];
467 wallInfo[patchPointi] = pointData
472 pointNormals[patchPointi]
477 List<pointData> edgeWallDist(mesh().nEdges());
478 PointEdgeWave<pointData> wallDistCalc
488 wallDistCalc.iterate(nMedialAxisIter);
492 wallDistCalc.getUnsetPoints(),
498 if (nMedialAxisIter > 0)
501 <<
" : Limited walk to " << nMedialAxisIter
502 <<
" steps. Not visited " << nUnvisit
503 <<
" out of " << mesh().globalData().nTotalPoints()
504 <<
" points" <<
endl;
509 <<
"Walking did not visit all points." <<
nl
510 <<
" Did not visit " << nUnvisit
511 <<
" out of " << mesh().globalData().nTotalPoints()
512 <<
" points. This is not necessarily a problem" <<
nl
513 <<
" and might be due to faceZones splitting of part"
514 <<
" of the domain." <<
nl <<
endl;
523 List<pointData> pointMedialDist(mesh().
nPoints());
524 List<pointData> edgeMedialDist(mesh().nEdges());
527 DynamicList<pointData> maxInfo(meshPoints.size());
528 DynamicList<label> maxPoints(meshPoints.size());
532 const edgeList& edges = mesh().edges();
536 const edge&
e = edges[edgei];
540 !pointWallDist[
e[0]].
valid(dummyTrackData)
541 || !pointWallDist[
e[1]].
valid(dummyTrackData)
549 if (!pointMedialDist[pointi].
valid(dummyTrackData))
551 maxPoints.append(pointi);
562 pointMedialDist[pointi] = maxInfo.last();
567 else if (isMaxEdge(pointWallDist, edgei, minMedialAxisAngleCos))
575 scalar eMag =
mag(eVec);
583 scalar dist0 = (p0-pointWallDist[
e[0]].origin()) & eVec;
584 scalar dist1 = (pointWallDist[
e[1]].origin()-p1) & eVec;
585 scalar
s = 0.5*(dist1+eMag+dist0);
592 else if (
s >= dist0+eMag)
598 medialAxisPt = p0+(
s-dist0)*eVec;
605 if (!pointMedialDist[pointi].
valid(dummyTrackData))
607 maxPoints.append(pointi);
618 pointMedialDist[pointi] = maxInfo.last();
627 const polyBoundaryMesh&
patches = mesh().boundaryMesh();
636 pointDisplacement().boundaryField()[
patchi];
641 && !isA<emptyPolyPatch>(pp)
642 && !adaptPatches.found(
patchi)
645 const labelList& meshPoints = pp.meshPoints();
652 if (pvf.fixesValue())
656 <<
" : Inserting all points on patch " << pp.name()
661 label pointi = meshPoints[i];
662 if (!pointMedialDist[pointi].
valid(dummyTrackData))
664 maxPoints.append(pointi);
675 pointMedialDist[pointi] = maxInfo.last();
687 <<
" : Inserting points on patch " << pp.name()
688 <<
" if angle to nearest layer patch > "
689 << slipFeatureAngle <<
" degrees." <<
endl;
702 label pointi = meshPoints[i];
706 pointWallDist[pointi].
valid(dummyTrackData)
707 && !pointMedialDist[pointi].
valid(dummyTrackData)
713 -pointWallDist[pointi].v()
716 if (cosAngle > slipFeatureAngleCos)
721 maxPoints.append(pointi);
732 pointMedialDist[pointi] = maxInfo.last();
749 PointEdgeWave<pointData> medialDistCalc
760 medialDistCalc.iterate(2*nMedialAxisIter);
764 forAll(pointMedialDist, pointi)
766 if (pointMedialDist[pointi].
valid(dummyTrackData))
770 pointMedialDist[pointi].distSqr()
772 medialVec_[pointi] = pointMedialDist[pointi].origin();
777 medialDist_[pointi] = 0.0;
778 medialVec_[pointi] =
point(1, 0, 0);
786 if (!pointWallDist[i].
valid(dummyTrackData))
788 dispVec_[i] =
vector(1, 0, 0);
792 dispVec_[i] = pointWallDist[i].v();
807 forAll(medialRatio_, pointi)
809 if (!pointWallDist[pointi].
valid(dummyTrackData))
811 medialRatio_[pointi] = 0.0;
815 scalar wDist2 = pointWallDist[pointi].distSqr();
816 scalar mDist = medialDist_[pointi];
818 if (wDist2 <
sqr(small) && mDist < small)
825 medialRatio_[pointi] = 0.0;
829 medialRatio_[pointi] = mDist / (
Foam::sqrt(wDist2) + mDist);
838 <<
" : Writing medial axis fields:" <<
nl
840 <<
"ratio of medial distance to wall distance : "
841 << medialRatio_.
name() <<
nl
842 <<
"distance to nearest medial axis : "
843 << medialDist_.name() <<
nl
844 <<
"nearest medial axis location : "
845 << medialVec_.name() <<
nl
846 <<
"normal at nearest wall : "
847 << dispVec_.name() <<
nl
852 medialRatio_.
write();
859 bool Foam::medialAxisMeshMover::unmarkExtrusion
861 const label patchPointi,
863 List<snappyLayerDriver::extrudeMode>& extrudeStatus
869 patchDisp[patchPointi] =
Zero;
875 patchDisp[patchPointi] =
Zero;
885 void Foam::medialAxisMeshMover::syncPatchDisplacement
889 List<snappyLayerDriver::extrudeMode>& extrudeStatus
893 const labelList& meshPoints = pp.meshPoints();
905 minMagSqrEqOp<vector>(),
912 if (
mag(patchDisp[i]) < minThickness[i])
914 if (unmarkExtrusion(i, patchDisp, extrudeStatus))
929 void Foam::medialAxisMeshMover::minSmoothField
931 const label nSmoothDisp,
932 const PackedBoolList& isPatchMasterPoint,
933 const PackedBoolList& isPatchMasterEdge,
940 const labelList& meshPoints = pp.meshPoints();
955 Info<< typeName <<
" : Smoothing field ..." <<
endl;
957 for (
label iter = 0; iter < nSmoothDisp; iter++)
981 average[pointi] < field[pointi]
982 &&
average[pointi] >= fieldMin[pointi]
985 field[pointi] =
average[pointi];
990 if ((iter % 10) == 0)
997 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
1004 void Foam::medialAxisMeshMover::
1005 handleFeatureAngleLayerTerminations
1007 const scalar minCos,
1008 const PackedBoolList& isPatchMasterPoint,
1010 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
1012 label& nPointCounter
1020 boolList extrudedFaces(pp.size(),
true);
1022 forAll(pp.localFaces(), facei)
1024 const face&
f = pp.localFaces()[facei];
1030 extrudedFaces[facei] =
false;
1045 List<List<point>> edgeFaceNormals(pp.nEdges());
1046 List<List<bool>> edgeFaceExtrude(pp.nEdges());
1049 const vectorField& faceNormals = pp.faceNormals();
1053 const labelList& eFaces = edgeFaces[edgei];
1055 edgeFaceNormals[edgei].
setSize(eFaces.size());
1056 edgeFaceExtrude[edgei].setSize(eFaces.size());
1059 label facei = eFaces[i];
1060 edgeFaceNormals[edgei][i] = faceNormals[facei];
1061 edgeFaceExtrude[edgei][i] = extrudedFaces[facei];
1070 globalMeshData::ListPlusEqOp<List<point>>(),
1079 globalMeshData::ListPlusEqOp<List<bool>>(),
1084 forAll(edgeFaceNormals, edgei)
1086 const List<point>& eFaceNormals = edgeFaceNormals[edgei];
1087 const List<bool>& eFaceExtrude = edgeFaceExtrude[edgei];
1089 if (eFaceNormals.size() == 2)
1091 const edge&
e = pp.edges()[edgei];
1101 if (!eFaceExtrude[0] || !eFaceExtrude[1])
1103 const vector& n0 = eFaceNormals[0];
1104 const vector& n1 = eFaceNormals[1];
1106 if ((n0 & n1) < minCos)
1108 if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
1110 if (isPatchMasterPoint[v0])
1115 if (unmarkExtrusion(v1, patchDisp, extrudeStatus))
1117 if (isPatchMasterPoint[v1])
1136 void Foam::medialAxisMeshMover::findIsolatedRegions
1138 const scalar minCosLayerTermination,
1139 const bool detectExtrusionIsland,
1140 const PackedBoolList& isPatchMasterPoint,
1141 const PackedBoolList& isPatchMasterEdge,
1144 List<snappyLayerDriver::extrudeMode>& extrudeStatus,
1150 const labelList& meshPoints = pp.meshPoints();
1152 Info<< typeName <<
" : Removing isolated regions ..." <<
endl;
1155 label nPointCounter = 0;
1158 autoPtr<OBJstream> str;
1165 mesh().time().path()
1166 /
"islandExcludePoints_"
1167 + mesh().time().
name()
1172 <<
" : Writing points surrounded by non-extruded points to "
1173 << str().name() <<
endl;
1180 handleFeatureAngleLayerTerminations
1182 minCosLayerTermination,
1191 syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
1202 boolList keptPoints(pp.nPoints(),
false);
1204 if (detectExtrusionIsland)
1213 const face&
f = pp.localFaces()[facei];
1219 if (islandPoint[facei] == -1)
1222 islandPoint[facei] =
f[fp];
1224 else if (islandPoint[facei] != -2)
1227 islandPoint[facei] = -2;
1239 forAll(pointFaces, patchPointi)
1248 if (islandPoint[facei] != patchPointi)
1250 keptPoints[patchPointi] =
true;
1263 boolList extrudedFaces(pp.size(),
true);
1264 forAll(pp.localFaces(), facei)
1266 const face&
f = pp.localFaces()[facei];
1271 extrudedFaces[facei] =
false;
1279 forAll(keptPoints, patchPointi)
1286 if (extrudedFaces[facei])
1288 keptPoints[patchPointi] =
true;
1307 forAll(keptPoints, patchPointi)
1309 if (!keptPoints[patchPointi])
1311 if (unmarkExtrusion(patchPointi, patchDisp, extrudeStatus))
1318 str().write(pp.points()[meshPoints[patchPointi]]);
1331 const edgeList& edges = pp.edges();
1337 labelList isolatedPoint(pp.nPoints(),0);
1341 if (isPatchMasterEdge[edgei])
1343 const edge&
e = edges[edgei];
1350 isolatedPoint[v0] += 1;
1354 isolatedPoint[v1] += 1;
1372 const face&
f = pp.localFaces()[facei];
1373 bool failed =
false;
1376 if (isolatedPoint[
f[fp]] > 2)
1382 bool allPointsExtruded =
true;
1389 allPointsExtruded =
false;
1394 if (allPointsExtruded)
1412 str().write(pp.points()[meshPoints[
f[fp]]]);
1420 reduce(nPointCounter, sumOp<label>());
1422 <<
" : Number of isolated points extrusion stopped : "<< nPointCounter
1427 void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
1429 const label nSmoothDisp,
1430 const PackedBoolList& isMeshMasterPoint,
1431 const PackedBoolList& isMeshMasterEdge,
1435 const edgeList& edges = mesh().edges();
1454 Info<< typeName <<
" : Smoothing displacement ..." <<
endl;
1456 const scalar
lambda = 0.33;
1457 const scalar
mu = -0.34;
1461 for (
label iter = 0; iter < nSmoothDisp; iter++)
1477 if (medialRatio_[i] > small && medialRatio_[i] < 1-small)
1498 if (medialRatio_[i] > small && medialRatio_[i] < 1-small)
1500 displacement[i] = (1-
mu)*displacement[i]+
mu*
average[i];
1506 if ((iter % 10) == 0)
1513 Info<<
" Iteration " << iter <<
" residual " << resid <<
endl;
1529 adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
1530 adaptPatchPtr_(getPatch(mesh(), adaptPatchIDs_)),
1536 pointDisplacement.time().
name(),
1537 pointDisplacement.db(),
1544 oldPoints_(mesh().
points()),
1561 pointDisplacement.time().
name(),
1562 pointDisplacement.db(),
1575 pointDisplacement.time().
name(),
1576 pointDisplacement.db(),
1589 pointDisplacement.time().
name(),
1590 pointDisplacement.db(),
1603 pointDisplacement.time().
name(),
1604 pointDisplacement.db(),
1625 void Foam::medialAxisMeshMover::calculateDisplacement
1633 Info<< typeName <<
" : Smoothing using Medial Axis ..." <<
endl;
1645 "nSmoothDisplacement",
1650 const scalar maxThicknessToMedialRatio =
1651 coeffDict.
lookup<scalar>(
"maxThicknessToMedialRatio");
1654 const scalar featureAngle = coeffDict.
lookup<scalar>(
"featureAngle");
1657 const scalar minCosLayerTermination =
Foam::cos
1663 const label nSmoothPatchThickness =
1670 mesh().globalData().nTotalPoints()
1676 "detectExtrusionIsland",
1715 thickness =
mag(patchDisp);
1717 forAll(thickness, patchPointi)
1721 thickness[patchPointi] = 0.0;
1725 label numThicknessRatioExclude = 0;
1730 autoPtr<OBJstream> str;
1737 mesh().time().path()
1738 /
"thicknessRatioExcludePoints_"
1739 + mesh().time().
name()
1744 <<
" : Writing points with too large an extrusion distance to "
1745 << str().name() <<
endl;
1748 autoPtr<OBJstream> medialVecStr;
1755 mesh().time().path()
1756 /
"thicknessRatioExcludeMedialVec_"
1757 + mesh().time().
name()
1762 <<
" : Writing medial axis vectors on points with too large"
1763 <<
" an extrusion distance to " << medialVecStr().name() <<
endl;
1766 forAll(meshPoints, patchPointi)
1770 label pointi = meshPoints[patchPointi];
1774 scalar mDist = medialDist_[pointi];
1775 scalar thicknessRatio = thickness[patchPointi]/(mDist+vSmall);
1780 patchDisp[patchPointi]
1781 / (
mag(patchDisp[patchPointi]) + vSmall);
1782 vector mVec = mesh().points()[pointi]-medialVec_[pointi];
1783 mVec /=
mag(mVec)+vSmall;
1784 thicknessRatio *= (
n&mVec);
1786 if (thicknessRatio > maxThicknessToMedialRatio)
1791 Pout<<
"truncating displacement at "
1792 << mesh().points()[pointi]
1793 <<
" from " << thickness[patchPointi]
1797 minThickness[patchPointi]
1798 +thickness[patchPointi]
1800 <<
" medial direction:" << mVec
1801 <<
" extrusion direction:" <<
n
1802 <<
" with thicknessRatio:" << thicknessRatio
1806 thickness[patchPointi] =
1807 0.5*(minThickness[patchPointi]+thickness[patchPointi]);
1809 patchDisp[patchPointi] = thickness[patchPointi]*
n;
1811 if (isPatchMasterPoint[patchPointi])
1813 numThicknessRatioExclude++;
1818 const point& pt = mesh().points()[pointi];
1819 str().write(
linePointRef(pt, pt+patchDisp[patchPointi]));
1821 if (medialVecStr.valid())
1823 const point& pt = mesh().points()[pointi];
1824 medialVecStr().write
1837 reduce(numThicknessRatioExclude, sumOp<label>());
1838 Info<< typeName <<
" : Reducing layer thickness at "
1839 << numThicknessRatioExclude
1840 <<
" nodes where thickness to medial axis distance is large " <<
endl;
1847 minCosLayerTermination,
1848 detectExtrusionIsland,
1860 forAll(thickness, patchPointi)
1864 thickness[patchPointi] = 0.0;
1872 nSmoothPatchThickness,
1882 int dummyTrackData = 0;
1884 List<pointData> pointWallDist(mesh().
nPoints());
1891 List<pointData> edgeWallDist(mesh().nEdges());
1895 List<pointData> wallInfo(meshPoints.
size());
1897 forAll(meshPoints, patchPointi)
1899 label pointi = meshPoints[patchPointi];
1900 wallPoints[patchPointi] = pointi;
1901 wallInfo[patchPointi] = pointData
1905 thickness[patchPointi],
1911 PointEdgeWave<pointData> wallDistCalc
1921 wallDistCalc.iterate(nMedialAxisIter);
1926 pointField& displacement = pointDisplacement_;
1928 forAll(displacement, pointi)
1930 if (!pointWallDist[pointi].
valid(dummyTrackData))
1932 displacement[pointi] =
Zero;
1941 displacement[pointi] =
1942 -medialRatio_[pointi]
1943 *pointWallDist[pointi].s()
1950 if (nSmoothDisplacement > 0)
1952 smoothLambdaMuDisplacement
1954 nSmoothDisplacement,
1963 bool Foam::medialAxisMeshMover::shrinkMesh
1965 const dictionary& meshQualityDict,
1966 const label nAllowableErrors,
1971 const label nSnap = meshQualityDict.lookup<
label>(
"nRelaxIter");
1975 meshMover_.setDisplacementPatchFields();
1977 Info<< typeName <<
" : Moving mesh ..." <<
endl;
1978 scalar oldErrorReduction = -1;
1980 bool meshOk =
false;
1982 for (
label iter = 0; iter < 2*nSnap ; iter++)
1985 <<
" : Iteration " << iter <<
endl;
1989 <<
" : Displacement scaling for error reduction set to 0."
1991 oldErrorReduction = meshMover_.setErrorReduction(0.0);
1996 meshMover_.scaleMesh
2000 meshMover_.paramDict(),
2007 Info<< typeName <<
" : Successfully moved mesh" <<
endl;
2013 if (oldErrorReduction >= 0)
2015 meshMover_.setErrorReduction(oldErrorReduction);
2018 Info<< typeName <<
" : Finished moving mesh ..." <<
endl;
2027 const label nAllowableErrors,
2035 const word minThicknessName =
word(moveDict.
lookup(
"minThicknessName"));
2040 movePoints(mesh().
points());
2050 if (minThicknessName ==
"none")
2056 (minThicknessName ==
"none")
2058 : mesh().lookupObject<
scalarField>(minThicknessName)
2064 pointDisplacement_.primitiveField(),
2073 forAll(extrudeStatus, pointi)
2075 if (
mag(patchDisp[pointi]) <= minThickness[pointi]+small)
2083 calculateDisplacement(moveDict, minThickness, extrudeStatus, patchDisp);
2100 adaptPatchPtr_().clearGeom();
2103 meshMover_.movePoints();
2106 meshMover_.correct();
#define forAll(list, i)
Loop across all elements in list.
Macros for easy insertion into run-time selection tables.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
virtual const fileName & name() const
Return the name of the stream.
void size(const label)
Override size to be inconsistent with allocated storage.
void setSize(const label)
Reset size of List.
virtual Ostream & write(const char)=0
Write character.
label nPoints() const
Return number of points supporting patch faces.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
label size() const
Return the number of elements in the UPtrList.
static const Form rootMax
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
A list of keyword definitions, which are a keyword followed by any number of values (e....
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Virtual base class for mesh movers with externally provided displacement field giving the boundary co...
virtual void movePoints(const pointField &)
Update local data for geometry changes.
static void calculateEdgeWeights(const polyMesh &mesh, const PackedBoolList &isMasterEdge, const labelList &meshPoints, const edgeList &edges, scalarField &edgeWeights, scalarField &invSumWeight)
Helper: calculate edge weights (1/length)
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
static PackedBoolList getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
static PackedBoolList getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
static void weightedSum(const polyMesh &mesh, const PackedBoolList &isMasterEdge, const labelList &meshPoints, const edgeList &edges, const scalarField &edgeWeights, const Field< Type > &data, Field< Type > &sum)
Helper: weighted sum (over all subset of mesh points) by.
Mesh representing a set of points created from polyMesh.
Mesh consisting of general polyhedral cells.
@ NOEXTRUDE
Do not extrude. No layers added.
A class for handling words, derived from string.
gmvFile<< "tracers "<< particles.size()<< nl;{ pointField positions(particles.size());label particlei=0;forAllConstIter(Cloud< passiveParticle >, particles, iter) { positions[particlei++]=iter().position(mesh);} for(i=0;i< pTraits< point >::nComponents;i++) { forAll(positions, particlei) { gmvFile<< component(positions[particlei], i)<< ' ';} gmvFile<< nl;}}forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
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.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const fvPatchList & patches
dimensionedScalar lambda(viscosity->lookup("lambda"))
#define WarningInFunction
Report a warning using Foam::Warning.
bool valid(const PtrList< ModelType > &l)
const dimensionedScalar mu
Atomic mass unit.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
line< point, const point & > linePointRef
Line using referred points.
word name(const bool)
Return a word representation of a bool.
const dimensionSet dimless
const dimensionSet dimLength
vectorField pointField
pointField is a vectorField.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
Ostream & incrIndent(Ostream &os)
Increment the indent level.
List< bool > boolList
Bool container classes.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Vector< scalar > vector
A scalar version of the templated Vector.
bool isA(const Type &t)
Check if a dynamic_cast to typeid is possible.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
List< labelList > labelListList
A List of labelList.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
PointField< vector > pointVectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
pointPatchField< vector > pointPatchVectorField
prefixOSstream Pout(cout, "Pout")
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar cos(const dimensionedScalar &ds)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
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, &mergedCyclicPolyPatch::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]
Unit conversion functions.