src/OpenFOAM/meshes/meshShapes/cellMatcher/prismMatcher.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "prismMatcher.H"
  27 #include "primitiveMesh.H"
  28 #include "ListOps.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 const Foam::label Foam::prismMatcher::vertPerCell = 6;
  33 const Foam::label Foam::prismMatcher::facePerCell = 5;
  34 const Foam::label Foam::prismMatcher::maxVertPerFace = 4;
  35
  36
  37 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  38
  39 Foam::prismMatcher::prismMatcher()
  40 :
  41     cellMatcher
  42     (
  43         vertPerCell,
  44         facePerCell,
  45         maxVertPerFace,
  46         "prism"
  47     )
  48 {}
  49
  50 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  51
  52 Foam::prismMatcher::~prismMatcher()
  53 {}
  54
  55 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  56
  57 bool Foam::prismMatcher::matchShape
  58 (
  59     const bool checkOnly,
  60     const faceList& faces,
  61     const labelList& owner,
  62     const label cellI,
  63     const labelList& myFaces
  64 )
  65 {
  66     if (!faceSizeMatch(faces, myFaces))
  67     {
  68         return false;
  69     }
  70
  71     // Calculate localFaces_ and mapping pointMap_, faceMap_
  72     label numVert = calcLocalFaces(faces, myFaces);
  73
  74     if (numVert != vertPerCell)
  75     {
  76         return false;
  77     }
  78
  79     // Set up 'edge' to face mapping.
  80     calcEdgeAddressing(numVert);
  81
  82     // Set up point on face to index-in-face mapping
  83     calcPointFaceIndex();
  84
  85     // Storage for maps -vertex to mesh and -face to mesh
  86     vertLabels_.setSize(vertPerCell);
  87     faceLabels_.setSize(facePerCell);
  88
  89     //
  90     // Try first triangular face.
  91     // Only need to try one orientation of this face since prism is
  92     // rotation symmetric
  93     //
  94
  95     label face0I = -1;
  96     forAll(faceSize_, faceI)
  97     {
  98         if (faceSize_[faceI] == 3)
  99         {
 100             face0I = faceI;
 101             break;
 102         }
 103     }
 104
 105     const face& face0 = localFaces_[face0I];
 106     label face0vert0 = 0;
 107
 108     //
 109     // Try to follow prespecified path on faces of cell,
 110     // starting at face0vert0
 111     //
 112
 113     vertLabels_[0] = pointMap_[face0[face0vert0]];
 114     faceLabels_[0] = faceMap_[face0I];
 115     //Info<< endl << "Prism vertex 0: vertex " <<  face0[face0vert0]
 116     //    << " at position " << face0vert0 << " in face " << face0
 117     //    << endl;
 118
 119     // Walk face 0 from vertex 0 to 1
 120     label face0vert1 =
 121         nextVert
 122         (
 123             face0vert0,
 124             faceSize_[face0I],
 125             !(owner[faceMap_[face0I]] == cellI)
 126         );
 127     vertLabels_[1] = pointMap_[face0[face0vert1]];
 128     //Info<< "Prism vertex 1: vertex " <<  face0[face0vert1]
 129     //    << " at position " << face0vert1 << " in face " << face0
 130     //    << endl;
 131
 132     // Jump edge from face0 to face4
 133     label face4I =
 134         otherFace
 135         (
 136             numVert,
 137             face0[face0vert0],
 138             face0[face0vert1],
 139             face0I
 140         );
 141     const face& face4 = localFaces_[face4I];
 142     //Info<< "Stepped to prism face 4 " << face4
 143     //    << " across edge " << face0[face0vert0] << " "
 144     //    << face0[face0vert1]
 145     //    << endl;
 146
 147     if (faceSize_[face4I] != 4)
 148     {
 149         //Info<< "Cannot be Prism Face 4 since size="
 150         //    << faceSize_[face4I] << endl;
 151         return false;
 152     }
 153     faceLabels_[4] = faceMap_[face4I];
 154
 155     label face4vert1 = pointFaceIndex_[face0[face0vert1]][face4I];
 156
 157     //Info<< "Prism vertex 1 also: vertex " <<  face4[face4vert1]
 158     //    << " at position " << face4vert1 << " in face " << face4
 159     //    << endl;
 160
 161     // Walk face 4 from vertex 1 to 4
 162     label face4vert4 =
 163         nextVert
 164         (
 165             face4vert1,
 166             faceSize_[face4I],
 167             (owner[faceMap_[face4I]] == cellI)
 168         );
 169     vertLabels_[4] = pointMap_[face4[face4vert4]];
 170     //Info<< "Prism vertex 4: vertex " <<  face4[face4vert4]
 171     //    << " at position " << face4vert4 << " in face " << face4
 172     //    << endl;
 173
 174     // Walk face 4 from vertex 1 to 3
 175     label face4vert3 =
 176         nextVert
 177         (
 178             face4vert4,
 179             faceSize_[face4I],
 180             (owner[faceMap_[face4I]] == cellI)
 181         );
 182     vertLabels_[3] = pointMap_[face4[face4vert3]];
 183     //Info<< "Prism vertex 3: vertex " <<  face4[face4vert3]
 184     //    << " at position " << face4vert3 << " in face " << face4
 185     //    << endl;
 186
 187     // Jump edge from face4 to face1
 188     label face1I =
 189         otherFace
 190         (
 191             numVert,
 192             face4[face4vert3],
 193             face4[face4vert4],
 194             face4I
 195         );
 196     //const face& face1 = localFaces_[face1I];
 197     //Info<< "Stepped to prism face 1 " << face1
 198     //    << " across edge " << face4[face4vert3] << " "
 199     //    << face4[face4vert4]
 200     //    << endl;
 201
 202     if (faceSize_[face1I] != 3)
 203     {
 204         //Info<< "Cannot be Prism Face 1 since size="
 205         //    << faceSize_[face1I] << endl;
 206         return false;
 207     }
 208
 209     // Is prism for sure now
 210     if (checkOnly)
 211     {
 212         return true;
 213     }
 214
 215     faceLabels_[1] = faceMap_[face1I];
 216
 217
 218     //
 219     // Walk to other faces and assign mapping.
 220     //
 221
 222
 223     // Walk face 0 from vertex 1 to 2
 224     label face0vert2 =
 225         nextVert
 226         (
 227             face0vert1,
 228             faceSize_[face0I],
 229             !(owner[faceMap_[face0I]] == cellI)
 230         );
 231     vertLabels_[2] = pointMap_[face0[face0vert2]];
 232     //Info<< "Prism vertex 2: vertex " <<  face0[face0vert2]
 233     //    << " at position " << face0vert2 << " in face " << face0
 234     //    << endl;
 235
 236     // Jump edge from face0 to face3
 237     label face3I =
 238         otherFace
 239         (
 240             numVert,
 241             face0[face0vert1],
 242             face0[face0vert2],
 243             face0I
 244         );
 245     faceLabels_[3] = faceMap_[face3I];
 246     const face& face3 = localFaces_[face3I];
 247     //Info<< "Stepped to prism face 3 " << face3
 248     //    << " across edge " << face0[face0vert1] << " "
 249     //    << face0[face0vert2]
 250     //    << endl;
 251
 252     label face3vert2 = pointFaceIndex_[face0[face0vert2]][face3I];
 253
 254     //Info<< "Prism vertex 2 also: vertex " <<  face3[face3vert2]
 255     //    << " at position " << face3vert2 << " in face " << face3
 256     //    << endl;
 257
 258     label face3vert5 =
 259         nextVert
 260         (
 261             face3vert2,
 262             faceSize_[face3I],
 263             (owner[faceMap_[face3I]] == cellI)
 264         );
 265     vertLabels_[5] = pointMap_[face3[face3vert5]];
 266     //Info<< "Prism vertex 5: vertex " <<  face3[face3vert5]
 267     //    << " at position " << face3vert5 << " in face " << face3
 268     //    << endl;
 269
 270     // Jump edge from face0 to face2
 271     label face2I =
 272         otherFace
 273         (
 274             numVert,
 275             face0[face0vert2],
 276             face0[face0vert0],
 277             face0I
 278         );
 279     faceLabels_[2] = faceMap_[face2I];
 280     //const face& face2 = localFaces_[face2I];
 281     //Info<< "Stepped to prism face 2 " << face2
 282     //    << " across edge " << face0[face0vert2] << " "
 283     //    << face0[face0vert0]
 284     //    << endl;
 285
 286     //label face2vert2 = pointFaceIndex_[face0[face0vert2]][face2I];
 287     //Info<< "Prism vertex 2 also: vertex " <<  face2[face2vert2]
 288     //    << " at position " << face2vert2 << " in face " << face2
 289     //    << endl;
 290
 291     return true;
 292 }
 293
 294
 295 Foam::label Foam::prismMatcher::faceHashValue() const
 296 {
 297     return 2*3 + 4*4;
 298 }
 299
 300
 301 bool Foam::prismMatcher::faceSizeMatch
 302 (
 303     const faceList& faces,
 304     const labelList& myFaces
 305 ) const
 306 {
 307     if (myFaces.size() != 5)
 308     {
 309         return false;
 310     }
 311
 312     label nTris = 0;
 313     label nQuads = 0;
 314
 315     forAll(myFaces, myFaceI)
 316     {
 317         label size = faces[myFaces[myFaceI]].size();
 318
 319         if (size == 3)
 320         {
 321             nTris++;
 322         }
 323         else if (size == 4)
 324         {
 325             nQuads++;
 326         }
 327         else
 328         {
 329             return false;
 330         }
 331     }
 332     if ((nTris == 2) && (nQuads == 3))
 333     {
 334         return true;
 335     }
 336     else
 337     {
 338         return false;
 339     }
 340 }
 341
 342
 343 bool Foam::prismMatcher::isA(const primitiveMesh& mesh, const label cellI)
 344 {
 345     return matchShape
 346     (
 347         true,
 348         mesh.faces(),
 349         mesh.faceOwner(),
 350         cellI,
 351         mesh.cells()[cellI]
 352     );
 353 }
 354
 355
 356 bool Foam::prismMatcher::isA(const faceList& faces)
 357 {
 358     // Do as if mesh with one cell only
 359     return matchShape
 360     (
 361         true,
 362         faces,                      // all faces in mesh
 363         labelList(faces.size(), 0), // cell 0 is owner of all faces
 364         0,                          // cell label
 365         identity(faces.size())      // faces of cell 0
 366     );
 367 }
 368
 369
 370 bool Foam::prismMatcher::matches
 371 (
 372     const primitiveMesh& mesh,
 373     const label cellI,
 374     cellShape& shape
 375 )
 376 {
 377     if
 378     (
 379         matchShape
 380         (
 381             false,
 382             mesh.faces(),
 383             mesh.faceOwner(),
 384             cellI,
 385             mesh.cells()[cellI]
 386         )
 387     )
 388     {
 389         shape = cellShape(model(), vertLabels());
 390
 391         return true;
 392     }
 393     else
 394     {
 395         return false;
 396     }
 397 }
 398
 399
 400 // ************************************************************************* //