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