Go to the documentation of this file.
74 void Foam::interfaceTrackingFvMesh::initializeData()
78 const word fsPatchName(
motion().get<word>(
"fsPatchName"));
85 <<
"Patch name " << fsPatchName <<
" not found."
89 fsPatchIndex_ =
patch.index();
96 for (
const word& patchName : pointNormalsCorrectionPatches_)
103 <<
"Patch name '" << patchName
104 <<
"' for point normals correction does not exist"
113 if (!normalMotionDir_)
123 "nonReflectingFreeSurfacePatches",
124 nonReflectingFreeSurfacePatches_
129 void Foam::interfaceTrackingFvMesh::makeUs()
const
132 <<
"making free-surface velocity field" <<
nl;
137 <<
"free-surface velocity field already exists"
144 zeroGradientFaPatchVectorField::typeName
152 == wedgeFaPatch::typeName
155 patchFieldTypes[patchI] =
156 wedgeFaPatchVectorField::typeName;
160 label ngbPolyPatchID =
163 if (ngbPolyPatchID != -1)
168 == wallFvPatch::typeName
171 patchFieldTypes[patchI] =
172 slipFaPatchVectorField::typeName;
178 for (
const word& patchName : fixedFreeSurfacePatches_)
180 const label fixedPatchID =
183 if (fixedPatchID == -1)
186 <<
"Wrong faPatch name '" << patchName
187 <<
"' in the fixedFreeSurfacePatches list"
188 <<
" defined in the dynamicMeshDict dictionary"
192 label ngbPolyPatchID =
195 if (ngbPolyPatchID != -1)
200 == wallFvPatch::typeName
203 patchFieldTypes[fixedPatchID] =
204 fixedValueFaPatchVectorField::typeName;
224 for (
const word& patchName : fixedFreeSurfacePatches_)
228 if (fixedPatchID == -1)
231 <<
"Wrong faPatch name '" << patchName
232 <<
"' in the fixedFreeSurfacePatches list"
233 <<
" defined in the dynamicMeshDict dictionary"
237 label ngbPolyPatchID =
240 if (ngbPolyPatchID != -1)
245 == wallFvPatch::typeName
248 UsPtr_->boundaryFieldRef()[fixedPatchID] ==
Zero;
255 void Foam::interfaceTrackingFvMesh::makeFsNetPhi()
const
258 <<
"making free-surface net flux" <<
nl;
263 <<
"free-surface net flux already exists"
283 void Foam::interfaceTrackingFvMesh::makeControlPoints()
286 <<
"making control points" <<
nl;
288 if (controlPointsPtr_)
291 <<
"control points already exists"
305 Info<<
"Reading control points" <<
endl;
321 Info<<
"Creating new control points" <<
endl;
333 aMesh().areaCentres().internalField()
336 initializeControlPointsPosition();
341 void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
344 <<
"making motion points mask" <<
nl;
346 if (motionPointsMaskPtr_)
349 <<
"motion points mask already exists"
366 == processorFaPatch::typeName
372 forAll(patchPoints, pointI)
374 motionPointsMask()[patchPoints[pointI]] = -1;
380 for (
const word& patchName : fixedFreeSurfacePatches_)
384 if (fixedPatchID == -1)
387 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
388 <<
" defined in the dynamicMeshDict dictionary"
395 forAll(patchPoints, pointI)
397 motionPointsMask()[patchPoints[pointI]] = 0;
403 void Foam::interfaceTrackingFvMesh::makeDirections()
406 <<
"make displacement directions for points and control points" <<
nl;
408 if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_)
411 <<
"points, control points displacement directions already exist"
415 pointsDisplacementDirPtr_ =
422 facesDisplacementDirPtr_ =
429 if (!normalMotionDir())
431 if (
mag(motionDir_) < SMALL)
434 <<
"Zero motion direction"
438 facesDisplacementDir() = motionDir_;
439 pointsDisplacementDir() = motionDir_;
442 updateDisplacementDirections();
446 void Foam::interfaceTrackingFvMesh::makePhis()
449 <<
"making free-surface flux" <<
nl;
454 <<
"free-surface flux already exists"
473 void Foam::interfaceTrackingFvMesh::makeSurfactConc()
const
476 <<
"making free-surface surfactant concentration field" <<
nl;
481 <<
"free-surface surfactant concentration field already exists"
504 void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc()
const
507 <<
"making volume surfactant concentration field" <<
nl;
509 if (bulkSurfactConcPtr_)
512 <<
"volume surfactant concentration field already exists"
537 bulkSurfactConc = surfactant().bulkConc();
543 void Foam::interfaceTrackingFvMesh::makeSurfaceTension()
const
546 <<
"making surface tension field" <<
nl;
548 if (surfaceTensionPtr_)
551 <<
"surface tension field already exists"
565 sigma() + surfactant().dSigma(surfactantConcentration())/rho_
570 void Foam::interfaceTrackingFvMesh::makeSurfactant()
const
573 <<
"making surfactant properties" <<
nl;
578 <<
"surfactant properties already exists"
583 motion().
subDict(
"surfactantProperties");
589 void Foam::interfaceTrackingFvMesh::makeContactAngle()
592 <<
"making contact angle field" <<
nl;
594 if (contactAnglePtr_)
597 <<
"contact angle already exists"
612 Info<<
"Reading contact angle field" <<
endl;
631 void Foam::interfaceTrackingFvMesh::updateDisplacementDirections()
633 if (normalMotionDir())
641 if (contactAnglePtr_)
643 label ngbPolyPatchID =
646 if (ngbPolyPatchID != -1)
651 == wallFvPatch::typeName
660 .ngbPolyPatchPointNormals()
663 forAll(patchPoints, pointI)
665 pointsDisplacementDir()[patchPoints[pointI]] -=
669 & pointsDisplacementDir()[patchPoints[pointI]]
672 pointsDisplacementDir()[patchPoints[pointI]] /=
675 pointsDisplacementDir()
687 facesDisplacementDir() =
694 facesDisplacementDir()
695 *(facesDisplacementDir()&(controlPoints() - Cf))
701 void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition()
709 correctPointDisplacement(sweptVolCorr, displacement);
717 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
724 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
731 deltaH[faceI] = sweptVol[faceI]/
732 ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL);
734 if (
mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL)
741 <<
"Something is wrong with specified motion direction"
746 for (
const word& patchName : fixedFreeSurfacePatches_)
750 if (fixedPatchID == -1)
753 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
754 <<
" defined in the freeSurfaceProperties dictionary"
763 deltaH[eFaces[edgeI]] *= 2.0;
767 controlPoints() += facesDisplacementDir()*deltaH;
772 void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh()
774 Info<<
"Smoothing free surface mesh" <<
endl;
788 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
794 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
799 sweptVol/(faceArea & facesDisplacementDir())
802 for (
const word& patchName : fixedFreeSurfacePatches_)
806 if (fixedPatchID == -1)
809 <<
"Wrong faPatch name fixedFreeSurfacePatches list"
810 <<
" defined in the dynamicMeshDict dictionary"
819 deltaHf[eFaces[edgeI]] *= 2.0;
823 controlPoints() += facesDisplacementDir()*deltaHf;
825 displacement = pointDisplacement();
828 refCast<velocityMotionSolver>
836 fixedValuePointPatchVectorField& fsPatchPointMeshU =
837 refCast<fixedValuePointPatchVectorField>
852 void Foam::interfaceTrackingFvMesh::updateSurfaceFlux()
858 void Foam::interfaceTrackingFvMesh::updateSurfactantConcentration()
860 if (!pureFreeSurface())
870 +
fam::div(Phis(), surfactantConcentration())
873 surfactant().diffusion(),
874 surfactantConcentration()
878 if (surfactant().soluble())
894 zeroGradientFaPatchScalarField::typeName
898 bulkSurfactantConcentration().boundaryField()[fsPatchIndex()];
899 Cb.correctBoundaryConditions();
904 surfactant().adsorptionCoeff()*
Cb
905 + surfactant().adsorptionCoeff()
906 *surfactant().desorptionCoeff(),
907 surfactantConcentration()
909 - surfactant().adsorptionCoeff()
910 *
Cb*surfactant().saturatedConc();
918 sigma() + surfactant().dSigma(surfactantConcentration())/rho_;
920 if (
neg(
min(surfaceTension().internalField().
field())))
923 <<
"Surface tension is negative"
930 Foam::vector Foam::interfaceTrackingFvMesh::totalPressureForce()
const
940 return gSum(pressureForces);
944 Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce()
const
964 U().boundaryField()[fsPatchIndex()].
snGrad()
977 return gSum(viscousForces);
981 Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce()
const
991 if (pureFreeSurface())
1002 surfTensionForces = surfaceTension().internalField().field()*
K*S*
n;
1005 return gSum(surfTensionForces);
1009 Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber()
1013 if (pureFreeSurface())
1019 mesh().time().deltaT().value()/
1038 mesh().time().deltaT().value()/
1050 void Foam::interfaceTrackingFvMesh::updateProperties()
1055 "transportProperties"
1064 void Foam::interfaceTrackingFvMesh::correctPointDisplacement
1079 for (
const word& patchName : fixedFreeSurfacePatches_)
1093 label curFace = eFaces[edgeI];
1095 const labelList& curPoints = faces[curFace];
1097 forAll(curPoints, pointI)
1099 label curPoint = curPoints[pointI];
1100 label index = pLabels.find(curPoint);
1116 forAll(corrPoints, pointI)
1118 label curPoint = corrPoints[pointI];
1124 label curFace =
pFaces[curPoint][faceI];
1126 label index = eFaces.find(curFace);
1140 forAll(corrPoints, pointI)
1142 label curPoint = corrPoints[pointI];
1146 const labelList& curPointFaces = corrPointFaces[pointI];
1148 forAll(curPointFaces, faceI)
1150 const face& curFace = faces[curPointFaces[faceI]];
1152 label ptInFace = curFace.
which(curPoint);
1153 label next = curFace.
nextLabel(ptInFace);
1154 label prev = curFace.
prevLabel(ptInFace);
1158 const vector&
c = pointsDisplacementDir()[curPoint];
1160 curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^
b)&
c);
1163 curDisp /= curPointFaces.size();
1165 displacement[curPoint] =
1166 curDisp*pointsDisplacementDir()[curPoint];
1171 for (
const word& patchName : nonReflectingFreeSurfacePatches_)
1173 label nonReflectingPatchID =
1186 forAll(corrPoints, pointI)
1188 label curPoint = corrPoints[pointI];
1194 label curFace =
pFaces[curPoint][faceI];
1196 label index = eFaces.find(curFace);
1211 forAll(corrPoints, pointI)
1213 label curPoint = corrPoints[pointI];
1217 const labelList& curPointFaces = corrPointFaces[pointI];
1219 forAll(curPointFaces, faceI)
1221 const face& curFace = faces[curPointFaces[faceI]];
1223 label ptInFace = curFace.
which(curPoint);
1224 label next = curFace.
nextLabel(ptInFace);
1225 label prev = curFace.
prevLabel(ptInFace);
1231 if (corrPoints.find(next) == -1)
1246 vector c0 = displacement[p1];
1248 scalar V0 =
mag((
a0^b0)&c0)/2;
1250 scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0;
1252 if (corrPoints.find(prev) != -1)
1256 (
points[next] + displacement[next])
1258 const vector&
c = pointsDisplacementDir()[curPoint];
1260 curDisp += 2*DV/((a^
b)&
c);
1265 - (
points[prev] + displacement[prev]);
1267 const vector&
c = pointsDisplacementDir()[curPoint];
1269 curDisp += 2*DV/((a^
b)&
c);
1273 curDisp /= curPointFaces.size();
1275 displacement[curPoint] =
1276 curDisp*pointsDisplacementDir()[curPoint];
1282 void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals()
1293 if (contactAnglePtr_ && correctContactLineNormals())
1295 Info<<
"Correcting contact line normals" <<
endl;
1334 == processorFaPatch::typeName
1338 refCast<const processorFaPatch>
1346 forAll(patchPointLabels, pointI)
1348 label curPoint = patchPointLabels[pointI];
1355 const labelList& curPointEdges = pointEdges[curPoint];
1357 forAll(curPointEdges, edgeI)
1359 label curEdge = curPointEdges[edgeI];
1361 if (edgeFaces[curEdge].size() == 1)
1368 label index = curEdges.find(curEdge);
1375 != processorFaPatch::typeName
1392 vector t = edges[curEdge].vec(oldPoints);
1393 t /=
mag(t) + SMALL;
1396 pointFaces[curPoint];
1398 forAll(curPointFaces, fI)
1400 tangent.ref().field()[curPointFaces[fI]] = t;
1407 tangent.correctBoundaryConditions();
1412 label ngbPolyPatchID =
1415 if (ngbPolyPatchID != -1)
1420 == wallFvPatch::typeName
1428 - contactAnglePtr_->boundaryField()[patchI]
1443 rotationAxis /=
mag(rotationAxis) + SMALL;
1453 forAll(pointEdges, pointI)
1457 forAll(pointEdges[pointI], eI)
1461 + pointEdges[pointI][eI];
1463 vector e = edges[curEdge].vec(oldPoints);
1465 e *= (
e&rotationAxis[pointI])
1466 /
mag(
e&rotationAxis[pointI]);
1468 e /=
mag(
e) + SMALL;
1473 if (pointEdges[pointI].size() == 1)
1475 label curPoint = patchPoints[pointI];
1482 label procPatchID = -1;
1488 forAll(curPointEdges, edgeI)
1490 label curEdge = curPointEdges[edgeI];
1492 if (edgeFaces[curEdge].size() == 1)
1500 curEdges.find(curEdge);
1507 == processorFaPatch::typeName
1519 if (procPatchID != -1)
1522 tangent.boundaryField()[procPatchID]
1523 .patchNeighbourField()()[
edgeID];
1525 t *= (t&rotationAxis[pointI])
1526 /(
mag(t&rotationAxis[pointI]) + SMALL);
1528 t /=
mag(t) + SMALL;
1534 rotationAxis[pointI] = rotAx/(
mag(rotAx) + SMALL);
1538 ngbN = ngbN*
cos(rotAngle)
1539 + rotationAxis*(rotationAxis & ngbN)*(1 -
cos(rotAngle))
1540 + (rotationAxis^ngbN)*
sin(rotAngle);
1543 forAll(patchPoints, pointI)
1545 N[patchPoints[pointI]] -=
1546 ngbN[pointI]*(ngbN[pointI]&
N[patchPoints[pointI]]);
1548 N[patchPoints[pointI]] /=
1549 mag(
N[patchPoints[pointI]]) + SMALL;
1562 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1569 aMeshPtr_(
new faMesh(*
this)),
1571 fixedFreeSurfacePatches_
1573 motion().get<wordList>(
"fixedFreeSurfacePatches")
1575 nonReflectingFreeSurfacePatches_(),
1576 pointNormalsCorrectionPatches_
1578 motion().get<wordList>(
"pointNormalsCorrectionPatches")
1582 motion().get<bool>(
"normalMotionDir")
1585 smoothing_(motion().getOrDefault(
"smoothing",
false)),
1586 pureFreeSurface_(motion().getOrDefault(
"pureFreeSurface",
true)),
1587 rigidFreeSurface_(motion().getOrDefault(
"rigidFreeSurface",
false)),
1588 correctContactLineNormals_
1590 motion().getOrDefault(
"correctContactLineNormals",
false)
1596 controlPointsPtr_(
nullptr),
1597 motionPointsMaskPtr_(
nullptr),
1598 pointsDisplacementDirPtr_(
nullptr),
1599 facesDisplacementDirPtr_(
nullptr),
1600 fsNetPhiPtr_(
nullptr),
1602 surfactConcPtr_(
nullptr),
1603 bulkSurfactConcPtr_(
nullptr),
1604 surfaceTensionPtr_(
nullptr),
1605 surfactantPtr_(
nullptr),
1606 contactAnglePtr_(
nullptr)
1612 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1627 std::move(allOwner),
1628 std::move(allNeighbour),
1631 aMeshPtr_(
new faMesh(*
this)),
1633 fixedFreeSurfacePatches_
1635 motion().get<wordList>(
"fixedFreeSurfacePatches")
1637 nonReflectingFreeSurfacePatches_(),
1638 pointNormalsCorrectionPatches_
1640 motion().get<wordList>(
"pointNormalsCorrectionPatches")
1644 motion().get<bool>(
"normalMotionDir")
1647 smoothing_(motion().getOrDefault(
"smoothing",
false)),
1648 pureFreeSurface_(motion().getOrDefault(
"pureFreeSurface",
true)),
1653 controlPointsPtr_(
nullptr),
1654 motionPointsMaskPtr_(
nullptr),
1655 pointsDisplacementDirPtr_(
nullptr),
1656 facesDisplacementDirPtr_(
nullptr),
1657 fsNetPhiPtr_(
nullptr),
1659 surfactConcPtr_(
nullptr),
1660 bulkSurfactConcPtr_(
nullptr),
1661 surfaceTensionPtr_(
nullptr),
1662 surfactantPtr_(
nullptr),
1663 contactAnglePtr_(
nullptr)
1719 return *fsNetPhiPtr_;
1730 return *fsNetPhiPtr_;
1736 forAll(
Us().boundaryField(), patchI)
1740 Us().boundaryField()[patchI].
type()
1741 == calculatedFaPatchVectorField::typeName
1746 pUs =
Us().boundaryField()[patchI].patchInternalField();
1748 label ngbPolyPatchID =
1751 if (ngbPolyPatchID != -1)
1756 U().boundaryField()[ngbPolyPatchID].
type()
1757 == slipFvPatchVectorField::typeName
1761 U().boundaryField()[ngbPolyPatchID].
type()
1762 == symmetryFvPatchVectorField::typeName
1777 Us().correctBoundaryConditions();
1785 Us().ref().field() =
U().boundaryField()[fsPatchIndex()];
1793 correctUsBoundaryConditions();
1799 return *getObjectPtr<const volVectorField>(
"U");
1805 return *getObjectPtr<const volScalarField>(
"p");
1811 return *getObjectPtr<const surfaceScalarField>(
"phi");
1819 auto& SnGradU = tSnGradU.ref();
1826 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1833 gradUs -=
n*(
n & gradUs);
1843 if (!pureFreeSurface() &&
max(
nu) > SMALL)
1845 tangentialSurfaceTensionForce =
1846 surfaceTensionGrad()().internalField();
1850 tangentialSurfaceTensionForce/(
nu + SMALL)
1862 auto& SnGradUn = tSnGradUn.ref();
1867 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1880 auto& pressureJump = tPressureJump.ref();
1894 + 2.0*
nu*freeSurfaceSnGradUn();
1896 if (pureFreeSurface())
1902 pressureJump -= surfaceTension().internalField()*
K;
1905 return tPressureJump;
1911 if (!controlPointsPtr_)
1913 makeControlPoints();
1916 return *controlPointsPtr_;
1922 if (!motionPointsMaskPtr_)
1924 makeMotionPointsMask();
1927 return *motionPointsMaskPtr_;
1933 if (!pointsDisplacementDirPtr_)
1938 return *pointsDisplacementDirPtr_;
1944 if (!facesDisplacementDirPtr_)
1949 return *facesDisplacementDirPtr_;
1966 if (!surfactConcPtr_)
1971 return *surfactConcPtr_;
1978 if (!surfactConcPtr_)
1983 return *surfactConcPtr_;
1990 if (!bulkSurfactConcPtr_)
1992 makeBulkSurfactConc();
1995 return *bulkSurfactConcPtr_;
2002 if (!bulkSurfactConcPtr_)
2004 makeBulkSurfactConc();
2007 return *bulkSurfactConcPtr_;
2014 if (!surfaceTensionPtr_)
2016 makeSurfaceTension();
2019 return *surfaceTensionPtr_;
2026 if (!surfaceTensionPtr_)
2028 makeSurfaceTension();
2031 return *surfaceTensionPtr_;
2038 if (!surfactantPtr_)
2043 return *surfactantPtr_;
2054 "surfaceTensionGrad",
2066 auto&
grad = tgrad.ref();
2071 grad.correctBoundaryConditions();
2081 if (smoothing_ && !rigidFreeSurface_)
2083 smoothFreeSurfaceMesh();
2084 clearControlPoints();
2087 updateDisplacementDirections();
2091 Info<<
"Maximal capillary Courant number: "
2092 << maxCourantNumber() <<
endl;
2096 Info<<
"Free surface curvature: min = " <<
gMin(
K)
2102 if (!rigidFreeSurface_)
2118 Info<<
"Free surface continuity error : sum local = "
2119 <<
gSum(
mag(sweptVolCorr)) <<
", global = " <<
gSum(sweptVolCorr)
2123 fsNetPhi().ref().field() = sweptVolCorr;
2129 "ddt(" +
U().
name() +
')'
2167 <<
"Unsupported temporal differencing scheme : "
2177 sweptVolCorr/(Sf*(Nf & facesDisplacementDir()))
2180 for (
const word& patchName : fixedFreeSurfacePatches_)
2182 label fixedPatchID =
2185 if (fixedPatchID == -1)
2188 <<
"Wrong faPatch name '" << patchName
2189 <<
"'in the fixedFreeSurfacePatches list"
2190 <<
" defined in dynamicMeshDict dictionary"
2199 deltaHf[eFaces[edgeI]] *= 2.0;
2203 controlPoints() += facesDisplacementDir()*deltaHf;
2205 pointField displacement(pointDisplacement());
2206 correctPointDisplacement(sweptVolCorr, displacement);
2209 refCast<velocityMotionSolver>
2217 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2218 refCast<fixedValuePointPatchVectorField>
2226 fsPatchPointMeshU ==
2231 correctContactLinePointNormals();
2242 refCast<velocityMotionSolver>
2250 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2251 refCast<fixedValuePointPatchVectorField>
2259 fsPatchPointMeshU ==
2267 updateSurfactantConcentration();
2276 OFstream mps(
mesh().time().timePath()/
"freeSurface.vtk");
2281 mps <<
"# vtk DataFile Version 2.0" <<
nl
2284 <<
"DATASET POLYDATA" <<
nl
2285 <<
"POINTS " <<
points.size() <<
" float" <<
nl;
2293 mlpBuffer[counter++] = float(
points[i].
x());
2294 mlpBuffer[counter++] = float(
points[i].
y());
2295 mlpBuffer[counter++] = float(
points[i].z());
2300 mps << mlpBuffer[i] <<
' ';
2302 if (i > 0 && (i % 10) == 0)
2309 label nFaceVerts = 0;
2313 nFaceVerts += faces[faceI].size() + 1;
2320 const face&
f = faces[faceI];
2322 mlfBuffer[counter++] =
f.size();
2326 mlfBuffer[counter++] =
f[fpI];
2331 mps <<
"POLYGONS " << faces.size() <<
' ' << nFaceVerts <<
endl;
2335 mps << mlfBuffer[i] <<
' ';
2337 if (i > 0 && (i % 10) == 0)
2359 Info<<
"Writing free surface control point to " <<
name <<
endl;
2361 mps <<
"# vtk DataFile Version 2.0" <<
nl
2364 <<
"DATASET POLYDATA" <<
nl
2365 <<
"POINTS " << controlPoints().size() <<
" float" <<
nl;
2367 forAll(controlPoints(), pointI)
2369 mps << controlPoints()[pointI].x() <<
' '
2370 << controlPoints()[pointI].y() <<
' '
2371 << controlPoints()[pointI].z() <<
nl;
2375 mps <<
"VERTICES " << controlPoints().size() <<
' '
2376 << controlPoints().size()*2 <<
nl;
2378 forAll(controlPoints(), pointI)
2380 mps << 1 <<
' ' << pointI <<
nl;
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
List< label > labelList
A List of labels.
IOField< vector > vectorIOField
vectorField with IO.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const edgeVectorField & Le() const
Return edge length vectors.
faMesh & aMesh()
Return reference to finite area mesh.
const labelListList & pointFaces() const
Return point-face addressing.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Dimensionless.
volScalarField & bulkSurfactantConcentration()
Return volume surfactant concentration field.
const labelListList & edgeFaces() const
Return edge-face addressing.
A class for handling words, derived from Foam::string.
const faceList & faces() const
Return faces.
A class for handling file names.
A primitive field of type <T> with automated input and output.
Info<< "Reading field U\n"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);volScalarField rho(IOobject("rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), thermo.rho());volVectorField rhoU(IOobject("rhoU", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *U);volScalarField rhoE(IOobject("rhoE", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *(e+0.5 *magSqr(U)));surfaceScalarField pos(IOobject("pos", runTime.timeName(), mesh), mesh, dimensionedScalar("pos", dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar("neg", dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\n"<< endl;autoPtr< compressible::turbulenceModel > turbulence(compressible::turbulenceModel::New(rho, U, phi, thermo))
dimensionedScalar deltaT() const
Return time step.
List< Key > toc() const
The table of contents (the keys) in unsorted order.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
areaScalarField & surfactantConcentration()
Return free-surface surfactant concentration field.
A class for managing temporary objects.
const indirectPrimitivePatch & patch() const
Return constant reference to primitive patch.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVelocity
tmp< faMatrix< Type > > laplacian(const GeometricField< Type, faPatchField, areaMesh > &vf)
const dimensionSet dimDensity
virtual bool update()
Update the mesh for both mesh motion and topology change.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
Type gAverage(const FieldField< Field, Type > &f)
Template functions to aid in the implementation of demand driven data.
dimensionedScalar sin(const dimensionedScalar &ds)
bool correctPatchPointNormals(const label patchID) const
Whether point normals should be corrected for a patch.
static bool & parRun()
Test if this a parallel run, or allow modify access.
The interfaceTrackingFvMesh.
Abstract base class for geometry and/or topology changing fvMesh.
void writeVTK() const
Write VTK freeSurface mesh.
tmp< GeometricField< typename outerProduct< vector, Type >::type, faPatchField, areaMesh >> grad(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Unit conversion functions.
const edgeList & edges() const
Return edges.
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0)
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
const Type & value() const
Return const reference to value.
const vectorField & pointAreaNormals() const
Return point area normals.
Type gSum(const FieldField< Field, Type > &f)
Template specialisation for scalar faMatrix.
A simple single-phase transport model based on viscosityModel.
const Internal & internalField() const
Return a const-reference to the dimensioned internal field.
const dimensionSet dimForce
Virtual base class for velocity motion solver.
const motionSolver & motion() const
Return the motionSolver.
tmp< GeometricField< Type, faPatchField, areaMesh > > div(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
tmp< faMatrix< Type > > ddt(const GeometricField< Type, faPatchField, areaMesh > &vf)
#define forAll(list, i)
Loop across all elements in list.
Field< vector > vectorField
Specialisation of Field<T> for vector.
tmp< scalarField > freeSurfacePressureJump()
Return free surface pressure jump.
CGAL::Exact_predicates_exact_constructions_kernel K
GeometricField< vector, faPatchField, areaMesh > areaVectorField
void deleteDemandDrivenData(DataPtr &dataPtr)
Ostream & flush(Ostream &os)
Flush stream.
tmp< faMatrix< Type > > Sp(const areaScalarField &sp, const GeometricField< Type, faPatchField, areaMesh > &vf)
void updateUs()
Update free-surface velocity field.
virtual bool found(const label id) const
Has the given index?
tmp< areaVectorField > surfaceTensionGrad()
Return surface tension gradient.
const dimensionSet dimArea(sqr(dimLength))
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]=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]
label which(const label pointLabel) const
Find local index on face for the point label,.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
const volScalarField & p() const
Return constant reference to pressure field.
const DimensionedField< scalar, areaMesh > & S() const
Return face areas.
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
messageStream Info
Information stream (uses stdout - output is on the master only)
#define DebugInFunction
Report an information message using Foam::Info.
word name(const complex &c)
Return string representation of complex.
fileName timePath() const
Return current time path.
const areaVectorField & faceAreaNormals() const
Return face area normals.
const edgeScalarField & lPN() const
Return reference to PN geodesic distance.
vectorField & pointsDisplacementDir()
Return reference to point displacement direction field.
const labelListList & pointEdges() const
Return point-edge addressing.
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
A List with indirect addressing.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const Type & lookupObject(const word &name, const bool recursive=false) const
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
virtual const pointField & oldPoints() const
Return old points (mesh motion)
tmp< GeometricField< Type, faPatchField, areaMesh > > edgeIntegrate(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Abstract base class for turbulence models (RAS, LES and laminar).
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
label findPatchID(const word &patchName) const
Find patch index given a name, return -1 if not found.
const Field< point_type > & points() const
Return reference to global points.
Macros for easy insertion into run-time selection tables.
const volVectorField & U() const
Return constant reference to velocity field.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
const areaScalarField & fsNetPhi() const
Return free-surface net flux.
tmp< vectorField > freeSurfaceSnGradU()
Return free surface normal derivative of velocity.
const pointField & points() const
Return mesh points.
const uniformDimensionedVectorField & g
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
errorManip< error > abort(error &err)
tmp< scalarField > freeSurfaceSnGradUn()
Return free surface normal derivative of normal velocity comp.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
areaScalarField & surfaceTension()
Return surface tension field.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
Output to file stream, using an OSstream.
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
void correctBoundaryConditions()
Correct boundary field.
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
tmp< GeometricField< Type, faePatchField, edgeMesh > > linearEdgeInterpolate(const GeometricField< Type, faPatchField, areaMesh > &vf)
const labelList & pointLabels() const
Return patch point labels.
labelList & motionPointsMask()
Return reference to motion points mask field.
label prevLabel(const label i) const
Previous vertex on face.
const dimensionedScalar a0
Bohr radius: default SI units: [m].
label nextLabel(const label i) const
Next vertex on face.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Virtual base class for mesh motion solver.
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep track of the Delaunay cells (tets)...
const areaScalarField & faceCurvatures() const
Return face curvatures.
Second-order backward-differencing ddt using the current and two previous time-step values.
const std::string patch
OpenFOAM patch number as a std::string.
dimensionedScalar sqrt(const dimensionedScalar &ds)
void writeVTKControlPoints()
Write VTK freeSurface control points.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
const dimensionedScalar e
Elementary charge.
vectorField & facesDisplacementDir()
Return reference to control points displacement direction field.
const faBoundaryMesh & boundary() const
Return constant reference to boundary mesh.
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
void correctUsBoundaryConditions()
Correct surface velocity boundary conditions.
pointVectorField & pointMotionU()
Return reference to the point motion velocity field.
const dimensionedScalar c
Speed of light in a vacuum.
const surfaceScalarField & phi() const
Return constant reference to velocity field.
Finite area mesh. Used for 2-D non-Euclidian finite area method.
const Time & time() const
Return the top-level database.
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
A face is a list of labels corresponding to mesh vertices.
CGAL::indexedCell< K > Cb
vectorField & controlPoints()
Return control points.
virtual bool update()
Update the mesh for both mesh motion and topology change.
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
const Vector< label > N(dict.get< Vector< label >>("N"))
const volScalarField & p0
Type gMin(const FieldField< Field, Type > &f)
label timeIndex() const
Return current time index.
const surfactantProperties & surfactant() const
Return surfactant properties.
const dimensionSet dimVolume(pow3(dimLength))
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
defineTypeNameAndDebug(combustionModel, 0)
Basic first-order Euler implicit/explicit ddt using only the current and previous time-step values.
dimensionedScalar neg(const dimensionedScalar &ds)
static tmp< GeometricField< Type, faePatchField, edgeMesh > > interpolate(const GeometricField< Type, faPatchField, areaMesh > &tvf, const edgeScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Type gMax(const FieldField< Field, Type > &f)
areaVectorField & Us()
Return free-surface velocity field.
static const gravity & New(const Time &runTime)
Construct on Time.
edgeScalarField & Phis()
Return free-surface fluid flux field.
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
const Boundary & boundaryField() const
Return const-reference to the boundary field.
~interfaceTrackingFvMesh()
Destructor.
The dynamicMotionSolverFvMesh.
tmp< faMatrix< Type > > div(const edgeScalarField &flux, const GeometricField< Type, faPatchField, areaMesh > &vf, const word &name)
const areaVectorField & areaCentres() const
Return face centres as areaVectorField.
dimensionedScalar cos(const dimensionedScalar &ds)