src/meshTools/twoDPointCorrector/twoDPointCorrector.C

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