tetrahedron.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) 2011-2012 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 "tetrahedron.H"
27 #include "triPointRef.H"
28 #include "scalarField.H"
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
32 template<class Point, class PointRef>
34 (
35  const tetrahedron<Point, PointRef>& tetB,
36  tetIntersectionList& insideTets,
37  label& nInside,
38  tetIntersectionList& outsideTets,
39  label& nOutside
40 ) const
41 {
42  // Work storage
43  tetIntersectionList cutInsideTets;
44  label nCutInside = 0;
45 
46  nInside = 0;
47  storeOp inside(insideTets, nInside);
48  storeOp cutInside(cutInsideTets, nCutInside);
49 
50  nOutside = 0;
51  storeOp outside(outsideTets, nOutside);
52 
53 
54  // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
55  // in aboveTets are inside tetB.
56 
57  {
58  // face0
59  plane pl0(tetB.b_, tetB.d_, tetB.c_);
60 
61  // Cut and insert subtets into cutInsideTets (either by getting
62  // an index from freeSlots or by appending to insideTets) or
63  // insert into outsideTets
64  sliceWithPlane(pl0, cutInside, outside);
65  }
66 
67  if (nCutInside == 0)
68  {
69  nInside = nCutInside;
70  return;
71  }
72 
73  {
74  // face1
75  plane pl1(tetB.a_, tetB.c_, tetB.d_);
76 
77  nInside = 0;
78 
79  for (label i = 0; i < nCutInside; i++)
80  {
81  cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
82  }
83 
84  if (nInside == 0)
85  {
86  return;
87  }
88  }
89 
90  {
91  // face2
92  plane pl2(tetB.a_, tetB.d_, tetB.b_);
93 
94  nCutInside = 0;
95 
96  for (label i = 0; i < nInside; i++)
97  {
98  insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
99  }
100 
101  if (nCutInside == 0)
102  {
103  nInside = nCutInside;
104  return;
105  }
106  }
107 
108  {
109  // face3
110  plane pl3(tetB.a_, tetB.b_, tetB.c_);
111 
112  nInside = 0;
113 
114  for (label i = 0; i < nCutInside; i++)
115  {
116  cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
117  }
118  }
119 }
120 
121 
122 // (Probably very inefficient) minimum containment sphere calculation.
123 // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
124 // Sphere ctr is smallest one of
125 // - tet circumcentre
126 // - triangle circumcentre
127 // - edge mids
128 template<class Point, class PointRef>
130 (
131  const scalar tol
132 ) const
133 {
134  const scalar fac = 1 + tol;
135 
136  // Halve order of tolerance for comparisons of sqr.
137  const scalar facSqr = Foam::sqrt(fac);
138 
139 
140  // 1. Circumcentre itself.
141 
142  pointHit pHit(circumCentre());
143  pHit.setHit();
144  scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
145 
146 
147  // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
148  // check if 4th point is inside.
149 
150  {
151  point ctr = triPointRef(a_, b_, c_).circumCentre();
152  scalar radiusSqr = magSqr(ctr - a_);
153 
154  if
155  (
156  radiusSqr < minRadiusSqr
157  && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
158  )
159  {
160  pHit.setMiss(false);
161  pHit.setPoint(ctr);
162  minRadiusSqr = radiusSqr;
163  }
164  }
165  {
166  point ctr = triPointRef(a_, b_, d_).circumCentre();
167  scalar radiusSqr = magSqr(ctr - a_);
168 
169  if
170  (
171  radiusSqr < minRadiusSqr
172  && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
173  )
174  {
175  pHit.setMiss(false);
176  pHit.setPoint(ctr);
177  minRadiusSqr = radiusSqr;
178  }
179  }
180  {
181  point ctr = triPointRef(a_, c_, d_).circumCentre();
182  scalar radiusSqr = magSqr(ctr - a_);
183 
184  if
185  (
186  radiusSqr < minRadiusSqr
187  && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
188  )
189  {
190  pHit.setMiss(false);
191  pHit.setPoint(ctr);
192  minRadiusSqr = radiusSqr;
193  }
194  }
195  {
196  point ctr = triPointRef(b_, c_, d_).circumCentre();
197  scalar radiusSqr = magSqr(ctr - b_);
198 
199  if
200  (
201  radiusSqr < minRadiusSqr
202  && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
203  )
204  {
205  pHit.setMiss(false);
206  pHit.setPoint(ctr);
207  minRadiusSqr = radiusSqr;
208  }
209  }
210 
211 
212  // 3. Try midpoints of edges
213 
214  // mid of edge A-B
215  {
216  point ctr = 0.5*(a_ + b_);
217  scalar radiusSqr = magSqr(a_ - ctr);
218  scalar testRadSrq = facSqr*radiusSqr;
219 
220  if
221  (
222  radiusSqr < minRadiusSqr
223  && magSqr(c_-ctr) <= testRadSrq
224  && magSqr(d_-ctr) <= testRadSrq)
225  {
226  pHit.setMiss(false);
227  pHit.setPoint(ctr);
228  minRadiusSqr = radiusSqr;
229  }
230  }
231 
232  // mid of edge A-C
233  {
234  point ctr = 0.5*(a_ + c_);
235  scalar radiusSqr = magSqr(a_ - ctr);
236  scalar testRadSrq = facSqr*radiusSqr;
237 
238  if
239  (
240  radiusSqr < minRadiusSqr
241  && magSqr(b_-ctr) <= testRadSrq
242  && magSqr(d_-ctr) <= testRadSrq
243  )
244  {
245  pHit.setMiss(false);
246  pHit.setPoint(ctr);
247  minRadiusSqr = radiusSqr;
248  }
249  }
250 
251  // mid of edge A-D
252  {
253  point ctr = 0.5*(a_ + d_);
254  scalar radiusSqr = magSqr(a_ - ctr);
255  scalar testRadSrq = facSqr*radiusSqr;
256 
257  if
258  (
259  radiusSqr < minRadiusSqr
260  && magSqr(b_-ctr) <= testRadSrq
261  && magSqr(c_-ctr) <= testRadSrq
262  )
263  {
264  pHit.setMiss(false);
265  pHit.setPoint(ctr);
266  minRadiusSqr = radiusSqr;
267  }
268  }
269 
270  // mid of edge B-C
271  {
272  point ctr = 0.5*(b_ + c_);
273  scalar radiusSqr = magSqr(b_ - ctr);
274  scalar testRadSrq = facSqr*radiusSqr;
275 
276  if
277  (
278  radiusSqr < minRadiusSqr
279  && magSqr(a_-ctr) <= testRadSrq
280  && magSqr(d_-ctr) <= testRadSrq
281  )
282  {
283  pHit.setMiss(false);
284  pHit.setPoint(ctr);
285  minRadiusSqr = radiusSqr;
286  }
287  }
288 
289  // mid of edge B-D
290  {
291  point ctr = 0.5*(b_ + d_);
292  scalar radiusSqr = magSqr(b_ - ctr);
293  scalar testRadSrq = facSqr*radiusSqr;
294 
295  if
296  (
297  radiusSqr < minRadiusSqr
298  && magSqr(a_-ctr) <= testRadSrq
299  && magSqr(c_-ctr) <= testRadSrq)
300  {
301  pHit.setMiss(false);
302  pHit.setPoint(ctr);
303  minRadiusSqr = radiusSqr;
304  }
305  }
306 
307  // mid of edge C-D
308  {
309  point ctr = 0.5*(c_ + d_);
310  scalar radiusSqr = magSqr(c_ - ctr);
311  scalar testRadSrq = facSqr*radiusSqr;
312 
313  if
314  (
315  radiusSqr < minRadiusSqr
316  && magSqr(a_-ctr) <= testRadSrq
317  && magSqr(b_-ctr) <= testRadSrq
318  )
319  {
320  pHit.setMiss(false);
321  pHit.setPoint(ctr);
322  minRadiusSqr = radiusSqr;
323  }
324  }
325 
326 
327  pHit.setDistance(sqrt(minRadiusSqr));
328 
329  return pHit;
330 }
331 
332 
333 template<class Point, class PointRef>
335 (
336  scalarField& buffer
337 ) const
338 {
339  // Change of sign between face area vector and gradient
340  // does not matter because of square
341 
342  // Warning: Added a mag to produce positive coefficients even if
343  // the tetrahedron is twisted inside out. This is pretty
344  // dangerous, but essential for mesh motion.
345  scalar magVol = Foam::mag(mag());
346 
347  buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
348  buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
349  buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
350  buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
351 }
352 
353 
354 template<class Point, class PointRef>
356 (
357  scalarField& buffer
358 ) const
359 {
360  // Warning. Ordering of edges needs to be the same for a tetrahedron
361  // class, a tetrahedron cell shape model and a tetCell
362 
363  // Warning: Added a mag to produce positive coefficients even if
364  // the tetrahedron is twisted inside out. This is pretty
365  // dangerous, but essential for mesh motion.
366 
367  // Double change of sign between face area vector and gradient
368 
369  scalar magVol = Foam::mag(mag());
370  vector sa = Sa();
371  vector sb = Sb();
372  vector sc = Sc();
373  vector sd = Sd();
374 
375  buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
376  buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
377  buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
378  buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
379  buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
380  buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
381 }
382 
383 
384 template<class Point, class PointRef>
386 (
387  tensorField& buffer
388 ) const
389 {
390  // Change of sign between face area vector and gradient
391  // does not matter because of square
392 
393  scalar magVol = Foam::mag(mag());
394 
395  buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
396  buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
397  buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
398  buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
399 }
400 
401 
402 template<class Point, class PointRef>
404 (
405  tensorField& buffer
406 ) const
407 {
408  // Warning. Ordering of edges needs to be the same for a tetrahedron
409  // class, a tetrahedron cell shape model and a tetCell
410 
411  // Double change of sign between face area vector and gradient
412 
413  scalar magVol = Foam::mag(mag());
414  vector sa = Sa();
415  vector sb = Sb();
416  vector sc = Sc();
417  vector sd = Sd();
418 
419  buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
420  buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
421  buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
422  buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
423  buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
424  buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
425 }
426 
427 
428 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
void gradNiSquared(scalarField &buffer) const
Fill buffer with shape function products.
Definition: tetrahedron.C:335
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void setPoint(const Point &p)
Definition: PointHit.H:181
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
void gradNiGradNj(tensorField &buffer) const
Definition: tetrahedron.C:404
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:53
pointHit containmentSphere(const scalar tol) const
Return (min)containment sphere, i.e. the smallest sphere with.
Definition: tetrahedron.C:130
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
void setDistance(const scalar d)
Definition: PointHit.H:186
A tetrahedron primitive.
Definition: tetrahedron.H:62
void setHit()
Definition: PointHit.H:169
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
void setMiss(const bool eligible)
Definition: PointHit.H:175
void gradNiGradNi(tensorField &buffer) const
Definition: tetrahedron.C:386
void gradNiDotGradNj(scalarField &buffer) const
Definition: tetrahedron.C:356
Store resulting tets.
Definition: tetrahedron.H:118
void tetOverlap(const tetrahedron< Point, PointRef > &tetB, tetIntersectionList &insideTets, label &nInside, tetIntersectionList &outsideTets, label &nOutside) const
Decompose tet into tets inside and outside other tet.
Definition: tetrahedron.C:34
dimensionedSymmTensor sqr(const dimensionedVector &dv)