src/mesh/cfMesh/utilities/meshes/triSurf/triSurf.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 \*---------------------------------------------------------------------------*/
  27
  28 #include "triSurf.H"
  29 #include "demandDrivenData.H"
  30 #include "IFstream.H"
  31 #include "OFstream.H"
  32 #include "gzstream.h"
  33 #include "triSurface.H"
  34
  35 #include "helperFunctions.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 void triSurf::readFromFTR(const fileName& fName)
  45 {
  46     IFstream fStream(fName);
  47
  48     fStream >> triSurfFacets::patches_;
  49
  50     fStream >> triSurfPoints::points_;
  51
  52     fStream >> triSurfFacets::triangles_;
  53 }
  54
  55 void triSurf::writeToFTR(const fileName& fName) const
  56 {
  57     OFstream fStream(fName);
  58
  59     fStream << triSurfFacets::patches_;
  60
  61     fStream << nl;
  62
  63     fStream << triSurfPoints::points_;
  64
  65     fStream << nl;
  66
  67     fStream << triSurfFacets::triangles_;
  68 }
  69
  70 void triSurf::readFromFMS(const fileName& fName)
  71 {
  72     IFstream fStream(fName);
  73
  74     //- read the list of patches defined on the surface mesh
  75     fStream >> triSurfFacets::patches_;
  76
  77     //- read points
  78     fStream >> triSurfPoints::points_;
  79
  80     //- read surface triangles
  81     fStream >> triSurfFacets::triangles_;
  82
  83     //- read feature edges
  84     fStream >> triSurfFeatureEdges::featureEdges_;
  85
  86     List<meshSubset> subsets;
  87
  88     //- read point subsets
  89     fStream >> subsets;
  90     forAll(subsets, subsetI)
  91         triSurfPoints::pointSubsets_.insert(subsetI, subsets[subsetI]);
  92
  93     subsets.clear();
  94
  95     //- read facet subsets
  96     fStream >> subsets;
  97     forAll(subsets, subsetI)
  98         triSurfFacets::facetSubsets_.insert(subsetI, subsets[subsetI]);
  99
 100     subsets.clear();
 101
 102     //- read subsets on feature edges
 103     fStream >> subsets;
 104     forAll(subsets, subsetI)
 105         triSurfFeatureEdges::featureEdgeSubsets_.insert
 106         (
 107             subsetI,
 108             subsets[subsetI]
 109         );
 110 }
 111
 112 void triSurf::writeToFMS(const fileName& fName) const
 113 {
 114     OFstream fStream(fName);
 115
 116     //- write patches
 117     fStream << triSurfFacets::patches_;
 118
 119     fStream << nl;
 120
 121     //- write points
 122     fStream << triSurfPoints::points_;
 123
 124     fStream << nl;
 125
 126     //- write triangles
 127     fStream << triSurfFacets::triangles_;
 128
 129     fStream << nl;
 130
 131     //- write feature edges
 132     fStream << triSurfFeatureEdges::featureEdges_;
 133
 134     fStream << nl;
 135
 136     //- write point subsets
 137     List<meshSubset> subsets;
 138     label i(0);
 139     subsets.setSize(pointSubsets_.size());
 140     forAllConstIter(Map<meshSubset>, pointSubsets_, it)
 141         subsets[i++] = it();
 142     fStream << subsets;
 143
 144     fStream << nl;
 145
 146     //- write subsets of facets
 147     subsets.setSize(triSurfFacets::facetSubsets_.size());
 148     i = 0;
 149     forAllConstIter(Map<meshSubset>, triSurfFacets::facetSubsets_, it)
 150         subsets[i++] = it();
 151     fStream << subsets;
 152
 153     fStream << nl;
 154
 155     //- write subets of feature edges
 156     subsets.setSize(triSurfFeatureEdges::featureEdgeSubsets_.size());
 157     i = 0;
 158     forAllConstIter
 159     (
 160         Map<meshSubset>,
 161         triSurfFeatureEdges::featureEdgeSubsets_,
 162         it
 163     )
 164         subsets[i++] = it();
 165     fStream << subsets;
 166 }
 167
 168 void triSurf::topologyCheck()
 169 {
 170     const pointField& pts = this->points();
 171     const LongList<labelledTri>& trias = this->facets();
 172
 173     //- check for inf and nan points
 174     # ifdef USE_OMP
 175     # pragma omp parallel for schedule(dynamic, 100)
 176     # endif
 177     forAll(pts, pointI)
 178     {
 179         const point& p = pts[pointI];
 180
 181         if( help::isnan(p) || help::isinf(p) )
 182         {
 183             # ifdef USE_OMP
 184             # pragma omp critical
 185             # endif
 186             {
 187                 FatalErrorIn
 188                 (
 189                     "void triSurf::topologyCheck()"
 190                 ) << "Point " << pointI << " has invalid coordinates "
 191                   << p << exit(FatalError);
 192             }
 193         }
 194     }
 195
 196     //- check whether the nodes are within the scope
 197     //- report duplicate nodes in the same triangle
 198     # ifdef USE_OMP
 199     # pragma omp parallel for schedule(dynamic, 100)
 200     # endif
 201     forAll(trias, triI)
 202     {
 203         const labelledTri& ltri = trias[triI];
 204
 205         forAll(ltri, pI)
 206         {
 207             if( ltri[pI] < 0 || (ltri[pI] >= pts.size()) )
 208             {
 209                 # ifdef USE_OMP
 210                 # pragma omp critical
 211                 # endif
 212                 FatalErrorIn
 213                 (
 214                     "void triSurf::topologyCheck()"
 215                 ) << "Point " << ltri[pI] << " in triangle " << ltri
 216                   << " is out of scope 0 " << pts.size() << exit(FatalError);
 217             }
 218
 219             if( ltri[pI] == ltri[(pI+1)%3] || ltri[pI] == ltri[(pI+2)%3] )
 220             {
 221                 # ifdef USE_OMP
 222                 # pragma omp critical
 223                 # endif
 224                 WarningIn
 225                 (
 226                     "void triSurf::topologyCheck()"
 227                 ) << "Triangle " << ltri << " has duplicated points. "
 228                   << "This may cause problems in the meshing process!" << endl;
 229             }
 230         }
 231     }
 232
 233     //- check feature edges
 234     const edgeLongList& featureEdges = this->featureEdges();
 235
 236     # ifdef USE_OMP
 237     # pragma omp parallel for schedule(dynamic, 100)
 238     # endif
 239     forAll(featureEdges, eI)
 240     {
 241         const edge& fe = featureEdges[eI];
 242
 243         forAll(fe, pI)
 244         {
 245             if( fe[pI] < 0 || (fe[pI] >= pts.size()) )
 246             {
 247                 # ifdef USE_OMP
 248                 # pragma omp critical
 249                 # endif
 250                 FatalErrorIn
 251                 (
 252                     "void triSurf::topologyCheck()"
 253                 ) << "Feature edge " << fe << " point " << fe[pI]
 254                   << " is out of scope 0 " << pts.size() << exit(FatalError);
 255             }
 256         }
 257
 258         if( fe.start() == fe.end() )
 259         {
 260             # ifdef USE_OMP
 261             # pragma omp critical
 262             # endif
 263             WarningIn
 264             (
 265                 "void triSurf::topologyCheck()"
 266             ) << "Feature edge " << fe << " has duplicated points. "
 267               << "This may cause problems in the meshing process!" << endl;
 268         }
 269     }
 270
 271     //- check point subsets
 272     DynList<label> subsetIds;
 273     this->pointSubsetIndices(subsetIds);
 274     forAll(subsetIds, i)
 275     {
 276         labelLongList elmts;
 277         this->pointsInSubset(subsetIds[i], elmts);
 278
 279         forAll(elmts, elmtI)
 280         {
 281             const label elI = elmts[elmtI];
 282
 283             if( elI < 0 || elI >= pts.size() )
 284             {
 285                 # ifdef USE_OMP
 286                 # pragma omp critical
 287                 # endif
 288                 FatalErrorIn
 289                 (
 290                     "void triSurf::topologyCheck()"
 291                 ) << "Point " << elI << " in point subset "
 292                   << this->pointSubsetName(subsetIds[i])
 293                   << " is out of scope 0 " << pts.size() << exit(FatalError);
 294             }
 295         }
 296     }
 297
 298     //- check face subsets
 299     subsetIds.clear();
 300     this->facetSubsetIndices(subsetIds);
 301     forAll(subsetIds, i)
 302     {
 303         labelLongList elmts;
 304         this->facetsInSubset(subsetIds[i], elmts);
 305
 306         forAll(elmts, elmtI)
 307         {
 308             const label elI = elmts[elmtI];
 309
 310             if( elI < 0 || elI >= trias.size() )
 311             {
 312                 # ifdef USE_OMP
 313                 # pragma omp critical
 314                 # endif
 315                 FatalErrorIn
 316                 (
 317                     "void triSurf::topologyCheck()"
 318                 ) << "Triangle " << elI << " in facet subset "
 319                   << this->facetSubsetName(subsetIds[i])
 320                   << " is out of scope 0 " << trias.size() << exit(FatalError);
 321             }
 322         }
 323     }
 324
 325     //- check feature edge subsets
 326     subsetIds.clear();
 327     this->edgeSubsetIndices(subsetIds);
 328     forAll(subsetIds, i)
 329     {
 330         labelLongList elmts;
 331         this->edgesInSubset(subsetIds[i], elmts);
 332
 333         forAll(elmts, elmtI)
 334         {
 335             const label elI = elmts[elmtI];
 336
 337             if( elI < 0 || elI >= featureEdges.size() )
 338             {
 339                 # ifdef USE_OMP
 340                 # pragma omp critical
 341                 # endif
 342                 FatalErrorIn
 343                 (
 344                     "void triSurf::topologyCheck()"
 345                 ) << "Feature edge " << elI << " in edge subset "
 346                   << this->edgeSubsetName(subsetIds[i])
 347                   << " is out of scope 0 " << featureEdges.size()
 348                   << exit(FatalError);
 349             }
 350         }
 351     }
 352 }
 353
 354 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 355
 356 triSurf::triSurf()
 357 :
 358     triSurfPoints(),
 359     triSurfFacets(),
 360     triSurfFeatureEdges(),
 361     triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
 362 {}
 363
 364 //- Construct from parts
 365 triSurf::triSurf
 366 (
 367     const LongList<labelledTri>& triangles,
 368     const geometricSurfacePatchList& patches,
 369     const edgeLongList& featureEdges,
 370     const pointField& points
 371 )
 372 :
 373     triSurfPoints(points),
 374     triSurfFacets(triangles, patches),
 375     triSurfFeatureEdges(featureEdges),
 376     triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
 377 {
 378     topologyCheck();
 379 }
 380
 381 //- Read from file
 382 triSurf::triSurf(const fileName& fName)
 383 :
 384     triSurfPoints(),
 385     triSurfFacets(),
 386     triSurfFeatureEdges(),
 387     triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
 388 {
 389     readSurface(fName);
 390
 391     topologyCheck();
 392 }
 393
 394 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 395
 396 triSurf::~triSurf()
 397 {}
 398
 399 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 400
 401 void triSurf::readSurface(const fileName& fName)
 402 {
 403     if( fName.ext() == "fms" || fName.ext() == "FMS" )
 404     {
 405         readFromFMS(fName);
 406     }
 407     else if( fName.ext() == "ftr" || fName.ext() == "FTR" )
 408     {
 409         readFromFTR(fName);
 410     }
 411     else
 412     {
 413         triSurface copySurface(fName);
 414
 415         //- copy the points
 416         triSurfPoints::points_.setSize(copySurface.points().size());
 417         forAll(copySurface.points(), pI)
 418             triSurfPoints::points_[pI] = copySurface.points()[pI];
 419
 420         //- copy the triangles
 421         triSurfFacets::triangles_.setSize(copySurface.size());
 422         forAll(copySurface, tI)
 423             triSurfFacets::triangles_[tI] = copySurface[tI];
 424
 425         //- copy patches
 426         triSurfFacets::patches_ = copySurface.patches();
 427     }
 428 }
 429
 430 void triSurf::writeSurface(const fileName& fName) const
 431 {
 432     if( fName.ext() == "fms" || fName.ext() == "FMS" )
 433     {
 434         writeToFMS(fName);
 435     }
 436     else if( fName.ext() == "ftr" || fName.ext() == "FTR" )
 437     {
 438         writeToFTR(fName);
 439     }
 440     else
 441     {
 442         const pointField& pts = this->points();
 443         const LongList<labelledTri>& facets = this->facets();
 444         const geometricSurfacePatchList& patches = this->patches();
 445
 446         List<labelledTri> newTrias(facets.size());
 447         forAll(facets, tI)
 448             newTrias[tI] = facets[tI];
 449
 450         triSurface newSurf(newTrias, patches, pts);
 451         newSurf.write(fName);
 452     }
 453 }
 454
 455 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 456
 457 } // End namespace Foam
 458
 459 // ************************************************************************* //