src/meshTools/twoDPointCorrector/twoDPointCorrector.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Class applies a two-dimensional correction to mesh motion point field.
  26     The correction guarantees that the mesh does not get twisted during motion
  27     and thus introduce a third dimension into a 2-D problem.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "twoDPointCorrector.H"
  32 #include "polyMesh.H"
  33 #include "wedgePolyPatch.H"
  34 #include "emptyPolyPatch.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 const scalar twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 void twoDPointCorrector::calcAddressing() const
  48 {
  49     // Find geometry normal
  50     planeNormalPtr_ = new vector(0, 0, 0);
  51     vector& pn = *planeNormalPtr_;
  52
  53     bool isWedge = false;
  54
  55     // Algorithm:
  56     // Attempt to find wedge patch and work out the normal from it.
  57     // If not found, find an empty patch with faces in it and use the
  58     // normal of the first face.  If neither is found, declare an
  59     // error.
  60
  61     // Try and find a wedge patch
  62     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  63
  64     forAll (patches, patchI)
  65     {
  66         if (isA<wedgePolyPatch>(patches[patchI]))
  67         {
  68             isWedge = true;
  69
  70             pn = refCast<const wedgePolyPatch>(patches[patchI]).centreNormal();
  71
  72             if (polyMesh::debug)
  73             {
  74                 Pout << "Found normal from wedge patch " << patchI;
  75             }
  76
  77             break;
  78         }
  79     }
  80
  81     // Try to find an empty patch with faces
  82     if (!isWedge)
  83     {
  84         forAll (patches, patchI)
  85         {
  86             if (isA<emptyPolyPatch>(patches[patchI]) && patches[patchI].size())
  87             {
  88                 pn = patches[patchI].faceAreas()[0];
  89
  90                 if (polyMesh::debug)
  91                 {
  92                     Pout << "Found normal from empty patch " << patchI;
  93                 }
  94
  95                 break;
  96             }
  97         }
  98     }
  99
 100
 101     if (mag(pn) < VSMALL)
 102     {
 103         FatalErrorIn
 104         (
 105             "twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh, "
 106             "const vector& n)"
 107         )   << "Cannot determine normal vector from patches."
 108             << abort(FatalError);
 109     }
 110     else
 111     {
 112         pn /= mag(pn);
 113     }
 114
 115     if (polyMesh::debug)
 116     {
 117         Pout << " twoDPointCorrector normal: " << pn << endl;
 118     }
 119
 120     // Select edges to be included in check.
 121     normalEdgeIndicesPtr_ = new labelList(mesh_.nEdges());
 122     labelList& neIndices = *normalEdgeIndicesPtr_;
 123
 124     const edgeList& meshEdges = mesh_.edges();
 125
 126     const pointField& meshPoints = mesh_.points();
 127
 128     label nNormalEdges = 0;
 129
 130     forAll (meshEdges, edgeI)
 131     {
 132         vector edgeVector =
 133             meshEdges[edgeI].vec(meshPoints)/
 134             (meshEdges[edgeI].mag(meshPoints) + VSMALL);
 135
 136         if (mag(edgeVector & pn) > edgeOrthogonalityTol)
 137         {
 138             // this edge is normal to the plane. Add it to the list
 139             neIndices[nNormalEdges++] = edgeI;
 140         }
 141     }
 142
 143     neIndices.setSize(nNormalEdges);
 144
 145     // Construction check: number of points in a read 2-D or wedge geometry
 146     // should be odd and the number of edges normal to the plane should be
 147     // exactly half the number of points
 148     if (!isWedge)
 149     {
 150         if (meshPoints.size() % 2 != 0)
 151         {
 152             WarningIn
 153             (
 154                 "twoDPointCorrector::twoDPointCorrector("
 155                 "const polyMesh& mesh, const vector& n)"
 156             )   << "the number of vertices in the geometry "
 157                 << "is odd - this should not be the case for a 2-D case. "
 158                 << "Please check the geometry."
 159                 << endl;
 160         }
 161
 162         if (2*nNormalEdges != meshPoints.size())
 163         {
 164             WarningIn
 165             (
 166                 "twoDPointCorrector::twoDPointCorrector("
 167                 "const polyMesh& mesh, const vector& n)"
 168             )   << "The number of points in the mesh is not equal "
 169                 << "to twice the number of edges normal to the plane "
 170                 << "- this may be OK only for wedge geometries.\n"
 171                 << "    Please check the geometry or adjust "
 172                 << "the orthogonality tolerance.\n" << endl
 173                 << "Number of normal edges: " << nNormalEdges
 174                 << " number of points: " << meshPoints.size()
 175                 << endl;
 176         }
 177     }
 178 }
 179
 180
 181 void twoDPointCorrector::clearAddressing() const
 182 {
 183     deleteDemandDrivenData(planeNormalPtr_);
 184     deleteDemandDrivenData(normalEdgeIndicesPtr_);
 185 }
 186
 187
 188 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 189
 190 twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
 191 :
 192     mesh_(mesh),
 193     required_(mesh_.nGeometricD() == 2),
 194     planeNormalPtr_(NULL),
 195     normalEdgeIndicesPtr_(NULL)
 196 {}
 197
 198
 199
 200 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 201
 202 Foam::twoDPointCorrector::~twoDPointCorrector()
 203 {
 204     clearAddressing();
 205 }
 206
 207
 208 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 209
 210 direction twoDPointCorrector::normalDir() const
 211 {
 212     const vector& pn = planeNormal();
 213
 214     if (mag(pn.x()) >= edgeOrthogonalityTol)
 215     {
 216         return vector::X;
 217     }
 218     else if (mag(pn.y()) >= edgeOrthogonalityTol)
 219     {
 220         return vector::Y;
 221     }
 222     else if (mag(pn.z()) >= edgeOrthogonalityTol)
 223     {
 224         return vector::Z;
 225     }
 226     else
 227     {
 228         FatalErrorIn("direction twoDPointCorrector::normalDir() const")
 229             << "Plane normal not aligned with the coordinate system" << nl
 230             << "    pn = " << pn
 231             << abort(FatalError);
 232
 233         return vector::Z;
 234     }
 235 }
 236
 237
 238 // Return plane normal
 239 const vector& twoDPointCorrector::planeNormal() const
 240 {
 241     if (!planeNormalPtr_)
 242     {
 243         calcAddressing();
 244     }
 245
 246     return *planeNormalPtr_;
 247 }
 248
 249
 250 // Return indices of normal edges.
 251 const labelList& twoDPointCorrector::normalEdgeIndices() const
 252 {
 253     if (!normalEdgeIndicesPtr_)
 254     {
 255         calcAddressing();
 256     }
 257
 258     return *normalEdgeIndicesPtr_;
 259 }
 260
 261
 262 void twoDPointCorrector::correctPoints(pointField& p) const
 263 {
 264     if (!required_) return;
 265
 266     // Algorithm:
 267     // Loop through all edges. Calculate the average point position A for
 268     // the front and the back. Correct the position of point P (in two planes)
 269     // such that vectors AP and planeNormal are parallel
 270
 271     // Get reference to edges
 272     const edgeList&  meshEdges = mesh_.edges();
 273
 274     const labelList& neIndices = normalEdgeIndices();
 275     const vector& pn = planeNormal();
 276
 277     forAll (neIndices, edgeI)
 278     {
 279         point& pStart = p[meshEdges[neIndices[edgeI]].start()];
 280
 281         point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
 282
 283         // calculate average point position
 284         const point A = 0.5*(pStart + pEnd);
 285
 286         // correct point locations
 287         pStart = A + pn*(pn & (pStart - A));
 288         pEnd = A + pn*(pn & (pEnd - A));
 289     }
 290 }
 291
 292
 293 void twoDPointCorrector::updateMesh()
 294 {
 295     clearAddressing();
 296 }
 297
 298
 299 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 300
 301 } // End namespace Foam
 302
 303 // ************************************************************************* //