src/foam/meshes/polyMesh/polyMeshInitMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMesh.H"
  27
  28 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  29
  30 void Foam::polyMesh::initMesh()
  31 {
  32     if (debug)
  33     {
  34         Info<< "void polyMesh::initMesh() : "
  35             << "initialising primitiveMesh" << endl;
  36     }
  37
  38     // For backward compatibility check if the owner array is the same
  39     // length as the number of allFaces and shrink to remove the -1s padding
  40     if (min(owner_) < 0)
  41     {
  42         label nActiveFaces = 0;
  43
  44         forAll(owner_, faceI)
  45         {
  46             if (owner_[faceI] == -1)
  47             {
  48                 break;
  49             }
  50             else
  51             {
  52                 nActiveFaces++;
  53             }
  54         }
  55
  56         InfoIn("void polyMesh::initMesh()")
  57             << "Truncating owner list at " << nActiveFaces
  58             << " for backward compatibility" << endl;
  59
  60         owner_.setSize(nActiveFaces);
  61     }
  62
  63     // For backward compatibility check if the neighbour array is the same
  64     // length as the owner and shrink to remove the -1s padding
  65     if (min(neighbour_) < 0)
  66     {
  67         label nIntFaces = 0;
  68
  69         forAll(neighbour_, faceI)
  70         {
  71             if (neighbour_[faceI] == -1)
  72             {
  73                 break;
  74             }
  75             else
  76             {
  77                 nIntFaces++;
  78             }
  79         }
  80
  81         InfoIn("void polyMesh::initMesh()")
  82             << "Truncating neighbour list at " << nIntFaces
  83             << " for backward compatibility" << endl;
  84
  85         neighbour_.setSize(nIntFaces);
  86     }
  87
  88     label nCells = -1;
  89
  90     forAll(owner_, facei)
  91     {
  92         nCells = max(nCells, owner_[facei]);
  93     }
  94
  95     // The neighbour array may or may not be the same length as the owner
  96     forAll(neighbour_, facei)
  97     {
  98         nCells = max(nCells, neighbour_[facei]);
  99     }
 100
 101     nCells++;
 102
 103     label nUsedFaces = 0;
 104
 105     // Use patch info if provided, use all faces otherwise
 106     if (boundary_.size())
 107     {
 108         nUsedFaces =
 109             boundary_[boundary_.size() - 1].start()
 110           + boundary_[boundary_.size() - 1].size();
 111     }
 112     else
 113     {
 114         // No patch info. Assume all faces are used.
 115         nUsedFaces = owner_.size();
 116     }
 117
 118
 119     label nUsedPoints = allPoints_.size();
 120
 121     // If not all faces are being used, check that all unused
 122     // faces are at the back of the list and check the number of
 123     // used points.
 124     if (nUsedFaces < allFaces_.size())
 125     {
 126         if (debug)
 127         {
 128             Info<< "void polyMesh::initMesh() : "
 129                 << "unused faces detected.  "
 130                 << "Number of used faces: " << nUsedFaces
 131                 << ".  Total number of faces: " << allFaces_.size() << endl;
 132         }
 133
 134         // Mark the points used by live faces.
 135         boolList usedPoints(allPoints_.size(), false);
 136
 137         for (label faceI = 0; faceI < nUsedFaces; faceI++)
 138         {
 139             const face& curFace = allFaces_[faceI];
 140
 141             forAll (curFace, pointI)
 142             {
 143                 usedPoints[curFace[pointI]] = true;
 144             }
 145         }
 146
 147         forAll (usedPoints, pointI)
 148         {
 149             if (!usedPoints[pointI])
 150             {
 151                 nUsedPoints = pointI;
 152                 break;
 153             }
 154         }
 155
 156         if (nUsedPoints < allPoints_.size())
 157         {
 158             if (debug)
 159             {
 160                 Info<< "void polyMesh::initMesh() : unused points "
 161                     << "detected.  Number of used points: "
 162                     << nUsedPoints << ". Total number of points: "
 163                     << allPoints_.size() << endl;
 164             }
 165
 166             for (label i = nUsedPoints; i < allPoints_.size(); i++)
 167             {
 168                 if (usedPoints[i])
 169                 {
 170                     FatalErrorIn("void polyMesh::initMesh()")
 171                         << "Error in point ordering: mixed used and unused "
 172                         << "points at the end of point list." << nl
 173                         << "Last used point: " << nUsedPoints
 174                         << " " << allPoints_[nUsedPoints] << nl
 175                         << "First unused point: " << nUsedPoints + 1
 176                         << " " << allPoints_[nUsedPoints + 1] << nl
 177                         << "and point " << i << " " << allPoints_[i]
 178                         << " is used by a live face." << endl;
 179
 180                     // Find out which face uses this point
 181                     // Problem point = i
 182                     for (label faceI = 0; faceI < nUsedFaces; faceI++)
 183                     {
 184                         const face& curFace = allFaces_[faceI];
 185
 186                         forAll (curFace, pointI)
 187                         {
 188                             if (curFace[pointI] == i)
 189                             {
 190                                 FatalError
 191                                     << "Face " << faceI << " " << curFace
 192                                     << " with points "
 193                                     << curFace.points(allPoints_)
 194                                     << endl;
 195                                 break;
 196                             }
 197                         }
 198                     }
 199
 200                     allPoints_.write();
 201                     allFaces_.write();
 202                     owner_.write();
 203                     neighbour_.write();
 204                     boundary_.write();
 205                     FatalError << "Done. " << abort(FatalError);
 206                 }
 207             }
 208         }
 209     }
 210
 211     // Set sliced lists
 212     points_.reset(allPoints_, nUsedPoints);
 213     faces_.reset(allFaces_, owner_.size());
 214
 215     // Reset the primitiveMesh with the sizes of the primitive arrays
 216     primitiveMesh::reset
 217     (
 218         nUsedPoints,
 219         neighbour_.size(),
 220         owner_.size(),
 221         nCells
 222     );
 223
 224     string meshInfo =
 225         "nPoints: " + Foam::name(nPoints())
 226       + " nAllPoints: " + Foam::name(allPoints().size())
 227       + " nInternalFaces: " + Foam::name(nInternalFaces())
 228       + " nFaces: " + Foam::name(nFaces())
 229       + " nAllFaces: " + Foam::name(allFaces().size())
 230       + " nCells: " + Foam::name(this->nCells());
 231
 232     owner_.note() = meshInfo;
 233     neighbour_.note() = meshInfo;
 234 }
 235
 236
 237 void Foam::polyMesh::initMesh(cellList& c)
 238 {
 239     if (debug)
 240     {
 241         Info<< "void polyMesh::initMesh(cellList& c) : "
 242             << "calculating owner-neighbour arrays" << endl;
 243     }
 244
 245     owner_.setSize(allFaces_.size(), -1);
 246     neighbour_.setSize(allFaces_.size(), -1);
 247
 248     boolList markedFaces(allFaces_.size(), false);
 249
 250     label nInternalFaces = 0;
 251
 252     forAll(c, cellI)
 253     {
 254         // get reference to face labels for current cell
 255         const labelList& cellfaces = c[cellI];
 256
 257         forAll (cellfaces, faceI)
 258         {
 259             if (!markedFaces[cellfaces[faceI]])
 260             {
 261                 // First visit: owner
 262                 owner_[cellfaces[faceI]] = cellI;
 263                 markedFaces[cellfaces[faceI]] = true;
 264             }
 265             else
 266             {
 267                 // Second visit: neighbour
 268                 neighbour_[cellfaces[faceI]] = cellI;
 269                 nInternalFaces++;
 270             }
 271         }
 272     }
 273
 274     // The neighbour array is initialised with the same length as the owner
 275     // padded with -1s and here it is truncated to the correct size of the
 276     // number of internal faces.
 277     neighbour_.setSize(nInternalFaces);
 278
 279     label nUsedPoints = allPoints_.size();
 280     label nUsedFaces = allFaces_.size();
 281
 282     forAll (owner_, faceI)
 283     {
 284         if (owner_[faceI] < 0)
 285         {
 286             nUsedFaces = faceI;
 287             break;
 288         }
 289     }
 290
 291     // If not all faces are being used, check that all unused
 292     // faces are at the back of the list and check the number of
 293     // used points.
 294     if (nUsedFaces < owner_.size())
 295     {
 296         if (debug)
 297         {
 298             Info<< "void polyMesh::initMesh(cellList& c) : "
 299                 << "unused faces detected.  "
 300                 << "Number of used faces: " << nUsedFaces
 301                 << ".  Total number of faces: " << owner_.size() << endl;
 302         }
 303
 304         for (label i = nUsedFaces; i < owner_.size(); i++)
 305         {
 306             if (owner_[i] >= 0)
 307             {
 308                 FatalErrorIn("void polyMesh::initMesh(cellList& c)")
 309                     << "Error in face ordering: mixed used and unused "
 310                     << "faces at the end of face list." << nl
 311                     << "Number of used faces: " << nUsedFaces
 312                     << "  and face " << i << " is owned by cell " << owner_[i]
 313                     << abort(FatalError);
 314             }
 315         }
 316
 317         // Reset the size of owner array to the number of live faces
 318         owner_.setSize(nUsedFaces);
 319
 320         // Mark the points used by live faces.
 321         boolList usedPoints(allPoints_.size(), false);
 322
 323         for (label faceI = 0; faceI < nUsedFaces; faceI++)
 324         {
 325             const face& curFace = allFaces_[faceI];
 326
 327             forAll (curFace, pointI)
 328             {
 329                 usedPoints[curFace[pointI]] = true;
 330             }
 331         }
 332
 333         forAll (usedPoints, pointI)
 334         {
 335             if (!usedPoints[pointI])
 336             {
 337                 nUsedPoints = pointI;
 338                 break;
 339             }
 340         }
 341
 342         if (nUsedPoints < allPoints_.size())
 343         {
 344             if (debug)
 345             {
 346                 Info<< "void polyMesh::initMesh(cellList& c) : unused points "
 347                     << "detected.  Number of used points: "
 348                     << nUsedPoints << ". Total number of points: "
 349                     << allPoints_.size() << endl;
 350             }
 351
 352             for (label i = nUsedPoints; i < allPoints_.size(); i++)
 353             {
 354                 if (usedPoints[i])
 355                 {
 356                     FatalErrorIn("void polyMesh::initMesh(cellList& c)")
 357                         << "Error in point ordering: mixed used and unused "
 358                         << "points at the end of point list." << nl
 359                         << "Number of used points: " << nUsedPoints
 360                         << "  and point " << i
 361                         << " is used by a live face." << endl;
 362
 363                     // Find out which face uses this point
 364                     // Problem point = i
 365                     for (label faceI = 0; faceI < nUsedFaces; faceI++)
 366                     {
 367                         const face& curFace = allFaces_[faceI];
 368
 369                         forAll (curFace, pointI)
 370                         {
 371                             if (curFace[pointI] == i)
 372                             {
 373                                 FatalError
 374                                     << "Face " << faceI << " " << curFace
 375                                     << endl;
 376                                 break;
 377                             }
 378                         }
 379                     }
 380
 381                     FatalError << "Done. " << abort(FatalError);
 382                 }
 383             }
 384         }
 385     }
 386
 387
 388     // Reset the primitiveMesh
 389     primitiveMesh::reset
 390     (
 391         nUsedPoints,
 392         neighbour_.size(),
 393         owner_.size(),
 394         c.size(),
 395         c
 396     );
 397
 398     string meshInfo =
 399         "nPoints: " + Foam::name(nPoints())
 400       + " nCells: " + Foam::name(nCells())
 401       + " nFaces: " + Foam::name(nFaces())
 402       + " nInternalFaces: " + Foam::name(this->nInternalFaces());
 403
 404     owner_.note() = meshInfo;
 405     neighbour_.note() = meshInfo;
 406 }
 407
 408
 409 // ************************************************************************* //