interfaceTrackingFvMesh.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb.
9 Copyright (C) 2020-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
31#include "motionSolver.H"
32#include "volFields.H"
33#include "processorFaPatch.H"
34#include "wedgeFaPatch.H"
35#include "wedgeFaPatchFields.H"
36#include "slipFaPatchFields.H"
38#include "slipFvPatchFields.H"
40#include "wallFvPatch.H"
41#include "polyPatchID.H"
42#include "fvcMeshPhi.H"
44#include "EulerDdtScheme.H"
46#include "backwardDdtScheme.H"
47#include "twoDPointCorrector.H"
48#include "gravityMeshObject.H"
50#include "demandDrivenData.H"
51#include "unitConversion.H"
53
54// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
55
56namespace Foam
57{
60 (
64 );
66 (
69 doInit
70 );
71}
72
73
74// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
75
76void Foam::interfaceTrackingFvMesh::initializeData()
77{
78 // Set free surface patch index
79 {
80 const word fsPatchName(motion().get<word>("fsPatchName"));
81
82 polyPatchID patch(fsPatchName, this->boundaryMesh());
83
84 if (!patch.active())
85 {
87 << "Patch name " << fsPatchName << " not found."
88 << abort(FatalError);
89 }
90
91 fsPatchIndex_ = patch.index();
92 }
93
94 // Set point normal correction for finite area mesh
95 {
97
98 for (const word& patchName : pointNormalsCorrectionPatches_)
99 {
100 label patchID = aMesh().boundary().findPatchID(patchName);
101
102 if (patchID == -1)
103 {
105 << "Patch name '" << patchName
106 << "' for point normals correction does not exist"
107 << abort(FatalError);
108 }
109
110 correction[patchID] = true;
111 }
112 }
113
114 // Read motion direction
115 if (!normalMotionDir_)
116 {
117 motionDir_ = normalised(motion().get<vector>("motionDir"));
118 }
119
120 // Check if contact angle is defined
121 makeContactAngle();
122
124 (
125 "nonReflectingFreeSurfacePatches",
126 nonReflectingFreeSurfacePatches_
127 );
128}
129
130
131void Foam::interfaceTrackingFvMesh::makeUs() const
132{
134 << "making free-surface velocity field" << nl;
135
136 if (UsPtr_)
137 {
139 << "free-surface velocity field already exists"
140 << abort(FatalError);
141 }
142
143 wordList patchFieldTypes
144 (
145 aMesh().boundary().size(),
146 zeroGradientFaPatchVectorField::typeName
147 );
148
149 forAll(aMesh().boundary(), patchI)
150 {
151 if
152 (
153 aMesh().boundary()[patchI].type()
155 )
156 {
157 patchFieldTypes[patchI] =
158 wedgeFaPatchVectorField::typeName;
159 }
160 else
161 {
162 label ngbPolyPatchID =
163 aMesh().boundary()[patchI].ngbPolyPatchIndex();
164
165 if (ngbPolyPatchID != -1)
166 {
167 if
168 (
169 mesh().boundary()[ngbPolyPatchID].type()
171 )
172 {
173 patchFieldTypes[patchI] =
174 slipFaPatchVectorField::typeName;
175 }
176 }
177 }
178 }
179
180 for (const word& patchName : fixedFreeSurfacePatches_)
181 {
182 const label fixedPatchID =
183 aMesh().boundary().findPatchID(patchName);
184
185 if (fixedPatchID == -1)
186 {
188 << "Wrong faPatch name '" << patchName
189 << "' in the fixedFreeSurfacePatches list"
190 << " defined in the dynamicMeshDict dictionary"
191 << abort(FatalError);
192 }
193
194 label ngbPolyPatchID =
195 aMesh().boundary()[fixedPatchID].ngbPolyPatchIndex();
196
197 if (ngbPolyPatchID != -1)
198 {
199 if
200 (
201 mesh().boundary()[ngbPolyPatchID].type()
203 )
204 {
205 patchFieldTypes[fixedPatchID] =
206 fixedValueFaPatchVectorField::typeName;
207 }
208 }
209 }
210
211 UsPtr_ = new areaVectorField
212 (
214 (
215 "Us",
216 mesh().time().timeName(),
217 mesh(),
220 ),
221 aMesh(),
223 patchFieldTypes
224 );
225
226 for (const word& patchName : fixedFreeSurfacePatches_)
227 {
228 const label fixedPatchID = aMesh().boundary().findPatchID(patchName);
229
230 if (fixedPatchID == -1)
231 {
233 << "Wrong faPatch name '" << patchName
234 << "' in the fixedFreeSurfacePatches list"
235 << " defined in the dynamicMeshDict dictionary"
236 << abort(FatalError);
237 }
238
239 label ngbPolyPatchID =
240 aMesh().boundary()[fixedPatchID].ngbPolyPatchIndex();
241
242 if (ngbPolyPatchID != -1)
243 {
244 if
245 (
246 mesh().boundary()[ngbPolyPatchID].type()
248 )
249 {
250 UsPtr_->boundaryFieldRef()[fixedPatchID] == Zero;
251 }
252 }
253 }
254}
255
256
257void Foam::interfaceTrackingFvMesh::makeFsNetPhi() const
258{
260 << "making free-surface net flux" << nl;
261
262 if (fsNetPhiPtr_)
263 {
265 << "free-surface net flux already exists"
266 << abort(FatalError);
267 }
268
269 fsNetPhiPtr_ = new areaScalarField
270 (
272 (
273 "fsNetPhi",
274 mesh().time().timeName(),
275 mesh(),
278 ),
279 aMesh(),
281 );
282}
283
284
285void Foam::interfaceTrackingFvMesh::makeControlPoints()
286{
288 << "making control points" << nl;
289
290 if (controlPointsPtr_)
291 {
293 << "control points already exists"
294 << abort(FatalError);
295 }
296
297 IOobject controlPointsHeader
298 (
299 "controlPoints",
300 mesh().time().timeName(),
301 mesh(),
303 );
304
305 if (controlPointsHeader.typeHeaderOk<vectorIOField>())
306 {
307 Info<< "Reading control points" << endl;
308 controlPointsPtr_ =
309 new vectorIOField
310 (
312 (
313 "controlPoints",
314 mesh().time().timeName(),
315 mesh(),
318 )
319 );
320 }
321 else
322 {
323 Info<< "Creating new control points" << endl;
324 controlPointsPtr_ =
325 new vectorIOField
326 (
328 (
329 "controlPoints",
330 mesh().time().timeName(),
331 mesh(),
334 ),
335 aMesh().areaCentres().internalField()
336 );
337
338 initializeControlPointsPosition();
339 }
340}
341
342
343void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
344{
346 << "making motion points mask" << nl;
347
348 if (motionPointsMaskPtr_)
349 {
351 << "motion points mask already exists"
352 << abort(FatalError);
353 }
354
355 motionPointsMaskPtr_ = new labelList
356 (
357 mesh().boundaryMesh()[fsPatchIndex()].nPoints(),
358 1
359 );
360
361 // Mark free surface boundary points
362 // that belong to processor patches
363 forAll(aMesh().boundary(), patchI)
364 {
365 if
366 (
367 aMesh().boundary()[patchI].type()
369 )
370 {
371 const labelList& patchPoints =
372 aMesh().boundary()[patchI].pointLabels();
373
374 forAll(patchPoints, pointI)
375 {
376 motionPointsMask()[patchPoints[pointI]] = -1;
377 }
378 }
379 }
380
381 // Mark fixed free surface boundary points
382 for (const word& patchName : fixedFreeSurfacePatches_)
383 {
384 const label fixedPatchID = aMesh().boundary().findPatchID(patchName);
385
386 if (fixedPatchID == -1)
387 {
389 << "Wrong faPatch name in the fixedFreeSurfacePatches list"
390 << " defined in the dynamicMeshDict dictionary"
391 << abort(FatalError);
392 }
393
394 const labelList& patchPoints =
395 aMesh().boundary()[fixedPatchID].pointLabels();
396
397 forAll(patchPoints, pointI)
398 {
399 motionPointsMask()[patchPoints[pointI]] = 0;
400 }
401 }
402}
403
404
405void Foam::interfaceTrackingFvMesh::makeDirections()
406{
408 << "make displacement directions for points and control points" << nl;
409
410 if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_)
411 {
413 << "points, control points displacement directions already exist"
414 << abort(FatalError);
415 }
416
417 pointsDisplacementDirPtr_ =
418 new vectorField
419 (
420 mesh().boundaryMesh()[fsPatchIndex()].nPoints(),
421 Zero
422 );
423
424 facesDisplacementDirPtr_ =
425 new vectorField
426 (
427 mesh().boundaryMesh()[fsPatchIndex()].size(),
428 Zero
429 );
430
431 if (!normalMotionDir())
432 {
433 if (mag(motionDir_) < SMALL)
434 {
436 << "Zero motion direction"
437 << abort(FatalError);
438 }
439
440 facesDisplacementDir() = motionDir_;
441 pointsDisplacementDir() = motionDir_;
442 }
443
444 updateDisplacementDirections();
445}
446
447
448void Foam::interfaceTrackingFvMesh::makePhis()
449{
451 << "making free-surface flux" << nl;
452
453 if (phisPtr_)
454 {
456 << "free-surface flux already exists"
457 << abort(FatalError);
458 }
459
460 phisPtr_ = new edgeScalarField
461 (
463 (
464 "phis",
465 mesh().time().timeName(),
466 mesh(),
469 ),
470 linearEdgeInterpolate(Us()) & aMesh().Le()
471 );
472}
473
474
475void Foam::interfaceTrackingFvMesh::makeSurfactConc() const
476{
478 << "making free-surface surfactant concentration field" << nl;
479
480 if (surfactConcPtr_)
481 {
483 << "free-surface surfactant concentration field already exists"
484 << abort(FatalError);
485 }
486
487 surfactConcPtr_ = new areaScalarField
488 (
490 (
491 "Cs",
492 mesh().time().timeName
493 (
494 mesh().time().startTime().value()
495 ),
496 // mesh().time().timeName(),
497 mesh(),
500 ),
501 aMesh()
502 );
503}
504
505
506void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc() const
507{
509 << "making volume surfactant concentration field" << nl;
510
511 if (bulkSurfactConcPtr_)
512 {
514 << "volume surfactant concentration field already exists"
515 << abort(FatalError);
516 }
517
518 bulkSurfactConcPtr_ = new volScalarField
519 (
521 (
522 "C",
523 mesh().time().timeName
524 (
525 mesh().time().startTime().value()
526 ),
527 // mesh().time().timeName(),
528 mesh(),
531 ),
532 mesh()
533 );
534 volScalarField& bulkSurfactConc = *bulkSurfactConcPtr_;
535
536 if (mesh().time().timeIndex()-1 == 0)
537 {
538 // Initialize uniform volume surfactant concentration
539 bulkSurfactConc = surfactant().bulkConc();
540 bulkSurfactConc.correctBoundaryConditions();
541 }
542}
543
544
545void Foam::interfaceTrackingFvMesh::makeSurfaceTension() const
546{
548 << "making surface tension field" << nl;
549
550 if (surfaceTensionPtr_)
551 {
553 << "surface tension field already exists"
554 << abort(FatalError);
555 }
556
557 surfaceTensionPtr_ = new areaScalarField
558 (
560 (
561 "surfaceTension",
562 mesh().time().timeName(),
563 mesh(),
566 ),
567 sigma() + surfactant().dSigma(surfactantConcentration())/rho_
568 );
569}
570
571
572void Foam::interfaceTrackingFvMesh::makeSurfactant() const
573{
575 << "making surfactant properties" << nl;
576
577 if (surfactantPtr_)
578 {
580 << "surfactant properties already exists"
581 << abort(FatalError);
582 }
583
584 const dictionary& surfactProp =
585 motion().subDict("surfactantProperties");
586
587 surfactantPtr_ = new surfactantProperties(surfactProp);
588}
589
590
591void Foam::interfaceTrackingFvMesh::makeContactAngle()
592{
594 << "making contact angle field" << nl;
595
596 if (contactAnglePtr_)
597 {
599 << "contact angle already exists"
600 << abort(FatalError);
601 }
602
603 // Check if contactAngle is defined
604 IOobject contactAngleHeader
605 (
606 "contactAngle",
607 mesh().time().timeName(),
608 mesh(),
610 );
611
612 if (contactAngleHeader.typeHeaderOk<areaScalarField>())
613 {
614 Info<< "Reading contact angle field" << endl;
615
616 contactAnglePtr_ =
618 (
620 (
621 "contactAngle",
622 mesh().time().timeName(),
623 mesh(),
626 ),
627 aMesh()
628 );
629 }
630}
631
632
633void Foam::interfaceTrackingFvMesh::updateDisplacementDirections()
634{
635 if (normalMotionDir())
636 {
637 // Update point displacement direction
638 pointsDisplacementDir() = aMesh().pointAreaNormals();
639
640 // Correct point displacement direction at contact line
641 forAll(aMesh().boundary(), patchI)
642 {
643 if (contactAnglePtr_)
644 {
645 label ngbPolyPatchID =
646 aMesh().boundary()[patchI].ngbPolyPatchIndex();
647
648 if (ngbPolyPatchID != -1)
649 {
650 if
651 (
652 mesh().boundary()[ngbPolyPatchID].type()
654 )
655 {
656 labelList patchPoints =
657 aMesh().boundary()[patchI].pointLabels();
658
660 (
661 aMesh().boundary()[patchI]
662 .ngbPolyPatchPointNormals()
663 );
664
665 forAll(patchPoints, pointI)
666 {
667 pointsDisplacementDir()[patchPoints[pointI]] -=
668 N[pointI]
669 *(
670 N[pointI]
671 & pointsDisplacementDir()[patchPoints[pointI]]
672 );
673
674 pointsDisplacementDir()[patchPoints[pointI]] /=
675 mag
676 (
677 pointsDisplacementDir()
678 [
679 patchPoints[pointI]
680 ]
681 ) + SMALL;
682 }
683 }
684 }
685 }
686 }
687
688 // Update face displacement direction
689 facesDisplacementDir() =
691
692 // Correction of control points position
693 const vectorField& Cf = aMesh().areaCentres().internalField();
694
695 controlPoints() =
696 facesDisplacementDir()
697 *(facesDisplacementDir()&(controlPoints() - Cf))
698 + Cf;
699 }
700}
701
702
703void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition()
704{
705 {
706 const faceList& faces = aMesh().faces();
707 const pointField& points = aMesh().points();
708
709 pointField displacement(pointDisplacement());
710 scalarField sweptVolCorr(faces.size(), Zero);
711 correctPointDisplacement(sweptVolCorr, displacement);
712
713 pointField newPoints(points + displacement);
714
715 scalarField sweptVol(faces.size(), Zero);
716
717 forAll(faces, faceI)
718 {
719 sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints);
720 }
721
722 vectorField faceArea(faces.size(), Zero);
723
724 forAll(faceArea, faceI)
725 {
726 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
727 }
728
729 scalarField deltaH = scalarField(aMesh().nFaces(), Zero);
730
731 forAll(deltaH, faceI)
732 {
733 deltaH[faceI] = sweptVol[faceI]/
734 ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL);
735
736 if (mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL)
737 {
738 // Info<< (faceArea[faceI] & facesDisplacementDir()[faceI])
739 // << ", " << faceArea[faceI]
740 // << ", " << facesDisplacementDir()[faceI] << endl;
741
743 << "Something is wrong with specified motion direction"
744 << abort(FatalError);
745 }
746 }
747
748 for (const word& patchName : fixedFreeSurfacePatches_)
749 {
750 label fixedPatchID = aMesh().boundary().findPatchID(patchName);
751
752 if (fixedPatchID == -1)
753 {
755 << "Wrong faPatch name in the fixedFreeSurfacePatches list"
756 << " defined in the freeSurfaceProperties dictionary"
757 << abort(FatalError);
758 }
759
760 const labelList& eFaces =
761 aMesh().boundary()[fixedPatchID].edgeFaces();
762
763 forAll(eFaces, edgeI)
764 {
765 deltaH[eFaces[edgeI]] *= 2.0;
766 }
767 }
768
769 controlPoints() += facesDisplacementDir()*deltaH;
770 }
771}
772
773
774void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh()
775{
776 Info<< "Smoothing free surface mesh" << endl;
777
778 controlPoints() = aMesh().areaCentres().internalField();
779
780 pointField displacement(pointDisplacement());
781
782 const faceList& faces = aMesh().faces();
783 const pointField& points = aMesh().points();
784
785 pointField newPoints(points + displacement);
786
787 scalarField sweptVol(faces.size(), Zero);
788 forAll(faces, faceI)
789 {
790 sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints);
791 }
792
793 vectorField faceArea(faces.size(), Zero);
794 forAll(faceArea, faceI)
795 {
796 faceArea[faceI] = faces[faceI].unitNormal(newPoints);
797 }
798
799 scalarField deltaHf
800 (
801 sweptVol/(faceArea & facesDisplacementDir())
802 );
803
804 for (const word& patchName : fixedFreeSurfacePatches_)
805 {
806 label fixedPatchID = aMesh().boundary().findPatchID(patchName);
807
808 if (fixedPatchID == -1)
809 {
811 << "Wrong faPatch name fixedFreeSurfacePatches list"
812 << " defined in the dynamicMeshDict dictionary"
813 << abort(FatalError);
814 }
815
816 const labelList& eFaces =
817 aMesh().boundary()[fixedPatchID].edgeFaces();
818
819 forAll(eFaces, edgeI)
820 {
821 deltaHf[eFaces[edgeI]] *= 2.0;
822 }
823 }
824
825 controlPoints() += facesDisplacementDir()*deltaHf;
826
827 displacement = pointDisplacement();
828
829 velocityMotionSolver& vMotion =
830 refCast<velocityMotionSolver>
831 (
832 const_cast<motionSolver&>(motion())
833 );
834
835 pointVectorField& pointMotionU = vMotion.pointMotionU();
836 pointMotionU.primitiveFieldRef() = Zero;
837
838 fixedValuePointPatchVectorField& fsPatchPointMeshU =
839 refCast<fixedValuePointPatchVectorField>
840 (
841 const_cast<pointPatchVectorField&>
842 (
843 pointMotionU.boundaryField()[fsPatchIndex()]
844 )
845 );
846
847 fsPatchPointMeshU ==
848 displacement/mesh().time().deltaT().value();
849
851}
852
853
854void Foam::interfaceTrackingFvMesh::updateSurfaceFlux()
855{
856 Phis() = fac::interpolate(Us()) & aMesh().Le();
857}
858
859
860void Foam::interfaceTrackingFvMesh::updateSurfactantConcentration()
861{
862 if (!pureFreeSurface())
863 {
864 Info<< "Correct surfactant concentration" << endl << flush;
865
866 updateSurfaceFlux();
867
868 // Crate and solve the surfactanta transport equation
869 faScalarMatrix CsEqn
870 (
871 fam::ddt(surfactantConcentration())
872 + fam::div(Phis(), surfactantConcentration())
874 (
875 surfactant().diffusion(),
876 surfactantConcentration()
877 )
878 );
879
880 if (surfactant().soluble())
881 {
882 #include "solveBulkSurfactant.H"
883
885 (
887 (
888 "Cb",
889 mesh().time().timeName(),
890 mesh(),
893 ),
894 aMesh(),
896 zeroGradientFaPatchScalarField::typeName
897 );
898
899 Cb.ref().field() =
900 bulkSurfactantConcentration().boundaryField()[fsPatchIndex()];
901 Cb.correctBoundaryConditions();
902
903 CsEqn +=
904 fam::Sp
905 (
906 surfactant().adsorptionCoeff()*Cb
907 + surfactant().adsorptionCoeff()
908 *surfactant().desorptionCoeff(),
909 surfactantConcentration()
910 )
911 - surfactant().adsorptionCoeff()
912 *Cb*surfactant().saturatedConc();
913 }
914
915 CsEqn.solve();
916
917 // Info<< "Correct surface tension" << endl;
918
919 surfaceTension() =
920 sigma() + surfactant().dSigma(surfactantConcentration())/rho_;
921
922 if (neg(min(surfaceTension().internalField().field())))
923 {
925 << "Surface tension is negative"
926 << abort(FatalError);
927 }
928 }
929}
930
931
932Foam::vector Foam::interfaceTrackingFvMesh::totalPressureForce() const
933{
934 const scalarField& S = aMesh().S();
935
937
938 const scalarField& P = p().boundaryField()[fsPatchIndex()];
939
940 vectorField pressureForces(S*P*n);
941
942 return gSum(pressureForces);
943}
944
945
946Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce() const
947{
948 const auto& turbulence =
949 mesh().lookupObject<turbulenceModel>("turbulenceProperties");
950
951 scalarField nu(turbulence.nuEff(fsPatchIndex()));
952
953 // const singlePhaseTransportModel& properties =
954 // mesh().lookupObject<singlePhaseTransportModel>
955 // (
956 // "transportProperties"
957 // );
958
959 // dimensionedScalar nu("nu", properties);
960
961 const scalarField& S = aMesh().S();
963
964 vectorField nGradU
965 (
966 U().boundaryField()[fsPatchIndex()].snGrad()
967 );
968
969 vectorField viscousForces
970 (
971 - nu*S
972 *(
973 nGradU
974 + (fac::grad(Us())().internalField()&n)
975 - (n*fac::div(Us())().internalField())
976 )
977 );
978
979 return gSum(viscousForces);
980}
981
982
983Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce() const
984{
985 const scalarField& S = aMesh().S();
986
988
990
991 vectorField surfTensionForces(n.size(), Zero);
992
993 if (pureFreeSurface())
994 {
995 surfTensionForces =
996 S*sigma().value()
998 (
999 aMesh().Le()*aMesh().edgeLengthCorrection()
1000 )().internalField();
1001 }
1002 else
1003 {
1004 surfTensionForces = surfaceTension().internalField().field()*K*S*n;
1005 }
1006
1007 return gSum(surfTensionForces);
1008}
1009
1010
1011Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber()
1012{
1013 scalar CoNum = 0;
1014
1015 if (pureFreeSurface())
1016 {
1017 const scalarField& dE = aMesh().lPN();
1018
1019 CoNum = gMax
1020 (
1021 mesh().time().deltaT().value()/
1022 sqrt
1023 (
1024 Foam::pow(dE, 3.0)/2.0/M_PI/(sigma().value() + SMALL)
1025 )
1026 );
1027 }
1028 else
1029 {
1030 scalarField sigmaE
1031 (
1032 linearEdgeInterpolate(surfaceTension())().internalField().field()
1033 + SMALL
1034 );
1035
1036 const scalarField& dE = aMesh().lPN();
1037
1038 CoNum = gMax
1039 (
1040 mesh().time().deltaT().value()/
1041 sqrt
1042 (
1043 Foam::pow(dE, 3.0)/2.0/M_PI/sigmaE
1044 )
1045 );
1046 }
1047
1048 return CoNum;
1049}
1050
1051
1052void Foam::interfaceTrackingFvMesh::updateProperties()
1053{
1054 const singlePhaseTransportModel& properties =
1056 (
1057 "transportProperties"
1058 );
1059
1060 rho_ = dimensionedScalar("rho", properties);
1061
1062 sigma0_ = dimensionedScalar("sigma", properties)/rho_;
1063}
1064
1065
1066void Foam::interfaceTrackingFvMesh::correctPointDisplacement
1067(
1068 const scalarField& sweptVolCorr,
1069 vectorField& displacement
1070)
1071{
1072 const labelListList& pFaces =
1073 aMesh().patch().pointFaces();
1074
1075 const faceList& faces =
1076 aMesh().patch().localFaces();
1077
1078 const pointField& points =
1079 aMesh().patch().localPoints();
1080
1081 for (const word& patchName : fixedFreeSurfacePatches_)
1082 {
1083 label fixedPatchID = aMesh().boundary().findPatchID(patchName);
1084
1085 const labelList& pLabels =
1086 aMesh().boundary()[fixedPatchID].pointLabels();
1087
1088 const labelList& eFaces =
1089 aMesh().boundary()[fixedPatchID].edgeFaces();
1090
1092
1093 forAll(eFaces, edgeI)
1094 {
1095 label curFace = eFaces[edgeI];
1096
1097 const labelList& curPoints = faces[curFace];
1098
1099 forAll(curPoints, pointI)
1100 {
1101 label curPoint = curPoints[pointI];
1102 label index = pLabels.find(curPoint);
1103
1104 if (index == -1)
1105 {
1106 pointSet.insert(curPoint);
1107 }
1108 }
1109 }
1110
1111 labelList corrPoints = pointSet.toc();
1112
1113 labelListList corrPointFaces(corrPoints.size());
1114
1115 forAll(corrPoints, pointI)
1116 {
1117 label curPoint = corrPoints[pointI];
1118
1120
1121 forAll(pFaces[curPoint], faceI)
1122 {
1123 label curFace = pFaces[curPoint][faceI];
1124
1125 label index = eFaces.find(curFace);
1126
1127 if (index != -1)
1128 {
1129 faceSet.insert(curFace);
1130 }
1131 }
1132
1133 corrPointFaces[pointI] = faceSet.toc();
1134 }
1135
1136 forAll(corrPoints, pointI)
1137 {
1138 label curPoint = corrPoints[pointI];
1139
1140 scalar curDisp = 0;
1141
1142 const labelList& curPointFaces = corrPointFaces[pointI];
1143
1144 forAll(curPointFaces, faceI)
1145 {
1146 const face& curFace = faces[curPointFaces[faceI]];
1147
1148 label ptInFace = curFace.which(curPoint);
1149 label next = curFace.nextLabel(ptInFace);
1150 label prev = curFace.prevLabel(ptInFace);
1151
1152 vector a = points[next] - points[curPoint];
1153 vector b = points[prev] - points[curPoint];
1154 const vector& c = pointsDisplacementDir()[curPoint];
1155
1156 curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^b)&c);
1157 }
1158
1159 curDisp /= curPointFaces.size();
1160
1161 displacement[curPoint] =
1162 curDisp*pointsDisplacementDir()[curPoint];
1163 }
1164 }
1165
1166
1167 for (const word& patchName : nonReflectingFreeSurfacePatches_)
1168 {
1169 label nonReflectingPatchID =
1170 aMesh().boundary().findPatchID(patchName);
1171
1172 const labelList& pLabels =
1173 aMesh().boundary()[nonReflectingPatchID].pointLabels();
1174
1175 const labelList& eFaces =
1176 aMesh().boundary()[nonReflectingPatchID].edgeFaces();
1177
1178 labelList corrPoints = pLabels;
1179
1180 labelListList corrPointFaces(corrPoints.size());
1181
1182 forAll(corrPoints, pointI)
1183 {
1184 label curPoint = corrPoints[pointI];
1185
1187
1188 forAll(pFaces[curPoint], faceI)
1189 {
1190 label curFace = pFaces[curPoint][faceI];
1191
1192 label index = eFaces.find(curFace);
1193
1194 if (index != -1)
1195 {
1196 faceSet.insert(curFace);
1197 }
1198 }
1199
1200 corrPointFaces[pointI] = faceSet.toc();
1201 }
1202
1203
1204 forAll(corrPoints, pointI)
1205 {
1206 label curPoint = corrPoints[pointI];
1207
1208 scalar curDisp = 0;
1209
1210 const labelList& curPointFaces = corrPointFaces[pointI];
1211
1212 forAll(curPointFaces, faceI)
1213 {
1214 const face& curFace = faces[curPointFaces[faceI]];
1215
1216 label ptInFace = curFace.which(curPoint);
1217 label next = curFace.nextLabel(ptInFace);
1218 label prev = curFace.prevLabel(ptInFace);
1219
1220 label p0 = -1;
1221 label p1 = -1;
1222 label p2 = -1;
1223
1224 if (corrPoints.find(next) == -1)
1225 {
1226 p0 = curPoint;
1227 p1 = next;
1228 p2 = curFace.nextLabel(curFace.which(next));
1229 }
1230 else
1231 {
1232 p0 = curFace.prevLabel(curFace.which(prev));
1233 p1 = prev;
1234 p2 = curPoint;
1235 }
1236
1237 vector a0 = points[p1] - points[p0];
1238 vector b0 = points[p2] - points[p1];
1239 vector c0 = displacement[p1];
1240
1241 scalar V0 = mag((a0^b0)&c0)/2;
1242
1243 scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0;
1244
1245 if (corrPoints.find(prev) != -1)
1246 {
1247 vector a = points[curPoint] - points[prev];
1248 vector b =
1249 (points[next] + displacement[next])
1250 - points[curPoint];
1251 const vector& c = pointsDisplacementDir()[curPoint];
1252
1253 curDisp += 2*DV/((a^b)&c);
1254 }
1255 else
1256 {
1257 vector a = points[curPoint]
1258 - (points[prev] + displacement[prev]);
1259 vector b = points[next] - points[curPoint];
1260 const vector& c = pointsDisplacementDir()[curPoint];
1261
1262 curDisp += 2*DV/((a^b)&c);
1263 }
1264 }
1265
1266 curDisp /= curPointFaces.size();
1267
1268 displacement[curPoint] =
1269 curDisp*pointsDisplacementDir()[curPoint];
1270 }
1271 }
1272}
1273
1274
1275void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals()
1276{
1277 // Correct normals for contact line points
1278 // according to specified contact angle
1279
1280 vectorField& N =
1281 const_cast<vectorField&>
1282 (
1284 );
1285
1286 if (contactAnglePtr_ && correctContactLineNormals())
1287 {
1288 Info<< "Correcting contact line normals" << endl;
1289
1290 vectorField oldPoints(aMesh().nPoints(), Zero);
1291
1292 const labelList& meshPoints = aMesh().patch().meshPoints();
1293
1294 forAll(oldPoints, ptI)
1295 {
1296 oldPoints[ptI] =
1297 mesh().oldPoints()[meshPoints[ptI]];
1298 }
1299
1300// # include "createTangentField.H"
1301 areaVectorField tangent
1302 (
1303 IOobject
1304 (
1305 "tangent",
1306 mesh().time().timeName(),
1307 mesh(),
1310 ),
1311 aMesh(),
1313 );
1314
1315 if (Pstream::parRun())
1316 {
1317 const labelListList& edgeFaces = aMesh().patch().edgeFaces();
1318 const labelListList& pointEdges = aMesh().patch().pointEdges();
1319 const labelListList& pointFaces = aMesh().patch().pointFaces();
1320 const edgeList& edges = aMesh().edges();
1321
1322 forAll(aMesh().boundary(), patchI)
1323 {
1324 if
1325 (
1326 aMesh().boundary()[patchI].type()
1328 )
1329 {
1330 const processorFaPatch& procPatch =
1331 refCast<const processorFaPatch>
1332 (
1333 aMesh().boundary()[patchI]
1334 );
1335
1336 const labelList& patchPointLabels =
1337 procPatch.pointLabels();
1338
1339 forAll(patchPointLabels, pointI)
1340 {
1341 label curPoint = patchPointLabels[pointI];
1342
1343 // Check if processor point is boundary point
1344
1345 label patchID = -1;
1346 label edgeID = -1;
1347
1348 const labelList& curPointEdges = pointEdges[curPoint];
1349
1350 forAll(curPointEdges, edgeI)
1351 {
1352 label curEdge = curPointEdges[edgeI];
1353
1354 if (edgeFaces[curEdge].size() == 1)
1355 {
1356 forAll(aMesh().boundary(), pI)
1357 {
1358 const labelList& curEdges =
1359 aMesh().boundary()[pI];
1360
1361 label index = curEdges.find(curEdge);
1362
1363 if (index != -1)
1364 {
1365 if
1366 (
1367 aMesh().boundary()[pI].type()
1369 )
1370 {
1371 patchID = pI;
1372 edgeID = index;
1373 break;
1374 }
1375 }
1376 }
1377 }
1378 }
1379
1380 if (patchID != -1)
1381 {
1382 label curEdge =
1384
1385 vector t = edges[curEdge].vec(oldPoints);
1386 t /= mag(t) + SMALL;
1387
1388 const labelList& curPointFaces =
1389 pointFaces[curPoint];
1390
1391 forAll(curPointFaces, fI)
1392 {
1393 tangent.ref().field()[curPointFaces[fI]] = t;
1394 }
1395 }
1396 }
1397 }
1398 }
1399
1400 tangent.correctBoundaryConditions();
1401 }
1402
1403 forAll(aMesh().boundary(), patchI)
1404 {
1405 label ngbPolyPatchID =
1406 aMesh().boundary()[patchI].ngbPolyPatchIndex();
1407
1408 if (ngbPolyPatchID != -1)
1409 {
1410 if
1411 (
1412 mesh().boundary()[ngbPolyPatchID].type()
1414 )
1415 {
1416 const scalar rotAngle = degToRad
1417 (
1418 gAverage
1419 (
1420 90
1421 - contactAnglePtr_->boundaryField()[patchI]
1422 )
1423 );
1424
1425 vectorField ngbN
1426 (
1427 aMesh().boundary()[patchI].ngbPolyPatchPointNormals()
1428 );
1429
1430 const labelList& patchPoints =
1431 aMesh().boundary()[patchI].pointLabels();
1432
1433 vectorField pN(N, patchPoints);
1434
1435 vectorField rotationAxis(ngbN^pN);
1436 rotationAxis /= mag(rotationAxis) + SMALL;
1437
1438
1439 // Calc rotation axis using edge vectors
1440
1441 const edgeList& edges = aMesh().edges();
1442
1443 const labelListList& pointEdges =
1444 aMesh().boundary()[patchI].pointEdges();
1445
1446 forAll(pointEdges, pointI)
1447 {
1448 vector rotAx = Zero;
1449
1450 forAll(pointEdges[pointI], eI)
1451 {
1452 label curEdge =
1453 aMesh().boundary()[patchI].start()
1454 + pointEdges[pointI][eI];
1455
1456 vector e = edges[curEdge].vec(oldPoints);
1457
1458 e *= (e&rotationAxis[pointI])
1459 /mag(e&rotationAxis[pointI]);
1460
1461 e /= mag(e) + SMALL;
1462
1463 rotAx += e;
1464 }
1465
1466 if (pointEdges[pointI].size() == 1)
1467 {
1468 label curPoint = patchPoints[pointI];
1469
1470 const labelListList& ptEdges =
1471 aMesh().patch().pointEdges();
1472 const labelList& curPointEdges =
1473 ptEdges[curPoint];
1474
1475 label procPatchID = -1;
1476 label edgeID = -1;
1477
1478 const labelListList& edgeFaces =
1479 aMesh().patch().edgeFaces();
1480
1481 forAll(curPointEdges, edgeI)
1482 {
1483 label curEdge = curPointEdges[edgeI];
1484
1485 if (edgeFaces[curEdge].size() == 1)
1486 {
1487 forAll(aMesh().boundary(), pI)
1488 {
1489 const labelList& curEdges =
1490 aMesh().boundary()[pI];
1491
1492 label index =
1493 curEdges.find(curEdge);
1494
1495 if (index != -1)
1496 {
1497 if
1498 (
1499 aMesh().boundary()[pI].type()
1501 )
1502 {
1503 procPatchID = pI;
1504 edgeID = index;
1505 break;
1506 }
1507 }
1508 }
1509 }
1510 }
1511
1512 if (procPatchID != -1)
1513 {
1514 vector t =
1515 tangent.boundaryField()[procPatchID]
1516 .patchNeighbourField()()[edgeID];
1517
1518 t *= (t&rotationAxis[pointI])
1519 /(mag(t&rotationAxis[pointI]) + SMALL);
1520
1521 t /= mag(t) + SMALL;
1522
1523 rotAx += t;
1524 }
1525 }
1526
1527 rotationAxis[pointI] = rotAx/(mag(rotAx) + SMALL);
1528 }
1529
1530 // Rodrigues' rotation formula
1531 ngbN = ngbN*cos(rotAngle)
1532 + rotationAxis*(rotationAxis & ngbN)*(1 - cos(rotAngle))
1533 + (rotationAxis^ngbN)*sin(rotAngle);
1534
1535 // Info<< aMesh().boundary()[patchI].name() << endl;
1536 forAll(patchPoints, pointI)
1537 {
1538 N[patchPoints[pointI]] -=
1539 ngbN[pointI]*(ngbN[pointI]&N[patchPoints[pointI]]);
1540
1541 N[patchPoints[pointI]] /=
1542 mag(N[patchPoints[pointI]]) + SMALL;
1543
1544 // Info<< N[patchPoints[pointI]] << endl;
1545 }
1546 }
1547 }
1548 }
1549 }
1550}
1551
1552
1553// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1554
1556(
1557 const IOobject& io,
1558 const bool doInit
1559)
1560:
1562 aMeshPtr_(nullptr),
1563 fsPatchIndex_(-1),
1564 fixedFreeSurfacePatches_(),
1565 nonReflectingFreeSurfacePatches_(),
1566 pointNormalsCorrectionPatches_(),
1567 normalMotionDir_(false),
1568 motionDir_(Zero),
1569 smoothing_(false),
1570 pureFreeSurface_(true),
1571 rigidFreeSurface_(false),
1572 correctContactLineNormals_(false),
1573 sigma0_("zero", dimForce/dimLength/dimDensity, Zero),
1574 rho_("one", dimDensity, 1.0),
1575 timeIndex_(-1),
1576 UsPtr_(nullptr),
1577 controlPointsPtr_(nullptr),
1578 motionPointsMaskPtr_(nullptr),
1579 pointsDisplacementDirPtr_(nullptr),
1580 facesDisplacementDirPtr_(nullptr),
1581 fsNetPhiPtr_(nullptr),
1582 phisPtr_(nullptr),
1583 surfactConcPtr_(nullptr),
1584 bulkSurfactConcPtr_(nullptr),
1585 surfaceTensionPtr_(nullptr),
1586 surfactantPtr_(nullptr),
1587 contactAnglePtr_(nullptr)
1588{
1589 if (doInit)
1590 {
1591 init(false); // do not initialise lower levels
1592 }
1593}
1594
1595/*
1596Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
1597(
1598 const IOobject& io,
1599 pointField&& points,
1600 faceList&& faces,
1601 labelList&& allOwner,
1602 labelList&& allNeighbour,
1603 const bool syncPar
1604)
1605:
1606 dynamicMotionSolverFvMesh
1607 (
1608 io,
1609 std::move(points),
1610 std::move(faces),
1611 std::move(allOwner),
1612 std::move(allNeighbour),
1613 syncPar
1614 ),
1615 aMeshPtr_(new faMesh(*this)),
1616 fsPatchIndex_(-1),
1617 fixedFreeSurfacePatches_(),
1618 nonReflectingFreeSurfacePatches_(),
1619 pointNormalsCorrectionPatches_(),
1620 normalMotionDir_(false),
1621 motionDir_(Zero),
1622 smoothing_(false),
1623 pureFreeSurface_(true),
1624 sigma0_("zero", dimForce/dimLength/dimDensity, Zero),
1625 rho_("one", dimDensity, 1.0),
1626 timeIndex_(-1),
1627 UsPtr_(nullptr),
1628 controlPointsPtr_(nullptr),
1629 motionPointsMaskPtr_(nullptr),
1630 pointsDisplacementDirPtr_(nullptr),
1631 facesDisplacementDirPtr_(nullptr),
1632 fsNetPhiPtr_(nullptr),
1633 phisPtr_(nullptr),
1634 surfactConcPtr_(nullptr),
1635 bulkSurfactConcPtr_(nullptr),
1636 surfaceTensionPtr_(nullptr),
1637 surfactantPtr_(nullptr),
1638 contactAnglePtr_(nullptr)
1639{}
1640*/
1641
1642// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1643
1645{
1646 deleteDemandDrivenData(UsPtr_);
1647 deleteDemandDrivenData(controlPointsPtr_);
1648 deleteDemandDrivenData(motionPointsMaskPtr_);
1649 deleteDemandDrivenData(pointsDisplacementDirPtr_);
1650 deleteDemandDrivenData(facesDisplacementDirPtr_);
1651 deleteDemandDrivenData(fsNetPhiPtr_);
1652 deleteDemandDrivenData(phisPtr_);
1653 deleteDemandDrivenData(surfactConcPtr_);
1654 deleteDemandDrivenData(bulkSurfactConcPtr_);
1655 deleteDemandDrivenData(surfaceTensionPtr_);
1656 deleteDemandDrivenData(surfactantPtr_);
1657 deleteDemandDrivenData(contactAnglePtr_);
1658}
1659
1660
1661// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1662
1664{
1665 if (doInit)
1666 {
1668 }
1669
1670 aMeshPtr_.reset(new faMesh(*this));
1671
1672 // Set motion-based data
1673 fixedFreeSurfacePatches_ =
1674 motion().get<wordList>("fixedFreeSurfacePatches");
1675
1676 pointNormalsCorrectionPatches_ =
1677 motion().get<wordList>("pointNormalsCorrectionPatches");
1678
1679 normalMotionDir_ = motion().get<bool>("normalMotionDir");
1680 smoothing_ = motion().getOrDefault("smoothing", false);
1681 pureFreeSurface_ = motion().getOrDefault("pureFreeSurface", true);
1682
1683 initializeData();
1684
1685 return true;
1686}
1687
1688
1690{
1691 if (!UsPtr_)
1692 {
1693 makeUs();
1694 }
1695
1696 return *UsPtr_;
1697}
1698
1699
1701{
1702 if (!UsPtr_)
1703 {
1704 makeUs();
1705 }
1706
1707 return *UsPtr_;
1708}
1709
1710
1712{
1713 if (!fsNetPhiPtr_)
1714 {
1715 makeFsNetPhi();
1716 }
1717
1718 return *fsNetPhiPtr_;
1719}
1720
1721
1723{
1724 if (!fsNetPhiPtr_)
1725 {
1726 makeFsNetPhi();
1727 }
1728
1729 return *fsNetPhiPtr_;
1730}
1731
1732
1734{
1735 forAll(Us().boundaryField(), patchI)
1736 {
1737 if
1738 (
1739 Us().boundaryField()[patchI].type()
1740 == calculatedFaPatchVectorField::typeName
1741 )
1742 {
1743 vectorField& pUs = Us().boundaryFieldRef()[patchI];
1744
1745 pUs = Us().boundaryField()[patchI].patchInternalField();
1746
1747 label ngbPolyPatchID =
1748 aMesh().boundary()[patchI].ngbPolyPatchIndex();
1749
1750 if (ngbPolyPatchID != -1)
1751 {
1752 if
1753 (
1754 (
1755 U().boundaryField()[ngbPolyPatchID].type()
1756 == slipFvPatchVectorField::typeName
1757 )
1758 ||
1759 (
1760 U().boundaryField()[ngbPolyPatchID].type()
1761 == symmetryFvPatchVectorField::typeName
1762 )
1763 )
1764 {
1766 (
1767 aMesh().boundary()[patchI].ngbPolyPatchFaceNormals()
1768 );
1769
1770 pUs -= N*(N&pUs);
1771 }
1772 }
1773 }
1774 }
1775
1776 Us().correctBoundaryConditions();
1777}
1778
1779
1781{
1782 // Info<< "Update Us" << endl;
1783
1784 Us().ref().field() = U().boundaryField()[fsPatchIndex()];
1785
1786 // // Correct normal component of free-surface velocity
1787 // const vectorField& nA = aMesh().faceAreaNormals().internalField();
1788 // vectorField UnFs = nA*phi().boundaryField()[fsPatchIndex()]
1789 // /mesh().boundary()[fsPatchIndex()].magSf();
1790 // Us().ref().field() += UnFs - nA*(nA&Us().internalField());
1791
1792 correctUsBoundaryConditions();
1793}
1794
1795
1797{
1798 return *getObjectPtr<const volVectorField>("U");
1799}
1800
1801
1803{
1804 return *getObjectPtr<const volScalarField>("p");
1805}
1806
1807
1809{
1810 return *getObjectPtr<const surfaceScalarField>("phi");
1811}
1812
1813
1816{
1817 auto tSnGradU = tmp<vectorField>::New(aMesh().nFaces(), Zero);
1818 auto& SnGradU = tSnGradU.ref();
1819
1821
1822 areaScalarField divUs
1823 (
1824 fac::div(Us())
1825 - aMesh().faceCurvatures()*(aMesh().faceAreaNormals()&Us())
1826 );
1827
1828 areaTensorField gradUs(fac::grad(Us()));
1829
1830 // Remove component of gradient normal to surface (area)
1832 gradUs -= n*(n & gradUs);
1834
1836 mesh().lookupObject<turbulenceModel>("turbulenceProperties");
1837
1838 scalarField nu(turbulence.nuEff(fsPatchIndex()));
1839
1840 vectorField tangentialSurfaceTensionForce(nA.size(), Zero);
1841
1842 if (!pureFreeSurface() && max(nu) > SMALL)
1843 {
1844 tangentialSurfaceTensionForce =
1845 surfaceTensionGrad()().internalField();
1846 }
1847
1848 SnGradU =
1849 tangentialSurfaceTensionForce/(nu + SMALL)
1850 - nA*divUs.internalField()
1851 - (gradUs.internalField()&nA);
1852
1853 return tSnGradU;
1854}
1855
1856
1859{
1860 auto tSnGradUn = tmp<scalarField>::New(aMesh().nFaces(), Zero);
1861 auto& SnGradUn = tSnGradUn.ref();
1862
1863 areaScalarField divUs
1864 (
1865 fac::div(Us())
1866 - aMesh().faceCurvatures()*(aMesh().faceAreaNormals()&Us())
1867 );
1868
1869 SnGradUn = -divUs.internalField();
1870
1871 return tSnGradUn;
1872}
1873
1874
1877{
1878 auto tPressureJump = tmp<scalarField>::New(aMesh().nFaces(), Zero);
1879 auto& pressureJump = tPressureJump.ref();
1880
1882
1885
1887 mesh().lookupObject<turbulenceModel>("turbulenceProperties");
1888
1889 scalarField nu(turbulence.nuEff(fsPatchIndex()));
1890
1891 pressureJump =
1892 - (g.value() & mesh().Cf().boundaryField()[fsPatchIndex()])
1893 + 2.0*nu*freeSurfaceSnGradUn();
1894
1895 if (pureFreeSurface())
1896 {
1897 pressureJump -= sigma().value()*K;
1898 }
1899 else
1900 {
1901 pressureJump -= surfaceTension().internalField()*K;
1902 }
1903
1904 return tPressureJump;
1905}
1906
1907
1909{
1910 if (!controlPointsPtr_)
1911 {
1912 makeControlPoints();
1913 }
1914
1915 return *controlPointsPtr_;
1916}
1917
1918
1920{
1921 if (!motionPointsMaskPtr_)
1922 {
1923 makeMotionPointsMask();
1924 }
1925
1926 return *motionPointsMaskPtr_;
1927}
1928
1929
1931{
1932 if (!pointsDisplacementDirPtr_)
1933 {
1934 makeDirections();
1935 }
1936
1937 return *pointsDisplacementDirPtr_;
1938}
1939
1940
1942{
1943 if (!facesDisplacementDirPtr_)
1944 {
1945 makeDirections();
1946 }
1947
1948 return *facesDisplacementDirPtr_;
1949}
1950
1951
1953{
1954 if (!phisPtr_)
1955 {
1956 makePhis();
1957 }
1958
1959 return *phisPtr_;
1960}
1961
1964{
1965 if (!surfactConcPtr_)
1966 {
1967 makeSurfactConc();
1968 }
1969
1970 return *surfactConcPtr_;
1971}
1972
1973
1976{
1977 if (!surfactConcPtr_)
1978 {
1979 makeSurfactConc();
1980 }
1981
1982 return *surfactConcPtr_;
1983}
1984
1985
1988{
1989 if (!bulkSurfactConcPtr_)
1990 {
1991 makeBulkSurfactConc();
1992 }
1993
1994 return *bulkSurfactConcPtr_;
1995}
1996
1997
2000{
2001 if (!bulkSurfactConcPtr_)
2002 {
2003 makeBulkSurfactConc();
2004 }
2005
2006 return *bulkSurfactConcPtr_;
2007}
2008
2009
2012{
2013 if (!surfaceTensionPtr_)
2014 {
2015 makeSurfaceTension();
2016 }
2017
2018 return *surfaceTensionPtr_;
2019}
2020
2021
2024{
2025 if (!surfaceTensionPtr_)
2026 {
2027 makeSurfaceTension();
2028 }
2029
2030 return *surfaceTensionPtr_;
2031}
2032
2033
2036{
2037 if (!surfactantPtr_)
2038 {
2039 makeSurfactant();
2040 }
2041
2042 return *surfactantPtr_;
2043}
2044
2045
2048{
2049 auto tgrad = tmp<areaVectorField>::New
2050 (
2051 IOobject
2052 (
2053 "surfaceTensionGrad",
2054 mesh().time().timeName(),
2055 mesh(),
2058 ),
2059 fac::grad(surfaceTension())
2060 // (-fac::grad(surfactantConcentration()/rho_)*
2061 // surfactant().surfactR()*surfactant().surfactT()/
2062 // (1.0 - surfactantConcentration()/
2063 // surfactant().surfactSaturatedConc()))()
2064 );
2065 auto& grad = tgrad.ref();
2066
2067 // Remove component of gradient normal to surface (area)
2069 grad -= n*(n & grad);
2070 grad.correctBoundaryConditions();
2071
2072 return tgrad;
2073}
2074
2075
2077{
2078 if (timeIndex_ != mesh().time().timeIndex())
2079 {
2080 if (smoothing_ && !rigidFreeSurface_)
2081 {
2082 smoothFreeSurfaceMesh();
2083 clearControlPoints();
2084 }
2085
2086 updateDisplacementDirections();
2087
2088 updateProperties();
2089
2090 Info<< "Maximal capillary Courant number: "
2091 << maxCourantNumber() << endl;
2092
2094
2095 Info<< "Free surface curvature: min = " << gMin(K)
2096 << ", max = " << gMax(K) << ", average = " << gAverage(K) << nl;
2097
2098 timeIndex_ = mesh().time().timeIndex();
2099 }
2100
2101 if (!rigidFreeSurface_)
2102 {
2103 // This is currently relaltive flux
2104 scalarField sweptVolCorr =
2105 phi().boundaryField()[fsPatchIndex()];
2106
2107 // Info<< "Free surface flux: sum local = "
2108 // << gSum(mag(sweptVolCorr))
2109 // << ", global = " << gSum(sweptVolCorr) << endl;
2110
2111 // if (mesh().moving())
2112 // {
2113 // sweptVolCorr -=
2114 // fvc::meshPhi(U())().boundaryField()[fsPatchIndex()];
2115 // }
2116
2117 Info<< "Free surface continuity error : sum local = "
2118 << gSum(mag(sweptVolCorr)) << ", global = " << gSum(sweptVolCorr)
2119 << endl;
2120
2121 // For postprocessing
2122 fsNetPhi().ref().field() = sweptVolCorr;
2123
2124 word ddtScheme
2125 (
2126 mesh().ddtScheme
2127 (
2128 "ddt(" + U().name() + ')'
2129 )
2130 );
2131
2132 if
2133 (
2134 ddtScheme
2136 )
2137 {
2138 sweptVolCorr *= (1.0/2.0)*mesh().time().deltaT().value();
2139 }
2140 else if
2141 (
2142 ddtScheme
2144 )
2145 {
2146 sweptVolCorr *= mesh().time().deltaT().value();
2147 }
2148 else if
2149 (
2150 ddtScheme
2152 )
2153 {
2154 if (mesh().time().timeIndex() == 1)
2155 {
2156 sweptVolCorr *= mesh().time().deltaT().value();
2157 }
2158 else
2159 {
2160 sweptVolCorr *= (2.0/3.0)*mesh().time().deltaT().value();
2161 }
2162 }
2163 else
2164 {
2166 << "Unsupported temporal differencing scheme : "
2167 << ddtScheme << nl
2168 << abort(FatalError);
2169 }
2170
2171 const scalarField& Sf = aMesh().S();
2173
2174 scalarField deltaHf
2175 (
2176 sweptVolCorr/(Sf*(Nf & facesDisplacementDir()))
2177 );
2178
2179 for (const word& patchName : fixedFreeSurfacePatches_)
2180 {
2181 label fixedPatchID =
2182 aMesh().boundary().findPatchID(patchName);
2183
2184 if (fixedPatchID == -1)
2185 {
2187 << "Wrong faPatch name '" << patchName
2188 << "'in the fixedFreeSurfacePatches list"
2189 << " defined in dynamicMeshDict dictionary"
2190 << abort(FatalError);
2191 }
2192
2193 const labelList& eFaces =
2194 aMesh().boundary()[fixedPatchID].edgeFaces();
2195
2196 forAll(eFaces, edgeI)
2197 {
2198 deltaHf[eFaces[edgeI]] *= 2.0;
2199 }
2200 }
2201
2202 controlPoints() += facesDisplacementDir()*deltaHf;
2203
2204 pointField displacement(pointDisplacement());
2205 correctPointDisplacement(sweptVolCorr, displacement);
2206
2207 velocityMotionSolver& vMotion =
2208 refCast<velocityMotionSolver>
2209 (
2210 const_cast<motionSolver&>(motion())
2211 );
2212
2213 pointVectorField& pointMotionU = vMotion.pointMotionU();
2214 pointMotionU.primitiveFieldRef() = Zero;
2215
2216 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2217 refCast<fixedValuePointPatchVectorField>
2218 (
2219 const_cast<pointPatchVectorField&>
2220 (
2221 pointMotionU.boundaryField()[fsPatchIndex()]
2222 )
2223 );
2224
2225 fsPatchPointMeshU ==
2226 displacement/mesh().time().deltaT().value();
2227
2229
2230 correctContactLinePointNormals();
2231 }
2232 else
2233 {
2234 vectorField displacement
2235 (
2236 mesh().boundaryMesh()[fsPatchIndex()].nPoints(),
2237 Zero
2238 );
2239
2240 velocityMotionSolver& vMotion =
2241 refCast<velocityMotionSolver>
2242 (
2243 const_cast<motionSolver&>(motion())
2244 );
2245
2246 pointVectorField& pointMotionU = vMotion.pointMotionU();
2247 pointMotionU.primitiveFieldRef() = Zero;
2248
2249 fixedValuePointPatchVectorField& fsPatchPointMeshU =
2250 refCast<fixedValuePointPatchVectorField>
2251 (
2252 const_cast<pointPatchVectorField&>
2253 (
2254 pointMotionU.boundaryField()[fsPatchIndex()]
2255 )
2256 );
2257
2258 fsPatchPointMeshU ==
2259 displacement/mesh().time().deltaT().value();
2260
2262 }
2263
2264 updateUs();
2265
2266 updateSurfactantConcentration();
2267
2268 return true;
2269}
2270
2271
2273{
2275 (
2276 aMesh().patch(),
2278 mesh().time().timePath()/"freeSurface",
2279 false // serial only
2280 );
2282}
2283
2284
2286{
2287 // Write control points into VTK
2288 OFstream os
2289 (
2290 mesh().time().timePath()/"freeSurfaceControlPoints.vtk"
2291 );
2292
2293 Info<< "Writing free surface control points to " << os.name() << nl;
2294
2295 os << "# vtk DataFile Version 2.0" << nl
2296 << "freeSurfaceControlPoints" << nl
2297 << "ASCII" << nl
2298 << "DATASET POLYDATA" << nl;
2299
2300 const label nPoints = controlPoints().size();
2301
2302 os << "POINTS " << nPoints << " float" << nl;
2303 for (const point& p : controlPoints())
2304 {
2305 os << float(p.x()) << ' '
2306 << float(p.y()) << ' '
2307 << float(p.z()) << nl;
2308 }
2309
2310 os << "VERTICES " << nPoints << ' ' << 2*nPoints << nl;
2311 for (label id = 0; id < nPoints; ++id)
2312 {
2313 os << 1 << ' ' << id << nl;
2314 }
2315}
2316
2317
2318// ************************************************************************* //
CGAL::Exact_predicates_exact_constructions_kernel K
CGAL::indexedCell< K > Cb
label n
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
Foam::label startTime
const uniformDimensionedVectorField & g
surfaceScalarField & phi
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep track of the Delaunay cells (tets)...
Definition: indexedCell.H:95
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
const Internal & internalField() const
Return a const-reference to the dimensioned internal field.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
List< Key > toc() const
The table of contents (the keys) in unsorted order.
Definition: HashTable.C:122
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Output to file stream, using an OSstream.
Definition: OFstream.H:57
virtual const fileName & name() const
Read/write access to the name of the stream.
Definition: OSstream.H:107
const labelListList & pointEdges() const
Return point-edge addressing.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const labelListList & pointFaces() const
Return point-face addressing.
const labelListList & edgeFaces() const
Return edge-face addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
label timeIndex() const noexcept
Return current time index.
Definition: TimeStateI.H:37
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:55
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type get(const label i) const
Definition: UList.H:528
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Generic dimensioned Type class.
const Type & value() const
Return const reference to value.
Abstract base class for geometry and/or topology changing fvMesh.
Definition: dynamicFvMesh.H:81
The dynamicMotionSolverFvMesh.
virtual bool update()
Update the mesh for both mesh motion and topology change.
const motionSolver & motion() const
Return the motionSolver.
const edgeScalarField & lPN() const
Return reference to PN geodesic distance.
UPtrList< const labelUList > edgeFaces() const
Return a list of edgeFaces for each patch.
label findPatchID(const word &patchName, const bool allowNotFound=true) const
Find patch index given a name, return -1 if not found.
label start() const
The start label of the edges in the faMesh edges list.
Finite area mesh (used for 2-D non-Euclidian finite area method) defined using a patch of faces on a ...
Definition: faMesh.H:100
const faBoundaryMesh & boundary() const noexcept
Return constant reference to boundary mesh.
Definition: faMeshI.H:38
const faceList & faces() const
Return local faces.
Definition: faMeshI.H:104
const DimensionedField< scalar, areaMesh > & S() const
Return face areas.
Definition: faMesh.C:780
bool correctPatchPointNormals(const label patchID) const
Whether point normals should be corrected for a patch.
Definition: faMesh.C:976
const edgeList & edges() const noexcept
Return local edges with reordered boundary.
Definition: faMeshI.H:92
const uindirectPrimitivePatch & patch() const
Return constant reference to primitive patch.
Definition: faMeshI.H:128
const edgeVectorField & Le() const
Return edge length vectors.
Definition: faMesh.C:735
const areaScalarField & faceCurvatures() const
Return face curvatures.
Definition: faMesh.C:870
const pointField & points() const
Return local points.
Definition: faMeshI.H:86
const vectorField & pointAreaNormals() const
Return point area normals.
Definition: faMesh.C:852
const areaVectorField & faceAreaNormals() const
Return face area normals.
Definition: faMesh.C:830
const areaVectorField & areaCentres() const
Return face centres as areaVectorField.
Definition: faMesh.C:757
const labelList & pointLabels() const
Return patch point labels.
Definition: faPatch.C:273
Template specialisation for scalar faMatrix.
A list of face labels.
Definition: faceSet.H:54
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
label which(const label pointLabel) const
Find local index on face for the point label, same as find()
Definition: faceI.H:163
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:175
label prevLabel(const label i) const
Previous vertex on face.
Definition: faceI.H:181
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
Second-oder Crank-Nicolson implicit ddt using the current and previous time-step fields as well as th...
Basic first-order Euler implicit/explicit ddt using only the current and previous time-step values.
Second-order backward-differencing ddt using the current and two previous time-step values.
The interfaceTrackingFvMesh.
vectorField & controlPoints()
Return control points.
void correctUsBoundaryConditions()
Correct surface velocity boundary conditions.
areaVectorField & Us()
Return free-surface velocity field.
const areaScalarField & fsNetPhi() const
Return free-surface net flux.
const surfaceScalarField & phi() const
Return constant reference to velocity field.
tmp< scalarField > freeSurfaceSnGradUn()
Return free surface normal derivative of normal velocity comp.
vectorField & pointsDisplacementDir()
Return reference to point displacement direction field.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
const volVectorField & U() const
Return constant reference to velocity field.
void writeVTK() const
Write VTK freeSurface mesh.
void updateUs()
Update free-surface velocity field.
areaScalarField & surfaceTension()
Return surface tension field.
vectorField & facesDisplacementDir()
Return reference to control points displacement direction field.
void writeVTKControlPoints()
Write VTK freeSurface control points.
faMesh & aMesh()
Return reference to finite area mesh.
virtual bool update()
Update the mesh for both mesh motion and topology change.
tmp< areaVectorField > surfaceTensionGrad()
Return surface tension gradient.
labelList & motionPointsMask()
Return reference to motion points mask field.
volScalarField & bulkSurfactantConcentration()
Return volume surfactant concentration field.
const volScalarField & p() const
Return constant reference to pressure field.
const surfactantProperties & surfactant() const
Return surfactant properties.
areaScalarField & surfactantConcentration()
Return free-surface surfactant concentration field.
edgeScalarField & Phis()
Return free-surface fluid flux field.
tmp< scalarField > freeSurfacePressureJump()
Return free surface pressure jump.
tmp< vectorField > freeSurfaceSnGradU()
Return free surface normal derivative of velocity.
Virtual base class for mesh motion solver.
Definition: motionSolver.H:61
const Type & lookupObject(const word &name, const bool recursive=false) const
A set of point labels.
Definition: pointSet.H:54
virtual const pointField & oldPoints() const
Return old points (mesh motion)
Definition: polyMesh.C:1133
A simple single-phase transport model based on viscosityModel.
A class for managing temporary objects.
Definition: tmp.H:65
Abstract base class for turbulence models (RAS, LES and laminar).
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Virtual base class for velocity motion solver.
pointVectorField & pointMotionU()
Return reference to the point motion velocity field.
Write concrete PrimitivePatch faces/points (optionally with fields) as a vtp file or a legacy vtk fil...
virtual bool writeGeometry()
Write mesh topology.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
volScalarField & p
rDeltaTY field()
const volScalarField & p0
Definition: EEqn.H:36
faceListList boundary
dynamicFvMesh & mesh
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
compressible::turbulenceModel & turbulence
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
const pointField & points
label nPoints
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
word timeName
Definition: getTimeIndex.H:3
#define DebugInFunction
Report an information message using Foam::Info.
const dimensionedScalar a0
Bohr radius: default SI units: [m].
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
const dimensionedScalar c
Speed of light in a vacuum.
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.
tmp< GeometricField< Type, faPatchField, areaMesh > > div(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Definition: facDiv.C:50
tmp< GeometricField< Type, faPatchField, areaMesh > > edgeIntegrate(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
tmp< GeometricField< typename outerProduct< vector, Type >::type, faPatchField, areaMesh > > grad(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Definition: facGrad.C:56
tmp< faMatrix< Type > > div(const edgeScalarField &flux, const GeometricField< Type, faPatchField, areaMesh > &vf, const word &name)
Definition: famDiv.C:48
tmp< faMatrix< Type > > laplacian(const GeometricField< Type, faPatchField, areaMesh > &vf)
Definition: famLaplacian.C:49
tmp< faMatrix< Type > > ddt(const GeometricField< Type, faPatchField, areaMesh > &vf)
Definition: famDdt.C:48
zeroField Sp(const Foam::zero, const GeometricField< Type, faPatchField, areaMesh > &)
A no-op source.
const std::string patch
OpenFOAM patch number as a std::string.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:47
@ LEGACY_ASCII
Legacy ASCII, legacyAsciiFormatter.
Namespace for OpenFOAM.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Type gSum(const FieldField< Field, Type > &f)
const dimensionSet dimless
Dimensionless.
List< label > labelList
A List of labels.
Definition: List.H:66
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0)
Definition: dimensionSets.H:55
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:59
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimVelocity
tmp< GeometricField< Type, faePatchField, edgeMesh > > linearEdgeInterpolate(const GeometricField< Type, faPatchField, areaMesh > &vf)
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
IOField< vector > vectorIOField
vectorField with IO.
Definition: vectorIOField.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
GeometricField< vector, faPatchField, areaMesh > areaVectorField
Definition: areaFieldsFwd.H:79
GeometricField< scalar, faePatchField, edgeMesh > edgeScalarField
Definition: edgeFieldsFwd.H:63
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
GeometricField< scalar, faPatchField, areaMesh > areaScalarField
Definition: areaFieldsFwd.H:78
errorManip< error > abort(error &err)
Definition: errorManip.H:144
const dimensionSet dimForce
Field< vector > vectorField
Specialisation of Field<T> for vector.
dimensionedScalar neg(const dimensionedScalar &ds)
Type gAverage(const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
error FatalError
Type gMin(const FieldField< Field, Type > &f)
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:60
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
const dimensionSet dimDensity
void deleteDemandDrivenData(DataPtr &dataPtr)
Type gMax(const FieldField< Field, Type > &f)
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:364
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label timeIndex
Definition: getTimeIndex.H:30
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]
Definition: readKivaGrid.H:235
volScalarField & nu
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
Us
faMesh aMesh(mesh)
scalar CoNum
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.
Unit conversion functions.
const Vector< label > N(dict.get< Vector< label > >("N"))