src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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 Description
  25     Calculation of shape function product for a tetrahedron
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "tetrahedron.H"
  30 #include "triPointRef.H"
  31 #include "scalarField.H"
  32
  33 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  34
  35 // (Probably very inefficient) minimum containment sphere calculation.
  36 // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
  37 // Sphere ctr is smallest one of
  38 // - tet circumcentre
  39 // - triangle circumcentre
  40 // - edge mids
  41 template<class Point, class PointRef>
  42 Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
  43 (
  44     const scalar tol
  45 ) const
  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>
 247 void Foam::tetrahedron<Point, PointRef>::gradNiSquared
 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>
 268 void Foam::tetrahedron<Point, PointRef>::gradNiDotGradNj
 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>
 298 void Foam::tetrahedron<Point, PointRef>::gradNiGradNi
 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>
 316 void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
 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 // ************************************************************************* //