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