levelSetTemplates.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 
33 template<class Type>
35 (
36  const fvMesh& mesh,
37  const scalarField& levelC,
38  const scalarField& levelP,
39  const Field<Type>& positiveC,
40  const Field<Type>& positiveP,
41  const Field<Type>& negativeC,
42  const Field<Type>& negativeP
43 )
44 {
45  tmp<Field<Type>> tResult(new Field<Type>(mesh.nCells(), Zero));
46  Field<Type>& result = tResult.ref();
47 
48  forAll(result, cI)
49  {
50  const List<tetIndices> cellTetIs =
51  polyMeshTetDecomposition::cellTetIndices(mesh, cI);
52 
53  scalar v = 0;
54  Type r = Zero;
55 
56  forAll(cellTetIs, cellTetI)
57  {
58  const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
59 
60  const FixedList<point, 4>
61  tet =
62  {
63  mesh.cellCentres()[cI],
64  mesh.points()[triIs[0]],
65  mesh.points()[triIs[1]],
66  mesh.points()[triIs[2]]
67  };
68  const FixedList<scalar, 4>
69  level =
70  {
71  levelC[cI],
72  levelP[triIs[0]],
73  levelP[triIs[1]],
74  levelP[triIs[2]]
75  };
76  const cutTriTet::volumeIntegrateOp<Type>
77  positive = FixedList<Type, 4>
78  ({
79  positiveC[cI],
80  positiveP[triIs[0]],
81  positiveP[triIs[1]],
82  positiveP[triIs[2]]
83  });
84  const cutTriTet::volumeIntegrateOp<Type>
85  negative = FixedList<Type, 4>
86  ({
87  negativeC[cI],
88  negativeP[triIs[0]],
89  negativeP[triIs[1]],
90  negativeP[triIs[2]]
91  });
92 
93  v += cutTriTet::volumeOp()(tet);
94 
95  r += tetCut(tet, level, positive, negative);
96  }
97 
98  result[cI] = r/v;
99  }
100 
101  return tResult;
102 }
103 
104 
105 template<class Type>
107 (
108  const fvPatch& patch,
109  const scalarField& levelF,
110  const scalarField& levelP,
111  const Field<Type>& positiveF,
112  const Field<Type>& positiveP,
113  const Field<Type>& negativeF,
114  const Field<Type>& negativeP
115 )
116 {
117  tmp<Field<Type>> tResult(new Field<Type>(patch.size(), Zero));
118  Field<Type>& result = tResult.ref();
119 
120  forAll(result, fI)
121  {
122  const face& f = patch.patch().localFaces()[fI];
123 
124  scalar a = 0;
125  Type r = Zero;
126 
127  for(label eI = 0; eI < f.size(); ++ eI)
128  {
129  const edge e = f.faceEdge(eI);
130 
131  const FixedList<point, 3>
132  tri =
133  {
134  patch.patch().faceCentres()[fI],
135  patch.patch().localPoints()[e[0]],
136  patch.patch().localPoints()[e[1]]
137  };
138  const FixedList<scalar, 3>
139  level =
140  {
141  levelF[fI],
142  levelP[e[0]],
143  levelP[e[1]]
144  };
145  const cutTriTet::areaMagIntegrateOp<Type>
146  positive = FixedList<Type, 3>
147  ({
148  positiveF[fI],
149  positiveP[e[0]],
150  positiveP[e[1]]
151  });
152  const cutTriTet::areaMagIntegrateOp<Type>
153  negative = FixedList<Type, 3>
154  ({
155  negativeF[fI],
156  negativeP[e[0]],
157  negativeP[e[1]]
158  });
159 
160  a += cutTriTet::areaMagOp()(tri);
161 
162  r += triCut(tri, level, positive, negative);
163  }
164 
165  result[fI] = r/a;
166  }
167 
168  return tResult;
169 }
170 
171 
172 template<class Type>
175 (
176  const volScalarField& levelC,
177  const pointScalarField& levelP,
178  const VolField<Type>& positiveC,
179  const PointField<Type>& positiveP,
180  const VolField<Type>& negativeC,
181  const PointField<Type>& negativeP
182 )
183 {
184  const fvMesh& mesh = levelC.mesh();
185 
186  tmp<VolField<Type>> tResult
187  (
189  (
190  positiveC.name() + ":levelSetAverage",
191  mesh,
192  dimensioned<Type>("0", positiveC.dimensions(), Zero)
193  )
194  );
195  VolField<Type>& result = tResult.ref();
196 
197  result.primitiveFieldRef() =
199  (
200  mesh,
201  levelC.primitiveField(),
202  levelP.primitiveField(),
203  positiveC.primitiveField(),
204  positiveP.primitiveField(),
205  negativeC.primitiveField(),
206  negativeP.primitiveField()
207  );
208 
209  forAll(mesh.boundary(), patchi)
210  {
211  result.boundaryFieldRef()[patchi] =
213  (
214  mesh.boundary()[patchi],
215  levelC.boundaryField()[patchi],
216  levelP.boundaryField()[patchi].patchInternalField()(),
217  positiveC.boundaryField()[patchi],
218  positiveP.boundaryField()[patchi].patchInternalField()(),
219  negativeC.boundaryField()[patchi],
220  negativeP.boundaryField()[patchi].patchInternalField()()
221  );
222  }
223 
224  return tResult;
225 }
226 
227 
228 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A class for managing temporary objects.
Definition: tmp.H:55
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.
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
PointField< scalar > pointScalarField
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
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.
tmp< Field< Type > > levelSetAverage(const fvMesh &mesh, const scalarField &levelC, const scalarField &levelP, const Field< Type > &positiveC, const Field< Type > &positiveP, const Field< Type > &negativeC, const Field< Type > &negativeP)
Calculate the average value of two fields, one on each side of a level set.
face triFace(3)
labelList f(nPoints)