Go to the documentation of this file.
65 void Foam::interfaceTrackingFvMesh::initializeData()
69 const word fsPatchName(
motion().get<word>(
"fsPatchName"));
76 <<
"Patch name " << fsPatchName <<
" not found."
80 fsPatchIndex_ =
patch.index();
87 for (
const word& patchName : pointNormalsCorrectionPatches_)
94 <<
"Patch name '" << patchName
95 <<
"' for point normals correction does not exist"
104 if (!normalMotionDir_)
114 "nonReflectingFreeSurfacePatches",
115 nonReflectingFreeSurfacePatches_
120 void Foam::interfaceTrackingFvMesh::makeUs()
const
123 <<
"making free-surface velocity field" <<
nl;
128 <<
"free-surface velocity field already exists"
135 zeroGradientFaPatchVectorField::typeName
143 == wedgeFaPatch::typeName
146 patchFieldTypes[patchI] =
147 wedgeFaPatchVectorField::typeName;
151 label ngbPolyPatchID =
154 if (ngbPolyPatchID != -1)
159 == wallFvPatch::typeName
162 patchFieldTypes[patchI] =
163 slipFaPatchVectorField::typeName;
169 for (
const word& patchName : fixedFreeSurfacePatches_)
171 const label fixedPatchID =
174 if (fixedPatchID == -1)
177 <<
"Wrong faPatch name '" << patchName
178 <<
"' in the fixedFreeSurfacePatches list"
179 <<
" defined in the dynamicMeshDict dictionary"
183 label ngbPolyPatchID =
186 if (ngbPolyPatchID != -1)
191 == wallFvPatch::typeName
194 patchFieldTypes[fixedPatchID] =
195 fixedValueFaPatchVectorField::typeName;
215 for (
const word& patchName : fixedFreeSurfacePatches_)
219 if (fixedPatchID == -1)
222 <<
"Wrong faPatch name '" << patchName
223 <<
"' in the fixedFreeSurfacePatches list"
224 <<
" defined in the dynamicMeshDict dictionary"
228 label ngbPolyPatchID =
231 if (ngbPolyPatchID != -1)
236 == wallFvPatch::typeName
239 UsPtr_->boundaryFieldRef()[fixedPatchID] ==
Zero;
246 void Foam::interfaceTrackingFvMesh::makeFsNetPhi()
const
249 <<
"making free-surface net flux" <<
nl;
254 <<
"free-surface net flux already exists"
274 void Foam::interfaceTrackingFvMesh::makeControlPoints()
277 <<
"making control points" <<
nl;
279 if (controlPointsPtr_)
282 <<
"control points already exists"
296 Info<<
"Reading control points" <<
endl;
312 Info<<
"Creating new control points" <<
endl;
324 aMesh().areaCentres().internalField()
327 initializeControlPointsPosition();
332 void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
335 <<
"making motion points mask" <<
nl;
337 if (motionPointsMaskPtr_)
340 <<
"motion points mask already exists"
357 == processorFaPatch::typeName
363 forAll(patchPoints, pointI)
365 motionPointsMask()[patchPoints[pointI]] = -1;
371 for (
const word& patchName : fixedFreeSurfacePatches_)
375 if (fixedPatchID == -1)
378 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
379 <<
" defined in the dynamicMeshDict dictionary"
386 forAll(patchPoints, pointI)
388 motionPointsMask()[patchPoints[pointI]] = 0;
394 void Foam::interfaceTrackingFvMesh::makeDirections()
397 <<
"make displacement directions for points and control points" <<
nl;
399 if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_)
402 <<
"points, control points displacement directions already exist"
406 pointsDisplacementDirPtr_ =
413 facesDisplacementDirPtr_ =
420 if (!normalMotionDir())
422 if (
mag(motionDir_) < SMALL)
425 <<
"Zero motion direction"
429 facesDisplacementDir() = motionDir_;
430 pointsDisplacementDir() = motionDir_;
433 updateDisplacementDirections();
437 void Foam::interfaceTrackingFvMesh::makePhis()
440 <<
"making free-surface flux" <<
nl;
445 <<
"free-surface flux already exists"
464 void Foam::interfaceTrackingFvMesh::makeSurfactConc()
const
467 <<
"making free-surface surfactant concentration field" <<
nl;
472 <<
"free-surface surfactant concentration field already exists"
495 void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc()
const
498 <<
"making volume surfactant concentration field" <<
nl;
500 if (bulkSurfactConcPtr_)
503 <<
"volume surfactant concentration field already exists"
528 bulkSurfactConc = surfactant().bulkConc();
534 void Foam::interfaceTrackingFvMesh::makeSurfaceTension()
const
537 <<
"making surface tension field" <<
nl;
539 if (surfaceTensionPtr_)
542 <<
"surface tension field already exists"
556 sigma() + surfactant().dSigma(surfactantConcentration())/rho_
561 void Foam::interfaceTrackingFvMesh::makeSurfactant()
const
564 <<
"making surfactant properties" <<
nl;
569 <<
"surfactant properties already exists"
574 motion().
subDict(
"surfactantProperties");
580 void Foam::interfaceTrackingFvMesh::makeContactAngle()
583 <<
"making contact angle field" <<
nl;
585 if (contactAnglePtr_)
588 <<
"contact angle already exists"
603 Info<<
"Reading contact angle field" <<
endl;
622 void Foam::interfaceTrackingFvMesh::updateDisplacementDirections()
624 if (normalMotionDir())
632 if (contactAnglePtr_)
634 label ngbPolyPatchID =
637 if (ngbPolyPatchID != -1)
642 == wallFvPatch::typeName
651 .ngbPolyPatchPointNormals()
654 forAll(patchPoints, pointI)
656 pointsDisplacementDir()[patchPoints[pointI]] -=
660 & pointsDisplacementDir()[patchPoints[pointI]]
663 pointsDisplacementDir()[patchPoints[pointI]] /=
666 pointsDisplacementDir()
678 facesDisplacementDir() =
685 facesDisplacementDir()
686 *(facesDisplacementDir()&(controlPoints() - Cf))
692 void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition()
700 correctPointDisplacement(sweptVolCorr, displacement);
708 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
715 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
722 deltaH[faceI] = sweptVol[faceI]/
723 ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL);
725 if (
mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL)
732 <<
"Something is wrong with specified motion direction"
737 for (
const word& patchName : fixedFreeSurfacePatches_)
741 if (fixedPatchID == -1)
744 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
745 <<
" defined in the freeSurfaceProperties dictionary"
754 deltaH[eFaces[edgeI]] *= 2.0;
758 controlPoints() += facesDisplacementDir()*deltaH;
763 void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh()
765 Info<<
"Smoothing free surface mesh" <<
endl;
779 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
785 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
790 sweptVol/(faceArea & facesDisplacementDir())
793 for (
const word& patchName : fixedFreeSurfacePatches_)
797 if (fixedPatchID == -1)
800 <<
"Wrong faPatch name fixedFreeSurfacePatches list"
801 <<
" defined in the dynamicMeshDict dictionary"
810 deltaHf[eFaces[edgeI]] *= 2.0;
814 controlPoints() += facesDisplacementDir()*deltaHf;
816 displacement = pointDisplacement();
819 refCast<velocityMotionSolver>
827 fixedValuePointPatchVectorField& fsPatchPointMeshU =
828 refCast<fixedValuePointPatchVectorField>
843 void Foam::interfaceTrackingFvMesh::updateSurfaceFlux()
849 void Foam::interfaceTrackingFvMesh::updateSurfactantConcentration()
851 if (!pureFreeSurface())
861 +
fam::div(Phis(), surfactantConcentration())
864 surfactant().diffusion(),
865 surfactantConcentration()
869 if (surfactant().soluble())
885 zeroGradientFaPatchScalarField::typeName
889 bulkSurfactantConcentration().boundaryField()[fsPatchIndex()];
890 Cb.correctBoundaryConditions();
895 surfactant().adsorptionCoeff()*
Cb
896 + surfactant().adsorptionCoeff()
897 *surfactant().desorptionCoeff(),
898 surfactantConcentration()
900 - surfactant().adsorptionCoeff()
901 *
Cb*surfactant().saturatedConc();
909 sigma() + surfactant().dSigma(surfactantConcentration())/rho_;
911 if (
neg(
min(surfaceTension().internalField().
field())))
914 <<
"Surface tension is negative"
921 Foam::vector Foam::interfaceTrackingFvMesh::totalPressureForce()
const
931 return gSum(pressureForces);
935 Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce()
const
955 U().boundaryField()[fsPatchIndex()].
snGrad()
968 return gSum(viscousForces);
972 Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce()
const
982 if (pureFreeSurface())
993 surfTensionForces = surfaceTension().internalField().field()*
K*S*
n;
996 return gSum(surfTensionForces);
1000 Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber()
1004 if (pureFreeSurface())
1010 mesh().time().deltaT().value()/
1029 mesh().time().deltaT().value()/
1041 void Foam::interfaceTrackingFvMesh::updateProperties()
1046 "transportProperties"
1055 void Foam::interfaceTrackingFvMesh::correctPointDisplacement
1070 for (
const word& patchName : fixedFreeSurfacePatches_)
1084 label curFace = eFaces[edgeI];
1086 const labelList& curPoints = faces[curFace];
1088 forAll(curPoints, pointI)
1090 label curPoint = curPoints[pointI];
1091 label index = pLabels.find(curPoint);
1107 forAll(corrPoints, pointI)
1109 label curPoint = corrPoints[pointI];
1117 label index = eFaces.find(curFace);
1131 forAll(corrPoints, pointI)
1133 label curPoint = corrPoints[pointI];
1137 const labelList& curPointFaces = corrPointFaces[pointI];
1139 forAll(curPointFaces, faceI)
1141 const face& curFace = faces[curPointFaces[faceI]];
1149 const vector&
c = pointsDisplacementDir()[curPoint];
1151 curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^
b)&
c);
1154 curDisp /= curPointFaces.size();
1156 displacement[curPoint] =
1157 curDisp*pointsDisplacementDir()[curPoint];
1162 for (
const word& patchName : nonReflectingFreeSurfacePatches_)
1164 label nonReflectingPatchID =
1177 forAll(corrPoints, pointI)
1179 label curPoint = corrPoints[pointI];
1187 label index = eFaces.find(curFace);
1202 forAll(corrPoints, pointI)
1204 label curPoint = corrPoints[pointI];
1208 const labelList& curPointFaces = corrPointFaces[pointI];
1210 forAll(curPointFaces, faceI)
1212 const face& curFace = faces[curPointFaces[faceI]];
1222 if (corrPoints.find(next) == -1)
1237 vector c0 = displacement[p1];
1239 scalar V0 =
mag((
a0^b0)&c0)/2;
1241 scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0;
1243 if (corrPoints.find(prev) != -1)
1247 (
points[next] + displacement[next])
1249 const vector&
c = pointsDisplacementDir()[curPoint];
1251 curDisp += 2*DV/((a^
b)&
c);
1256 - (
points[prev] + displacement[prev]);
1258 const vector&
c = pointsDisplacementDir()[curPoint];
1260 curDisp += 2*DV/((a^
b)&
c);
1264 curDisp /= curPointFaces.size();
1266 displacement[curPoint] =
1267 curDisp*pointsDisplacementDir()[curPoint];
1273 void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals()
1284 if (contactAnglePtr_ && correctContactLineNormals())
1286 Info<<
"Correcting contact line normals" <<
endl;
1325 == processorFaPatch::typeName
1329 refCast<const processorFaPatch>
1337 forAll(patchPointLabels, pointI)
1339 label curPoint = patchPointLabels[pointI];
1346 const labelList& curPointEdges = pointEdges[curPoint];
1348 forAll(curPointEdges, edgeI)
1350 label curEdge = curPointEdges[edgeI];
1352 if (edgeFaces[curEdge].size() == 1)
1359 label index = curEdges.find(curEdge);
1366 != processorFaPatch::typeName
1383 vector t = edges[curEdge].vec(oldPoints);
1384 t /=
mag(t) + SMALL;
1387 pointFaces[curPoint];
1389 forAll(curPointFaces, fI)
1391 tangent.ref().field()[curPointFaces[fI]] = t;
1398 tangent.correctBoundaryConditions();
1403 label ngbPolyPatchID =
1406 if (ngbPolyPatchID != -1)
1411 == wallFvPatch::typeName
1419 - contactAnglePtr_->boundaryField()[patchI]
1434 rotationAxis /=
mag(rotationAxis) + SMALL;
1444 forAll(pointEdges, pointI)
1448 forAll(pointEdges[pointI], eI)
1452 + pointEdges[pointI][eI];
1454 vector e = edges[curEdge].vec(oldPoints);
1456 e *= (
e&rotationAxis[pointI])
1457 /
mag(
e&rotationAxis[pointI]);
1459 e /=
mag(
e) + SMALL;
1464 if (pointEdges[pointI].size() == 1)
1466 label curPoint = patchPoints[pointI];
1473 label procPatchID = -1;
1479 forAll(curPointEdges, edgeI)
1481 label curEdge = curPointEdges[edgeI];
1483 if (edgeFaces[curEdge].size() == 1)
1491 curEdges.find(curEdge);
1498 == processorFaPatch::typeName
1510 if (procPatchID != -1)
1513 tangent.boundaryField()[procPatchID]
1514 .patchNeighbourField()()[
edgeID];
1516 t *= (t&rotationAxis[pointI])
1517 /(
mag(t&rotationAxis[pointI]) + SMALL);
1519 t /=
mag(t) + SMALL;
1525 rotationAxis[pointI] = rotAx/(
mag(rotAx) + SMALL);
1529 ngbN = ngbN*
cos(rotAngle)
1530 + rotationAxis*(rotationAxis & ngbN)*(1 -
cos(rotAngle))
1531 + (rotationAxis^ngbN)*
sin(rotAngle);
1534 forAll(patchPoints, pointI)
1536 N[patchPoints[pointI]] -=
1537 ngbN[pointI]*(ngbN[pointI]&
N[patchPoints[pointI]]);
1539 N[patchPoints[pointI]] /=
1540 mag(
N[patchPoints[pointI]]) + SMALL;
1553 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh(
const IOobject& io)
1556 aMeshPtr_(new
faMesh(*this)),
1558 fixedFreeSurfacePatches_
1560 motion().
lookup(
"fixedFreeSurfacePatches")
1562 nonReflectingFreeSurfacePatches_(),
1563 pointNormalsCorrectionPatches_
1565 motion().
lookup(
"pointNormalsCorrectionPatches")
1569 motion().
lookup(
"normalMotionDir")
1572 smoothing_(motion().lookupOrDefault(
"smoothing", false)),
1573 pureFreeSurface_(motion().lookupOrDefault(
"pureFreeSurface", true)),
1574 rigidFreeSurface_(motion().lookupOrDefault(
"rigidFreeSurface", false)),
1575 correctContactLineNormals_
1577 motion().lookupOrDefault(
"correctContactLineNormals", false)
1583 controlPointsPtr_(nullptr),
1584 motionPointsMaskPtr_(nullptr),
1585 pointsDisplacementDirPtr_(nullptr),
1586 facesDisplacementDirPtr_(nullptr),
1587 fsNetPhiPtr_(nullptr),
1589 surfactConcPtr_(nullptr),
1590 bulkSurfactConcPtr_(nullptr),
1591 surfaceTensionPtr_(nullptr),
1592 surfactantPtr_(nullptr),
1593 contactAnglePtr_(nullptr)
1599 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1614 std::move(allOwner),
1615 std::move(allNeighbour),
1618 aMeshPtr_(
new faMesh(*
this)),
1620 fixedFreeSurfacePatches_
1622 motion().
lookup(
"fixedFreeSurfacePatches")
1624 nonReflectingFreeSurfacePatches_(),
1625 pointNormalsCorrectionPatches_
1627 motion().
lookup(
"pointNormalsCorrectionPatches")
1631 motion().
lookup(
"normalMotionDir")
1634 smoothing_(motion().lookupOrDefault(
"smoothing",
false)),
1635 pureFreeSurface_(motion().lookupOrDefault(
"pureFreeSurface",
true)),
1640 controlPointsPtr_(
nullptr),
1641 motionPointsMaskPtr_(
nullptr),
1642 pointsDisplacementDirPtr_(
nullptr),
1643 facesDisplacementDirPtr_(
nullptr),
1644 fsNetPhiPtr_(
nullptr),
1646 surfactConcPtr_(
nullptr),
1647 bulkSurfactConcPtr_(
nullptr),
1648 surfaceTensionPtr_(
nullptr),
1649 surfactantPtr_(
nullptr),
1650 contactAnglePtr_(
nullptr)
1706 return *fsNetPhiPtr_;
1717 return *fsNetPhiPtr_;
1723 forAll(
Us().boundaryField(), patchI)
1727 Us().boundaryField()[patchI].
type()
1728 == calculatedFaPatchVectorField::typeName
1733 pUs =
Us().boundaryField()[patchI].patchInternalField();
1735 label ngbPolyPatchID =
1738 if (ngbPolyPatchID != -1)
1743 U().boundaryField()[ngbPolyPatchID].
type()
1744 == slipFvPatchVectorField::typeName
1748 U().boundaryField()[ngbPolyPatchID].
type()
1749 == symmetryFvPatchVectorField::typeName
1764 Us().correctBoundaryConditions();
1772 Us().ref().field() =
U().boundaryField()[fsPatchIndex()];
1780 correctUsBoundaryConditions();
1786 return *getObjectPtr<const volVectorField>(
"U");
1792 return *getObjectPtr<const volScalarField>(
"p");
1798 return *getObjectPtr<const surfaceScalarField>(
"phi");
1806 auto& SnGradU = tSnGradU.ref();
1813 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1820 gradUs -=
n*(
n & gradUs);
1830 if (!pureFreeSurface() &&
max(
nu) > SMALL)
1832 tangentialSurfaceTensionForce =
1833 surfaceTensionGrad()().internalField();
1837 tangentialSurfaceTensionForce/(
nu + SMALL)
1849 auto& SnGradUn = tSnGradUn.ref();
1854 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1867 auto& pressureJump = tPressureJump.ref();
1881 + 2.0*
nu*freeSurfaceSnGradUn();
1883 if (pureFreeSurface())
1889 pressureJump -= surfaceTension().internalField()*
K;
1892 return tPressureJump;
1898 if (!controlPointsPtr_)
1900 makeControlPoints();
1903 return *controlPointsPtr_;
1909 if (!motionPointsMaskPtr_)
1911 makeMotionPointsMask();
1914 return *motionPointsMaskPtr_;
1920 if (!pointsDisplacementDirPtr_)
1925 return *pointsDisplacementDirPtr_;
1931 if (!facesDisplacementDirPtr_)
1936 return *facesDisplacementDirPtr_;
1953 if (!surfactConcPtr_)
1958 return *surfactConcPtr_;
1965 if (!surfactConcPtr_)
1970 return *surfactConcPtr_;
1977 if (!bulkSurfactConcPtr_)
1979 makeBulkSurfactConc();
1982 return *bulkSurfactConcPtr_;
1989 if (!bulkSurfactConcPtr_)
1991 makeBulkSurfactConc();
1994 return *bulkSurfactConcPtr_;
2001 if (!surfaceTensionPtr_)
2003 makeSurfaceTension();
2006 return *surfaceTensionPtr_;
2013 if (!surfaceTensionPtr_)
2015 makeSurfaceTension();
2018 return *surfaceTensionPtr_;
2025 if (!surfactantPtr_)
2030 return *surfactantPtr_;
2041 "surfaceTensionGrad",
2053 auto&
grad = tgrad.ref();
2058 grad.correctBoundaryConditions();
2068 if (smoothing_ && !rigidFreeSurface_)
2070 smoothFreeSurfaceMesh();
2071 clearControlPoints();
2074 updateDisplacementDirections();
2078 Info<<
"Maximal capillary Courant number: "
2079 << maxCourantNumber() <<
endl;
2083 Info<<
"Free surface curvature: min = " <<
gMin(
K)
2089 if (!rigidFreeSurface_)
2105 Info<<
"Free surface continuity error : sum local = "
2106 <<
gSum(
mag(sweptVolCorr)) <<
", global = " <<
gSum(sweptVolCorr)
2110 fsNetPhi().ref().field() = sweptVolCorr;
2116 "ddt(" +
U().
name() +
')'
2154 <<
"Unsupported temporal differencing scheme : "
2164 sweptVolCorr/(Sf*(Nf & facesDisplacementDir()))
2167 for (
const word& patchName : fixedFreeSurfacePatches_)
2169 label fixedPatchID =
2172 if (fixedPatchID == -1)
2175 <<
"Wrong faPatch name '" << patchName
2176 <<
"'in the fixedFreeSurfacePatches list"
2177 <<
" defined in dynamicMeshDict dictionary"
2186 deltaHf[eFaces[edgeI]] *= 2.0;
2190 controlPoints() += facesDisplacementDir()*deltaHf;
2192 pointField displacement(pointDisplacement());
2193 correctPointDisplacement(sweptVolCorr, displacement);
2196 refCast<velocityMotionSolver>
2204 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2205 refCast<fixedValuePointPatchVectorField>
2213 fsPatchPointMeshU ==
2218 correctContactLinePointNormals();
2229 refCast<velocityMotionSolver>
2237 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2238 refCast<fixedValuePointPatchVectorField>
2246 fsPatchPointMeshU ==
2254 updateSurfactantConcentration();
2263 OFstream mps(
mesh().time().timePath()/
"freeSurface.vtk");
2268 mps <<
"# vtk DataFile Version 2.0" <<
nl
2271 <<
"DATASET POLYDATA" <<
nl
2272 <<
"POINTS " <<
points.size() <<
" float" <<
nl;
2280 mlpBuffer[counter++] = float(
points[i].
x());
2281 mlpBuffer[counter++] = float(
points[i].
y());
2282 mlpBuffer[counter++] = float(
points[i].z());
2287 mps << mlpBuffer[i] <<
' ';
2289 if (i > 0 && (i % 10) == 0)
2296 label nFaceVerts = 0;
2300 nFaceVerts += faces[faceI].size() + 1;
2307 const face&
f = faces[faceI];
2309 mlfBuffer[counter++] =
f.size();
2313 mlfBuffer[counter++] =
f[fpI];
2318 mps <<
"POLYGONS " << faces.size() <<
' ' << nFaceVerts <<
endl;
2322 mps << mlfBuffer[i] <<
' ';
2324 if (i > 0 && (i % 10) == 0)
2346 Info<<
"Writing free surface control point to " <<
name <<
endl;
2348 mps <<
"# vtk DataFile Version 2.0" <<
nl
2351 <<
"DATASET POLYDATA" <<
nl
2352 <<
"POINTS " << controlPoints().size() <<
" float" <<
nl;
2354 forAll(controlPoints(), pointI)
2356 mps << controlPoints()[pointI].x() <<
' '
2357 << controlPoints()[pointI].y() <<
' '
2358 << controlPoints()[pointI].z() <<
nl;
2362 mps <<
"VERTICES " << controlPoints().size() <<
' '
2363 << controlPoints().size()*2 <<
nl;
2365 forAll(controlPoints(), pointI)
2367 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 Field< PointType > & points() const
Return reference to global points.
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.
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)
const Field< PointType > & localPoints() const
Return pointField of points in patch.
dimensionedScalar sin(const dimensionedScalar &ds)
bool correctPatchPointNormals(const label patchID) const
Whether point normals should be corrected for a patch.
static bool & parRun()
Is this a parallel run?
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)
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 > &)
Return the correction form of the given matrix.
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))
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
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 cellModel * lookup(const word &modelName)
Deprecated(2017-11) equivalent to cellModel::ptr static method.
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.
Lookup type of boundary radiation properties.
virtual const pointField & oldPoints() const
Return old points for 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.
Macros for easy insertion into run-time selection tables.
const volVectorField & U() const
Return constant reference to velocity field.
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.
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.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
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)