76 void Foam::interfaceTrackingFvMesh::initializeData()
80 const word fsPatchName(
motion().get<word>(
"fsPatchName"));
87 <<
"Patch name " << fsPatchName <<
" not found." 91 fsPatchIndex_ =
patch.index();
95 if (!pointNormalsCorrectionPatches_.empty())
99 for (
const word& patchName : pointNormalsCorrectionPatches_)
106 <<
"Patch name '" << patchName
107 <<
"' for point normals correction does not exist" 116 if (!normalMotionDir_)
126 "nonReflectingFreeSurfacePatches",
127 nonReflectingFreeSurfacePatches_
132 void Foam::interfaceTrackingFvMesh::makeUs()
const 135 <<
"making free-surface velocity field" <<
nl;
140 <<
"free-surface velocity field already exists" 155 == wedgeFaPatch::typeName
158 patchFieldTypes[patchI] =
159 wedgeFaPatchVectorField::typeName;
163 label ngbPolyPatchID =
164 aMesh().boundary()[patchI].ngbPolyPatchIndex();
166 if (ngbPolyPatchID != -1)
171 == wallFvPatch::typeName
174 patchFieldTypes[patchI] =
175 slipFaPatchVectorField::typeName;
181 for (
const word& patchName : fixedFreeSurfacePatches_)
183 const label fixedPatchID =
184 aMesh().boundary().findPatchID(patchName);
186 if (fixedPatchID == -1)
189 <<
"Wrong faPatch name '" << patchName
190 <<
"' in the fixedFreeSurfacePatches list" 191 <<
" defined in the dynamicMeshDict dictionary" 195 label ngbPolyPatchID =
196 aMesh().boundary()[fixedPatchID].ngbPolyPatchIndex();
198 if (ngbPolyPatchID != -1)
203 == wallFvPatch::typeName
206 patchFieldTypes[fixedPatchID] =
207 fixedValueFaPatchVectorField::typeName;
227 for (
const word& patchName : fixedFreeSurfacePatches_)
229 const label fixedPatchID = aMesh().boundary().findPatchID(patchName);
231 if (fixedPatchID == -1)
234 <<
"Wrong faPatch name '" << patchName
235 <<
"' in the fixedFreeSurfacePatches list" 236 <<
" defined in the dynamicMeshDict dictionary" 240 label ngbPolyPatchID =
241 aMesh().boundary()[fixedPatchID].ngbPolyPatchIndex();
243 if (ngbPolyPatchID != -1)
248 == wallFvPatch::typeName
251 UsPtr_->boundaryFieldRef()[fixedPatchID] ==
Zero;
258 void Foam::interfaceTrackingFvMesh::makeFsNetPhi()
const 261 <<
"making free-surface net flux" <<
nl;
266 <<
"free-surface net flux already exists" 286 void Foam::interfaceTrackingFvMesh::makeControlPoints()
289 <<
"making control points" <<
nl;
291 if (controlPointsPtr_)
294 <<
"control points already exists" 308 Info<<
"Reading control points" <<
endl;
324 Info<<
"Creating new control points" <<
endl;
336 aMesh().areaCentres().internalField()
339 initializeControlPointsPosition();
344 void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
347 <<
"making motion points mask" <<
nl;
349 if (motionPointsMaskPtr_)
352 <<
"motion points mask already exists" 369 == processorFaPatch::typeName
373 aMesh().boundary()[patchI].pointLabels();
375 forAll(patchPoints, pointI)
377 motionPointsMask()[patchPoints[pointI]] = -1;
383 for (
const word& patchName : fixedFreeSurfacePatches_)
385 const label fixedPatchID = aMesh().boundary().findPatchID(patchName);
387 if (fixedPatchID == -1)
390 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list" 391 <<
" defined in the dynamicMeshDict dictionary" 396 aMesh().boundary()[fixedPatchID].pointLabels();
398 forAll(patchPoints, pointI)
400 motionPointsMask()[patchPoints[pointI]] = 0;
406 void Foam::interfaceTrackingFvMesh::makeDirections()
409 <<
"make displacement directions for points and control points" <<
nl;
411 if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_)
414 <<
"points, control points displacement directions already exist" 418 pointsDisplacementDirPtr_ =
425 facesDisplacementDirPtr_ =
432 if (!normalMotionDir())
434 if (
mag(motionDir_) < SMALL)
437 <<
"Zero motion direction" 441 facesDisplacementDir() = motionDir_;
442 pointsDisplacementDir() = motionDir_;
445 updateDisplacementDirections();
449 void Foam::interfaceTrackingFvMesh::makePhis()
452 <<
"making free-surface flux" <<
nl;
457 <<
"free-surface flux already exists" 476 void Foam::interfaceTrackingFvMesh::makeSurfactConc()
const 479 <<
"making free-surface surfactant concentration field" <<
nl;
484 <<
"free-surface surfactant concentration field already exists" 507 void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc()
const 510 <<
"making volume surfactant concentration field" <<
nl;
512 if (bulkSurfactConcPtr_)
515 <<
"volume surfactant concentration field already exists" 540 bulkSurfactConc = surfactant().bulkConc();
546 void Foam::interfaceTrackingFvMesh::makeSurfaceTension()
const 549 <<
"making surface tension field" <<
nl;
551 if (surfaceTensionPtr_)
554 <<
"surface tension field already exists" 568 sigma() + surfactant().dSigma(surfactantConcentration())/rho_
573 void Foam::interfaceTrackingFvMesh::makeSurfactant()
const 576 <<
"making surfactant properties" <<
nl;
581 <<
"surfactant properties already exists" 586 motion().
subDict(
"surfactantProperties");
592 void Foam::interfaceTrackingFvMesh::makeContactAngle()
595 <<
"making contact angle field" <<
nl;
597 if (contactAnglePtr_)
600 <<
"contact angle already exists" 616 Info<<
"Reading contact angle field" <<
endl;
623 void Foam::interfaceTrackingFvMesh::updateDisplacementDirections()
625 if (normalMotionDir())
628 pointsDisplacementDir() = aMesh().pointAreaNormals();
633 if (contactAnglePtr_)
635 label ngbPolyPatchID =
636 aMesh().boundary()[patchI].ngbPolyPatchIndex();
638 if (ngbPolyPatchID != -1)
643 == wallFvPatch::typeName
647 aMesh().boundary()[patchI].pointLabels();
652 .ngbPolyPatchPointNormals()
655 forAll(patchPoints, pointI)
657 pointsDisplacementDir()[patchPoints[pointI]] -=
661 & pointsDisplacementDir()[patchPoints[pointI]]
664 pointsDisplacementDir()[patchPoints[pointI]] /=
667 pointsDisplacementDir()
679 facesDisplacementDir() =
680 aMesh().faceAreaNormals().internalField();
683 const vectorField& Cf = aMesh().areaCentres().internalField();
686 facesDisplacementDir()
687 *(facesDisplacementDir()&(controlPoints() - Cf))
693 void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition()
696 const faceList& faces = aMesh().faces();
701 correctPointDisplacement(sweptVolCorr, displacement);
709 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
716 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
723 deltaH[faceI] = sweptVol[faceI]/
724 ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL);
726 if (
mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL)
733 <<
"Something is wrong with specified motion direction" 738 for (
const word& patchName : fixedFreeSurfacePatches_)
740 label fixedPatchID = aMesh().boundary().findPatchID(patchName);
742 if (fixedPatchID == -1)
745 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list" 746 <<
" defined in the freeSurfaceProperties dictionary" 751 aMesh().boundary()[fixedPatchID].edgeFaces();
755 deltaH[eFaces[edgeI]] *= 2.0;
759 controlPoints() += facesDisplacementDir()*deltaH;
764 void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh()
766 Info<<
"Smoothing free surface mesh" <<
endl;
768 controlPoints() = aMesh().areaCentres().internalField();
772 const faceList& faces = aMesh().faces();
780 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
786 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
791 sweptVol/(faceArea & facesDisplacementDir())
794 for (
const word& patchName : fixedFreeSurfacePatches_)
796 label fixedPatchID = aMesh().boundary().findPatchID(patchName);
798 if (fixedPatchID == -1)
801 <<
"Wrong faPatch name fixedFreeSurfacePatches list" 802 <<
" defined in the dynamicMeshDict dictionary" 807 aMesh().boundary()[fixedPatchID].edgeFaces();
811 deltaHf[eFaces[edgeI]] *= 2.0;
815 controlPoints() += facesDisplacementDir()*deltaHf;
817 displacement = pointDisplacement();
820 refCast<velocityMotionSolver>
828 fixedValuePointPatchVectorField& fsPatchPointMeshU =
829 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())
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 925 const vectorField&
n = aMesh().faceAreaNormals().internalField();
931 return gSum(pressureForces);
935 Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce()
const 951 const vectorField&
n = aMesh().faceAreaNormals().internalField();
955 U().boundaryField()[fsPatchIndex()].
snGrad()
968 return gSum(viscousForces);
972 Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce()
const 976 const vectorField&
n = aMesh().faceAreaNormals().internalField();
978 const scalarField&
K = aMesh().faceCurvatures().internalField();
982 if (pureFreeSurface())
988 aMesh().Le()*aMesh().edgeLengthCorrection()
993 surfTensionForces = surfaceTension().internalField().field()*
K*S*
n;
996 return gSum(surfTensionForces);
1000 Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber()
1004 if (pureFreeSurface())
1010 mesh().time().deltaTValue()/
1029 mesh().time().deltaTValue()/
1041 void Foam::interfaceTrackingFvMesh::updateProperties()
1046 "transportProperties" 1055 void Foam::interfaceTrackingFvMesh::correctPointDisplacement
1062 aMesh().patch().pointFaces();
1065 aMesh().patch().localFaces();
1068 aMesh().patch().localPoints();
1070 for (
const word& patchName : fixedFreeSurfacePatches_)
1072 label fixedPatchID = aMesh().boundary().findPatchID(patchName);
1075 aMesh().boundary()[fixedPatchID].pointLabels();
1078 aMesh().boundary()[fixedPatchID].edgeFaces();
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);
1104 forAll(corrPoints, pointI)
1106 label curPoint = corrPoints[pointI];
1112 label curFace =
pFaces[curPoint][faceI];
1114 label index = eFaces.
find(curFace);
1125 forAll(corrPoints, pointI)
1127 label curPoint = corrPoints[pointI];
1131 const labelList& curPointFaces = corrPointFaces[pointI];
1133 forAll(curPointFaces, faceI)
1135 const face& curFace = faces[curPointFaces[faceI]];
1137 label ptInFace = curFace.
which(curPoint);
1138 label next = curFace.
nextLabel(ptInFace);
1139 label prev = curFace.
prevLabel(ptInFace);
1143 const vector&
c = pointsDisplacementDir()[curPoint];
1145 curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^
b)&
c);
1148 curDisp /= curPointFaces.
size();
1150 displacement[curPoint] =
1151 curDisp*pointsDisplacementDir()[curPoint];
1156 for (
const word& patchName : nonReflectingFreeSurfacePatches_)
1158 label nonReflectingPatchID =
1159 aMesh().boundary().findPatchID(patchName);
1162 aMesh().boundary()[nonReflectingPatchID].pointLabels();
1165 aMesh().boundary()[nonReflectingPatchID].edgeFaces();
1171 forAll(corrPoints, pointI)
1173 label curPoint = corrPoints[pointI];
1179 label curFace =
pFaces[curPoint][faceI];
1181 label index = eFaces.
find(curFace);
1193 forAll(corrPoints, pointI)
1195 label curPoint = corrPoints[pointI];
1199 const labelList& curPointFaces = corrPointFaces[pointI];
1201 forAll(curPointFaces, faceI)
1203 const face& curFace = faces[curPointFaces[faceI]];
1205 label ptInFace = curFace.
which(curPoint);
1206 label next = curFace.
nextLabel(ptInFace);
1207 label prev = curFace.
prevLabel(ptInFace);
1213 if (corrPoints.
find(next) == -1)
1228 vector c0 = displacement[p1];
1230 scalar V0 =
mag((
a0^b0)&c0)/2;
1232 scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0;
1234 if (corrPoints.
find(prev) != -1)
1238 (
points[next] + displacement[next])
1240 const vector&
c = pointsDisplacementDir()[curPoint];
1242 curDisp += 2*DV/((a^
b)&
c);
1247 - (
points[prev] + displacement[prev]);
1249 const vector&
c = pointsDisplacementDir()[curPoint];
1251 curDisp += 2*DV/((a^
b)&
c);
1255 curDisp /= curPointFaces.size();
1257 displacement[curPoint] =
1258 curDisp*pointsDisplacementDir()[curPoint];
1264 void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals()
1272 aMesh().pointAreaNormals()
1275 if (contactAnglePtr_ && correctContactLineNormals())
1277 Info<<
"Correcting contact line normals" <<
endl;
1281 const labelList& meshPoints = aMesh().patch().meshPoints();
1306 const labelListList& edgeFaces = aMesh().patch().edgeFaces();
1307 const labelListList& pointEdges = aMesh().patch().pointEdges();
1308 const labelListList& pointFaces = aMesh().patch().pointFaces();
1309 const edgeList& edges = aMesh().edges();
1316 == processorFaPatch::typeName
1320 refCast<const processorFaPatch>
1322 aMesh().boundary()[patchI]
1328 forAll(patchPointLabels, pointI)
1330 label curPoint = patchPointLabels[pointI];
1337 const labelList& curPointEdges = pointEdges[curPoint];
1339 forAll(curPointEdges, edgeI)
1341 label curEdge = curPointEdges[edgeI];
1343 if (edgeFaces[curEdge].size() == 1)
1348 aMesh().boundary()[pI];
1350 label index = curEdges.
find(curEdge);
1357 != processorFaPatch::typeName
1374 vector t = edges[curEdge].vec(oldPoints);
1375 t /=
mag(t) + SMALL;
1378 pointFaces[curPoint];
1380 forAll(curPointFaces, fI)
1382 tangent.ref().field()[curPointFaces[fI]] = t;
1389 tangent.correctBoundaryConditions();
1394 label ngbPolyPatchID =
1395 aMesh().boundary()[patchI].ngbPolyPatchIndex();
1397 if (ngbPolyPatchID != -1)
1402 == wallFvPatch::typeName
1410 - contactAnglePtr_->boundaryField()[patchI]
1416 aMesh().
boundary()[patchI].ngbPolyPatchPointNormals()
1420 aMesh().boundary()[patchI].pointLabels();
1425 rotationAxis /=
mag(rotationAxis) + SMALL;
1430 const edgeList& edges = aMesh().edges();
1433 aMesh().boundary()[patchI].pointEdges();
1435 forAll(pointEdges, pointI)
1439 forAll(pointEdges[pointI], eI)
1442 aMesh().boundary()[patchI].start()
1443 + pointEdges[pointI][eI];
1445 vector e = edges[curEdge].vec(oldPoints);
1447 e *= (
e&rotationAxis[pointI])
1448 /
mag(
e&rotationAxis[pointI]);
1450 e /=
mag(
e) + SMALL;
1455 if (pointEdges[pointI].size() == 1)
1457 label curPoint = patchPoints[pointI];
1460 aMesh().patch().pointEdges();
1464 label procPatchID = -1;
1468 aMesh().patch().edgeFaces();
1470 forAll(curPointEdges, edgeI)
1472 label curEdge = curPointEdges[edgeI];
1474 if (edgeFaces[curEdge].size() == 1)
1479 aMesh().boundary()[pI];
1482 curEdges.
find(curEdge);
1489 == processorFaPatch::typeName
1501 if (procPatchID != -1)
1504 tangent.boundaryField()[procPatchID]
1505 .patchNeighbourField()()[
edgeID];
1507 t *= (t&rotationAxis[pointI])
1508 /(
mag(t&rotationAxis[pointI]) + SMALL);
1510 t /=
mag(t) + SMALL;
1516 rotationAxis[pointI] = rotAx/(
mag(rotAx) + SMALL);
1520 ngbN = ngbN*
cos(rotAngle)
1521 + rotationAxis*(rotationAxis & ngbN)*(1 -
cos(rotAngle))
1522 + (rotationAxis^ngbN)*
sin(rotAngle);
1525 forAll(patchPoints, pointI)
1527 N[patchPoints[pointI]] -=
1528 ngbN[pointI]*(ngbN[pointI]&
N[patchPoints[pointI]]);
1530 N[patchPoints[pointI]] /=
1531 mag(
N[patchPoints[pointI]]) + SMALL;
1544 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1553 fixedFreeSurfacePatches_(),
1554 nonReflectingFreeSurfacePatches_(),
1555 pointNormalsCorrectionPatches_(),
1556 normalMotionDir_(false),
1559 pureFreeSurface_(true),
1560 rigidFreeSurface_(false),
1561 correctContactLineNormals_(false),
1566 controlPointsPtr_(nullptr),
1567 motionPointsMaskPtr_(nullptr),
1568 pointsDisplacementDirPtr_(nullptr),
1569 facesDisplacementDirPtr_(nullptr),
1570 fsNetPhiPtr_(nullptr),
1572 surfactConcPtr_(nullptr),
1573 bulkSurfactConcPtr_(nullptr),
1574 surfaceTensionPtr_(nullptr),
1575 surfactantPtr_(nullptr),
1576 contactAnglePtr_(nullptr)
1659 aMeshPtr_.reset(
new faMesh(*
this));
1662 fixedFreeSurfacePatches_ =
1663 motion().get<
wordList>(
"fixedFreeSurfacePatches");
1665 pointNormalsCorrectionPatches_ =
1666 motion().
get<
wordList>(
"pointNormalsCorrectionPatches");
1668 normalMotionDir_ = motion().
get<
bool>(
"normalMotionDir");
1669 smoothing_ = motion().getOrDefault(
"smoothing",
false);
1670 pureFreeSurface_ = motion().getOrDefault(
"pureFreeSurface",
true);
1707 return *fsNetPhiPtr_;
1718 return *fsNetPhiPtr_;
1724 forAll(
Us().boundaryField(), patchI)
1728 Us().boundaryField()[patchI].
type()
1729 == calculatedFaPatchVectorField::typeName
1734 pUs =
Us().boundaryField()[patchI].patchInternalField();
1736 label ngbPolyPatchID =
1737 aMesh().boundary()[patchI].ngbPolyPatchIndex();
1739 if (ngbPolyPatchID != -1)
1744 U().boundaryField()[ngbPolyPatchID].
type()
1745 == slipFvPatchVectorField::typeName
1749 U().boundaryField()[ngbPolyPatchID].
type()
1750 == symmetryFvPatchVectorField::typeName
1756 aMesh().
boundary()[patchI].ngbPolyPatchFaceNormals()
1765 Us().correctBoundaryConditions();
1773 Us().ref().field() =
U().boundaryField()[fsPatchIndex()];
1781 correctUsBoundaryConditions();
1787 return *getObjectPtr<const volVectorField>(
"U");
1793 return *getObjectPtr<const volScalarField>(
"p");
1799 return *getObjectPtr<const surfaceScalarField>(
"phi");
1807 auto& SnGradU = tSnGradU.ref();
1809 const vectorField& nA = aMesh().faceAreaNormals().internalField();
1814 - aMesh().faceCurvatures()*(aMesh().faceAreaNormals()&
Us())
1821 gradUs -=
n*(
n & gradUs);
1822 gradUs.correctBoundaryConditions();
1831 if (!pureFreeSurface() &&
max(
nu) > SMALL)
1833 tangentialSurfaceTensionForce =
1834 surfaceTensionGrad()().internalField();
1838 tangentialSurfaceTensionForce/(
nu + SMALL)
1839 - nA*divUs.internalField()
1840 - (gradUs.internalField()&nA);
1850 auto& SnGradUn = tSnGradUn.ref();
1855 - aMesh().faceCurvatures()*(aMesh().faceAreaNormals()&
Us())
1868 auto& pressureJump = tPressureJump.ref();
1870 const scalarField&
K = aMesh().faceCurvatures().internalField();
1882 + 2.0*
nu*freeSurfaceSnGradUn();
1884 if (pureFreeSurface())
1890 pressureJump -= surfaceTension().internalField()*
K;
1893 return tPressureJump;
1899 if (!controlPointsPtr_)
1901 makeControlPoints();
1904 return *controlPointsPtr_;
1910 if (!motionPointsMaskPtr_)
1912 makeMotionPointsMask();
1915 return *motionPointsMaskPtr_;
1921 if (!pointsDisplacementDirPtr_)
1926 return *pointsDisplacementDirPtr_;
1932 if (!facesDisplacementDirPtr_)
1937 return *facesDisplacementDirPtr_;
1954 if (!surfactConcPtr_)
1959 return *surfactConcPtr_;
1966 if (!surfactConcPtr_)
1971 return *surfactConcPtr_;
1978 if (!bulkSurfactConcPtr_)
1980 makeBulkSurfactConc();
1983 return *bulkSurfactConcPtr_;
1990 if (!bulkSurfactConcPtr_)
1992 makeBulkSurfactConc();
1995 return *bulkSurfactConcPtr_;
2002 if (!surfaceTensionPtr_)
2004 makeSurfaceTension();
2007 return *surfaceTensionPtr_;
2014 if (!surfaceTensionPtr_)
2016 makeSurfaceTension();
2019 return *surfaceTensionPtr_;
2026 if (!surfactantPtr_)
2031 return *surfactantPtr_;
2042 "surfaceTensionGrad",
2054 auto&
grad = tgrad.ref();
2059 grad.correctBoundaryConditions();
2069 if (smoothing_ && !rigidFreeSurface_)
2071 smoothFreeSurfaceMesh();
2072 clearControlPoints();
2075 updateDisplacementDirections();
2079 Info<<
"Maximal capillary Courant number: " 2080 << maxCourantNumber() <<
endl;
2082 const scalarField&
K = aMesh().faceCurvatures().internalField();
2084 Info<<
"Free surface curvature: min = " <<
gMin(
K)
2090 if (!rigidFreeSurface_)
2094 phi().boundaryField()[fsPatchIndex()];
2106 Info<<
"Free surface continuity error : sum local = " 2107 <<
gSum(
mag(sweptVolCorr)) <<
", global = " <<
gSum(sweptVolCorr)
2111 fsNetPhi().ref().field() = sweptVolCorr;
2117 "ddt(" +
U().
name() +
')' 2155 <<
"Unsupported temporal differencing scheme : " 2161 const vectorField& Nf = aMesh().faceAreaNormals().internalField();
2165 sweptVolCorr/(Sf*(Nf & facesDisplacementDir()))
2168 for (
const word& patchName : fixedFreeSurfacePatches_)
2170 label fixedPatchID =
2171 aMesh().boundary().findPatchID(patchName);
2173 if (fixedPatchID == -1)
2176 <<
"Wrong faPatch name '" << patchName
2177 <<
"'in the fixedFreeSurfacePatches list" 2178 <<
" defined in dynamicMeshDict dictionary" 2183 aMesh().boundary()[fixedPatchID].edgeFaces();
2187 deltaHf[eFaces[edgeI]] *= 2.0;
2191 controlPoints() += facesDisplacementDir()*deltaHf;
2193 pointField displacement(pointDisplacement());
2194 correctPointDisplacement(sweptVolCorr, displacement);
2197 refCast<velocityMotionSolver>
2205 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2206 refCast<fixedValuePointPatchVectorField>
2218 correctContactLinePointNormals();
2229 refCast<velocityMotionSolver>
2237 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2238 refCast<fixedValuePointPatchVectorField>
2253 updateSurfactantConcentration();
2265 mesh().time().timePath()/
"freeSurface",
2277 mesh().time().timePath()/
"freeSurfaceControlPoints.vtk" 2280 Info<<
"Writing free surface control points to " <<
os.
name() <<
nl;
2282 os <<
"# vtk DataFile Version 2.0" <<
nl 2283 <<
"freeSurfaceControlPoints" <<
nl 2285 <<
"DATASET POLYDATA" <<
nl;
2287 const label
nPoints = controlPoints().size();
2290 for (
const point&
p : controlPoints())
2292 os << float(
p.x()) <<
' ' 2293 <<
float(
p.y()) <<
' ' 2294 <<
float(
p.z()) <<
nl;
2298 for (label
id = 0;
id <
nPoints; ++id)
2300 os << 1 <<
' ' <<
id <<
nl;
edgeScalarField & Phis()
Return free-surface fluid flux field.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
Finite area mesh (used for 2-D non-Euclidian finite area method) defined using a patch of faces on a ...
Write concrete PrimitivePatch faces/points (optionally with fields) as a vtp file or a legacy vtk fil...
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix by subtracting the matrix multiplied by the current fi...
const Type & value() const noexcept
Return const reference to value.
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
void size(const label n)
Older name for setAddressableSize.
const Type & lookupObject(const word &name, const bool recursive=false) const
Lookup and return const reference to the object of the given Type. Fatal if not found or the wrong ty...
vectorField & facesDisplacementDir()
Return reference to control points displacement direction field.
scalar deltaTValue() const noexcept
Return time step value.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
CGAL::indexedCell< K > Cb
vtk::lineWriter writer(edgeCentres, edgeList::null(), fileName(aMesh.time().globalPath()/"finiteArea-edgesCentres"))
const areaScalarField & fsNetPhi() const
Return free-surface net flux.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
Type gMin(const FieldField< Field, Type > &f)
virtual const fileName & name() const override
Read/write access to the name of the stream.
Virtual base class for velocity motion solver.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
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.
Unit conversion functions.
const Internal & internalField() const noexcept
Return a const-reference to the dimensioned internal field.
Output to file stream, using an OSstream.
tmp< scalarField > freeSurfacePressureJump()
Return free surface pressure jump.
constexpr char nl
The newline '\n' character (0x0a)
vectorField & pointsDisplacementDir()
Return reference to point displacement direction field.
tmp< GeometricField< Type, faPatchField, areaMesh > > div(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
virtual bool update()
Update the mesh for both mesh motion and topology change.
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
void updateUs()
Update free-surface velocity field.
const labelList & pointLabels() const
Return patch point labels.
static bool & parRun() noexcept
Test if this a parallel run.
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep track of the Delaunay cells (tets)...
label nextLabel(const label i) const
Next vertex on face.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type get(const label i) const
Return bool value at specified position, always false for out-of-range access.
const volVectorField & U() const
Return constant reference to velocity field.
Virtual base class for mesh motion solver.
void correctUsBoundaryConditions()
Correct surface velocity boundary conditions.
Ignore writing from objectRegistry::writeObject()
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
const dimensionSet dimless
Dimensionless.
const Time & time() const
Return the top-level database.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
CGAL::Exact_predicates_exact_constructions_kernel K
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Abstract base class for turbulence models (RAS, LES and laminar).
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
faMesh & aMesh()
Return reference to finite area mesh.
dimensionedScalar neg(const dimensionedScalar &ds)
static const word & zeroGradientType() noexcept
The type name for zeroGradient patch fields.
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
tmp< faMatrix< Type > > laplacian(const GeometricField< Type, faPatchField, areaMesh > &vf)
label prevLabel(const label i) const
Previous vertex on face.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
volScalarField & bulkSurfactantConcentration()
Return volume surfactant concentration field.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
const dimensionSet dimVolume(pow3(dimLength))
zeroField Sp(const Foam::zero, const GeometricField< Type, faPatchField, areaMesh > &)
A no-op source.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
label findPatchID(const word &patchName, const bool allowNotFound=true) const
Find patch index given a name, return -1 if not found.
tmp< GeometricField< Type, faePatchField, edgeMesh > > linearEdgeInterpolate(const GeometricField< Type, faPatchField, areaMesh > &vf)
const dimensionedScalar e
Elementary charge.
Type gSum(const FieldField< Field, Type > &f)
dimensionedScalar cos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
* ~interfaceTrackingFvMesh()
Destructor.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
#define DebugInFunction
Report an information message using Foam::Info.
areaVectorField & Us()
Return free-surface velocity field.
virtual const pointField & oldPoints() const
Return old points (mesh motion)
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
label which(const label vertex) const
Find local vertex on face for the vertex label, same as find()
label timeIndex() const noexcept
Return the current time index.
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Second-order backward-differencing ddt using the current and two previous time-step values...
errorManip< error > abort(error &err)
vectorField & controlPoints()
Return control points.
label find(const T &val) const
Find index of the first occurrence of the value.
virtual bool update()
Update the mesh for both mesh motion and topology change.
areaScalarField & surfactantConcentration()
Return free-surface surfactant concentration field.
tmp< faMatrix< Type > > ddt(const GeometricField< Type, faPatchField, areaMesh > &vf)
const uniformDimensionedVectorField & g
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimForce
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
tmp< faMatrix< Type > > div(const edgeScalarField &flux, const GeometricField< Type, faPatchField, areaMesh > &vf, const word &name)
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val. FatalIOError if it is found and the number of tokens i...
Type gMax(const FieldField< Field, Type > &f)
OBJstream os(runTime.globalPath()/outputName)
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0)
Info<< "Reading field U\"<< 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\"<< endl;autoPtr< compressible::turbulenceModel > turbulence(compressible::turbulenceModel::New(rho, U, phi, thermo))
defineTypeNameAndDebug(combustionModel, 0)
const volScalarField & p() const
Return constant reference to pressure field.
static const gravity & New(const word &name, const Time &runTime)
Return named gravity field cached or construct on Time.
tmp< GeometricField< Type, faPatchField, areaMesh > > edgeIntegrate(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
const faBoundaryMesh & boundary() const noexcept
Return constant reference to boundary mesh.
const surfaceScalarField & phi() const
Return constant reference to velocity field.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
labelList & motionPointsMask()
Return reference to motion points mask field.
const dimensionSet dimDensity
const Vector< label > N(dict.get< Vector< label >>("N"))
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
const surfactantProperties & surfactant() const
Return surfactant properties.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Template functions to aid in the implementation of demand driven data.
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field values.
Ostream & flush(Ostream &os)
Flush stream.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Type gAverage(const FieldField< Field, Type > &f)
Abstract base class for geometry and/or topology changing fvMesh.
const dimensionedScalar c
Speed of light in a vacuum.
Automatically write from objectRegistry::writeObject()
void writeVTKControlPoints()
Write VTK freeSurface control points.
const dimensionedScalar a0
Bohr radius: default SI units: [m].
Template specialisation for scalar faMatrix.
const std::string patch
OpenFOAM patch number as a std::string.
const motionSolver & motion() const
Return the motionSolver.
tmp< scalarField > freeSurfaceSnGradUn()
Return free surface normal derivative of normal velocity comp.
void correctBoundaryConditions()
Correct boundary field.
Legacy ASCII, legacyAsciiFormatter.
bool correctPatchPointNormals(const label patchID) const
Whether point normals should be corrected for a patch.
messageStream Info
Information stream (stdout output on master, null elsewhere)
virtual bool writeGeometry()
Write patch topology.
Field< vector > vectorField
Specialisation of Field<T> for vector.
tmp< areaVectorField > surfaceTensionGrad()
Return surface tension gradient.
The dynamicMotionSolverFvMesh.
List< Key > toc() const
The table of contents (the keys) in unsorted order.
void writeVTK() const
Write VTK freeSurface mesh.
List< label > labelList
A List of labels.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, IOobject::NO_REGISTER)
A class for managing temporary objects.
tmp< vectorField > freeSurfaceSnGradU()
Return free surface normal derivative of velocity.
Basic first-order Euler implicit/explicit ddt using only the current and previous time-step values...
A simple single-phase transport model based on viscosityModel.
void deleteDemandDrivenData(DataPtr &dataPtr)
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Defines the attributes of an object for which implicit objectRegistry management is supported...
GeometricField< vector, faPatchField, areaMesh > areaVectorField
A primitive field of type <T> with automated input and output.
IOField< vector > vectorIOField
IO for a Field of vector.
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
const volScalarField & p0
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
const dimensionSet dimArea(sqr(dimLength))
tmp< GeometricField< typename outerProduct< vector, Type >::type, faPatchField, areaMesh >> grad(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
areaScalarField & surfaceTension()
Return surface tension field.
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
The interfaceTrackingFvMesh.
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVelocity