src/dynamicMesh/layerAdditionRemoval/setLayerPairing.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 Description
  26     Remove a layer of cells and prepare addressing data
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "layerAdditionRemoval.H"
  31 #include "polyMesh.H"
  32 #include "primitiveMesh.H"
  33 #include "polyTopoChange.H"
  34 #include "oppositeFace.H"
  35 #include "polyTopoChanger.H"
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 bool Foam::layerAdditionRemoval::setLayerPairing() const
  40 {
  41     // Note:
  42     // This is also the most complex part of the topological change.
  43     // Therefore it will be calculated here and stored as temporary
  44     // data until the actual topological change, after which it will
  45     // be cleared.
  46
  47     // Algorithm for point collapse
  48     // 1)  Go through the master cell layer and for every face of
  49     //     the face zone find the opposite face in the master cell.
  50     //     Check the direction of the opposite face and adjust as
  51     //     necessary.  Check other faces to find an edge defining
  52     //     relative orientation of the two faces and adjust the face
  53     //     as necessary.  Once the face is adjusted, record the
  54     //     addressing between the master and slave vertex layer.
  55
  56     const polyMesh& mesh = topoChanger().mesh();
  57
  58     const labelList& mc =
  59         mesh.faceZones()[faceZoneID_.index()].masterCells();
  60
  61     const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
  62
  63     const boolList& mfFlip =
  64         mesh.faceZones()[faceZoneID_.index()].flipMap();
  65
  66     const faceList& faces = mesh.faces();
  67     const cellList& cells = mesh.cells();
  68
  69     // Grab the local faces from the master zone
  70     const faceList& mlf =
  71         mesh.faceZones()[faceZoneID_.index()]().localFaces();
  72
  73     const labelList& meshPoints =
  74         mesh.faceZones()[faceZoneID_.index()]().meshPoints();
  75
  76     // Create a list of points to collapse for every point of
  77     // the master patch
  78     if (pointsPairingPtr_ || facesPairingPtr_)
  79     {
  80         FatalErrorIn
  81         (
  82             "void Foam::layerAdditionRemoval::setLayerPairing() const"
  83         )   << "Problem with layer pairing data"
  84             << abort(FatalError);
  85     }
  86
  87     pointsPairingPtr_ = new labelList(meshPoints.size(), -1);
  88     labelList& ptc = *pointsPairingPtr_;
  89
  90     facesPairingPtr_ = new labelList(mf.size(), -1);
  91     labelList& ftc = *facesPairingPtr_;
  92 //     Pout << "meshPoints: " << meshPoints << nl
  93 //          << "localPoints: " << mesh.faceZones()[faceZoneID_.index()]().localPoints() << endl;
  94
  95     // For all faces, create the mapping
  96     label nPointErrors = 0;
  97     label nFaceErrors = 0;
  98
  99     forAll (mf, faceI)
 100     {
 101         // Get the local master face
 102         face curLocalFace = mlf[faceI];
 103
 104         // Flip face based on flip index to recover original orientation
 105         if (mfFlip[faceI])
 106         {
 107             curLocalFace = curLocalFace.reverseFace();
 108         }
 109
 110         // Get the opposing face from the master cell
 111         oppositeFace lidFace =
 112             cells[mc[faceI]].opposingFace(mf[faceI], faces);
 113
 114         if (!lidFace.found())
 115         {
 116             // This is not a valid layer; cannot continue
 117             nFaceErrors++;
 118             continue;
 119         }
 120
 121 // Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
 122 //     << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
 123 //     << "curLocalFace: " << curLocalFace << nl
 124 //     << "lidFace: " << lidFace
 125 //     << " master index: " << lidFace.masterIndex()
 126 //     << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
 127
 128         // Grab the opposite face for face collapse addressing
 129         ftc[faceI] = lidFace.oppositeIndex();
 130
 131         // Using the local face insert the points into the lid list
 132         forAll (curLocalFace, pointI)
 133         {
 134             const label clp = curLocalFace[pointI];
 135
 136             if (ptc[clp] == -1)
 137             {
 138                 // Point not mapped yet.  Insert the label
 139                 ptc[clp] = lidFace[pointI];
 140             }
 141             else
 142             {
 143                 // Point mapped from some other face.  Check the label
 144                 if (ptc[clp] != lidFace[pointI])
 145                 {
 146                     nPointErrors++;
 147 //                     Pout<< "Topological error in cell layer pairing.  "
 148 //                         << "This mesh is either topologically incorrect "
 149 //                         << "or the master afce layer is not defined "
 150 //                         << "consistently.  Please check the "
 151 //                         << "face zone flip map." << nl
 152 //                         << "First index: " << ptc[clp]
 153 //                         << " new index: " << lidFace[pointI] << endl;
 154                 }
 155             }
 156         }
 157 //         Pout << "ptc: " << ptc << endl;
 158     }
 159
 160     reduce(nPointErrors, sumOp<label>());
 161     reduce(nFaceErrors, sumOp<label>());
 162
 163     if (nPointErrors > 0 || nFaceErrors > 0)
 164     {
 165         clearAddressing();
 166
 167         return false;
 168     }
 169     else
 170     {
 171         // Valid layer
 172         return true;
 173     }
 174 }
 175
 176
 177 const Foam::labelList& Foam::layerAdditionRemoval::pointsPairing() const
 178 {
 179     if (!pointsPairingPtr_)
 180     {
 181         FatalErrorIn
 182         (
 183             "const labelList& layerAdditionRemoval::pointsPairing() const"
 184         )   << "Problem with layer pairing data for object " << name()
 185             << abort(FatalError);
 186     }
 187
 188     return *pointsPairingPtr_;
 189 }
 190
 191 const Foam::labelList& Foam::layerAdditionRemoval::facesPairing() const
 192 {
 193     if (!facesPairingPtr_)
 194     {
 195         FatalErrorIn
 196         (
 197             "const labelList& layerAdditionRemoval::facesPairing() const"
 198         )   << "Problem with layer pairing data for object " << name()
 199             << abort(FatalError);
 200     }
 201
 202     return *facesPairingPtr_;
 203 }
 204
 205
 206 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 207
 208 void Foam::layerAdditionRemoval::modifyMotionPoints
 209 (
 210     pointField& motionPoints
 211 ) const
 212 {
 213     if (debug)
 214     {
 215         Pout<< "void layerAdditionRemoval::modifyMotionPoints("
 216             << "pointField& motionPoints) const for object "
 217             << name() << " : ";
 218     }
 219
 220     if (debug)
 221     {
 222         Pout << "No motion point adjustment" << endl;
 223     }
 224 }
 225
 226
 227 // ************************************************************************* //