src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 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 \*---------------------------------------------------------------------------*/
  26
  27 #include "cellPointWeight.H"
  28 #include "polyMesh.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 int Foam::cellPointWeight::debug(debug::debugSwitch("cellPointWeight", 0));
  33 Foam::scalar Foam::cellPointWeight::tol(SMALL);
  34
  35 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  36
  37 void Foam::cellPointWeight::findTetrahedron
  38 (
  39     const polyMesh& mesh,
  40     const vector& position,
  41     const label cellIndex
  42 )
  43 {
  44     if (debug)
  45     {
  46         Pout<< "\nFoam::cellPointWeight::findTetrahedron" << nl
  47             << "position = " << position << nl
  48             << "cellIndex = " << cellIndex << endl;
  49     }
  50
  51     // Initialise closest triangle variables
  52     scalar minUVWClose = VGREAT;
  53     label pointIClose = 0;
  54     label faceClose = 0;
  55
  56     const vector& P0 = mesh.cellCentres()[cellIndex];
  57     const labelList& cellFaces = mesh.cells()[cellIndex];
  58     const scalar cellVolume = mesh.cellVolumes()[cellIndex];
  59
  60     // Find the tet that the point occupies
  61     forAll(cellFaces, faceI)
  62     {
  63         // Decompose each face into triangles, making a tet when
  64         // augmented by the cell centre
  65         const labelList& facePoints = mesh.faces()[cellFaces[faceI]];
  66
  67         label pointI = 1;
  68         while ((pointI + 1) < facePoints.size())
  69         {
  70             // Cartesian co-ordinates of the triangle vertices
  71             const vector& P1 = mesh.points()[facePoints[0]];
  72             const vector& P2 = mesh.points()[facePoints[pointI]];
  73             const vector& P3 = mesh.points()[facePoints[pointI + 1]];
  74
  75             // Edge vectors
  76             const vector e1 = P1 - P0;
  77             const vector e2 = P2 - P0;
  78             const vector e3 = P3 - P0;
  79
  80             // Solve for interpolation weighting factors
  81
  82             // Source term
  83             const vector rhs = position - P0;
  84
  85             // Determinant of coefficients matrix
  86             // Note: if det(A) = 0 the tet is degenerate
  87             const scalar detA =
  88                 e1.x()*e2.y()*e3.z() + e2.x()*e3.y()*e1.z()
  89               + e3.x()*e1.y()*e2.z() - e1.x()*e3.y()*e2.z()
  90               - e2.x()*e1.y()*e3.z() - e3.x()*e2.y()*e1.z();
  91
  92             if (mag(detA/cellVolume) > tol)
  93             {
  94                 // Solve using Cramers' rule
  95                 const scalar u =
  96                 (
  97                    rhs.x()*e2.y()*e3.z() + e2.x()*e3.y()*rhs.z()
  98                   + e3.x()*rhs.y()*e2.z() - rhs.x()*e3.y()*e2.z()
  99                   - e2.x()*rhs.y()*e3.z() - e3.x()*e2.y()*rhs.z()
 100                 )/detA;
 101
 102                 const scalar v =
 103                 (
 104                     e1.x()*rhs.y()*e3.z() + rhs.x()*e3.y()*e1.z()
 105                   + e3.x()*e1.y()*rhs.z() - e1.x()*e3.y()*rhs.z()
 106                   - rhs.x()*e1.y()*e3.z() - e3.x()*rhs.y()*e1.z()
 107                 )/detA;
 108
 109                 const scalar w =
 110                 (
 111                     e1.x()*e2.y()*rhs.z() + e2.x()*rhs.y()*e1.z()
 112                   + rhs.x()*e1.y()*e2.z() - e1.x()*rhs.y()*e2.z()
 113                   - e2.x()*e1.y()*rhs.z() - rhs.x()*e2.y()*e1.z()
 114                 )/detA;
 115
 116                 // Check if point is in tet
 117                 // value = 0 indicates position lies on a tet face
 118                 if
 119                 (
 120                    (u + tol > 0) && (v + tol > 0) && (w + tol > 0)
 121                 && (u + v + w < 1 + tol)
 122                 )
 123                 {
 124                     faceVertices_[0] = facePoints[0];
 125                     faceVertices_[1] = facePoints[pointI];
 126                     faceVertices_[2] = facePoints[pointI + 1];
 127
 128                     weights_[0] = u;
 129                     weights_[1] = v;
 130                     weights_[2] = w;
 131                     weights_[3] = 1.0 - (u + v + w);
 132
 133                     return;
 134                 }
 135                 else
 136                 {
 137                     scalar minU = mag(u);
 138                     scalar minV = mag(v);
 139                     scalar minW = mag(w);
 140                     if (minU > 1.0)
 141                     {
 142                         minU -= 1.0;
 143                     }
 144                     if (minV > 1.0)
 145                     {
 146                         minV -= 1.0;
 147                     }
 148                     if (minW > 1.0)
 149                     {
 150                         minW -= 1.0;
 151                     }
 152                     const scalar minUVW = mag(minU + minV + minW);
 153
 154                     if (minUVW < minUVWClose)
 155                     {
 156                         minUVWClose = minUVW;
 157                         pointIClose = pointI;
 158                         faceClose = faceI;
 159                     }
 160                 }
 161             }
 162
 163             pointI++;
 164         }
 165     }
 166
 167     if (debug)
 168     {
 169         Pout<< "cellPointWeight::findTetrahedron" << nl
 170             << "    Tetrahedron search failed; using closest tet values to "
 171             << "point " << nl << "    cell: " << cellIndex << nl << endl;
 172     }
 173
 174     const labelList& facePointsClose = mesh.faces()[cellFaces[faceClose]];
 175     faceVertices_[0] = facePointsClose[0];
 176     faceVertices_[1] = facePointsClose[pointIClose];
 177     faceVertices_[2] = facePointsClose[pointIClose + 1];
 178
 179     weights_[0] = 0.25;
 180     weights_[1] = 0.25;
 181     weights_[2] = 0.25;
 182     weights_[3] = 0.25;
 183 }
 184
 185
 186 void Foam::cellPointWeight::findTriangle
 187 (
 188     const polyMesh& mesh,
 189     const vector& position,
 190     const label faceIndex
 191 )
 192 {
 193     if (debug)
 194     {
 195         Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
 196             << "position = " << position << nl
 197             << "faceIndex = " << faceIndex << endl;
 198     }
 199
 200     // Initialise closest triangle variables
 201     scalar minUVClose = VGREAT;
 202     label pointIClose = 0;
 203
 204     // Decompose each face into triangles, making a tet when
 205     // augmented by the cell centre
 206     const labelList& facePoints = mesh.faces()[faceIndex];
 207
 208     const scalar faceArea2 = magSqr(mesh.faceAreas()[faceIndex]);
 209
 210     label pointI = 1;
 211     while ((pointI + 1) < facePoints.size())
 212     {
 213         // Cartesian co-ordinates of the triangle vertices
 214         const vector& P1 = mesh.points()[facePoints[0]];
 215         const vector& P2 = mesh.points()[facePoints[pointI]];
 216         const vector& P3 = mesh.points()[facePoints[pointI + 1]];
 217
 218         // Direction vectors
 219         vector v1 = position - P1;
 220         const vector v2 = P2 - P1;
 221         const vector v3 = P3 - P1;
 222
 223         // Plane normal
 224         vector n = v2 ^ v3;
 225         n /= mag(n);
 226
 227         // Remove any offset to plane
 228         v1 -= (n & v1)*v1;
 229
 230         // Helper variables
 231         const scalar d12 = v1 & v2;
 232         const scalar d13 = v1 & v3;
 233         const scalar d22 = v2 & v2;
 234         const scalar d23 = v2 & v3;
 235         const scalar d33 = v3 & v3;
 236
 237         // Determinant of coefficients matrix
 238         // Note: if det(A) = 0 the triangle is degenerate
 239         const scalar detA = d22*d33 - d23*d23;
 240
 241         if (0.25*detA/faceArea2 > tol)
 242         {
 243             // Solve using Cramers' rule
 244             const scalar u = (d12*d33 - d23*d13)/detA;
 245             const scalar v = (d22*d13 - d12*d23)/detA;
 246
 247             // Check if point is in triangle
 248             if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
 249             {
 250                 // Indices of the cell vertices making up the triangle
 251                 faceVertices_[0] = facePoints[0];
 252                 faceVertices_[1] = facePoints[pointI];
 253                 faceVertices_[2] = facePoints[pointI + 1];
 254
 255                 weights_[0] = u;
 256                 weights_[1] = v;
 257                 weights_[2] = 1.0 - (u + v);
 258                 weights_[3] = 0.0;
 259
 260                 return;
 261             }
 262             else
 263             {
 264                 scalar minU = mag(u);
 265                 scalar minV = mag(v);
 266                 if (minU > 1.0)
 267                 {
 268                     minU -= 1.0;
 269                 }
 270                 if (minV > 1.0)
 271                 {
 272                     minV -= 1.0;
 273                 }
 274                 const scalar minUV = mag(minU + minV);
 275
 276                 if (minUV < minUVClose)
 277                 {
 278                     minUVClose = minUV;
 279                     pointIClose = pointI;
 280                 }
 281             }
 282         }
 283
 284         pointI++;
 285     }
 286
 287     if (debug)
 288     {
 289         Pout<< "Foam::cellPointWeight::findTriangle"
 290             << "Triangle search failed; using closest triangle to point" << nl
 291             << "    cell face: " << faceIndex << nl << endl;
 292     }
 293
 294     // Indices of the cell vertices making up the triangle
 295     faceVertices_[0] = facePoints[0];
 296     faceVertices_[1] = facePoints[pointIClose];
 297     faceVertices_[2] = facePoints[pointIClose + 1];
 298
 299     weights_[0] = 1.0/3.0;
 300     weights_[1] = 1.0/3.0;
 301     weights_[2] = 1.0/3.0;
 302     weights_[3] = 0.0;
 303 }
 304
 305
 306 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 307
 308 Foam::cellPointWeight::cellPointWeight
 309 (
 310     const polyMesh& mesh,
 311     const vector& position,
 312     const label cellIndex,
 313     const label faceIndex
 314 )
 315 :
 316     cellIndex_(cellIndex)
 317 {
 318     if (faceIndex < 0)
 319     {
 320         // Face data not supplied
 321         findTetrahedron(mesh, position, cellIndex);
 322     }
 323     else
 324     {
 325         // Face data supplied
 326         findTriangle(mesh, position, faceIndex);
 327     }
 328 }
 329
 330
 331 // ************************************************************************* //