levelSet.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2017-2022 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "levelSet.H"
27 #include "cutTriTet.H"
29 #include "tetIndices.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
34 (
35  const fvMesh& mesh,
36  const scalarField& levelC,
37  const scalarField& levelP,
38  const bool above
39 )
40 {
41  typedef cutTriTet::noOp noOp;
42  typedef cutTriTet::volumeOp volumeOp;
43 
44  tmp<scalarField> tResult(new scalarField(mesh.nCells(), Zero));
45  scalarField& result = tResult.ref();
46 
47  forAll(result, cI)
48  {
49  const List<tetIndices> cellTetIs =
50  polyMeshTetDecomposition::cellTetIndices(mesh, cI);
51 
52  scalar v = 0, r = 0;
53 
54  forAll(cellTetIs, cellTetI)
55  {
56  const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
57 
59  tet =
60  {
61  mesh.cellCentres()[cI],
62  mesh.points()[triIs[0]],
63  mesh.points()[triIs[1]],
64  mesh.points()[triIs[2]]
65  };
67  level =
68  {
69  levelC[cI],
70  levelP[triIs[0]],
71  levelP[triIs[1]],
72  levelP[triIs[2]]
73  };
74 
75  v += volumeOp()(tet);
76 
77  if (above)
78  {
79  r += tetCut(tet, level, volumeOp(), noOp());
80  }
81  else
82  {
83  r += tetCut(tet, level, noOp(), volumeOp());
84  }
85  }
86 
87  result[cI] = r/v;
88  }
89 
90  return tResult;
91 }
92 
93 
95 (
96  const fvPatch& patch,
97  const scalarField& levelF,
98  const scalarField& levelP,
99  const bool above
100 )
101 {
102  typedef cutTriTet::noOp noOp;
103  typedef cutTriTet::areaMagOp areaMagOp;
104 
105  tmp<scalarField> tResult(new scalarField(patch.size(), 0));
106  scalarField& result = tResult.ref();
107 
108  forAll(result, fI)
109  {
110  const face& f = patch.patch().localFaces()[fI];
111 
112  scalar a = 0, r = 0;
113 
114  for(label eI = 0; eI < f.size(); ++ eI)
115  {
116  const edge e = f.faceEdge(eI);
117 
118  const FixedList<point, 3>
119  tri =
120  {
121  patch.patch().faceCentres()[fI],
122  patch.patch().localPoints()[e[0]],
123  patch.patch().localPoints()[e[1]]
124  };
126  level =
127  {
128  levelF[fI],
129  levelP[e[0]],
130  levelP[e[1]]
131  };
132 
133  a += areaMagOp()(tri);
134 
135  if (above)
136  {
137  r += triCut(tri, level, areaMagOp(), noOp());
138  }
139  else
140  {
141  r += triCut(tri, level, noOp(), areaMagOp());
142  }
143  }
144 
145  result[fI] = r/a;
146  }
147 
148  return tResult;
149 }
150 
151 
153 (
154  const volScalarField& levelC,
155  const pointScalarField& levelP,
156  const bool above
157 )
158 {
159  const fvMesh& mesh = levelC.mesh();
160 
161  tmp<volScalarField> tResult
162  (
164  (
165  "levelSetFraction",
166  mesh,
168  )
169  );
170  volScalarField& result = tResult.ref();
171 
172  result.primitiveFieldRef() =
174  (
175  mesh,
176  levelC.primitiveField(),
177  levelP.primitiveField(),
178  above
179  );
180 
181  forAll(mesh.boundary(), patchi)
182  {
183  result.boundaryFieldRef()[patchi] =
185  (
186  mesh.boundary()[patchi],
187  levelC.boundaryField()[patchi],
188  levelP.boundaryField()[patchi].patchInternalField()(),
189  above
190  );
191  }
192 
193  return tResult;
194 }
195 
196 
197 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const Mesh & mesh() const
Return mesh.
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:78
Generic GeometricField class.
Internal::FieldType & primitiveFieldRef()
Return a reference to the primitive field.
const Internal::FieldType & primitiveField() const
Return a const-reference to the primitive field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:91
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:857
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:64
virtual label size() const
Return size.
Definition: fvPatch.H:138
const polyPatch & patch() const
Return the polyPatch.
Definition: fvPatch.H:120
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:276
label nCells() const
const vectorField & cellCentres() const
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:71
volScalarField scalarField(fieldObject, mesh)
label patchi
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
const dimensionedScalar e
Elementary charge.
tmp< scalarField > levelSetFraction(const fvMesh &mesh, const scalarField &levelC, const scalarField &levelP, const bool above)
Calculate the volume-fraction that a level set occupies. This gives the.
Definition: levelSet.C:34
static const zero Zero
Definition: zero.H:97
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const dimensionSet dimless
cutTriTet::opAddResult< AboveOp, BelowOp >::type triCut(const FixedList< Point, 3 > &tri, const FixedList< scalar, 3 > &level, const AboveOp &aboveOp, const BelowOp &belowOp)
Cut a triangle along the zero plane defined by the given levels. Templated.
cutTriTet::opAddResult< AboveOp, BelowOp >::type tetCut(const FixedList< Point, 4 > &tet, const FixedList< scalar, 4 > &level, const AboveOp &aboveOp, const BelowOp &belowOp)
As triCut, but for a tetrahedron.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
labelList f(nPoints)