src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/tetPolyMeshLduAddressingCellDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Class
  26     tetPolyMeshLduAddressingCellDecomp
  27
  28 Description
  29
  30 Author
  31     Hrvoje Jasak.  All rights reserved.
  32
  33 \*----------------------------------------------------------------------------*/
  34
  35 #include "tetPolyMeshLduAddressingCellDecomp.H"
  36 #include "globalMeshData.H"
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 Foam::tetPolyMeshLduAddressingCellDecomp::tetPolyMeshLduAddressingCellDecomp
  41 (
  42     const tetPolyMeshCellDecomp& mesh
  43 )
  44 :
  45     lduAddressing(mesh.nPoints()),
  46     lowerAddr_(mesh.nEdges(), -1),
  47     upperAddr_(mesh.nEdges(), -1),
  48     patchAddr_(mesh.boundary().size()),
  49     patchSchedule_(mesh.globalData().patchSchedule())
  50 {
  51     // Create owner and neighbour addressing list.
  52     // At the same time fill in the owner start lookup list
  53
  54     const faceList& faces = mesh().faces();
  55
  56     const labelListList& pointFaces = mesh().pointFaces();
  57     const labelListList& pointCells = mesh().pointCells();
  58
  59     // Count the added owners and neighbours
  60     label nCreatedEdges = 0;
  61
  62     label pointInFace, prev, next, f0;
  63
  64     // Loop through all points
  65     forAll (pointFaces, pointI)
  66     {
  67         const labelList& curFaces = pointFaces[pointI];
  68
  69         // Create a list of labels to keep the neighbours that
  70         // have already been added.  Size is estimated
  71         labelHashSet addedNeighbours
  72         (
  73             2*curFaces.size()*primitiveMesh::pointsPerFace_
  74         );
  75
  76         forAll (curFaces, faceI)
  77         {
  78             const face& f = faces[curFaces[faceI]];
  79
  80             // Grab zeroth label of face
  81             f0 = f[0];
  82
  83             labelList neighbourPointsFromFace(f.size() - 1, -1);
  84
  85             if (f0 == pointI)
  86             {
  87                 // If the current point is the zero point of the face,
  88                 // it is connected to all other points
  89                 for (label nbrI = 1; nbrI < f.size(); nbrI++)
  90                 {
  91                     if (f[nbrI] > pointI)
  92                     {
  93                         addedNeighbours.insert(f[nbrI]);
  94                     }
  95                 }
  96             }
  97             else if (f[1] == pointI)
  98             {
  99                 // If the current point is the first point of the face,
 100                 // it is connected to all other points
 101                 // if it is the last point, it is connected to point zero
 102                 // and the penultimate point
 103                 if (f0 > pointI)
 104                 {
 105                     addedNeighbours.insert(f0);
 106                 }
 107
 108                 if (f[2] > pointI)
 109                 {
 110                     addedNeighbours.insert(f[2]);
 111                 }
 112             }
 113             else if (f[f.size() - 1] == pointI)
 114             {
 115                 // If it is the last point, it is connected to point zero
 116                 // and the penultimate point
 117                 if (f[f.size() - 2] > pointI)
 118                 {
 119                     addedNeighbours.insert(f[f.size() - 2]);
 120                 }
 121
 122                 if (f0 > pointI)
 123                 {
 124                     addedNeighbours.insert(f0);
 125                 }
 126             }
 127             else
 128             {
 129                 // Otherwise, it is connected to the previous and the next
 130                 // point and additionally to point zero
 131                 pointInFace = f.which(pointI);
 132                 prev = f.prevLabel(pointInFace);
 133                 next = f.nextLabel(pointInFace);
 134
 135                 if (prev > pointI)
 136                 {
 137                     addedNeighbours.insert(prev);
 138                 }
 139
 140                 if (next > pointI)
 141                 {
 142                     addedNeighbours.insert(next);
 143                 }
 144
 145                 if (f0 > pointI)
 146                 {
 147                     addedNeighbours.insert(f0);
 148                 }
 149             }
 150         }
 151
 152         // All neighbours for the current point found. Before adding
 153         // them to the list, it is necessary to sort them in the
 154         // increasing order of the neighbouring point.
 155
 156         // Make real list out of SLList to simplify the manipulation.
 157         labelList an(addedNeighbours.toc());
 158
 159         // Use a simple sort to sort the an list.
 160         sort(an);
 161
 162         // Adding the neighbours
 163         forAll (an, edgeI)
 164         {
 165             lowerAddr_[nCreatedEdges] = pointI;
 166             upperAddr_[nCreatedEdges] = an[edgeI];
 167             nCreatedEdges++;
 168         }
 169
 170         // Now add cell neighbours
 171         const labelList& curPointCells = pointCells[pointI];
 172
 173         forAll (curPointCells, cellI)
 174         {
 175             // Add as neighbour
 176             lowerAddr_[nCreatedEdges] = pointI;
 177
 178             upperAddr_[nCreatedEdges] =
 179                 mesh.cellOffset() + curPointCells[cellI];
 180
 181             nCreatedEdges++;
 182         }
 183     }
 184
 185     // Add dummy boundary addressing
 186     forAll (patchAddr_, patchI)
 187     {
 188         patchAddr_[patchI].setSize(0);
 189     }
 190
 191     if (nCreatedEdges != mesh.nEdges())
 192     {
 193         FatalErrorIn
 194         (
 195             "tetPolyMeshLduAddressingCellDecomp::"
 196             "tetPolyMeshLduAddressingCellDecomp\n"
 197             "(\n"
 198             "    const tetPolyMeshCellDecomp& mesh\n"
 199             ")"
 200         )   << "Problem with edge counting in lduAddressing: "
 201             << "the cell decomposition is multiply connected or otherwise "
 202             << "invalid.  Please Use face decomposition instead.  "
 203             << "nCreatedEdges: " << nCreatedEdges
 204             << " nEdges: " << mesh.nEdges()
 205             << abort(FatalError);
 206     }
 207 }
 208
 209
 210 // ************************************************************************* //