tetrahedron.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) 2011-2018 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 scalar tol
36 ) const
37 {
38  // (Probably very inefficient) minimum containment sphere calculation.
39  // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
40  // Sphere ctr is smallest one of
41  // - tet circumcentre
42  // - triangle circumcentre
43  // - edge mids
44 
45  const scalar fac = 1 + tol;
46 
47  // Halve order of tolerance for comparisons of sqr.
48  const scalar facSqr = Foam::sqrt(fac);
49 
50 
51  // 1. Circumcentre itself.
52 
53  pointHit pHit(circumCentre());
54  pHit.setHit();
55  scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
56 
57 
58  // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
59  // check if 4th point is inside.
60 
61  {
62  point ctr = triPointRef(a_, b_, c_).circumCentre();
63  scalar radiusSqr = magSqr(ctr - a_);
64 
65  if
66  (
67  radiusSqr < minRadiusSqr
68  && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
69  )
70  {
71  pHit.setMiss(false);
72  pHit.setPoint(ctr);
73  minRadiusSqr = radiusSqr;
74  }
75  }
76  {
77  point ctr = triPointRef(a_, b_, d_).circumCentre();
78  scalar radiusSqr = magSqr(ctr - a_);
79 
80  if
81  (
82  radiusSqr < minRadiusSqr
83  && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
84  )
85  {
86  pHit.setMiss(false);
87  pHit.setPoint(ctr);
88  minRadiusSqr = radiusSqr;
89  }
90  }
91  {
92  point ctr = triPointRef(a_, c_, d_).circumCentre();
93  scalar radiusSqr = magSqr(ctr - a_);
94 
95  if
96  (
97  radiusSqr < minRadiusSqr
98  && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
99  )
100  {
101  pHit.setMiss(false);
102  pHit.setPoint(ctr);
103  minRadiusSqr = radiusSqr;
104  }
105  }
106  {
107  point ctr = triPointRef(b_, c_, d_).circumCentre();
108  scalar radiusSqr = magSqr(ctr - b_);
109 
110  if
111  (
112  radiusSqr < minRadiusSqr
113  && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
114  )
115  {
116  pHit.setMiss(false);
117  pHit.setPoint(ctr);
118  minRadiusSqr = radiusSqr;
119  }
120  }
121 
122 
123  // 3. Try midpoints of edges
124 
125  // mid of edge A-B
126  {
127  point ctr = 0.5*(a_ + b_);
128  scalar radiusSqr = magSqr(a_ - ctr);
129  scalar testRadSrq = facSqr*radiusSqr;
130 
131  if
132  (
133  radiusSqr < minRadiusSqr
134  && magSqr(c_-ctr) <= testRadSrq
135  && magSqr(d_-ctr) <= testRadSrq)
136  {
137  pHit.setMiss(false);
138  pHit.setPoint(ctr);
139  minRadiusSqr = radiusSqr;
140  }
141  }
142 
143  // mid of edge A-C
144  {
145  point ctr = 0.5*(a_ + c_);
146  scalar radiusSqr = magSqr(a_ - ctr);
147  scalar testRadSrq = facSqr*radiusSqr;
148 
149  if
150  (
151  radiusSqr < minRadiusSqr
152  && magSqr(b_-ctr) <= testRadSrq
153  && magSqr(d_-ctr) <= testRadSrq
154  )
155  {
156  pHit.setMiss(false);
157  pHit.setPoint(ctr);
158  minRadiusSqr = radiusSqr;
159  }
160  }
161 
162  // mid of edge A-D
163  {
164  point ctr = 0.5*(a_ + d_);
165  scalar radiusSqr = magSqr(a_ - ctr);
166  scalar testRadSrq = facSqr*radiusSqr;
167 
168  if
169  (
170  radiusSqr < minRadiusSqr
171  && magSqr(b_-ctr) <= testRadSrq
172  && magSqr(c_-ctr) <= testRadSrq
173  )
174  {
175  pHit.setMiss(false);
176  pHit.setPoint(ctr);
177  minRadiusSqr = radiusSqr;
178  }
179  }
180 
181  // mid of edge B-C
182  {
183  point ctr = 0.5*(b_ + c_);
184  scalar radiusSqr = magSqr(b_ - ctr);
185  scalar testRadSrq = facSqr*radiusSqr;
186 
187  if
188  (
189  radiusSqr < minRadiusSqr
190  && magSqr(a_-ctr) <= testRadSrq
191  && magSqr(d_-ctr) <= testRadSrq
192  )
193  {
194  pHit.setMiss(false);
195  pHit.setPoint(ctr);
196  minRadiusSqr = radiusSqr;
197  }
198  }
199 
200  // mid of edge B-D
201  {
202  point ctr = 0.5*(b_ + d_);
203  scalar radiusSqr = magSqr(b_ - ctr);
204  scalar testRadSrq = facSqr*radiusSqr;
205 
206  if
207  (
208  radiusSqr < minRadiusSqr
209  && magSqr(a_-ctr) <= testRadSrq
210  && magSqr(c_-ctr) <= testRadSrq)
211  {
212  pHit.setMiss(false);
213  pHit.setPoint(ctr);
214  minRadiusSqr = radiusSqr;
215  }
216  }
217 
218  // mid of edge C-D
219  {
220  point ctr = 0.5*(c_ + d_);
221  scalar radiusSqr = magSqr(c_ - ctr);
222  scalar testRadSrq = facSqr*radiusSqr;
223 
224  if
225  (
226  radiusSqr < minRadiusSqr
227  && magSqr(a_-ctr) <= testRadSrq
228  && magSqr(b_-ctr) <= testRadSrq
229  )
230  {
231  pHit.setMiss(false);
232  pHit.setPoint(ctr);
233  minRadiusSqr = radiusSqr;
234  }
235  }
236 
237 
238  pHit.setDistance(sqrt(minRadiusSqr));
239 
240  return pHit;
241 }
242 
243 
244 template<class Point, class PointRef>
246 (
247  scalarField& buffer
248 ) const
249 {
250  // Change of sign between face area vector and gradient
251  // does not matter because of square
252 
253  // Warning: Added a mag to produce positive coefficients even if
254  // the tetrahedron is twisted inside out. This is pretty
255  // dangerous, but essential for mesh motion.
256  scalar magVol = Foam::mag(mag());
257 
258  buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
259  buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
260  buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
261  buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
262 }
263 
264 
265 template<class Point, class PointRef>
267 (
268  scalarField& buffer
269 ) const
270 {
271  // Warning. Ordering of edges needs to be the same for a tetrahedron
272  // class, a tetrahedron cell shape model and a tetCell
273 
274  // Warning: Added a mag to produce positive coefficients even if
275  // the tetrahedron is twisted inside out. This is pretty
276  // dangerous, but essential for mesh motion.
277 
278  // Double change of sign between face area vector and gradient
279 
280  scalar magVol = Foam::mag(mag());
281  vector sa = Sa();
282  vector sb = Sb();
283  vector sc = Sc();
284  vector sd = Sd();
285 
286  buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
287  buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
288  buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
289  buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
290  buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
291  buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
292 }
293 
294 
295 template<class Point, class PointRef>
297 (
298  tensorField& buffer
299 ) const
300 {
301  // Change of sign between face area vector and gradient
302  // does not matter because of square
303 
304  scalar magVol = Foam::mag(mag());
305 
306  buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
307  buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
308  buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
309  buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
310 }
311 
312 
313 template<class Point, class PointRef>
315 (
316  tensorField& buffer
317 ) const
318 {
319  // Warning. Ordering of edges needs to be the same for a tetrahedron
320  // class, a tetrahedron cell shape model and a tetCell
321 
322  // Double change of sign between face area vector and gradient
323 
324  scalar magVol = Foam::mag(mag());
325  vector sa = Sa();
326  vector sb = Sb();
327  vector sc = Sc();
328  vector sd = Sd();
329 
330  buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
331  buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
332  buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
333  buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
334  buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
335  buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
336 }
337 
338 
339 template<class Point, class PointRef>
341 {
342  return
343  boundBox
344  (
345  min(a(), min(b(), min(c(), d()))),
346  max(a(), max(b(), max(c(), d())))
347  );
348 }
349 
350 
351 // ************************************************************************* //
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
pointHit containmentSphere(const scalar tol) const
Return (min)containment sphere, i.e. the smallest sphere with.
Definition: tetrahedron.C:34
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void gradNiSquared(scalarField &buffer) const
Fill buffer with shape function products.
Definition: tetrahedron.C:246
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
dimensionedScalar sqrt(const dimensionedScalar &ds)
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
const dimensionedScalar c
Speed of light in a vacuum.
void setHit()
Definition: PointHit.H:169
void setDistance(const scalar d)
Definition: PointHit.H:186
void gradNiGradNi(tensorField &buffer) const
Definition: tetrahedron.C:297
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
This class describes the interaction of a face and a point. It carries the info of a successful hit a...
Definition: PointHit.H:51
void setMiss(const bool eligible)
Definition: PointHit.H:175
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
boundBox bounds() const
Calculate the bounding box.
Definition: tetrahedron.C:340
void gradNiGradNj(tensorField &buffer) const
Definition: tetrahedron.C:315
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
dimensioned< scalar > mag(const dimensioned< Type > &)
void gradNiDotGradNj(scalarField &buffer) const
Definition: tetrahedron.C:267
void setPoint(const Point &p)
Definition: PointHit.H:181