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