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-2010 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
  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: " << mesh.faceZones()[faceZoneID_.index()]().localPoints() << endl;
  93
  94     // For all faces, create the mapping
  95     label nPointErrors = 0;
  96     label nFaceErrors = 0;
  97
  98     forAll (mf, faceI)
  99     {
 100         // Get the local master face
 101         face curLocalFace = mlf[faceI];
 102
 103         // Flip face based on flip index to recover original orientation
 104         if (mfFlip[faceI])
 105         {
 106             curLocalFace = curLocalFace.reverseFace();
 107         }
 108
 109         // Get the opposing face from the master cell
 110         oppositeFace lidFace =
 111             cells[mc[faceI]].opposingFace(mf[faceI], faces);
 112
 113         if (!lidFace.found())
 114         {
 115             // This is not a valid layer; cannot continue
 116             nFaceErrors++;
 117             continue;
 118         }
 119
 120 // Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
 121 //     << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
 122 //     << "curLocalFace: " << curLocalFace << nl
 123 //     << "lidFace: " << lidFace
 124 //     << " master index: " << lidFace.masterIndex()
 125 //     << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
 126
 127         // Grab the opposite face for face collapse addressing
 128         ftc[faceI] = lidFace.oppositeIndex();
 129
 130         // Using the local face insert the points into the lid list
 131         forAll (curLocalFace, pointI)
 132         {
 133             const label clp = curLocalFace[pointI];
 134
 135             if (ptc[clp] == -1)
 136             {
 137                 // Point not mapped yet.  Insert the label
 138                 ptc[clp] = lidFace[pointI];
 139             }
 140             else
 141             {
 142                 // Point mapped from some other face.  Check the label
 143                 if (ptc[clp] != lidFace[pointI])
 144                 {
 145                     nPointErrors++;
 146 //                     Pout<< "Topological error in cell layer pairing.  "
 147 //                         << "This mesh is either topologically incorrect "
 148 //                         << "or the master afce layer is not defined "
 149 //                         << "consistently.  Please check the "
 150 //                         << "face zone flip map." << nl
 151 //                         << "First index: " << ptc[clp]
 152 //                         << " new index: " << lidFace[pointI] << endl;
 153                 }
 154             }
 155         }
 156 //         Pout << "ptc: " << ptc << endl;
 157     }
 158
 159     reduce(nPointErrors, sumOp<label>());
 160     reduce(nFaceErrors, sumOp<label>());
 161
 162     if (nPointErrors > 0 || nFaceErrors > 0)
 163     {
 164         clearAddressing();
 165
 166         return false;
 167     }
 168     else
 169     {
 170         // Valid layer
 171         return true;
 172     }
 173 }
 174
 175
 176 const Foam::labelList& Foam::layerAdditionRemoval::pointsPairing() const
 177 {
 178     if (!pointsPairingPtr_)
 179     {
 180         FatalErrorIn
 181         (
 182             "const labelList& layerAdditionRemoval::pointsPairing() const"
 183         )   << "Problem with layer pairing data for object " << name()
 184             << abort(FatalError);
 185     }
 186
 187     return *pointsPairingPtr_;
 188 }
 189
 190 const Foam::labelList& Foam::layerAdditionRemoval::facesPairing() const
 191 {
 192     if (!facesPairingPtr_)
 193     {
 194         FatalErrorIn
 195         (
 196             "const labelList& layerAdditionRemoval::facesPairing() const"
 197         )   << "Problem with layer pairing data for object " << name()
 198             << abort(FatalError);
 199     }
 200
 201     return *facesPairingPtr_;
 202 }
 203
 204
 205 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 206
 207 void Foam::layerAdditionRemoval::modifyMotionPoints
 208 (
 209     pointField& motionPoints
 210 ) const
 211 {
 212     if (debug)
 213     {
 214         Pout<< "void layerAdditionRemoval::modifyMotionPoints("
 215             << "pointField& motionPoints) const for object "
 216             << name() << " : ";
 217     }
 218
 219     if (debug)
 220     {
 221         Pout << "No motion point adjustment" << endl;
 222     }
 223 }
 224
 225
 226 // ************************************************************************* //