src/mesh/cfMesh/utilities/decomposeCells/decomposeCellsPyramids.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*----------------------p-----------------------------------------------------*/
  27
  28 #include "decomposeCells.H"
  29 #include "helperFunctions.H"
  30 #include "triFace.H"
  31
  32 //#define DEBUGDecompose
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 void decomposeCells::findAddressingForCell
  40 (
  41     const label cellI,
  42     DynList<label, 32>& vrt,
  43     DynList<edge, 64>& edges,
  44     DynList<DynList<label, 8> >& faceEdges,
  45     DynList<DynList<label, 2>, 64>& edgeFaces
  46 ) const
  47 {
  48     const cell& c = mesh_.cells()[cellI];
  49
  50     vrt.clear();
  51     edges.clear();
  52     edgeFaces.clear();
  53     faceEdges.setSize(c.size());
  54
  55     const faceListPMG& faces = mesh_.faces();
  56     forAll(faceEdges, feI)
  57     {
  58         faceEdges[feI].setSize(faces[c[feI]].size());
  59         faceEdges[feI] = -1;
  60     }
  61
  62     forAll(c, fI)
  63     {
  64         const face& f = faces[c[fI]];
  65
  66         forAll(f, eI)
  67         {
  68             const edge e = f.faceEdge(eI);
  69
  70             bool store(true);
  71             forAll(vrt, vI)
  72                 if( vrt[vI] == f[eI] )
  73                 {
  74                     store = false;
  75                     break;
  76                 }
  77             if( store )
  78             {
  79                 vrt.append(f[eI]);
  80             }
  81
  82             //- check if the edge alreready exists
  83             store = true;
  84
  85             forAll(edges, eJ)
  86                 if( e == edges[eJ] )
  87                 {
  88                     store = false;
  89                     faceEdges[fI][eI] = eJ;
  90                     edgeFaces[eJ].append(fI);
  91                     break;
  92                 }
  93
  94             if( store )
  95             {
  96                 faceEdges[fI][eI] = edges.size();
  97                 DynList<label, 2> ef;
  98                 ef.append(fI);
  99                 edgeFaces.append(ef);
 100                 edges.append(e);
 101             }
 102         }
 103     }
 104
 105     //- check if the cell is topologically closed
 106     forAll(edgeFaces, efI)
 107         if( edgeFaces[efI].size() != 2 )
 108         {
 109             forAll(c, fI)
 110                 Info << "Face " << c[fI] << " is " << faces[c[fI]] << endl;
 111
 112             Info << "Edges " << edges << endl;
 113             Info << "faceEdges " << faceEdges << endl;
 114             Info << "edgeFaces " << edgeFaces << endl;
 115             mesh_.write();
 116             FatalErrorIn
 117             (
 118                 "void decomposeCells::findAddressingForCell"
 119                 "(const label, DynList<label, 32>&, DynList<edge, 64>&"
 120                 ", DynList<DynList<label, 8> >&"
 121                 ", DynList<DynList<label, 2>, 64>&) const"
 122             ) << "Cell " << cellI << " is not topologically closed!"
 123                 << abort(FatalError);
 124         }
 125 }
 126
 127 label decomposeCells::findTopVertex
 128 (
 129     const label cellI,
 130     const DynList<label, 32>& vrt,
 131     const DynList<edge, 64>& edges,
 132     const DynList<DynList<label, 2>, 64>& edgeFaces
 133 )
 134 {
 135     const cell& c = mesh_.cells()[cellI];
 136     const faceListPMG& faces = mesh_.faces();
 137
 138     pointFieldPMG& pointsAccess = mesh_.points();
 139
 140     //- there is no vertex in 3 or more patches
 141     //- find boundary faces
 142     label topVertex(-1);
 143
 144     const labelList cp = c.labels(faces);
 145     point p(vector::zero);
 146     forAll(cp, cpI)
 147         p += pointsAccess[cp[cpI]];
 148     p /= cp.size();
 149
 150     topVertex = pointsAccess.size();
 151     pointsAccess.append(p);
 152
 153     # ifdef DEBUGDecompose
 154     Info << "Top vertex is " << topVertex << endl;
 155     # endif
 156
 157     return topVertex;
 158 }
 159
 160 void decomposeCells::decomposeCellIntoPyramids(const label cellI)
 161 {
 162     const cellListPMG& cells = mesh_.cells();
 163     const faceListPMG& faces = mesh_.faces();
 164     const labelList& owner = mesh_.owner();
 165
 166     const cell& c = cells[cellI];
 167
 168     # ifdef DEBUGDecompose
 169     Info << "Starting decomposing cell " << cellI << endl;
 170     Info << "Cell consists of faces " << c << endl;
 171     forAll(c, fI)
 172         Info << "Face " << c[fI] << " is " << faces[c[fI]] << endl;
 173     # endif
 174
 175     //- calculate edges, faceEdges and edgeFaces addressing
 176     DynList<label, 32> vrt;
 177     DynList<edge, 64> edges;
 178     DynList<DynList<label, 8> > faceEdges;
 179     faceEdges.setSize(c.size());
 180     DynList<DynList<label, 2>, 64> edgeFaces;
 181     findAddressingForCell(cellI, vrt, edges, faceEdges, edgeFaces);
 182
 183     // find a vertex which will be the top of the pyramids
 184     //- if there exist a corner vertex which is in 3 or more patches then
 185     //- it is selected as the top vertex
 186     label topVertex = findTopVertex(cellI, vrt, edges, edgeFaces);
 187
 188     //- start generating pyramids
 189     forAll(c, fI)
 190     {
 191         # ifdef DEBUGDecompose
 192         Info << "Face " << faces[c[fI]] << " is a base face" << endl;
 193         #endif
 194         const face& f = faces[c[fI]];
 195         DynList<DynList<label, 8> > cellFaces;
 196         cellFaces.setSize(f.size() + 1);
 197
 198         DynList<triFace> triFaces;
 199         triFaces.setSize(f.size());
 200         forAll(triFaces, pI)
 201         {
 202             triFaces[pI][0] = f.nextLabel(pI);
 203             triFaces[pI][1] = f[pI];
 204             triFaces[pI][2] = topVertex;
 205         }
 206
 207         label cfI(0);
 208         if( owner[c[fI]] == cellI )
 209         {
 210             cellFaces[cfI++] = faces[c[fI]];
 211
 212             forAll(triFaces, tfI)
 213             {
 214                 cellFaces[cfI++] = triFaces[tfI];
 215             }
 216         }
 217         else
 218         {
 219             cellFaces[cfI++] = faces[c[fI]].reverseFace();
 220
 221             forAll(triFaces, tfI)
 222             {
 223                 triFace rf;
 224                 rf[0] = triFaces[tfI][0];
 225                 rf[1] = triFaces[tfI][2];
 226                 rf[2] = triFaces[tfI][1];
 227                 cellFaces[cfI++] = rf;
 228             }
 229         }
 230
 231         # ifdef DEBUGDecompose
 232         Info << "Cell for face is " << cellFaces << endl;
 233
 234         DynList<edge, 64> cEdges;
 235         DynList<DynList<label, 2>, 64> eFaces;
 236         forAll(cellFaces, fI)
 237         {
 238             const DynList<label, 8>& f = cellFaces[fI];
 239             forAll(f, eI)
 240             {
 241                 const edge e(f[eI], f[(eI+1)%f.size()]);
 242
 243                 const label pos = cEdges.contains(e);
 244
 245                 if( pos < 0 )
 246                 {
 247                     cEdges.append(e);
 248                     DynList<label, 2> ef;
 249                     ef.append(fI);
 250                     eFaces.append(ef);
 251                 }
 252                 else
 253                 {
 254                     eFaces[pos].append(fI);
 255                 }
 256             }
 257         }
 258
 259         forAll(eFaces, eI)
 260             if( eFaces[eI].size() != 2 )
 261                 Pout << "This pyrmid is not topologically closed" << endl;
 262         # endif
 263
 264         facesOfNewCells_.appendGraph(cellFaces);
 265     }
 266 }
 267
 268 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 269
 270 } // End namespace Foam
 271
 272 // ************************************************************************* //