applications/utilities/mesh/manipulation/mirrorMesh/mirrorFvMesh.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 "mirrorFvMesh.H"
  28 #include "Time.H"
  29 #include "plane.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 const Foam::label Foam::mirrorFvMesh::cellRenumber[8][8] =
  34 {
  35     {-1, -1, -1, -1, -1, -1, -1, -1},    // unknown
  36     {-1, -1, -1, -1, -1, -1, -1, -1},    //
  37     {-1, -1, -1, -1, -1, -1, -1, -1},    //
  38     { 0,  3,  2,  1,  4,  7,  6,  5},    // hex
  39     { 2,  1,  0,  5,  4,  3,  6, -1},    // wedge
  40     { 0,  2,  1,  3,  5,  4, -1, -1},    // prism
  41     { 0,  3,  2,  1,  4, -1, -1, -1},    // pyramid
  42     { 2,  1,  0,  3, -1, -1, -1, -1},    // tet
  43 };
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 Foam::mirrorFvMesh::mirrorFvMesh(const IOobject& io)
  48 :
  49     fvMesh(io),
  50     mirrorMeshDict_
  51     (
  52         IOobject
  53         (
  54             "mirrorMeshDict",
  55             time().system(),
  56             *this,
  57             IOobject::MUST_READ,
  58             IOobject::NO_WRITE
  59         )
  60     ),
  61     mirrorMeshPtr_(NULL)
  62 {
  63     plane mirrorPlane(mirrorMeshDict_);
  64
  65     scalar planeTolerance
  66     (
  67         readScalar(mirrorMeshDict_.lookup("planeTolerance"))
  68     );
  69
  70     const pointField& oldPoints = points();
  71     const faceList& oldFaces = faces();
  72     const cellList& oldCells = cells();
  73     const label nOldInternalFaces = nInternalFaces();
  74     const polyPatchList& oldPatches = boundaryMesh();
  75
  76     // Mirror the points
  77     Info << "Mirroring points. Old points: " << oldPoints.size();
  78
  79     pointField newPoints(2*oldPoints.size());
  80     label nNewPoints = 0;
  81
  82     labelList mirrorPointLookup(oldPoints.size(), -1);
  83
  84     // Grab the old points
  85     forAll (oldPoints, pointI)
  86     {
  87         newPoints[nNewPoints] = oldPoints[pointI];
  88         nNewPoints++;
  89     }
  90
  91     forAll (oldPoints, pointI)
  92     {
  93         scalar alpha =
  94             mirrorPlane.normalIntersect
  95             (
  96                 oldPoints[pointI],
  97                 mirrorPlane.normal()
  98             );
  99
 100         // Check plane on tolerance
 101         if (mag(alpha) > planeTolerance)
 102         {
 103             // The point gets mirrored
 104             newPoints[nNewPoints] =
 105                 oldPoints[pointI] + 2.0*alpha*mirrorPlane.normal();
 106
 107             // remember the point correspondence
 108             mirrorPointLookup[pointI] = nNewPoints;
 109             nNewPoints++;
 110         }
 111         else
 112         {
 113             // The point is on the plane and does not get mirrored
 114             // Adjust plane location
 115             newPoints[nNewPoints] =
 116                 oldPoints[pointI] + alpha*mirrorPlane.normal();
 117
 118             mirrorPointLookup[pointI] = pointI;
 119         }
 120     }
 121
 122     // Reset the size of the point list
 123     Info << " New points: " << nNewPoints << endl;
 124     newPoints.setSize(nNewPoints);
 125
 126     Info << "Mirroring faces. Old faces: " << oldFaces.size();
 127
 128     // Algorithm:
 129     // During mirroring, the faces that were previously boundary faces
 130     // in the mirror plane may become ineternal faces. In order to
 131     // deal with the ordering of the faces, the algorithm is split
 132     // into two parts.  For original faces, the internal faces are
 133     // distributed to their owner cells.  Once all internal faces are
 134     // distributed, the boundary faces are visited and if they are in
 135     // the mirror plane they are added to the master cells (the future
 136     // boundary faces are not touched).  After the first phase, the
 137     // internal faces are collected in the cell order and numbering
 138     // information is added.  Then, the internal faces are mirrored
 139     // and the face numbering data is stored for the mirrored section.
 140     // Once all the internal faces are mirrored, the boundary faces
 141     // are added by mirroring the faces patch by patch.
 142
 143     // Distribute internal faces
 144     labelListList newCellFaces(oldCells.size());
 145
 146     const unallocLabelList& oldOwnerStart = lduAddr().ownerStartAddr();
 147
 148     forAll (newCellFaces, cellI)
 149     {
 150         labelList& curFaces = newCellFaces[cellI];
 151
 152         const label s = oldOwnerStart[cellI];
 153         const label e = oldOwnerStart[cellI + 1];
 154
 155         curFaces.setSize(e - s);
 156
 157         forAll (curFaces, i)
 158         {
 159             curFaces[i] = s + i;
 160         }
 161     }
 162
 163     // Distribute boundary faces.  Remember the faces that have been inserted
 164     // as internal
 165     boolListList insertedBouFace(oldPatches.size());
 166
 167     forAll (oldPatches, patchI)
 168     {
 169         const polyPatch& curPatch = oldPatches[patchI];
 170         boolList& curInsBouFace = insertedBouFace[patchI];
 171
 172         curInsBouFace.setSize(curPatch.size());
 173         curInsBouFace = false;
 174
 175         // Get faceCells for face insertion
 176         const unallocLabelList& curFaceCells = curPatch.faceCells();
 177         const label curStart = curPatch.start();
 178
 179         forAll (curPatch, faceI)
 180         {
 181             // Find out if the mirrored face is identical to the
 182             // original.  If so, the face needs to become internal and
 183             // added to its owner cell
 184             const face& origFace = curPatch[faceI];
 185
 186             face mirrorFace(origFace.size());
 187             forAll (mirrorFace, pointI)
 188             {
 189                 mirrorFace[pointI] = mirrorPointLookup[origFace[pointI]];
 190             }
 191
 192             if (origFace == mirrorFace)
 193             {
 194                 // The mirror is identical to current face.  This will
 195                 // become an internal face
 196                 const label oldSize = newCellFaces[curFaceCells[faceI]].size();
 197
 198                 newCellFaces[curFaceCells[faceI]].setSize(oldSize + 1);
 199                 newCellFaces[curFaceCells[faceI]][oldSize] = curStart + faceI;
 200
 201                 curInsBouFace[faceI] = true;
 202             }
 203         }
 204     }
 205
 206     // Construct the new list of faces.  Boundary faces are added
 207     // last, cush that each patch is mirrored separately.  The
 208     // addressing is stored in two separate arrays: first for the
 209     // original cells (face order has changed) and then for the
 210     // mirrored cells.
 211     labelList masterFaceLookup(oldFaces.size(), -1);
 212     labelList mirrorFaceLookup(oldFaces.size(), -1);
 213
 214     faceList newFaces(2*oldFaces.size());
 215     label nNewFaces = 0;
 216
 217     // Insert original (internal) faces
 218     forAll (newCellFaces, cellI)
 219     {
 220         const labelList& curCellFaces = newCellFaces[cellI];
 221
 222         forAll (curCellFaces, cfI)
 223         {
 224             newFaces[nNewFaces] = oldFaces[curCellFaces[cfI]];
 225             masterFaceLookup[curCellFaces[cfI]] = nNewFaces;
 226
 227             nNewFaces++;
 228         }
 229     }
 230
 231     // Mirror internal faces
 232     for (label faceI = 0; faceI < nOldInternalFaces; faceI++)
 233     {
 234         const face& oldFace = oldFaces[faceI];
 235         face& nf = newFaces[nNewFaces];
 236         nf.setSize(oldFace.size());
 237
 238         nf[0] = mirrorPointLookup[oldFace[0]];
 239
 240         for (label i = 1; i < oldFace.size(); i++)
 241         {
 242             nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
 243         }
 244
 245         mirrorFaceLookup[faceI] = nNewFaces;
 246         nNewFaces++;
 247     }
 248
 249     // Mirror boundary faces patch by patch
 250
 251     wordList newPatchTypes(boundary().size());
 252     wordList newPatchNames(boundary().size());
 253     labelList newPatchSizes(boundary().size(), -1);
 254     labelList newPatchStarts(boundary().size(), -1);
 255     label nNewPatches = 0;
 256
 257     forAll (boundaryMesh(), patchI)
 258     {
 259         const label curPatchSize = boundaryMesh()[patchI].size();
 260         const label curPatchStart = boundaryMesh()[patchI].start();
 261         const boolList& curInserted = insertedBouFace[patchI];
 262
 263         newPatchStarts[nNewPatches] = nNewFaces;
 264
 265         // Master side
 266         for (label faceI = 0; faceI < curPatchSize; faceI++)
 267         {
 268             // Check if the face has already been added.  If not, add it and
 269             // insert the numbering details.
 270             if (!curInserted[faceI])
 271             {
 272                 newFaces[nNewFaces] = oldFaces[curPatchStart + faceI];
 273
 274                 masterFaceLookup[curPatchStart + faceI] = nNewFaces;
 275                 nNewFaces++;
 276             }
 277         }
 278
 279         // Mirror side
 280         for (label faceI = 0; faceI < curPatchSize; faceI++)
 281         {
 282             // Check if the face has already been added.  If not, add it and
 283             // insert the numbering details.
 284             if (!curInserted[faceI])
 285             {
 286                 const face& oldFace = oldFaces[curPatchStart + faceI];
 287                 face& nf = newFaces[nNewFaces];
 288                 nf.setSize(oldFace.size());
 289
 290                 nf[0] = mirrorPointLookup[oldFace[0]];
 291
 292                 for (label i = 1; i < oldFace.size(); i++)
 293                 {
 294                     nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
 295                 }
 296
 297                 mirrorFaceLookup[curPatchStart + faceI] = nNewFaces;
 298                 nNewFaces++;
 299             }
 300             else
 301             {
 302                 // Grab the index of the master face for the mirror side
 303                 mirrorFaceLookup[curPatchStart + faceI] =
 304                     masterFaceLookup[curPatchStart + faceI];
 305             }
 306         }
 307
 308         // If patch exists, grab the name and type of the original patch
 309         if (nNewFaces > newPatchStarts[nNewPatches])
 310         {
 311             newPatchTypes[nNewPatches] = boundaryMesh()[patchI].type();
 312             newPatchNames[nNewPatches] = boundaryMesh()[patchI].name();
 313             newPatchSizes[nNewPatches] =
 314                 nNewFaces - newPatchStarts[nNewPatches];
 315
 316             nNewPatches++;
 317         }
 318     }
 319
 320     // Tidy up the lists
 321     newFaces.setSize(nNewFaces);
 322     Info << " New faces: " << nNewFaces << endl;
 323
 324     newPatchTypes.setSize(nNewPatches);
 325     newPatchNames.setSize(nNewPatches);
 326     newPatchSizes.setSize(nNewPatches);
 327     newPatchStarts.setSize(nNewPatches);
 328
 329     Info << "Mirroring patches. Old patches: " << boundary().size()
 330         << " New patches: " << nNewPatches << endl;
 331
 332     Info<< "Mirroring cells.  Old cells: " << oldCells.size()
 333         << " New cells: " << 2*oldCells.size() << endl;
 334
 335     cellList newCells(2*oldCells.size());
 336     label nNewCells = 0;
 337
 338     // Grab the original cells.  Take care of face renumbering.
 339     forAll (oldCells, cellI)
 340     {
 341         const cell& oc = oldCells[cellI];
 342
 343         cell& nc = newCells[nNewCells];
 344         nc.setSize(oc.size());
 345
 346         forAll (oc, i)
 347         {
 348             nc[i] = masterFaceLookup[oc[i]];
 349         }
 350
 351         nNewCells++;
 352     }
 353
 354     // Mirror the cells
 355     forAll (oldCells, cellI)
 356     {
 357         const cell& oc = oldCells[cellI];
 358
 359         cell& nc = newCells[nNewCells];
 360         nc.setSize(oc.size());
 361
 362         forAll (oc, i)
 363         {
 364             nc[i] = mirrorFaceLookup[oc[i]];
 365         }
 366
 367         nNewCells++;
 368     }
 369
 370     // Mirror the cell shapes
 371     Info << "Mirroring cell shapes." << endl;
 372
 373     Info << nl << "Creating new mesh" << endl;
 374     mirrorMeshPtr_ = new fvMesh
 375     (
 376         io,
 377         xferMove(newPoints),
 378         xferMove(newFaces),
 379         xferMove(newCells)
 380     );
 381
 382     fvMesh& pMesh = *mirrorMeshPtr_;
 383
 384     // Add the boundary patches
 385     List<polyPatch*> p(newPatchTypes.size());
 386
 387     forAll (p, patchI)
 388     {
 389         p[patchI] = polyPatch::New
 390         (
 391             newPatchTypes[patchI],
 392             newPatchNames[patchI],
 393             newPatchSizes[patchI],
 394             newPatchStarts[patchI],
 395             patchI,
 396             pMesh.boundaryMesh()
 397         ).ptr();
 398     }
 399
 400     pMesh.addPatches(p);
 401 }
 402
 403
 404 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 405
 406 Foam::mirrorFvMesh::~mirrorFvMesh()
 407 {}
 408
 409
 410 // ************************************************************************* //