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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "tetrahedron.H"
29 #include "triPointRef.H"
30 #include "scalarField.H"
31 
32 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
33 
34 template<class Point, class PointRef>
36 (
37  const scalar tol
38 ) const
39 {
40  // (Probably very inefficient) minimum containment sphere calculation.
41  // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
42  // Sphere ctr is smallest one of
43  // - tet circumcentre
44  // - triangle circumcentre
45  // - edge mids
46 
47  const scalar fac = 1 + tol;
48 
49  // Halve order of tolerance for comparisons of sqr.
50  const scalar facSqr = Foam::sqrt(fac);
51 
52 
53  // 1. Circumcentre itself.
54 
55  pointHit pHit(circumCentre());
56  pHit.setHit();
57  scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
58 
59 
60  // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
61  // check if 4th point is inside.
62 
63  {
64  point ctr = triPointRef(a_, b_, c_).circumCentre();
65  scalar radiusSqr = magSqr(ctr - a_);
66 
67  if
68  (
69  radiusSqr < minRadiusSqr
70  && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
71  )
72  {
73  pHit.setMiss(false);
74  pHit.setPoint(ctr);
75  minRadiusSqr = radiusSqr;
76  }
77  }
78  {
79  point ctr = triPointRef(a_, b_, d_).circumCentre();
80  scalar radiusSqr = magSqr(ctr - a_);
81 
82  if
83  (
84  radiusSqr < minRadiusSqr
85  && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
86  )
87  {
88  pHit.setMiss(false);
89  pHit.setPoint(ctr);
90  minRadiusSqr = radiusSqr;
91  }
92  }
93  {
94  point ctr = triPointRef(a_, c_, d_).circumCentre();
95  scalar radiusSqr = magSqr(ctr - a_);
96 
97  if
98  (
99  radiusSqr < minRadiusSqr
100  && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
101  )
102  {
103  pHit.setMiss(false);
104  pHit.setPoint(ctr);
105  minRadiusSqr = radiusSqr;
106  }
107  }
108  {
109  point ctr = triPointRef(b_, c_, d_).circumCentre();
110  scalar radiusSqr = magSqr(ctr - b_);
111 
112  if
113  (
114  radiusSqr < minRadiusSqr
115  && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
116  )
117  {
118  pHit.setMiss(false);
119  pHit.setPoint(ctr);
120  minRadiusSqr = radiusSqr;
121  }
122  }
123 
124 
125  // 3. Try midpoints of edges
126 
127  // mid of edge A-B
128  {
129  point ctr = 0.5*(a_ + b_);
130  scalar radiusSqr = magSqr(a_ - ctr);
131  scalar testRadSrq = facSqr*radiusSqr;
132 
133  if
134  (
135  radiusSqr < minRadiusSqr
136  && magSqr(c_-ctr) <= testRadSrq
137  && magSqr(d_-ctr) <= testRadSrq)
138  {
139  pHit.setMiss(false);
140  pHit.setPoint(ctr);
141  minRadiusSqr = radiusSqr;
142  }
143  }
144 
145  // mid of edge A-C
146  {
147  point ctr = 0.5*(a_ + c_);
148  scalar radiusSqr = magSqr(a_ - ctr);
149  scalar testRadSrq = facSqr*radiusSqr;
150 
151  if
152  (
153  radiusSqr < minRadiusSqr
154  && magSqr(b_-ctr) <= testRadSrq
155  && magSqr(d_-ctr) <= testRadSrq
156  )
157  {
158  pHit.setMiss(false);
159  pHit.setPoint(ctr);
160  minRadiusSqr = radiusSqr;
161  }
162  }
163 
164  // mid of edge A-D
165  {
166  point ctr = 0.5*(a_ + d_);
167  scalar radiusSqr = magSqr(a_ - ctr);
168  scalar testRadSrq = facSqr*radiusSqr;
169 
170  if
171  (
172  radiusSqr < minRadiusSqr
173  && magSqr(b_-ctr) <= testRadSrq
174  && magSqr(c_-ctr) <= testRadSrq
175  )
176  {
177  pHit.setMiss(false);
178  pHit.setPoint(ctr);
179  minRadiusSqr = radiusSqr;
180  }
181  }
182 
183  // mid of edge B-C
184  {
185  point ctr = 0.5*(b_ + c_);
186  scalar radiusSqr = magSqr(b_ - ctr);
187  scalar testRadSrq = facSqr*radiusSqr;
188 
189  if
190  (
191  radiusSqr < minRadiusSqr
192  && magSqr(a_-ctr) <= testRadSrq
193  && magSqr(d_-ctr) <= testRadSrq
194  )
195  {
196  pHit.setMiss(false);
197  pHit.setPoint(ctr);
198  minRadiusSqr = radiusSqr;
199  }
200  }
201 
202  // mid of edge B-D
203  {
204  point ctr = 0.5*(b_ + d_);
205  scalar radiusSqr = magSqr(b_ - ctr);
206  scalar testRadSrq = facSqr*radiusSqr;
207 
208  if
209  (
210  radiusSqr < minRadiusSqr
211  && magSqr(a_-ctr) <= testRadSrq
212  && magSqr(c_-ctr) <= testRadSrq)
213  {
214  pHit.setMiss(false);
215  pHit.setPoint(ctr);
216  minRadiusSqr = radiusSqr;
217  }
218  }
219 
220  // mid of edge C-D
221  {
222  point ctr = 0.5*(c_ + d_);
223  scalar radiusSqr = magSqr(c_ - ctr);
224  scalar testRadSrq = facSqr*radiusSqr;
225 
226  if
227  (
228  radiusSqr < minRadiusSqr
229  && magSqr(a_-ctr) <= testRadSrq
230  && magSqr(b_-ctr) <= testRadSrq
231  )
232  {
233  pHit.setMiss(false);
234  pHit.setPoint(ctr);
235  minRadiusSqr = radiusSqr;
236  }
237  }
238 
239 
240  pHit.setDistance(sqrt(minRadiusSqr));
241 
242  return pHit;
243 }
244 
245 
246 template<class Point, class PointRef>
248 (
249  scalarField& buffer
250 ) const
251 {
252  // Change of sign between face area vector and gradient
253  // does not matter because of square
254 
255  // Warning: Added a mag to produce positive coefficients even if
256  // the tetrahedron is twisted inside out. This is pretty
257  // dangerous, but essential for mesh motion.
258  scalar magVol = Foam::mag(mag());
259 
260  buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
261  buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
262  buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
263  buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
264 }
265 
266 
267 template<class Point, class PointRef>
269 (
270  scalarField& buffer
271 ) const
272 {
273  // Warning. Ordering of edges needs to be the same for a tetrahedron
274  // class, a tetrahedron cell shape model and a tetCell
275 
276  // Warning: Added a mag to produce positive coefficients even if
277  // the tetrahedron is twisted inside out. This is pretty
278  // dangerous, but essential for mesh motion.
279 
280  // Double change of sign between face area vector and gradient
281 
282  scalar magVol = Foam::mag(mag());
283  vector sa = Sa();
284  vector sb = Sb();
285  vector sc = Sc();
286  vector sd = Sd();
287 
288  buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
289  buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
290  buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
291  buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
292  buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
293  buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
294 }
295 
296 
297 template<class Point, class PointRef>
299 (
300  tensorField& buffer
301 ) const
302 {
303  // Change of sign between face area vector and gradient
304  // does not matter because of square
305 
306  scalar magVol = Foam::mag(mag());
307 
308  buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
309  buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
310  buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
311  buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
312 }
313 
314 
315 template<class Point, class PointRef>
317 (
318  tensorField& buffer
319 ) const
320 {
321  // Warning. Ordering of edges needs to be the same for a tetrahedron
322  // class, a tetrahedron cell shape model and a tetCell
323 
324  // Double change of sign between face area vector and gradient
325 
326  scalar magVol = Foam::mag(mag());
327  vector sa = Sa();
328  vector sb = Sb();
329  vector sc = Sc();
330  vector sd = Sd();
331 
332  buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
333  buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
334  buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
335  buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
336  buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
337  buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
338 }
339 
340 
341 template<class Point, class PointRef>
343 {
344  return
345  boundBox
346  (
347  min(a(), min(b(), min(c(), d()))),
348  max(a(), max(b(), max(c(), d())))
349  );
350 }
351 
352 
353 // ************************************************************************* //
Foam::PointHit::setHit
void setHit() noexcept
Set the hit status on.
Definition: PointHit.H:181
Foam::PointHit::setMiss
void setMiss(const bool eligible) noexcept
Set the hit status off and set the eligible miss status.
Definition: PointHit.H:188
Foam::PointHit::setDistance
void setDistance(const scalar d) noexcept
Set the distance.
Definition: PointHit.H:201
Foam::PointHit
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:53
scalarField.H
triPointRef.H
Foam::tetrahedron::gradNiGradNi
void gradNiGradNi(tensorField &buffer) const
Definition: tetrahedron.C:299
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::tetrahedron::gradNiDotGradNj
void gradNiDotGradNj(scalarField &buffer) const
Definition: tetrahedron.C:269
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::Field< scalar >
Foam::tetrahedron::containmentSphere
pointHit containmentSphere(const scalar tol) const
Return (min)containment sphere, i.e. the smallest sphere with.
Definition: tetrahedron.C:36
Foam::PointHit::rawPoint
const point_type & rawPoint() const noexcept
The point, no checks.
Definition: PointHit.H:172
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::tetrahedron::gradNiSquared
void gradNiSquared(scalarField &buffer) const
Fill buffer with shape function products.
Definition: tetrahedron.C:248
fac
Calculate the second temporal derivative.
Foam::PointHit::setPoint
void setPoint(const point_type &p)
Set the point.
Definition: PointHit.H:195
Foam::tetrahedron::gradNiGradNj
void gradNiGradNj(tensorField &buffer) const
Definition: tetrahedron.C:317
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::Vector< scalar >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::boundBox
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:63
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::tetrahedron::bounds
boundBox bounds() const
Calculate the bounding box.
Definition: tetrahedron.C:342
Foam::triPointRef
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
tetrahedron.H