applications/utilities/mesh/conversion/fluentMeshToFoam/create3DCellShape.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     Construct a cell shape from face information
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "cellShapeRecognition.H"
  30 #include "labelList.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 cellShape create3DCellShape
  40 (
  41     const label cellIndex,
  42     const labelList& faceLabels,
  43     const faceList& faces,
  44     const labelList& owner,
  45     const labelList& neighbour,
  46     const label fluentCellModelID
  47 )
  48 {
  49     // List of pointers to shape models for 3-D shape recognition
  50     static List<const cellModel*> fluentCellModelLookup
  51     (
  52         7,
  53         reinterpret_cast<const cellModel*>(0)
  54     );
  55
  56     fluentCellModelLookup[2] = cellModeller::lookup("tet");
  57     fluentCellModelLookup[4] = cellModeller::lookup("hex");
  58     fluentCellModelLookup[5] = cellModeller::lookup("pyr");
  59     fluentCellModelLookup[6] = cellModeller::lookup("prism");
  60
  61     static label faceMatchingOrder[7][6] =
  62     {
  63         {-1, -1, -1, -1, -1, -1},
  64         {-1, -1, -1, -1, -1, -1},
  65         { 0,  1,  2,  3, -1, -1},    // tet
  66         {-1, -1, -1, -1, -1, -1},
  67         { 0,  2,  4,  3,  5,  1},    // hex
  68         { 0,  1,  2,  3,  4, -1},    // pyr
  69         { 0,  2,  3,  4,  1, -1},    // prism
  70     };
  71
  72     const cellModel& curModel = *fluentCellModelLookup[fluentCellModelID];
  73
  74     // Checking
  75     if (faceLabels.size() != curModel.nFaces())
  76     {
  77         FatalErrorIn
  78         (
  79             "create3DCellShape(const label cellIndex, "
  80             "const labelList& faceLabels, const labelListList& faces, "
  81             "const labelList& owner, const labelList& neighbour, "
  82             "const label fluentCellModelID)"
  83         )   << "Number of face labels not equal to"
  84             << "number of face in the model. "
  85             << "Number of face labels: " << faceLabels.size()
  86             << " number of faces in model: " << curModel.nFaces()
  87             << abort(FatalError);
  88     }
  89
  90     // make a list of outward-pointing faces
  91     labelListList localFaces(faceLabels.size());
  92
  93     forAll(faceLabels, faceI)
  94     {
  95         const label curFaceLabel = faceLabels[faceI];
  96
  97         const labelList& curFace = faces[curFaceLabel];
  98
  99         if (owner[curFaceLabel] == cellIndex)
 100         {
 101             localFaces[faceI] = curFace;
 102         }
 103         else if (neighbour[curFaceLabel] == cellIndex)
 104         {
 105             // Reverse the face
 106             localFaces[faceI].setSize(curFace.size());
 107
 108             forAllReverse(curFace, i)
 109             {
 110                 localFaces[faceI][curFace.size() - i - 1] =
 111                     curFace[i];
 112             }
 113         }
 114         else
 115         {
 116             FatalErrorIn
 117             (
 118                 "create3DCellShape(const label cellIndex, "
 119                 "const labelList& faceLabels, const labelListList& faces, "
 120                 "const labelList& owner, const labelList& neighbour, "
 121                 "const label fluentCellModelID)"
 122             )   << "face " << curFaceLabel
 123                 << " does not belong to cell " << cellIndex
 124                 << ". Face owner: " << owner[curFaceLabel] << " neighbour: "
 125                 << neighbour[curFaceLabel]
 126                 << abort(FatalError);
 127         }
 128     }
 129
 130     // Algorithm:
 131     // Make an empty list of pointLabels and initialise it with -1. Pick the
 132     // first face from modelFaces and look through the faces to find one with
 133     // the same number of labels. Insert face by copying its labels into
 134     // pointLabels. Mark the face as used. Loop through all model faces.
 135     // For each model face loop through faces. If the face is unused and the
 136     // numbers of labels fit, try to match the face onto the point labels. If
 137     // at least one edge is matched, insert the face into pointLabels. If at
 138     // any stage the matching algorithm reaches the end of faces, the matching
 139     // algorithm has failed. Once all the faces are matched, the list of
 140     // pointLabels defines the model.
 141
 142     // Make a list of empty pointLabels
 143     labelList pointLabels(curModel.nPoints(), -1);
 144
 145     // Follow the used mesh faces
 146     List<bool> meshFaceUsed(localFaces.size(), false);
 147
 148     // Get the raw model faces
 149     const faceList& modelFaces = curModel.modelFaces();
 150
 151     // Insert the first face into the list
 152     const labelList& firstModelFace =
 153         modelFaces[faceMatchingOrder[fluentCellModelID][0]];
 154
 155     bool found = false;
 156
 157     forAll(localFaces, meshFaceI)
 158     {
 159         if (localFaces[meshFaceI].size() == firstModelFace.size())
 160         {
 161             // Match. Insert points into the pointLabels
 162             found = true;
 163
 164             const labelList& curMeshFace = localFaces[meshFaceI];
 165
 166             meshFaceUsed[meshFaceI] = true;
 167
 168             forAll(curMeshFace, pointI)
 169             {
 170                 pointLabels[firstModelFace[pointI]] = curMeshFace[pointI];
 171             }
 172
 173             break;
 174         }
 175     }
 176
 177     if (!found)
 178     {
 179         FatalErrorIn
 180         (
 181             "create3DCellShape(const label cellIndex, "
 182             "const labelList& faceLabels, const labelListList& faces, "
 183             "const labelList& owner, const labelList& neighbour, "
 184             "const label fluentCellModelID)"
 185         )   << "Cannot find match for first face. "
 186             << "cell model: " << curModel.name() << " first model face: "
 187             << firstModelFace << " Mesh faces: " << localFaces
 188             << abort(FatalError);
 189     }
 190
 191     for (label modelFaceI = 1; modelFaceI < modelFaces.size(); modelFaceI++)
 192     {
 193         // get the next model face
 194         const labelList& curModelFace =
 195             modelFaces
 196             [faceMatchingOrder[fluentCellModelID][modelFaceI]];
 197
 198         found = false;
 199
 200         // Loop through mesh faces until a match is found
 201         forAll(localFaces, meshFaceI)
 202         {
 203             if
 204             (
 205                 !meshFaceUsed[meshFaceI]
 206              && localFaces[meshFaceI].size() == curModelFace.size()
 207             )
 208             {
 209                 // A possible match. A mesh face will be rotated, so make a copy
 210                 labelList meshFaceLabels = localFaces[meshFaceI];
 211
 212                 for
 213                 (
 214                     label rotation = 0;
 215                     rotation < meshFaceLabels.size();
 216                     rotation++
 217                 )
 218                 {
 219                     // try matching the face
 220                     label nMatchedLabels = 0;
 221
 222                     forAll(meshFaceLabels, pointI)
 223                     {
 224                         if
 225                         (
 226                             pointLabels[curModelFace[pointI]]
 227                          == meshFaceLabels[pointI]
 228                         )
 229                         {
 230                             nMatchedLabels++;
 231                         }
 232                     }
 233
 234                     if (nMatchedLabels >= 2)
 235                     {
 236                         // match!
 237                         found = true;
 238                     }
 239
 240                     if (found)
 241                     {
 242                         // match found. Insert mesh face
 243                         forAll(meshFaceLabels, pointI)
 244                         {
 245                             pointLabels[curModelFace[pointI]] =
 246                                 meshFaceLabels[pointI];
 247                         }
 248
 249                         meshFaceUsed[meshFaceI] = true;
 250
 251                         break;
 252                     }
 253                     else
 254                     {
 255                         // No match found. Rotate face
 256                         label firstLabel = meshFaceLabels[0];
 257
 258                         for (label i = 1; i < meshFaceLabels.size(); i++)
 259                         {
 260                             meshFaceLabels[i - 1] = meshFaceLabels[i];
 261                         }
 262
 263                         meshFaceLabels.last() = firstLabel;
 264                     }
 265                 }
 266
 267                 if (found) break;
 268             }
 269         }
 270
 271         if (!found)
 272         {
 273             // A model face is not matched. Shape detection failed
 274             FatalErrorIn
 275             (
 276                 "create3DCellShape(const label cellIndex, "
 277                 "const labelList& faceLabels, const labelListList& faces, "
 278                 "const labelList& owner, const labelList& neighbour, "
 279                 "const label fluentCellModelID)"
 280             )   << "Cannot find match for face "
 281                 << modelFaceI
 282                 << ".\nModel: " << curModel.name() << " model face: "
 283                 << curModelFace << " Mesh faces: " << localFaces
 284                 << "Matched points: " << pointLabels
 285                 << abort(FatalError);
 286         }
 287     }
 288
 289     return cellShape(curModel, pointLabels);
 290 }
 291
 292
 293 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 294
 295 } // End namespace Foam
 296
 297 // ************************************************************************* //