levelSetTemplates.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 | Copyright (C) 2017 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 
33 template<class Type>
35 (
36  const fvMesh& mesh,
37  const scalarField& levelC,
38  const scalarField& levelP,
39  const DimensionedField<Type, volMesh>& positiveC,
40  const DimensionedField<Type, pointMesh>& positiveP,
41  const DimensionedField<Type, volMesh>& negativeC,
42  const DimensionedField<Type, pointMesh>& negativeP
43 )
44 {
45  tmp<DimensionedField<Type, volMesh>> tResult
46  (
47  new DimensionedField<Type, volMesh>
48  (
49  IOobject
50  (
51  positiveC.name() + ":levelSetAverage",
52  mesh.time().timeName(),
53  mesh
54  ),
55  mesh,
56  dimensioned<Type>("0", positiveC.dimensions(), Zero)
57  )
58  );
59  DimensionedField<Type, volMesh>& result = tResult.ref();
60 
61  forAll(result, cI)
62  {
63  const List<tetIndices> cellTetIs =
64  polyMeshTetDecomposition::cellTetIndices(mesh, cI);
65 
66  scalar v = 0;
67  Type r = Zero;
68 
69  forAll(cellTetIs, cellTetI)
70  {
71  const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
72 
73  const FixedList<point, 4>
74  tet =
75  {
76  mesh.cellCentres()[cI],
77  mesh.points()[triIs[0]],
78  mesh.points()[triIs[1]],
79  mesh.points()[triIs[2]]
80  };
81  const FixedList<scalar, 4>
82  level =
83  {
84  levelC[cI],
85  levelP[triIs[0]],
86  levelP[triIs[1]],
87  levelP[triIs[2]]
88  };
89  const cut::volumeIntegrateOp<Type>
90  positive = FixedList<Type, 4>
91  ({
92  positiveC[cI],
93  positiveP[triIs[0]],
94  positiveP[triIs[1]],
95  positiveP[triIs[2]]
96  });
97  const cut::volumeIntegrateOp<Type>
98  negative = FixedList<Type, 4>
99  ({
100  negativeC[cI],
101  negativeP[triIs[0]],
102  negativeP[triIs[1]],
103  negativeP[triIs[2]]
104  });
105 
106  v += cut::volumeOp()(tet);
107 
108  r += tetCut(tet, level, positive, negative);
109  }
110 
111  result[cI] = r/v;
112  }
113 
114  return tResult;
115 }
116 
117 
118 template<class Type>
120 (
121  const fvPatch& patch,
122  const scalarField& levelF,
123  const scalarField& levelP,
124  const Field<Type>& positiveF,
125  const Field<Type>& positiveP,
126  const Field<Type>& negativeF,
127  const Field<Type>& negativeP
128 )
129 {
130  typedef typename outerProduct<Type, vector>::type sumType;
131 
132  tmp<Field<Type>> tResult(new Field<Type>(patch.size(), Zero));
133  Field<Type>& result = tResult.ref();
134 
135  forAll(result, fI)
136  {
137  const face& f = patch.patch().localFaces()[fI];
138 
139  vector a = vector::zero;
140  sumType r = Zero;
141 
142  for(label eI = 0; eI < f.size(); ++ eI)
143  {
144  const edge e = f.faceEdge(eI);
145 
146  const FixedList<point, 3>
147  tri =
148  {
149  patch.patch().faceCentres()[fI],
150  patch.patch().localPoints()[e[0]],
151  patch.patch().localPoints()[e[1]]
152  };
153  const FixedList<scalar, 3>
154  level =
155  {
156  levelF[fI],
157  levelP[e[0]],
158  levelP[e[1]]
159  };
160  const cut::areaIntegrateOp<Type>
161  positive = FixedList<Type, 3>
162  ({
163  positiveF[fI],
164  positiveP[e[0]],
165  positiveP[e[1]]
166  });
167  const cut::areaIntegrateOp<Type>
168  negative = FixedList<Type, 3>
169  ({
170  negativeF[fI],
171  negativeP[e[0]],
172  negativeP[e[1]]
173  });
174 
175  a += cut::areaOp()(tri);
176 
177  r += triCut(tri, level, positive, negative);
178  }
179 
180  result[fI] = a/magSqr(a) & r;
181  }
182 
183  return tResult;
184 }
185 
186 
187 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
type
Types of root.
Definition: Roots.H:52
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
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.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
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.
tmp< DimensionedField< Type, volMesh > > levelSetAverage(const fvMesh &mesh, const scalarField &levelC, const scalarField &levelP, const DimensionedField< Type, volMesh > &positiveC, const DimensionedField< Type, pointMesh > &positiveP, const DimensionedField< Type, volMesh > &negativeC, const DimensionedField< Type, pointMesh > &negativeP)
Calculate the average value of two fields, one on each side of a level set.
dynamicFvMesh & mesh
face triFace(3)
static const zero Zero
Definition: zero.H:91
dimensioned< scalar > magSqr(const dimensioned< Type > &)
volScalarField scalarField(fieldObject, mesh)
A class for managing temporary objects.
Definition: PtrList.H:53