src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/tetMetrics/tetMetrics.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Class
  26     tetMetrics
  27
  28 Description
  29     Implementation of tetrahedral mesh-quality metrics.
  30
  31     Most measures are referenced from:
  32     Jonathan R. Shewchuk. What is a good linear finite element?
  33     Interpolation, Conditioning, Anisotropy and Quality Measures.
  34     Eleventh International Meshing Roundtable (Ithaca, NY), pp 115-126, 2002.
  35
  36     NOTE:
  37     For metrics involving volume the function assumes points (p0-p1-p2)
  38     are in counter-clockwise fashion when viewed from p3.
  39
  40 Author
  41     Sandeep Menon
  42     University of Massachusetts Amherst
  43     All rights reserved
  44
  45 \*----------------------------------------------------------------------------*/
  46
  47 #include "tetMetrics.H"
  48 #include "addToMemberFunctionSelectionTable.H"
  49
  50 namespace Foam
  51 {
  52
  53 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  54
  55 defineTypeNameAndDebug(Knupp,0);
  56 defineTypeNameAndDebug(Dihedral,0);
  57 defineTypeNameAndDebug(cubicMeanRatio,0);
  58 defineTypeNameAndDebug(Frobenius,0);
  59 defineTypeNameAndDebug(PGH,0);
  60 defineTypeNameAndDebug(CSG,0);
  61
  62
  63 addToMemberFunctionSelectionTable(tetMetric, Knupp, metric, Point);
  64 addToMemberFunctionSelectionTable(tetMetric, Dihedral, metric, Point);
  65 addToMemberFunctionSelectionTable(tetMetric, cubicMeanRatio, metric, Point);
  66 addToMemberFunctionSelectionTable(tetMetric, Frobenius, metric, Point);
  67 addToMemberFunctionSelectionTable(tetMetric, PGH, metric, Point);
  68 addToMemberFunctionSelectionTable(tetMetric, CSG, metric, Point);
  69
  70
  71 // Enumeration for tets
  72 label Dihedral::tetEnum[6][4] =
  73 {
  74     {0,1,2,3},
  75     {0,2,3,1},
  76     {0,3,1,2},
  77     {1,2,0,3},
  78     {1,3,0,2},
  79     {2,3,0,1}
  80 };
  81
  82 // * * * * * * * * * * * * * Static Members Functions * * * * * * * * * * *  //
  83
  84 // Tetrahedral mesh-quality metric suggested by Knupp [2003].
  85 scalar Knupp::metric
  86 (
  87     const point& p0,
  88     const point& p1,
  89     const point& p2,
  90     const point& p3
  91 )
  92 {
  93     // Obtain signed tet volume
  94     scalar V = (1.0/6.0)*(((p1 - p0) ^ (p2 - p0)) & (p3 - p0));
  95
  96     // Obtain the magSqr edge-lengths
  97     scalar Le = ((p1-p0) & (p1-p0))
  98               + ((p2-p0) & (p2-p0))
  99               + ((p3-p0) & (p3-p0))
 100               + ((p2-p1) & (p2-p1))
 101               + ((p3-p1) & (p3-p1))
 102               + ((p3-p2) & (p3-p2));
 103
 104     // Return signed quality
 105     return sign(V)*((24.96100588*::cbrt(V*V))/Le);
 106 }
 107
 108
 109 // Minimum dihedral angle among six edges of the tetrahedron. Normalized
 110 // by 70.529 degrees (equilateral tet) and signed by volume.
 111 scalar Dihedral::metric
 112 (
 113     const point& p0,
 114     const point& p1,
 115     const point& p2,
 116     const point& p3
 117 )
 118 {
 119     scalar minAngle = 0.0;
 120     FixedList<scalar,6> cosAngles(1.0);
 121     FixedList<vector,4> pts(vector::zero);
 122
 123     // Assign point-positions
 124     pts[0] = p0;
 125     pts[1] = p1;
 126     pts[2] = p2;
 127     pts[3] = p3;
 128
 129     // Permute over all six edges
 130     for (label i = 0; i < 6; i++)
 131     {
 132         // Normalize the axis
 133         vector v0 = (pts[tetEnum[i][1]] - pts[tetEnum[i][0]])
 134                     /mag(pts[tetEnum[i][1]] - pts[tetEnum[i][0]]);
 135
 136         // Obtain plane-vectors
 137         vector v1 = (pts[tetEnum[i][2]] - pts[tetEnum[i][0]]);
 138         vector v2 = (pts[tetEnum[i][3]] - pts[tetEnum[i][0]]);
 139
 140         v1 -= (v1 & v0)*v0;
 141         v2 -= (v2 & v0)*v0;
 142
 143         cosAngles[i] = acos((v1/mag(v1)) & (v2/mag(v2)));
 144
 145         // Compute minimum angle on-the-fly
 146         if (i > 0)
 147         {
 148             minAngle = minAngle < cosAngles[i] ? minAngle : cosAngles[i];
 149         }
 150         else
 151         {
 152             minAngle = cosAngles[0];
 153         }
 154     }
 155
 156     // Compute signed volume and multiply by the normalized angle
 157     return sign(((p1 - p0) ^ (p2 - p0)) & (p3 - p0))*(minAngle/1.2309632);
 158 }
 159
 160
 161 // Cubic Mean Ratio Tetrahedral mesh metric
 162 // Liu,A. and Joe, B., “On the shape of tetrahedra from bisection”
 163 // Mathematics of Computation, Vol. 63, 1994, pp. 141–154.
 164 scalar cubicMeanRatio::metric
 165 (
 166     const point& p0,
 167     const point& p1,
 168     const point& p2,
 169     const point& p3
 170 )
 171 {
 172     // Obtain signed tet volume
 173     scalar V = (1.0/6.0)*(((p1 - p0) ^ (p2 - p0)) & (p3 - p0));
 174
 175     // Obtain the magSqr edge-lengths
 176     scalar Le = ((p1-p0) & (p1-p0))
 177               + ((p2-p0) & (p2-p0))
 178               + ((p3-p0) & (p3-p0))
 179               + ((p2-p1) & (p2-p1))
 180               + ((p3-p1) & (p3-p1))
 181               + ((p3-p2) & (p3-p2));
 182
 183     // Return signed quality
 184     return  sign(V)*((15552.0*V*V)/(Le*Le*Le));
 185 }
 186
 187
 188 // Tetrahedral mesh-metric based on the Frobenius Condition Number
 189 // Patrick M. Knupp. Matrix Norms & the Condition Number: A General Framework
 190 // to Improve Mesh Quality via Node-Movement. Eighth International Meshing
 191 // Roundtable (Lake Tahoe, California), pages 13–22, October 1999.
 192 scalar Frobenius::metric
 193 (
 194     const point& p0,
 195     const point& p1,
 196     const point& p2,
 197     const point& p3
 198 )
 199 {
 200     // Obtain signed tet volume
 201     scalar V = (1.0/6.0)*(((p1 - p0) ^ (p2 - p0)) & (p3 - p0));
 202
 203     // Obtain the magSqr edge-lengths
 204     scalar Le = ((p1-p0) & (p1-p0))
 205               + ((p2-p0) & (p2-p0))
 206               + ((p3-p0) & (p3-p0))
 207               + ((p2-p1) & (p2-p1))
 208               + ((p3-p1) & (p3-p1))
 209               + ((p3-p2) & (p3-p2));
 210
 211     // Compute magSqr of face-areas
 212     scalar A = magSqr(0.5*((p1-p0) ^ (p2-p0)))
 213              + magSqr(0.5*((p1-p0) ^ (p3-p0)))
 214              + magSqr(0.5*((p2-p0) ^ (p3-p0)))
 215              + magSqr(0.5*((p3-p1) ^ (p2-p1)));
 216
 217     // Return signed quality
 218     return 3.67423461*(V/sqrt((Le/6.0)*(A/4.0)));
 219 }
 220
 221
 222 // Tetrahedral mesh-metric suggested by:
 223 // V. N. Parthasarathy, C. M. Graichen, and A. F. Hathaway.
 224 // Fast Evaluation & Improvement of Tetrahedral 3-D Grid Quality. [1991]
 225 scalar PGH::metric
 226 (
 227     const point& p0,
 228     const point& p1,
 229     const point& p2,
 230     const point& p3
 231 )
 232 {
 233     // Obtain signed tet volume
 234     scalar V = (1.0/6.0)*(((p1 - p0) ^ (p2 - p0)) & (p3 - p0));
 235
 236     // Obtain the magSqr edge-lengths
 237     scalar Le = ((p1-p0) & (p1-p0))
 238               + ((p2-p0) & (p2-p0))
 239               + ((p3-p0) & (p3-p0))
 240               + ((p2-p1) & (p2-p1))
 241               + ((p3-p1) & (p3-p1))
 242               + ((p3-p2) & (p3-p2));
 243
 244     // Return signed quality
 245     return 8.48528137*(V/pow(Le/4.0,1.5));
 246 }
 247
 248
 249 // Metric suggested by:
 250 // Hugues L. de Cougny, Mark S. Shephard, and Marcel K. Georges.
 251 // Explicit Node Point Smoothing Within Octree. Technical Report 10-1990,
 252 // Scientific Computation Research Center, Rensselaer Polytechnic Institute,
 253 // Troy, New York. [1990]
 254 scalar CSG::metric
 255 (
 256     const point& p0,
 257     const point& p1,
 258     const point& p2,
 259     const point& p3
 260 )
 261 {
 262     // Obtain signed tet volume
 263     scalar V = (1.0/6.0)*(((p1 - p0) ^ (p2 - p0)) & (p3 - p0));
 264
 265     // Compute magSqr of face-areas
 266     scalar A = magSqr(0.5*((p1-p0) ^ (p2-p0)))
 267              + magSqr(0.5*((p1-p0) ^ (p3-p0)))
 268              + magSqr(0.5*((p2-p0) ^ (p3-p0)))
 269              + magSqr(0.5*((p3-p1) ^ (p2-p1)));
 270
 271     // Return signed quality
 272     return 6.83852117*(V/pow(A,0.75));
 273 }
 274
 275
 276 } // End namespace Foam
 277
 278 // ************************************************************************* //