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