src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchCheck.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 Description
  25     Checks topology of the patch.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "PrimitivePatchTemplate.H"
  30 #include "Map.H"
  31 #include "ListOps.H"
  32 #include "OFstream.H"
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 template
  37 <
  38     class Face,
  39     template<class> class FaceList,
  40     class PointField,
  41     class PointType
  42 >
  43 void
  44 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
  45 visitPointRegion
  46 (
  47     const label pointI,
  48     const labelList& pFaces,
  49     const label startFaceI,
  50     const label startEdgeI,
  51     boolList& pFacesHad
  52 ) const
  53 {
  54     label index = findIndex(pFaces, startFaceI);
  55
  56     if (!pFacesHad[index])
  57     {
  58         // Mark face as been visited.
  59         pFacesHad[index] = true;
  60
  61         // Step to next edge on face which is still using pointI
  62         const labelList& fEdges = faceEdges()[startFaceI];
  63
  64         label nextEdgeI = -1;
  65
  66         forAll(fEdges, i)
  67         {
  68             label edgeI = fEdges[i];
  69
  70             const edge& e = edges()[edgeI];
  71
  72             if (edgeI != startEdgeI && (e[0] == pointI || e[1] == pointI))
  73             {
  74                 nextEdgeI = edgeI;
  75
  76                 break;
  77             }
  78         }
  79
  80         if (nextEdgeI == -1)
  81         {
  82             FatalErrorIn
  83             (
  84                 "PrimitivePatch<Face, FaceList, PointField, PointType>::"
  85                 "visitPointRegion"
  86             )   << "Problem: cannot find edge out of " << fEdges
  87                 << "on face " << startFaceI << " that uses point " << pointI
  88                 << " and is not edge " << startEdgeI << abort(FatalError);
  89         }
  90
  91         // Walk to next face(s) across edge.
  92         const labelList& eFaces = edgeFaces()[nextEdgeI];
  93
  94         forAll(eFaces, i)
  95         {
  96             if (eFaces[i] != startFaceI)
  97             {
  98                 visitPointRegion
  99                 (
 100                     pointI,
 101                     pFaces,
 102                     eFaces[i],
 103                     nextEdgeI,
 104                     pFacesHad
 105                 );
 106             }
 107         }
 108     }
 109 }
 110
 111
 112 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 113
 114 template
 115 <
 116     class Face,
 117     template<class> class FaceList,
 118     class PointField,
 119     class PointType
 120 >
 121 typename Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::surfaceTopo
 122 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 123 surfaceType() const
 124 {
 125     if (debug)
 126     {
 127         Info<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
 128                "surfaceType() : "
 129                "calculating patch topology"
 130             << endl;
 131     }
 132
 133     const labelListList& edgeFcs = edgeFaces();
 134
 135     surfaceTopo pType = MANIFOLD;
 136
 137     forAll(edgeFcs, edgeI)
 138     {
 139         label nNbrs = edgeFcs[edgeI].size();
 140
 141         if (nNbrs < 1 || nNbrs > 2)
 142         {
 143             pType = ILLEGAL;
 144
 145             // Can exit now. Surface is illegal.
 146             return pType;
 147         }
 148         else if (nNbrs == 1)
 149         {
 150             // Surface might be open or illegal so keep looping.
 151             pType = OPEN;
 152         }
 153     }
 154
 155     if (debug)
 156     {
 157         Info<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
 158                "surfaceType() : "
 159                "finished calculating patch topology"
 160             << endl;
 161     }
 162
 163     return pType;
 164 }
 165
 166
 167 template
 168 <
 169     class Face,
 170     template<class> class FaceList,
 171     class PointField,
 172     class PointType
 173 >
 174 bool
 175 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 176 checkTopology
 177 (
 178     const bool report,
 179     labelHashSet* setPtr
 180 ) const
 181 {
 182     if (debug)
 183     {
 184         Info<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
 185                "checkTopology(const bool, labelHashSet&) : "
 186                "checking patch topology"
 187             << endl;
 188     }
 189
 190     // Check edgeFaces
 191
 192     const labelListList& edgeFcs = edgeFaces();
 193
 194     surfaceTopo surfaceType = MANIFOLD;
 195
 196     forAll(edgeFcs, edgeI)
 197     {
 198         label nNbrs = edgeFcs[edgeI].size();
 199
 200         if (nNbrs < 1 || nNbrs > 2)
 201         {
 202             surfaceType = ILLEGAL;
 203
 204             if (report)
 205             {
 206                 Info<< "Edge " << edgeI << " with vertices:" << edges()[edgeI]
 207                     << " has " << nNbrs << " face neighbours"
 208                     << endl;
 209             }
 210
 211             if (setPtr)
 212             {
 213                 const edge& e = edges()[edgeI];
 214
 215                 setPtr->insert(meshPoints()[e.start()]);
 216                 setPtr->insert(meshPoints()[e.end()]);
 217             }
 218         }
 219         else if (nNbrs == 1)
 220         {
 221             surfaceType = OPEN;
 222         }
 223     }
 224
 225     if (debug)
 226     {
 227         Info<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
 228                "checkTopology(const bool, labelHashSet&) : "
 229                "finished checking patch topology"
 230             << endl;
 231     }
 232
 233     return surfaceType == ILLEGAL;
 234 }
 235
 236
 237 template
 238 <
 239     class Face,
 240     template<class> class FaceList,
 241     class PointField,
 242     class PointType
 243 >
 244 bool
 245 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 246 checkPointManifold
 247 (
 248     const bool report,
 249     labelHashSet* setPtr
 250 ) const
 251 {
 252     const labelListList& pf = pointFaces();
 253     const labelListList& pe = pointEdges();
 254     const labelListList& ef = edgeFaces();
 255     const labelList& mp = meshPoints();
 256
 257     bool foundError = false;
 258
 259     forAll(pf, pointI)
 260     {
 261         const labelList& pFaces = pf[pointI];
 262
 263         // Visited faces (as indices into pFaces)
 264         boolList pFacesHad(pFaces.size(), false);
 265
 266         // Starting edge
 267         const labelList& pEdges = pe[pointI];
 268         label startEdgeI = pEdges[0];
 269
 270         const labelList& eFaces = ef[startEdgeI];
 271
 272         forAll(eFaces, i)
 273         {
 274             // Visit all faces using pointI, starting from eFaces[i] and
 275             // startEdgeI. Mark off all faces visited in pFacesHad.
 276             this->visitPointRegion
 277             (
 278                 pointI,
 279                 pFaces,
 280                 eFaces[i],  // starting face for walk
 281                 startEdgeI, // starting edge for walk
 282                 pFacesHad
 283             );
 284         }
 285
 286         // After this all faces using pointI should have been visited and
 287         // marked off in pFacesHad.
 288
 289         label unset = findIndex(pFacesHad, false);
 290
 291         if (unset != -1)
 292         {
 293             foundError = true;
 294
 295             label meshPointI = mp[pointI];
 296
 297             if (setPtr)
 298             {
 299                 setPtr->insert(meshPointI);
 300             }
 301
 302             if (report)
 303             {
 304                 Info<< "Point " << meshPointI
 305                     << " uses faces which are not connected through an edge"
 306                     << nl
 307                     << "This means that the surface formed by this patched"
 308                     << " is multiply connected at this point" << nl
 309                     << "Connected (patch) faces:" << nl;
 310
 311                 forAll(pFacesHad, i)
 312                 {
 313                     if (pFacesHad[i])
 314                     {
 315                         Info<< "    " << pFaces[i] << endl;
 316                     }
 317                 }
 318
 319                 Info<< nl << "Unconnected (patch) faces:" << nl;
 320                 forAll(pFacesHad, i)
 321                 {
 322                     if (!pFacesHad[i])
 323                     {
 324                         Info<< "    " << pFaces[i] << endl;
 325                     }
 326                 }
 327             }
 328         }
 329     }
 330
 331     return foundError;
 332 }
 333
 334
 335 template
 336 <
 337     class Face,
 338     template<class> class FaceList,
 339     class PointField,
 340     class PointType
 341 >
 342 void
 343 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 344 writeVTK
 345 (
 346     const fileName& name,
 347     const FaceListType& faces,
 348     const Field<PointType>& points
 349 )
 350 {
 351     // Write patch and points into VTK
 352     OFstream mps(name + ".vtk");
 353
 354     mps << "# vtk DataFile Version 2.0" << nl
 355         << name << ".vtk" << nl
 356         << "ASCII" << nl
 357         << "DATASET POLYDATA" << nl
 358         << "POINTS " << points.size() << " float" << nl;
 359
 360     // Write points
 361     List<float> mlpBuffer(3*points.size());
 362
 363     label counter = 0;
 364     forAll (points, i)
 365     {
 366         mlpBuffer[counter++] = float(points[i].x());
 367         mlpBuffer[counter++] = float(points[i].y());
 368         mlpBuffer[counter++] = float(points[i].z());
 369     }
 370
 371     forAll (mlpBuffer, i)
 372     {
 373         mps << mlpBuffer[i] << ' ';
 374
 375         if (i > 0 && (i % 10) == 0)
 376         {
 377             mps << nl;
 378         }
 379     }
 380
 381     // Write faces
 382     label nFaceVerts = 0;
 383
 384     forAll (faces, faceI)
 385     {
 386         nFaceVerts += faces[faceI].size() + 1;
 387     }
 388     labelList mlfBuffer(nFaceVerts);
 389
 390     counter = 0;
 391     forAll (faces, faceI)
 392     {
 393         const Face& f = faces[faceI];
 394
 395         mlfBuffer[counter++] = f.size();
 396
 397         forAll (f, fpI)
 398         {
 399             mlfBuffer[counter++] = f[fpI];
 400         }
 401     }
 402     mps << nl;
 403
 404     mps << "POLYGONS " << faces.size() << ' ' << nFaceVerts << endl;
 405
 406     forAll (mlfBuffer, i)
 407     {
 408         mps << mlfBuffer[i] << ' ';
 409
 410         if (i > 0 && (i % 10) == 0)
 411         {
 412             mps << nl;
 413         }
 414     }
 415     mps << nl;
 416 }
 417
 418
 419 template
 420 <
 421     class Face,
 422     template<class> class FaceList,
 423     class PointField,
 424     class PointType
 425 >
 426 void
 427 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 428 writeVTKNormals
 429 (
 430     const fileName& name,
 431     const FaceListType& faces,
 432     const Field<PointType>& points
 433 )
 434 {
 435     // Write patch and points into VTK
 436     OFstream mps(name + ".vtk");
 437
 438     mps << "# vtk DataFile Version 2.0" << nl
 439         << name << ".vtk" << nl
 440         << "ASCII" << nl
 441         << "DATASET POLYDATA" << nl
 442         << "POINTS " << faces.size() << " float" << nl;
 443
 444     // Write points
 445     List<float> mlPointBuffer(3*faces.size());
 446
 447     label counter = 0;
 448     forAll (faces, i)
 449     {
 450         const vector c = faces[i].centre(points);
 451
 452         mlPointBuffer[counter++] = float(c.x());
 453         mlPointBuffer[counter++] = float(c.y());
 454         mlPointBuffer[counter++] = float(c.z());
 455     }
 456
 457     forAll (mlPointBuffer, i)
 458     {
 459         mps << mlPointBuffer[i] << ' ';
 460
 461         if (i > 0 && (i % 10) == 0)
 462         {
 463             mps << nl;
 464         }
 465     }
 466     mps << nl;
 467
 468     // Write normals
 469     mps << "POINT_DATA " << faces.size() << nl
 470         << "FIELD attributes " << 1 << nl
 471         << "normals" << " 3 "
 472         << faces.size() << " float" << nl;
 473
 474     List<float> mlNormalBuffer(3*faces.size());
 475
 476     counter = 0;
 477     forAll (faces, i)
 478     {
 479         const vector n = faces[i].normal(points);
 480
 481         mlNormalBuffer[counter++] = float(n.x());
 482         mlNormalBuffer[counter++] = float(n.y());
 483         mlNormalBuffer[counter++] = float(n.z());
 484     }
 485
 486     forAll (mlNormalBuffer, i)
 487     {
 488         mps << mlNormalBuffer[i] << ' ';
 489
 490         if (i > 0 && (i % 10) == 0)
 491         {
 492             mps << nl;
 493         }
 494     }
 495     mps << nl;
 496 }
 497
 498
 499 template
 500 <
 501     class Face,
 502     template<class> class FaceList,
 503     class PointField,
 504     class PointType
 505 >
 506 void
 507 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 508 writeVTK
 509 (
 510     const fileName& name
 511 ) const
 512 {
 513     PrimitivePatch<Face, FaceList, PointField, PointType>::writeVTK
 514     (
 515         name,
 516         this->localFaces(),
 517         this->localPoints()
 518     );
 519 }
 520
 521
 522 template
 523 <
 524     class Face,
 525     template<class> class FaceList,
 526     class PointField,
 527     class PointType
 528 >
 529 void
 530 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 531 writeVTKNormals
 532 (
 533     const fileName& name
 534 ) const
 535 {
 536     PrimitivePatch<Face, FaceList, PointField, PointType>::writeVTKNormals
 537     (
 538         name,
 539         this->localFaces(),
 540         this->localPoints()
 541     );
 542 }
 543
 544
 545 // ************************************************************************* //