All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cutTemplates.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 "cut.H"
27 
28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
29 
30 template<class Point, class AboveOp, class BelowOp>
32 (
33  const FixedList<Point, 3>& tri,
34  const FixedList<scalar, 3>& level,
35  const AboveOp& aboveOp,
36  const BelowOp& belowOp
37 )
38 {
39  // Quick return if all levels are zero or have the same sign
40  if (level[0] == 0 && level[1] == 0 && level[2] == 0)
41  {
42  return aboveOp() + belowOp();
43  }
44  if (level[0] >= 0 && level[1] >= 0 && level[2] >= 0)
45  {
46  return aboveOp(tri) + belowOp();
47  }
48  if (level[0] <= 0 && level[1] <= 0 && level[2] <= 0)
49  {
50  return aboveOp() + belowOp(tri);
51  }
52 
53  // There will be just one edge without a sign change. Find it, and put it
54  // opposite the first vertex. This may change the sign of the tri.
55  FixedList<label, 3> indices({0, 1, 2});
56  label i;
57  for (i = 0; i < 3; ++ i)
58  {
59  if (level[(i + 1)%3]*level[(i + 2)%3] >= 0)
60  {
61  Swap(indices[0], indices[i]);
62  break;
63  }
64  }
65  if (i == 3)
66  {
68  << "The number of tri vertices above the level set should always "
69  << "be 1" << exit(FatalError);
70  }
71 
72  // Correct the sign
73  if (indices[0] != 0)
74  {
75  Swap(indices[1], indices[2]);
76  }
77 
78  // Permute the data
79  const FixedList<Point, 3> p = triReorder(tri, indices);
80  const FixedList<scalar, 3> l = triReorder(level, indices);
81  AboveOp a = triReorder(aboveOp, indices);
82  BelowOp b = triReorder(belowOp, indices);
83 
84  // Slice off one corner to form a tri and a quad
85  Pair<scalar> f;
86  for (label i = 0; i < 2; ++ i)
87  {
88  f[i] = l[0]/(l[0] - l[i+1]);
89  }
90  if (l[0] > 0)
91  {
92  return triCutTri(a, p, f) + triCutQuad(b, p, f);
93  }
94  else
95  {
96  return triCutQuad(a, p, f) + triCutTri(b, p, f);
97  }
98 }
99 
100 
101 template<class AboveOp, class BelowOp>
103 (
104  const FixedList<point, 3>& tri,
105  const plane& p,
106  const AboveOp& aboveOp,
107  const BelowOp& belowOp
108 )
109 {
110  // Set the level set to the signed distance from the plane
111  FixedList<scalar, 3> level;
112  for (label i = 0; i < 3; ++ i)
113  {
114  level[i] = (tri[i] - p.refPoint()) & p.normal();
115  }
116 
117  // Run the level set function
118  return triCut(tri, level, aboveOp, belowOp);
119 }
120 
121 
122 template<class Point, class AboveOp, class BelowOp>
124 (
125  const FixedList<Point, 4>& tet,
126  const FixedList<scalar, 4>& level,
127  const AboveOp& aboveOp,
128  const BelowOp& belowOp
129 )
130 {
131  // Get the min and max over all four vertices and quick return if
132  // all levels are zero or have the same sign
133  scalar levelMin = vGreat, levelMax = - vGreat;
134  for (label i = 0; i < 4; ++ i)
135  {
136  levelMin = min(levelMin, level[i]);
137  levelMax = max(levelMax, level[i]);
138  }
139  if (levelMin == 0 && levelMax == 0)
140  {
141  return aboveOp() + belowOp();
142  }
143  if (levelMin >= 0)
144  {
145  return aboveOp(tet) + belowOp();
146  }
147  if (levelMax <= 0)
148  {
149  return aboveOp() + belowOp(tet);
150  }
151 
152  // Partition the level so that positive values are at the start. This is
153  // like a single iteration of quick-sort, except that the pivot is a hard-
154  // coded zero, rather than an element of the array. This can change the sign
155  // of the tet.
156  FixedList<label, 4> indices({0, 1, 2, 3});
157  bool signChange = false;
158  label i = 0, j = 3;
159  while (true)
160  {
161  while (i < j && level[indices[i]] > 0)
162  {
163  i ++;
164  }
165  while (j > i && level[indices[j]] <= 0)
166  {
167  j --;
168  }
169  if (i == j)
170  {
171  break;
172  }
173  Swap(indices[i], indices[j]);
174  signChange = !signChange;
175  }
176 
177  // The number of vertices above the slice
178  label n = i;
179 
180  // If there are more positives than negatives then reverse the order so that
181  // the negatives are at the start
182  if (n > 2)
183  {
184  n = 4 - n;
185  for (label i = 0; i < 2; ++ i)
186  {
187  Swap(indices[i], indices[3-i]);
188  }
189  }
190 
191  // Correct the sign
192  if (signChange)
193  {
194  Swap(indices[2], indices[3]);
195  }
196 
197  // Permute the data
198  const FixedList<Point, 4> p = tetReorder(tet, indices);
199  const FixedList<scalar, 4> l = tetReorder(level, indices);
200  AboveOp a = tetReorder(aboveOp, indices);
201  BelowOp b = tetReorder(belowOp, indices);
202 
203  // Calculate the integrals above and below the level set
204  if (n == 1)
205  {
206  // Slice off one corner to form a tet and a prism
207  FixedList<scalar, 3> f;
208  for (label i = 0; i < 3; ++ i)
209  {
210  f[i] = l[0]/(l[0] - l[i+1]);
211  }
212  if (l[0] > 0)
213  {
214  return tetCutTet(a, p, f) + tetCutPrism0(b, p, f);
215  }
216  else
217  {
218  return tetCutPrism0(a, p, f) + tetCutTet(b, p, f);
219  }
220  }
221  else if (n == 2)
222  {
223  // Slice off two corners to form two prisms
224  FixedList<scalar, 4> f;
225  for (label i = 0; i < 2; ++ i)
226  {
227  for (label j = 0; j < 2; ++ j)
228  {
229  f[2*i+j] = l[i]/(l[i] - l[j+2]);
230  }
231  }
232  if (l[0] > 0)
233  {
234  return tetCutPrism01(a, p, f) + tetCutPrism23(b, p, f);
235  }
236  else
237  {
238  return tetCutPrism23(a, p, f) + tetCutPrism01(b, p, f);
239  }
240  }
241 
243  << "The number of tet vertices above the level set should always be "
244  << "either 1 or 2" << exit(FatalError);
245 
246  return aboveOp() + belowOp();
247 }
248 
249 
250 template<class AboveOp, class BelowOp>
252 (
253  const FixedList<point, 4>& tet,
254  const plane& p,
255  const AboveOp& aboveOp,
256  const BelowOp& belowOp
257 )
258 {
259  // Set the level set to the signed distance from the plane
260  FixedList<scalar, 4> level;
261  for (label i = 0; i < 4; ++ i)
262  {
263  level[i] = (tet[i] - p.refPoint()) & p.normal();
264  }
265 
266  // Run the level set function
267  return tetCut(tet, level, aboveOp, belowOp);
268 }
269 
270 // ************************************************************************* //
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
Trait to determine the result of the addition of two operations.
Definition: cut.H:527
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const cut::uniformOp< Type > & tetCutPrism0(const cut::uniformOp< Type > &x, const FixedList< scalar, 3 > &)
Modify a uniform operation for cutting prism0 from a tet (does nothing)
Definition: cutI.H:98
const cut::uniformOp< Type > & triCutTri(const cut::uniformOp< Type > &x, const Pair< scalar > &)
Modify a uniform operation for cutting a tri from a tri (does nothing)
Definition: cutI.H:50
const cut::uniformOp< Type > & tetCutPrism23(const cut::uniformOp< Type > &x, const FixedList< scalar, 4 > &)
Modify a uniform operation for cutting prism23 from a tet (does nothing)
Definition: cutI.H:122
const cut::uniformOp< Type > & triCutQuad(const cut::uniformOp< Type > &x, const Pair< scalar > &)
Modify a uniform operation for cutting a quad from a tri (does nothing)
Definition: cutI.H:62
void Swap(T &a, T &b)
Definition: Swap.H:43
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.
const cut::uniformOp< Type > & tetCutPrism01(const cut::uniformOp< Type > &x, const FixedList< scalar, 4 > &)
Modify a uniform operation for cutting prism01 from a tet (does nothing)
Definition: cutI.H:110
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const cut::uniformOp< Type > & tetReorder(const cut::uniformOp< Type > &x, const FixedList< label, 4 > &)
Modify a uniform operation for reordering a tet (does nothing)
Definition: cutI.H:74
const cut::uniformOp< Type > & triReorder(const cut::uniformOp< Type > &x, const FixedList< label, 3 > &)
Modify a uniform operation for reordering a tri (does nothing)
Definition: cutI.H:38
const cut::uniformOp< Type > & tetCutTet(const cut::uniformOp< Type > &x, const FixedList< scalar, 3 > &)
Modify a uniform operation for cutting a tet from a tet (does nothing)
Definition: cutI.H:86
label n
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.