Go to the documentation of this file.
68 void Foam::interfaceTrackingFvMesh::initializeData()
72 const word fsPatchName(
motion().get<word>(
"fsPatchName"));
79 <<
"Patch name " << fsPatchName <<
" not found."
83 fsPatchIndex_ =
patch.index();
90 for (
const word& patchName : pointNormalsCorrectionPatches_)
97 <<
"Patch name '" << patchName
98 <<
"' for point normals correction does not exist"
107 if (!normalMotionDir_)
117 "nonReflectingFreeSurfacePatches",
118 nonReflectingFreeSurfacePatches_
123 void Foam::interfaceTrackingFvMesh::makeUs()
const
126 <<
"making free-surface velocity field" <<
nl;
131 <<
"free-surface velocity field already exists"
138 zeroGradientFaPatchVectorField::typeName
146 == wedgeFaPatch::typeName
149 patchFieldTypes[patchI] =
150 wedgeFaPatchVectorField::typeName;
154 label ngbPolyPatchID =
157 if (ngbPolyPatchID != -1)
162 == wallFvPatch::typeName
165 patchFieldTypes[patchI] =
166 slipFaPatchVectorField::typeName;
172 for (
const word& patchName : fixedFreeSurfacePatches_)
174 const label fixedPatchID =
177 if (fixedPatchID == -1)
180 <<
"Wrong faPatch name '" << patchName
181 <<
"' in the fixedFreeSurfacePatches list"
182 <<
" defined in the dynamicMeshDict dictionary"
186 label ngbPolyPatchID =
189 if (ngbPolyPatchID != -1)
194 == wallFvPatch::typeName
197 patchFieldTypes[fixedPatchID] =
198 fixedValueFaPatchVectorField::typeName;
218 for (
const word& patchName : fixedFreeSurfacePatches_)
222 if (fixedPatchID == -1)
225 <<
"Wrong faPatch name '" << patchName
226 <<
"' in the fixedFreeSurfacePatches list"
227 <<
" defined in the dynamicMeshDict dictionary"
231 label ngbPolyPatchID =
234 if (ngbPolyPatchID != -1)
239 == wallFvPatch::typeName
242 UsPtr_->boundaryFieldRef()[fixedPatchID] ==
Zero;
249 void Foam::interfaceTrackingFvMesh::makeFsNetPhi()
const
252 <<
"making free-surface net flux" <<
nl;
257 <<
"free-surface net flux already exists"
277 void Foam::interfaceTrackingFvMesh::makeControlPoints()
280 <<
"making control points" <<
nl;
282 if (controlPointsPtr_)
285 <<
"control points already exists"
299 Info<<
"Reading control points" <<
endl;
315 Info<<
"Creating new control points" <<
endl;
327 aMesh().areaCentres().internalField()
330 initializeControlPointsPosition();
335 void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
338 <<
"making motion points mask" <<
nl;
340 if (motionPointsMaskPtr_)
343 <<
"motion points mask already exists"
360 == processorFaPatch::typeName
366 forAll(patchPoints, pointI)
368 motionPointsMask()[patchPoints[pointI]] = -1;
374 for (
const word& patchName : fixedFreeSurfacePatches_)
378 if (fixedPatchID == -1)
381 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
382 <<
" defined in the dynamicMeshDict dictionary"
389 forAll(patchPoints, pointI)
391 motionPointsMask()[patchPoints[pointI]] = 0;
397 void Foam::interfaceTrackingFvMesh::makeDirections()
400 <<
"make displacement directions for points and control points" <<
nl;
402 if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_)
405 <<
"points, control points displacement directions already exist"
409 pointsDisplacementDirPtr_ =
416 facesDisplacementDirPtr_ =
423 if (!normalMotionDir())
425 if (
mag(motionDir_) < SMALL)
428 <<
"Zero motion direction"
432 facesDisplacementDir() = motionDir_;
433 pointsDisplacementDir() = motionDir_;
436 updateDisplacementDirections();
440 void Foam::interfaceTrackingFvMesh::makePhis()
443 <<
"making free-surface flux" <<
nl;
448 <<
"free-surface flux already exists"
467 void Foam::interfaceTrackingFvMesh::makeSurfactConc()
const
470 <<
"making free-surface surfactant concentration field" <<
nl;
475 <<
"free-surface surfactant concentration field already exists"
498 void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc()
const
501 <<
"making volume surfactant concentration field" <<
nl;
503 if (bulkSurfactConcPtr_)
506 <<
"volume surfactant concentration field already exists"
531 bulkSurfactConc = surfactant().bulkConc();
537 void Foam::interfaceTrackingFvMesh::makeSurfaceTension()
const
540 <<
"making surface tension field" <<
nl;
542 if (surfaceTensionPtr_)
545 <<
"surface tension field already exists"
559 sigma() + surfactant().dSigma(surfactantConcentration())/rho_
564 void Foam::interfaceTrackingFvMesh::makeSurfactant()
const
567 <<
"making surfactant properties" <<
nl;
572 <<
"surfactant properties already exists"
577 motion().
subDict(
"surfactantProperties");
583 void Foam::interfaceTrackingFvMesh::makeContactAngle()
586 <<
"making contact angle field" <<
nl;
588 if (contactAnglePtr_)
591 <<
"contact angle already exists"
606 Info<<
"Reading contact angle field" <<
endl;
625 void Foam::interfaceTrackingFvMesh::updateDisplacementDirections()
627 if (normalMotionDir())
635 if (contactAnglePtr_)
637 label ngbPolyPatchID =
640 if (ngbPolyPatchID != -1)
645 == wallFvPatch::typeName
654 .ngbPolyPatchPointNormals()
657 forAll(patchPoints, pointI)
659 pointsDisplacementDir()[patchPoints[pointI]] -=
663 & pointsDisplacementDir()[patchPoints[pointI]]
666 pointsDisplacementDir()[patchPoints[pointI]] /=
669 pointsDisplacementDir()
681 facesDisplacementDir() =
688 facesDisplacementDir()
689 *(facesDisplacementDir()&(controlPoints() - Cf))
695 void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition()
703 correctPointDisplacement(sweptVolCorr, displacement);
711 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
718 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
725 deltaH[faceI] = sweptVol[faceI]/
726 ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL);
728 if (
mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL)
735 <<
"Something is wrong with specified motion direction"
740 for (
const word& patchName : fixedFreeSurfacePatches_)
744 if (fixedPatchID == -1)
747 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
748 <<
" defined in the freeSurfaceProperties dictionary"
757 deltaH[eFaces[edgeI]] *= 2.0;
761 controlPoints() += facesDisplacementDir()*deltaH;
766 void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh()
768 Info<<
"Smoothing free surface mesh" <<
endl;
782 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
788 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
793 sweptVol/(faceArea & facesDisplacementDir())
796 for (
const word& patchName : fixedFreeSurfacePatches_)
800 if (fixedPatchID == -1)
803 <<
"Wrong faPatch name fixedFreeSurfacePatches list"
804 <<
" defined in the dynamicMeshDict dictionary"
813 deltaHf[eFaces[edgeI]] *= 2.0;
817 controlPoints() += facesDisplacementDir()*deltaHf;
819 displacement = pointDisplacement();
822 refCast<velocityMotionSolver>
830 fixedValuePointPatchVectorField& fsPatchPointMeshU =
831 refCast<fixedValuePointPatchVectorField>
846 void Foam::interfaceTrackingFvMesh::updateSurfaceFlux()
852 void Foam::interfaceTrackingFvMesh::updateSurfactantConcentration()
854 if (!pureFreeSurface())
864 +
fam::div(Phis(), surfactantConcentration())
867 surfactant().diffusion(),
868 surfactantConcentration()
872 if (surfactant().soluble())
888 zeroGradientFaPatchScalarField::typeName
892 bulkSurfactantConcentration().boundaryField()[fsPatchIndex()];
893 Cb.correctBoundaryConditions();
898 surfactant().adsorptionCoeff()*
Cb
899 + surfactant().adsorptionCoeff()
900 *surfactant().desorptionCoeff(),
901 surfactantConcentration()
903 - surfactant().adsorptionCoeff()
904 *
Cb*surfactant().saturatedConc();
912 sigma() + surfactant().dSigma(surfactantConcentration())/rho_;
914 if (
neg(
min(surfaceTension().internalField().
field())))
917 <<
"Surface tension is negative"
924 Foam::vector Foam::interfaceTrackingFvMesh::totalPressureForce()
const
934 return gSum(pressureForces);
938 Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce()
const
958 U().boundaryField()[fsPatchIndex()].
snGrad()
971 return gSum(viscousForces);
975 Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce()
const
985 if (pureFreeSurface())
996 surfTensionForces = surfaceTension().internalField().field()*
K*S*
n;
999 return gSum(surfTensionForces);
1003 Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber()
1007 if (pureFreeSurface())
1013 mesh().time().deltaT().value()/
1032 mesh().time().deltaT().value()/
1044 void Foam::interfaceTrackingFvMesh::updateProperties()
1049 "transportProperties"
1058 void Foam::interfaceTrackingFvMesh::correctPointDisplacement
1073 for (
const word& patchName : fixedFreeSurfacePatches_)
1087 label curFace = eFaces[edgeI];
1089 const labelList& curPoints = faces[curFace];
1091 forAll(curPoints, pointI)
1093 label curPoint = curPoints[pointI];
1094 label index = pLabels.find(curPoint);
1110 forAll(corrPoints, pointI)
1112 label curPoint = corrPoints[pointI];
1118 label curFace =
pFaces[curPoint][faceI];
1120 label index = eFaces.find(curFace);
1134 forAll(corrPoints, pointI)
1136 label curPoint = corrPoints[pointI];
1140 const labelList& curPointFaces = corrPointFaces[pointI];
1142 forAll(curPointFaces, faceI)
1144 const face& curFace = faces[curPointFaces[faceI]];
1146 label ptInFace = curFace.
which(curPoint);
1147 label next = curFace.
nextLabel(ptInFace);
1148 label prev = curFace.
prevLabel(ptInFace);
1152 const vector&
c = pointsDisplacementDir()[curPoint];
1154 curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^
b)&
c);
1157 curDisp /= curPointFaces.size();
1159 displacement[curPoint] =
1160 curDisp*pointsDisplacementDir()[curPoint];
1165 for (
const word& patchName : nonReflectingFreeSurfacePatches_)
1167 label nonReflectingPatchID =
1180 forAll(corrPoints, pointI)
1182 label curPoint = corrPoints[pointI];
1188 label curFace =
pFaces[curPoint][faceI];
1190 label index = eFaces.find(curFace);
1205 forAll(corrPoints, pointI)
1207 label curPoint = corrPoints[pointI];
1211 const labelList& curPointFaces = corrPointFaces[pointI];
1213 forAll(curPointFaces, faceI)
1215 const face& curFace = faces[curPointFaces[faceI]];
1217 label ptInFace = curFace.
which(curPoint);
1218 label next = curFace.
nextLabel(ptInFace);
1219 label prev = curFace.
prevLabel(ptInFace);
1225 if (corrPoints.find(next) == -1)
1240 vector c0 = displacement[p1];
1242 scalar V0 =
mag((
a0^b0)&c0)/2;
1244 scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0;
1246 if (corrPoints.find(prev) != -1)
1250 (
points[next] + displacement[next])
1252 const vector&
c = pointsDisplacementDir()[curPoint];
1254 curDisp += 2*DV/((a^
b)&
c);
1259 - (
points[prev] + displacement[prev]);
1261 const vector&
c = pointsDisplacementDir()[curPoint];
1263 curDisp += 2*DV/((a^
b)&
c);
1267 curDisp /= curPointFaces.size();
1269 displacement[curPoint] =
1270 curDisp*pointsDisplacementDir()[curPoint];
1276 void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals()
1287 if (contactAnglePtr_ && correctContactLineNormals())
1289 Info<<
"Correcting contact line normals" <<
endl;
1328 == processorFaPatch::typeName
1332 refCast<const processorFaPatch>
1340 forAll(patchPointLabels, pointI)
1342 label curPoint = patchPointLabels[pointI];
1349 const labelList& curPointEdges = pointEdges[curPoint];
1351 forAll(curPointEdges, edgeI)
1353 label curEdge = curPointEdges[edgeI];
1355 if (edgeFaces[curEdge].size() == 1)
1362 label index = curEdges.find(curEdge);
1369 != processorFaPatch::typeName
1386 vector t = edges[curEdge].vec(oldPoints);
1387 t /=
mag(t) + SMALL;
1390 pointFaces[curPoint];
1392 forAll(curPointFaces, fI)
1394 tangent.ref().field()[curPointFaces[fI]] = t;
1401 tangent.correctBoundaryConditions();
1406 label ngbPolyPatchID =
1409 if (ngbPolyPatchID != -1)
1414 == wallFvPatch::typeName
1422 - contactAnglePtr_->boundaryField()[patchI]
1437 rotationAxis /=
mag(rotationAxis) + SMALL;
1447 forAll(pointEdges, pointI)
1451 forAll(pointEdges[pointI], eI)
1455 + pointEdges[pointI][eI];
1457 vector e = edges[curEdge].vec(oldPoints);
1459 e *= (
e&rotationAxis[pointI])
1460 /
mag(
e&rotationAxis[pointI]);
1462 e /=
mag(
e) + SMALL;
1467 if (pointEdges[pointI].size() == 1)
1469 label curPoint = patchPoints[pointI];
1476 label procPatchID = -1;
1482 forAll(curPointEdges, edgeI)
1484 label curEdge = curPointEdges[edgeI];
1486 if (edgeFaces[curEdge].size() == 1)
1494 curEdges.find(curEdge);
1501 == processorFaPatch::typeName
1513 if (procPatchID != -1)
1516 tangent.boundaryField()[procPatchID]
1517 .patchNeighbourField()()[
edgeID];
1519 t *= (t&rotationAxis[pointI])
1520 /(
mag(t&rotationAxis[pointI]) + SMALL);
1522 t /=
mag(t) + SMALL;
1528 rotationAxis[pointI] = rotAx/(
mag(rotAx) + SMALL);
1532 ngbN = ngbN*
cos(rotAngle)
1533 + rotationAxis*(rotationAxis & ngbN)*(1 -
cos(rotAngle))
1534 + (rotationAxis^ngbN)*
sin(rotAngle);
1537 forAll(patchPoints, pointI)
1539 N[patchPoints[pointI]] -=
1540 ngbN[pointI]*(ngbN[pointI]&
N[patchPoints[pointI]]);
1542 N[patchPoints[pointI]] /=
1543 mag(
N[patchPoints[pointI]]) + SMALL;
1556 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh(
const IOobject& io)
1559 aMeshPtr_(new
faMesh(*this)),
1561 fixedFreeSurfacePatches_
1563 motion().get<
wordList>(
"fixedFreeSurfacePatches")
1565 nonReflectingFreeSurfacePatches_(),
1566 pointNormalsCorrectionPatches_
1568 motion().get<
wordList>(
"pointNormalsCorrectionPatches")
1572 motion().get<
bool>(
"normalMotionDir")
1575 smoothing_(motion().getOrDefault(
"smoothing", false)),
1576 pureFreeSurface_(motion().getOrDefault(
"pureFreeSurface", true)),
1577 rigidFreeSurface_(motion().getOrDefault(
"rigidFreeSurface", false)),
1578 correctContactLineNormals_
1580 motion().getOrDefault(
"correctContactLineNormals", false)
1586 controlPointsPtr_(nullptr),
1587 motionPointsMaskPtr_(nullptr),
1588 pointsDisplacementDirPtr_(nullptr),
1589 facesDisplacementDirPtr_(nullptr),
1590 fsNetPhiPtr_(nullptr),
1592 surfactConcPtr_(nullptr),
1593 bulkSurfactConcPtr_(nullptr),
1594 surfaceTensionPtr_(nullptr),
1595 surfactantPtr_(nullptr),
1596 contactAnglePtr_(nullptr)
1602 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1617 std::move(allOwner),
1618 std::move(allNeighbour),
1621 aMeshPtr_(
new faMesh(*
this)),
1623 fixedFreeSurfacePatches_
1625 motion().get<wordList>(
"fixedFreeSurfacePatches")
1627 nonReflectingFreeSurfacePatches_(),
1628 pointNormalsCorrectionPatches_
1630 motion().get<wordList>(
"pointNormalsCorrectionPatches")
1634 motion().get<bool>(
"normalMotionDir")
1637 smoothing_(motion().getOrDefault(
"smoothing",
false)),
1638 pureFreeSurface_(motion().getOrDefault(
"pureFreeSurface",
true)),
1643 controlPointsPtr_(
nullptr),
1644 motionPointsMaskPtr_(
nullptr),
1645 pointsDisplacementDirPtr_(
nullptr),
1646 facesDisplacementDirPtr_(
nullptr),
1647 fsNetPhiPtr_(
nullptr),
1649 surfactConcPtr_(
nullptr),
1650 bulkSurfactConcPtr_(
nullptr),
1651 surfaceTensionPtr_(
nullptr),
1652 surfactantPtr_(
nullptr),
1653 contactAnglePtr_(
nullptr)
1709 return *fsNetPhiPtr_;
1720 return *fsNetPhiPtr_;
1726 forAll(
Us().boundaryField(), patchI)
1730 Us().boundaryField()[patchI].
type()
1731 == calculatedFaPatchVectorField::typeName
1736 pUs =
Us().boundaryField()[patchI].patchInternalField();
1738 label ngbPolyPatchID =
1741 if (ngbPolyPatchID != -1)
1746 U().boundaryField()[ngbPolyPatchID].
type()
1747 == slipFvPatchVectorField::typeName
1751 U().boundaryField()[ngbPolyPatchID].
type()
1752 == symmetryFvPatchVectorField::typeName
1767 Us().correctBoundaryConditions();
1775 Us().ref().field() =
U().boundaryField()[fsPatchIndex()];
1783 correctUsBoundaryConditions();
1789 return *getObjectPtr<const volVectorField>(
"U");
1795 return *getObjectPtr<const volScalarField>(
"p");
1801 return *getObjectPtr<const surfaceScalarField>(
"phi");
1809 auto& SnGradU = tSnGradU.ref();
1816 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1823 gradUs -=
n*(
n & gradUs);
1833 if (!pureFreeSurface() &&
max(
nu) > SMALL)
1835 tangentialSurfaceTensionForce =
1836 surfaceTensionGrad()().internalField();
1840 tangentialSurfaceTensionForce/(
nu + SMALL)
1852 auto& SnGradUn = tSnGradUn.ref();
1857 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1870 auto& pressureJump = tPressureJump.ref();
1884 + 2.0*
nu*freeSurfaceSnGradUn();
1886 if (pureFreeSurface())
1892 pressureJump -= surfaceTension().internalField()*
K;
1895 return tPressureJump;
1901 if (!controlPointsPtr_)
1903 makeControlPoints();
1906 return *controlPointsPtr_;
1912 if (!motionPointsMaskPtr_)
1914 makeMotionPointsMask();
1917 return *motionPointsMaskPtr_;
1923 if (!pointsDisplacementDirPtr_)
1928 return *pointsDisplacementDirPtr_;
1934 if (!facesDisplacementDirPtr_)
1939 return *facesDisplacementDirPtr_;
1956 if (!surfactConcPtr_)
1961 return *surfactConcPtr_;
1968 if (!surfactConcPtr_)
1973 return *surfactConcPtr_;
1980 if (!bulkSurfactConcPtr_)
1982 makeBulkSurfactConc();
1985 return *bulkSurfactConcPtr_;
1992 if (!bulkSurfactConcPtr_)
1994 makeBulkSurfactConc();
1997 return *bulkSurfactConcPtr_;
2004 if (!surfaceTensionPtr_)
2006 makeSurfaceTension();
2009 return *surfaceTensionPtr_;
2016 if (!surfaceTensionPtr_)
2018 makeSurfaceTension();
2021 return *surfaceTensionPtr_;
2028 if (!surfactantPtr_)
2033 return *surfactantPtr_;
2044 "surfaceTensionGrad",
2056 auto&
grad = tgrad.ref();
2061 grad.correctBoundaryConditions();
2071 if (smoothing_ && !rigidFreeSurface_)
2073 smoothFreeSurfaceMesh();
2074 clearControlPoints();
2077 updateDisplacementDirections();
2081 Info<<
"Maximal capillary Courant number: "
2082 << maxCourantNumber() <<
endl;
2086 Info<<
"Free surface curvature: min = " <<
gMin(
K)
2092 if (!rigidFreeSurface_)
2108 Info<<
"Free surface continuity error : sum local = "
2109 <<
gSum(
mag(sweptVolCorr)) <<
", global = " <<
gSum(sweptVolCorr)
2113 fsNetPhi().ref().field() = sweptVolCorr;
2119 "ddt(" +
U().
name() +
')'
2157 <<
"Unsupported temporal differencing scheme : "
2167 sweptVolCorr/(Sf*(Nf & facesDisplacementDir()))
2170 for (
const word& patchName : fixedFreeSurfacePatches_)
2172 label fixedPatchID =
2175 if (fixedPatchID == -1)
2178 <<
"Wrong faPatch name '" << patchName
2179 <<
"'in the fixedFreeSurfacePatches list"
2180 <<
" defined in dynamicMeshDict dictionary"
2189 deltaHf[eFaces[edgeI]] *= 2.0;
2193 controlPoints() += facesDisplacementDir()*deltaHf;
2195 pointField displacement(pointDisplacement());
2196 correctPointDisplacement(sweptVolCorr, displacement);
2199 refCast<velocityMotionSolver>
2207 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2208 refCast<fixedValuePointPatchVectorField>
2216 fsPatchPointMeshU ==
2221 correctContactLinePointNormals();
2232 refCast<velocityMotionSolver>
2240 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2241 refCast<fixedValuePointPatchVectorField>
2249 fsPatchPointMeshU ==
2257 updateSurfactantConcentration();
2266 OFstream mps(
mesh().time().timePath()/
"freeSurface.vtk");
2271 mps <<
"# vtk DataFile Version 2.0" <<
nl
2274 <<
"DATASET POLYDATA" <<
nl
2275 <<
"POINTS " <<
points.size() <<
" float" <<
nl;
2283 mlpBuffer[counter++] = float(
points[i].
x());
2284 mlpBuffer[counter++] = float(
points[i].
y());
2285 mlpBuffer[counter++] = float(
points[i].z());
2290 mps << mlpBuffer[i] <<
' ';
2292 if (i > 0 && (i % 10) == 0)
2299 label nFaceVerts = 0;
2303 nFaceVerts += faces[faceI].size() + 1;
2310 const face&
f = faces[faceI];
2312 mlfBuffer[counter++] =
f.size();
2316 mlfBuffer[counter++] =
f[fpI];
2321 mps <<
"POLYGONS " << faces.size() <<
' ' << nFaceVerts <<
endl;
2325 mps << mlfBuffer[i] <<
' ';
2327 if (i > 0 && (i % 10) == 0)
2349 Info<<
"Writing free surface control point to " <<
name <<
endl;
2351 mps <<
"# vtk DataFile Version 2.0" <<
nl
2354 <<
"DATASET POLYDATA" <<
nl
2355 <<
"POINTS " << controlPoints().size() <<
" float" <<
nl;
2357 forAll(controlPoints(), pointI)
2359 mps << controlPoints()[pointI].x() <<
' '
2360 << controlPoints()[pointI].y() <<
' '
2361 << controlPoints()[pointI].z() <<
nl;
2365 mps <<
"VERTICES " << controlPoints().size() <<
' '
2366 << controlPoints().size()*2 <<
nl;
2368 forAll(controlPoints(), pointI)
2370 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()
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)
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 > &)
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))
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 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.
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)