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-2019 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 "cut.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  tmp<scalarField> tResult(new scalarField(mesh.nCells(), Zero));
42  scalarField& result = tResult.ref();
43 
44  forAll(result, cI)
45  {
46  const List<tetIndices> cellTetIs =
47  polyMeshTetDecomposition::cellTetIndices(mesh, cI);
48 
49  scalar v = 0, r = 0;
50 
51  forAll(cellTetIs, cellTetI)
52  {
53  const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
54 
56  tet =
57  {
58  mesh.cellCentres()[cI],
59  mesh.points()[triIs[0]],
60  mesh.points()[triIs[1]],
61  mesh.points()[triIs[2]]
62  };
64  level =
65  {
66  levelC[cI],
67  levelP[triIs[0]],
68  levelP[triIs[1]],
69  levelP[triIs[2]]
70  };
71 
72  v += cut::volumeOp()(tet);
73 
74  if (above)
75  {
76  r += tetCut(tet, level, cut::volumeOp(), cut::noOp());
77  }
78  else
79  {
80  r += tetCut(tet, level, cut::noOp(), cut::volumeOp());
81  }
82  }
83 
84  result[cI] = r/v;
85  }
86 
87  return tResult;
88 }
89 
90 
92 (
93  const fvPatch& patch,
94  const scalarField& levelF,
95  const scalarField& levelP,
96  const bool above
97 )
98 {
99  tmp<scalarField> tResult(new scalarField(patch.size(), 0));
100  scalarField& result = tResult.ref();
101 
102  forAll(result, fI)
103  {
104  const face& f = patch.patch().localFaces()[fI];
105 
106  scalar a = 0, r = 0;
107 
108  for(label eI = 0; eI < f.size(); ++ eI)
109  {
110  const edge e = f.faceEdge(eI);
111 
112  const FixedList<point, 3>
113  tri =
114  {
115  patch.patch().faceCentres()[fI],
116  patch.patch().localPoints()[e[0]],
117  patch.patch().localPoints()[e[1]]
118  };
120  level =
121  {
122  levelF[fI],
123  levelP[e[0]],
124  levelP[e[1]]
125  };
126 
127  a += cut::areaMagOp()(tri);
128 
129  if (above)
130  {
131  r += triCut(tri, level, cut::areaMagOp(), cut::noOp());
132  }
133  else
134  {
135  r += triCut(tri, level, cut::noOp(), cut::areaMagOp());
136  }
137  }
138 
139  result[fI] = r/a;
140  }
141 
142  return tResult;
143 }
144 
145 
147 (
148  const volScalarField& levelC,
149  const pointScalarField& levelP,
150  const bool above
151 )
152 {
153  const fvMesh& mesh = levelC.mesh();
154 
155  tmp<volScalarField> tResult
156  (
158  (
159  "levelSetFraction",
160  mesh,
162  )
163  );
164  volScalarField& result = tResult.ref();
165 
166  result.primitiveFieldRef() =
168  (
169  mesh,
170  levelC.primitiveField(),
171  levelP.primitiveField(),
172  above
173  );
174 
175  forAll(mesh.boundary(), patchi)
176  {
177  result.boundaryFieldRef()[patchi] =
179  (
180  mesh.boundary()[patchi],
181  levelC.boundaryField()[patchi],
182  levelP.boundaryField()[patchi].patchInternalField()(),
183  above
184  );
185  }
186 
187  return tResult;
188 }
189 
190 
191 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
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
const Boundary & boundaryField() const
Return const-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: HashTable.H:59
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label nCells() const
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:63
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
const dimensionSet dimless
const Field< PointType > & localPoints() const
Return pointField of points in patch.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
cut::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.
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
const polyPatch & patch() const
Return the polyPatch.
Definition: fvPatch.H:139
static const zero Zero
Definition: zero.H:97
const vectorField & cellCentres() const
virtual label size() const
Return size.
Definition: fvPatch.H:157
volScalarField scalarField(fieldObject, mesh)
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
const Mesh & mesh() const
Return mesh.
labelList f(nPoints)
label patchi
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
A class for managing temporary objects.
Definition: PtrList.H:53
edge faceEdge(const label n) const
Return n-th face edge.
Definition: faceI.H:94
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
cut::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.
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:800
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:303