Go to the documentation of this file.
75 void Foam::interfaceTrackingFvMesh::initializeData()
79 const word fsPatchName(
motion().get<word>(
"fsPatchName"));
86 <<
"Patch name " << fsPatchName <<
" not found."
90 fsPatchIndex_ =
patch.index();
97 for (
const word& patchName : pointNormalsCorrectionPatches_)
104 <<
"Patch name '" << patchName
105 <<
"' for point normals correction does not exist"
114 if (!normalMotionDir_)
124 "nonReflectingFreeSurfacePatches",
125 nonReflectingFreeSurfacePatches_
130 void Foam::interfaceTrackingFvMesh::makeUs()
const
133 <<
"making free-surface velocity field" <<
nl;
138 <<
"free-surface velocity field already exists"
145 zeroGradientFaPatchVectorField::typeName
153 == wedgeFaPatch::typeName
156 patchFieldTypes[patchI] =
157 wedgeFaPatchVectorField::typeName;
161 label ngbPolyPatchID =
164 if (ngbPolyPatchID != -1)
169 == wallFvPatch::typeName
172 patchFieldTypes[patchI] =
173 slipFaPatchVectorField::typeName;
179 for (
const word& patchName : fixedFreeSurfacePatches_)
181 const label fixedPatchID =
184 if (fixedPatchID == -1)
187 <<
"Wrong faPatch name '" << patchName
188 <<
"' in the fixedFreeSurfacePatches list"
189 <<
" defined in the dynamicMeshDict dictionary"
193 label ngbPolyPatchID =
196 if (ngbPolyPatchID != -1)
201 == wallFvPatch::typeName
204 patchFieldTypes[fixedPatchID] =
205 fixedValueFaPatchVectorField::typeName;
225 for (
const word& patchName : fixedFreeSurfacePatches_)
229 if (fixedPatchID == -1)
232 <<
"Wrong faPatch name '" << patchName
233 <<
"' in the fixedFreeSurfacePatches list"
234 <<
" defined in the dynamicMeshDict dictionary"
238 label ngbPolyPatchID =
241 if (ngbPolyPatchID != -1)
246 == wallFvPatch::typeName
249 UsPtr_->boundaryFieldRef()[fixedPatchID] ==
Zero;
256 void Foam::interfaceTrackingFvMesh::makeFsNetPhi()
const
259 <<
"making free-surface net flux" <<
nl;
264 <<
"free-surface net flux already exists"
284 void Foam::interfaceTrackingFvMesh::makeControlPoints()
287 <<
"making control points" <<
nl;
289 if (controlPointsPtr_)
292 <<
"control points already exists"
306 Info<<
"Reading control points" <<
endl;
322 Info<<
"Creating new control points" <<
endl;
334 aMesh().areaCentres().internalField()
337 initializeControlPointsPosition();
342 void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
345 <<
"making motion points mask" <<
nl;
347 if (motionPointsMaskPtr_)
350 <<
"motion points mask already exists"
367 == processorFaPatch::typeName
373 forAll(patchPoints, pointI)
375 motionPointsMask()[patchPoints[pointI]] = -1;
381 for (
const word& patchName : fixedFreeSurfacePatches_)
385 if (fixedPatchID == -1)
388 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
389 <<
" defined in the dynamicMeshDict dictionary"
396 forAll(patchPoints, pointI)
398 motionPointsMask()[patchPoints[pointI]] = 0;
404 void Foam::interfaceTrackingFvMesh::makeDirections()
407 <<
"make displacement directions for points and control points" <<
nl;
409 if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_)
412 <<
"points, control points displacement directions already exist"
416 pointsDisplacementDirPtr_ =
423 facesDisplacementDirPtr_ =
430 if (!normalMotionDir())
432 if (
mag(motionDir_) < SMALL)
435 <<
"Zero motion direction"
439 facesDisplacementDir() = motionDir_;
440 pointsDisplacementDir() = motionDir_;
443 updateDisplacementDirections();
447 void Foam::interfaceTrackingFvMesh::makePhis()
450 <<
"making free-surface flux" <<
nl;
455 <<
"free-surface flux already exists"
474 void Foam::interfaceTrackingFvMesh::makeSurfactConc()
const
477 <<
"making free-surface surfactant concentration field" <<
nl;
482 <<
"free-surface surfactant concentration field already exists"
505 void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc()
const
508 <<
"making volume surfactant concentration field" <<
nl;
510 if (bulkSurfactConcPtr_)
513 <<
"volume surfactant concentration field already exists"
538 bulkSurfactConc = surfactant().bulkConc();
544 void Foam::interfaceTrackingFvMesh::makeSurfaceTension()
const
547 <<
"making surface tension field" <<
nl;
549 if (surfaceTensionPtr_)
552 <<
"surface tension field already exists"
566 sigma() + surfactant().dSigma(surfactantConcentration())/rho_
571 void Foam::interfaceTrackingFvMesh::makeSurfactant()
const
574 <<
"making surfactant properties" <<
nl;
579 <<
"surfactant properties already exists"
584 motion().
subDict(
"surfactantProperties");
590 void Foam::interfaceTrackingFvMesh::makeContactAngle()
593 <<
"making contact angle field" <<
nl;
595 if (contactAnglePtr_)
598 <<
"contact angle already exists"
613 Info<<
"Reading contact angle field" <<
endl;
632 void Foam::interfaceTrackingFvMesh::updateDisplacementDirections()
634 if (normalMotionDir())
642 if (contactAnglePtr_)
644 label ngbPolyPatchID =
647 if (ngbPolyPatchID != -1)
652 == wallFvPatch::typeName
661 .ngbPolyPatchPointNormals()
664 forAll(patchPoints, pointI)
666 pointsDisplacementDir()[patchPoints[pointI]] -=
670 & pointsDisplacementDir()[patchPoints[pointI]]
673 pointsDisplacementDir()[patchPoints[pointI]] /=
676 pointsDisplacementDir()
688 facesDisplacementDir() =
695 facesDisplacementDir()
696 *(facesDisplacementDir()&(controlPoints() - Cf))
702 void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition()
710 correctPointDisplacement(sweptVolCorr, displacement);
718 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
725 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
732 deltaH[faceI] = sweptVol[faceI]/
733 ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL);
735 if (
mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL)
742 <<
"Something is wrong with specified motion direction"
747 for (
const word& patchName : fixedFreeSurfacePatches_)
751 if (fixedPatchID == -1)
754 <<
"Wrong faPatch name in the fixedFreeSurfacePatches list"
755 <<
" defined in the freeSurfaceProperties dictionary"
764 deltaH[eFaces[edgeI]] *= 2.0;
768 controlPoints() += facesDisplacementDir()*deltaH;
773 void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh()
775 Info<<
"Smoothing free surface mesh" <<
endl;
789 sweptVol[faceI] = -faces[faceI].sweptVol(
points, newPoints);
795 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
800 sweptVol/(faceArea & facesDisplacementDir())
803 for (
const word& patchName : fixedFreeSurfacePatches_)
807 if (fixedPatchID == -1)
810 <<
"Wrong faPatch name fixedFreeSurfacePatches list"
811 <<
" defined in the dynamicMeshDict dictionary"
820 deltaHf[eFaces[edgeI]] *= 2.0;
824 controlPoints() += facesDisplacementDir()*deltaHf;
826 displacement = pointDisplacement();
829 refCast<velocityMotionSolver>
837 fixedValuePointPatchVectorField& fsPatchPointMeshU =
838 refCast<fixedValuePointPatchVectorField>
853 void Foam::interfaceTrackingFvMesh::updateSurfaceFlux()
859 void Foam::interfaceTrackingFvMesh::updateSurfactantConcentration()
861 if (!pureFreeSurface())
871 +
fam::div(Phis(), surfactantConcentration())
874 surfactant().diffusion(),
875 surfactantConcentration()
879 if (surfactant().soluble())
895 zeroGradientFaPatchScalarField::typeName
899 bulkSurfactantConcentration().boundaryField()[fsPatchIndex()];
900 Cb.correctBoundaryConditions();
905 surfactant().adsorptionCoeff()*
Cb
906 + surfactant().adsorptionCoeff()
907 *surfactant().desorptionCoeff(),
908 surfactantConcentration()
910 - surfactant().adsorptionCoeff()
911 *
Cb*surfactant().saturatedConc();
919 sigma() + surfactant().dSigma(surfactantConcentration())/rho_;
921 if (
neg(
min(surfaceTension().internalField().
field())))
924 <<
"Surface tension is negative"
931 Foam::vector Foam::interfaceTrackingFvMesh::totalPressureForce()
const
941 return gSum(pressureForces);
945 Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce()
const
965 U().boundaryField()[fsPatchIndex()].
snGrad()
978 return gSum(viscousForces);
982 Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce()
const
992 if (pureFreeSurface())
1003 surfTensionForces = surfaceTension().internalField().field()*
K*S*
n;
1006 return gSum(surfTensionForces);
1010 Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber()
1014 if (pureFreeSurface())
1020 mesh().time().deltaT().value()/
1039 mesh().time().deltaT().value()/
1051 void Foam::interfaceTrackingFvMesh::updateProperties()
1056 "transportProperties"
1065 void Foam::interfaceTrackingFvMesh::correctPointDisplacement
1080 for (
const word& patchName : fixedFreeSurfacePatches_)
1094 label curFace = eFaces[edgeI];
1096 const labelList& curPoints = faces[curFace];
1098 forAll(curPoints, pointI)
1100 label curPoint = curPoints[pointI];
1101 label index = pLabels.find(curPoint);
1114 forAll(corrPoints, pointI)
1116 label curPoint = corrPoints[pointI];
1122 label curFace =
pFaces[curPoint][faceI];
1124 label index = eFaces.find(curFace);
1132 corrPointFaces[pointI] =
faceSet.toc();
1135 forAll(corrPoints, pointI)
1137 label curPoint = corrPoints[pointI];
1141 const labelList& curPointFaces = corrPointFaces[pointI];
1143 forAll(curPointFaces, faceI)
1145 const face& curFace = faces[curPointFaces[faceI]];
1147 label ptInFace = curFace.
which(curPoint);
1148 label next = curFace.
nextLabel(ptInFace);
1149 label prev = curFace.
prevLabel(ptInFace);
1153 const vector&
c = pointsDisplacementDir()[curPoint];
1155 curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^
b)&
c);
1158 curDisp /= curPointFaces.size();
1160 displacement[curPoint] =
1161 curDisp*pointsDisplacementDir()[curPoint];
1166 for (
const word& patchName : nonReflectingFreeSurfacePatches_)
1168 label nonReflectingPatchID =
1181 forAll(corrPoints, pointI)
1183 label curPoint = corrPoints[pointI];
1189 label curFace =
pFaces[curPoint][faceI];
1191 label index = eFaces.find(curFace);
1199 corrPointFaces[pointI] =
faceSet.toc();
1203 forAll(corrPoints, pointI)
1205 label curPoint = corrPoints[pointI];
1209 const labelList& curPointFaces = corrPointFaces[pointI];
1211 forAll(curPointFaces, faceI)
1213 const face& curFace = faces[curPointFaces[faceI]];
1215 label ptInFace = curFace.
which(curPoint);
1216 label next = curFace.
nextLabel(ptInFace);
1217 label prev = curFace.
prevLabel(ptInFace);
1223 if (corrPoints.find(next) == -1)
1238 vector c0 = displacement[p1];
1240 scalar V0 =
mag((
a0^b0)&c0)/2;
1242 scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0;
1244 if (corrPoints.find(prev) != -1)
1248 (
points[next] + displacement[next])
1250 const vector&
c = pointsDisplacementDir()[curPoint];
1252 curDisp += 2*DV/((a^
b)&
c);
1257 - (
points[prev] + displacement[prev]);
1259 const vector&
c = pointsDisplacementDir()[curPoint];
1261 curDisp += 2*DV/((a^
b)&
c);
1265 curDisp /= curPointFaces.size();
1267 displacement[curPoint] =
1268 curDisp*pointsDisplacementDir()[curPoint];
1274 void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals()
1285 if (contactAnglePtr_ && correctContactLineNormals())
1287 Info<<
"Correcting contact line normals" <<
endl;
1326 == processorFaPatch::typeName
1330 refCast<const processorFaPatch>
1338 forAll(patchPointLabels, pointI)
1340 label curPoint = patchPointLabels[pointI];
1347 const labelList& curPointEdges = pointEdges[curPoint];
1349 forAll(curPointEdges, edgeI)
1351 label curEdge = curPointEdges[edgeI];
1353 if (edgeFaces[curEdge].size() == 1)
1360 label index = curEdges.find(curEdge);
1367 != processorFaPatch::typeName
1384 vector t = edges[curEdge].vec(oldPoints);
1385 t /=
mag(t) + SMALL;
1388 pointFaces[curPoint];
1390 forAll(curPointFaces, fI)
1392 tangent.ref().field()[curPointFaces[fI]] = t;
1399 tangent.correctBoundaryConditions();
1404 label ngbPolyPatchID =
1407 if (ngbPolyPatchID != -1)
1412 == wallFvPatch::typeName
1420 - contactAnglePtr_->boundaryField()[patchI]
1435 rotationAxis /=
mag(rotationAxis) + SMALL;
1445 forAll(pointEdges, pointI)
1449 forAll(pointEdges[pointI], eI)
1453 + pointEdges[pointI][eI];
1455 vector e = edges[curEdge].vec(oldPoints);
1457 e *= (
e&rotationAxis[pointI])
1458 /
mag(
e&rotationAxis[pointI]);
1460 e /=
mag(
e) + SMALL;
1465 if (pointEdges[pointI].size() == 1)
1467 label curPoint = patchPoints[pointI];
1474 label procPatchID = -1;
1480 forAll(curPointEdges, edgeI)
1482 label curEdge = curPointEdges[edgeI];
1484 if (edgeFaces[curEdge].size() == 1)
1492 curEdges.find(curEdge);
1499 == processorFaPatch::typeName
1511 if (procPatchID != -1)
1514 tangent.boundaryField()[procPatchID]
1515 .patchNeighbourField()()[
edgeID];
1517 t *= (t&rotationAxis[pointI])
1518 /(
mag(t&rotationAxis[pointI]) + SMALL);
1520 t /=
mag(t) + SMALL;
1526 rotationAxis[pointI] = rotAx/(
mag(rotAx) + SMALL);
1530 ngbN = ngbN*
cos(rotAngle)
1531 + rotationAxis*(rotationAxis & ngbN)*(1 -
cos(rotAngle))
1532 + (rotationAxis^ngbN)*
sin(rotAngle);
1535 forAll(patchPoints, pointI)
1537 N[patchPoints[pointI]] -=
1538 ngbN[pointI]*(ngbN[pointI]&
N[patchPoints[pointI]]);
1540 N[patchPoints[pointI]] /=
1541 mag(
N[patchPoints[pointI]]) + SMALL;
1554 Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1563 fixedFreeSurfacePatches_(),
1564 nonReflectingFreeSurfacePatches_(),
1565 pointNormalsCorrectionPatches_(),
1566 normalMotionDir_(
false),
1569 pureFreeSurface_(
true),
1570 rigidFreeSurface_(
false),
1571 correctContactLineNormals_(
false),
1576 controlPointsPtr_(
nullptr),
1577 motionPointsMaskPtr_(
nullptr),
1578 pointsDisplacementDirPtr_(
nullptr),
1579 facesDisplacementDirPtr_(
nullptr),
1580 fsNetPhiPtr_(
nullptr),
1582 surfactConcPtr_(
nullptr),
1583 bulkSurfactConcPtr_(
nullptr),
1584 surfaceTensionPtr_(
nullptr),
1585 surfactantPtr_(
nullptr),
1586 contactAnglePtr_(
nullptr)
1669 aMeshPtr_.reset(
new faMesh(*
this));
1672 fixedFreeSurfacePatches_ =
1673 motion().get<
wordList>(
"fixedFreeSurfacePatches");
1675 pointNormalsCorrectionPatches_ =
1676 motion().get<
wordList>(
"pointNormalsCorrectionPatches");
1678 normalMotionDir_ = motion().get<
bool>(
"normalMotionDir");
1679 smoothing_ = motion().getOrDefault(
"smoothing",
false);
1680 pureFreeSurface_ = motion().getOrDefault(
"pureFreeSurface",
true);
1717 return *fsNetPhiPtr_;
1728 return *fsNetPhiPtr_;
1734 forAll(
Us().boundaryField(), patchI)
1738 Us().boundaryField()[patchI].
type()
1739 == calculatedFaPatchVectorField::typeName
1744 pUs =
Us().boundaryField()[patchI].patchInternalField();
1746 label ngbPolyPatchID =
1749 if (ngbPolyPatchID != -1)
1754 U().boundaryField()[ngbPolyPatchID].
type()
1755 == slipFvPatchVectorField::typeName
1759 U().boundaryField()[ngbPolyPatchID].
type()
1760 == symmetryFvPatchVectorField::typeName
1775 Us().correctBoundaryConditions();
1783 Us().ref().field() =
U().boundaryField()[fsPatchIndex()];
1791 correctUsBoundaryConditions();
1797 return *getObjectPtr<const volVectorField>(
"U");
1803 return *getObjectPtr<const volScalarField>(
"p");
1809 return *getObjectPtr<const surfaceScalarField>(
"phi");
1817 auto& SnGradU = tSnGradU.ref();
1824 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1831 gradUs -=
n*(
n & gradUs);
1841 if (!pureFreeSurface() &&
max(
nu) > SMALL)
1843 tangentialSurfaceTensionForce =
1844 surfaceTensionGrad()().internalField();
1848 tangentialSurfaceTensionForce/(
nu + SMALL)
1860 auto& SnGradUn = tSnGradUn.ref();
1865 -
aMesh().faceCurvatures()*(
aMesh().faceAreaNormals()&
Us())
1878 auto& pressureJump = tPressureJump.ref();
1892 + 2.0*
nu*freeSurfaceSnGradUn();
1894 if (pureFreeSurface())
1900 pressureJump -= surfaceTension().internalField()*
K;
1903 return tPressureJump;
1909 if (!controlPointsPtr_)
1911 makeControlPoints();
1914 return *controlPointsPtr_;
1920 if (!motionPointsMaskPtr_)
1922 makeMotionPointsMask();
1925 return *motionPointsMaskPtr_;
1931 if (!pointsDisplacementDirPtr_)
1936 return *pointsDisplacementDirPtr_;
1942 if (!facesDisplacementDirPtr_)
1947 return *facesDisplacementDirPtr_;
1964 if (!surfactConcPtr_)
1969 return *surfactConcPtr_;
1976 if (!surfactConcPtr_)
1981 return *surfactConcPtr_;
1988 if (!bulkSurfactConcPtr_)
1990 makeBulkSurfactConc();
1993 return *bulkSurfactConcPtr_;
2000 if (!bulkSurfactConcPtr_)
2002 makeBulkSurfactConc();
2005 return *bulkSurfactConcPtr_;
2012 if (!surfaceTensionPtr_)
2014 makeSurfaceTension();
2017 return *surfaceTensionPtr_;
2024 if (!surfaceTensionPtr_)
2026 makeSurfaceTension();
2029 return *surfaceTensionPtr_;
2036 if (!surfactantPtr_)
2041 return *surfactantPtr_;
2052 "surfaceTensionGrad",
2064 auto&
grad = tgrad.ref();
2069 grad.correctBoundaryConditions();
2079 if (smoothing_ && !rigidFreeSurface_)
2081 smoothFreeSurfaceMesh();
2082 clearControlPoints();
2085 updateDisplacementDirections();
2089 Info<<
"Maximal capillary Courant number: "
2090 << maxCourantNumber() <<
endl;
2094 Info<<
"Free surface curvature: min = " <<
gMin(
K)
2100 if (!rigidFreeSurface_)
2116 Info<<
"Free surface continuity error : sum local = "
2117 <<
gSum(
mag(sweptVolCorr)) <<
", global = " <<
gSum(sweptVolCorr)
2121 fsNetPhi().ref().field() = sweptVolCorr;
2127 "ddt(" +
U().
name() +
')'
2165 <<
"Unsupported temporal differencing scheme : "
2175 sweptVolCorr/(Sf*(Nf & facesDisplacementDir()))
2178 for (
const word& patchName : fixedFreeSurfacePatches_)
2180 label fixedPatchID =
2183 if (fixedPatchID == -1)
2186 <<
"Wrong faPatch name '" << patchName
2187 <<
"'in the fixedFreeSurfacePatches list"
2188 <<
" defined in dynamicMeshDict dictionary"
2197 deltaHf[eFaces[edgeI]] *= 2.0;
2201 controlPoints() += facesDisplacementDir()*deltaHf;
2203 pointField displacement(pointDisplacement());
2204 correctPointDisplacement(sweptVolCorr, displacement);
2207 refCast<velocityMotionSolver>
2215 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2216 refCast<fixedValuePointPatchVectorField>
2224 fsPatchPointMeshU ==
2229 correctContactLinePointNormals();
2240 refCast<velocityMotionSolver>
2248 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2249 refCast<fixedValuePointPatchVectorField>
2257 fsPatchPointMeshU ==
2265 updateSurfactantConcentration();
2277 mesh().time().timePath()/
"freeSurface",
2289 mesh().time().timePath()/
"freeSurfaceControlPoints.vtk"
2292 Info<<
"Writing free surface control points to " <<
os.
name() <<
nl;
2294 os <<
"# vtk DataFile Version 2.0" <<
nl
2295 <<
"freeSurfaceControlPoints" <<
nl
2297 <<
"DATASET POLYDATA" <<
nl;
2299 const label
nPoints = controlPoints().size();
2302 for (
const point&
p : controlPoints())
2304 os << float(
p.x()) <<
' '
2305 <<
float(
p.y()) <<
' '
2306 <<
float(
p.z()) <<
nl;
2310 for (label
id = 0;
id <
nPoints; ++id)
2312 os << 1 <<
' ' <<
id <<
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.
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 local patch faces.
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.
virtual const fileName & name() const
Read/write access to the name of the stream.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
areaScalarField & surfactantConcentration()
Return free-surface surfactant concentration field.
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVelocity
Write concrete PrimitivePatch faces/points (optionally with fields) as a vtp file or a legacy vtk fil...
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.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
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 dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0)
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
UPtrList< const labelUList > edgeFaces() const
Return a list of edgeFaces for each patch.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
const Type & value() const
Return const reference to value.
const vectorField & pointAreaNormals() const
Return point area normals.
Type gSum(const FieldField< Field, Type > &f)
Template specialisation for scalar faMatrix.
A simple single-phase transport model based on viscosityModel.
const Internal & internalField() const
Return a const-reference to the dimensioned internal field.
const dimensionSet dimForce
Virtual base class for velocity motion solver.
const motionSolver & motion() const
Return the motionSolver.
tmp< GeometricField< Type, faPatchField, areaMesh > > div(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Legacy ASCII, legacyAsciiFormatter.
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.
label timeIndex() const noexcept
Return current time index.
const edgeList & edges() const noexcept
Return local patch edges with reordered boundary.
tmp< areaVectorField > surfaceTensionGrad()
Return surface tension gradient.
const dimensionSet dimArea(sqr(dimLength))
label which(const label pointLabel) const
Find local index on face for the point label, same as find()
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 (stdout output on master, null elsewhere)
#define DebugInFunction
Report an information message using Foam::Info.
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.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const Type & lookupObject(const word &name, const bool recursive=false) const
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
virtual const pointField & oldPoints() const
Return old points (mesh motion)
tmp< GeometricField< Type, faPatchField, areaMesh > > edgeIntegrate(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Abstract base class for turbulence models (RAS, LES and laminar).
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
label findPatchID(const word &patchName) const
Find patch index given a name, return -1 if not found.
OBJstream os(runTime.globalPath()/outputName)
Macros for easy insertion into run-time selection tables.
const volVectorField & U() const
Return constant reference to velocity field.
Base class for graphics format writing. Entry points are.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
const faBoundaryMesh & boundary() const noexcept
Return constant reference to boundary mesh.
const areaScalarField & fsNetPhi() const
Return free-surface net flux.
tmp< vectorField > freeSurfaceSnGradU()
Return free surface normal derivative of velocity.
const pointField & points() const
Return local patch 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.
label start() const
The start label of the edges in the faMesh edges list.
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.
static bool & parRun() noexcept
Test if this a parallel run.
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.
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.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
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)
const surfactantProperties & surfactant() const
Return surfactant properties.
const dimensionSet dimVolume(pow3(dimLength))
const labelList & meshPoints() const
Return labelList of mesh points in patch.
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
const uindirectPrimitivePatch & patch() const
Return constant reference to primitive patch.
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.
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]
Type gMax(const FieldField< Field, Type > &f)
areaVectorField & Us()
Return free-surface velocity field.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
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.
const dimensionSet dimless
Dimensionless.
* ~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)