src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMesh.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 "demandDrivenData.H"
  29 #include "partTriMesh.H"
  30 #include "meshSurfacePartitioner.H"
  31 #include "VRWGraphList.H"
  32 #include "triSurfModifier.H"
  33 #include "helperFunctions.H"
  34
  35 #include <map>
  36
  37 # ifdef USE_OMP
  38 #include <omp.h>
  39 # endif
  40
  41 //#define DEBUGSmooth
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 partTriMesh::partTriMesh(const meshSurfacePartitioner& mPart)
  51 :
  52     mPart_(mPart),
  53     surf_(),
  54     pointLabelInMeshSurface_(),
  55     meshSurfacePointLabelInTriMesh_(),
  56     pointType_(),
  57     globalPointLabelPtr_(NULL),
  58     pAtProcsPtr_(NULL),
  59     globalToLocalPointAddressingPtr_(NULL),
  60     neiProcsPtr_(NULL),
  61     pAtParallelBoundariesPtr_(NULL),
  62     pAtBufferLayersPtr_(NULL)
  63 {
  64     const meshSurfaceEngine& meshSurface = mPart.surfaceEngine();
  65     List<direction> useFace(meshSurface.boundaryFaces().size(), direction(1));
  66
  67     createPointsAndTrias(useFace);
  68 }
  69
  70 partTriMesh::partTriMesh
  71 (
  72     const meshSurfacePartitioner& mPart,
  73     const labelHashSet& invertedPoints,
  74     const label additionalLayers
  75 )
  76 :
  77     mPart_(mPart),
  78     surf_(),
  79     pointLabelInMeshSurface_(),
  80     meshSurfacePointLabelInTriMesh_(),
  81     pointType_(),
  82     globalPointLabelPtr_(NULL),
  83     pAtProcsPtr_(NULL),
  84     globalToLocalPointAddressingPtr_(NULL),
  85     neiProcsPtr_(NULL),
  86     pAtParallelBoundariesPtr_(NULL),
  87     pAtBufferLayersPtr_(NULL)
  88 {
  89     const meshSurfaceEngine& meshSurface = mPart.surfaceEngine();
  90     const VRWGraph& pointFaces = meshSurface.pointFaces();
  91     const faceList::subList& bFaces = meshSurface.boundaryFaces();
  92     const labelList& bPoints = meshSurface.boundaryPoints();
  93     const labelList& bp = meshSurface.bp();
  94
  95     List<direction> useFace(bFaces.size(), direction(0));
  96
  97     //- select cells containing at least one vertex of the bad faces
  98     forAll(pointFaces, bpI)
  99         if( invertedPoints.found(bPoints[bpI]) )
 100         {
 101             forAllRow(pointFaces, bpI, pfI)
 102                 useFace[pointFaces(bpI, pfI)] = 1;
 103         }
 104
 105     //- add additional layer of cells
 106     for(direction layerI=1;layerI<(additionalLayers+direction(1));++layerI)
 107     {
 108         forAll(useFace, bfI)
 109             if( useFace[bfI] == layerI )
 110             {
 111                 const face& bf = bFaces[bfI];
 112
 113                 forAll(bf, pI)
 114                 {
 115                     const label bpI = bp[bf[pI]];
 116
 117                     forAllRow(pointFaces, bpI, pfI)
 118                     {
 119                         const label fLabel = pointFaces(bpI, pfI);
 120
 121                         if( !useFace[fLabel] )
 122                             useFace[fLabel] = layerI + 1;
 123                     }
 124                 }
 125             }
 126
 127         if( Pstream::parRun() )
 128         {
 129             const labelList& globalPointLabel =
 130                 meshSurface.globalBoundaryPointLabel();
 131             const VRWGraph& pProcs = meshSurface.bpAtProcs();
 132             const Map<label>& globalToLocal =
 133                 meshSurface.globalToLocalBndPointAddressing();
 134
 135             std::map<label, labelLongList> eData;
 136             forAllConstIter(Map<label>, globalToLocal, iter)
 137             {
 138                 const label bpI = iter();
 139
 140                 forAllRow(pProcs, bpI, procI)
 141                 {
 142                     const label neiProc = pProcs(bpI, procI);
 143                     if( neiProc == Pstream::myProcNo() )
 144                         continue;
 145
 146                     if( eData.find(neiProc) == eData.end() )
 147                     {
 148                         eData.insert
 149                         (
 150                             std::make_pair(neiProc, labelLongList())
 151                         );
 152                     }
 153
 154                     forAllRow(pointFaces, bpI, pfI)
 155                         if( useFace[pointFaces(bpI, pfI)] == layerI )
 156                         {
 157                             eData[neiProc].append(globalPointLabel[bpI]);
 158                             break;
 159                         }
 160                 }
 161             }
 162
 163             //- exchange data with other processors
 164             labelLongList receivedData;
 165             help::exchangeMap(eData, receivedData);
 166
 167             forAll(receivedData, i)
 168             {
 169                 const label bpI = globalToLocal[receivedData[i]];
 170
 171                 forAllRow(pointFaces, bpI, pfI)
 172                 {
 173                     const label fLabel = pointFaces(bpI, pfI);
 174                     if( !useFace[fLabel] )
 175                         useFace[fLabel] = layerI + 1;
 176                 }
 177             }
 178         }
 179     }
 180
 181     createPointsAndTrias(useFace);
 182 }
 183
 184 partTriMesh::~partTriMesh()
 185 {
 186     deleteDemandDrivenData(globalPointLabelPtr_);
 187     deleteDemandDrivenData(pAtProcsPtr_);
 188     deleteDemandDrivenData(globalToLocalPointAddressingPtr_);
 189     deleteDemandDrivenData(neiProcsPtr_);
 190     deleteDemandDrivenData(pAtParallelBoundariesPtr_);
 191     deleteDemandDrivenData(pAtBufferLayersPtr_);
 192 }
 193
 194 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 195
 196 void partTriMesh::updateVertex(const label pI, const point& newP)
 197 {
 198     triSurfModifier sMod(surf_);
 199     pointField& pts = sMod.pointsAccess();
 200     const VRWGraph& pointFacets = surf_.pointFacets();
 201
 202     pts[pI] = newP;
 203
 204     if( pointType_[pI] & FACECENTRE )
 205     {
 206         Warning << "Smoothing auxiliary vertex."
 207             << " This has no effect on the original mesh" << endl;
 208         return;
 209     }
 210
 211     //- find face centres attached
 212     DynList<label> helper;
 213     forAllRow(pointFacets, pI, ptI)
 214     {
 215         const label centreI = surf_[pointFacets(pI, ptI)][2];
 216         if( pointType_[centreI] & FACECENTRE )
 217             helper.appendIfNotIn(centreI);
 218     }
 219
 220     //- update coordinates of FACECENTRE vertices
 221     forAll(helper, i)
 222     {
 223         const label centreI = helper[i];
 224
 225         point centre(vector::zero);
 226         scalar faceArea(0.0);
 227         forAllRow(pointFacets, centreI, ptI)
 228         {
 229             const labelledTri& tri = surf_[pointFacets(centreI, ptI)];
 230             point c(vector::zero);
 231             for(label i=0;i<3;++i)
 232                 c += pts[tri[i]];
 233             c /= 3;
 234             const scalar area = tri.mag(pts) + VSMALL;
 235
 236             centre += c * area;
 237             faceArea += area;
 238         }
 239
 240         pts[centreI] = centre / faceArea;
 241     }
 242 }
 243
 244 void partTriMesh::updateVerticesSMP(const List<LongList<labelledPoint> >& np)
 245 {
 246     triSurfModifier sMod(surf_);
 247     pointField& pts = sMod.pointsAccess();
 248     const VRWGraph& pointFacets = surf_.pointFacets();
 249
 250     List<direction> updateType(pts.size(), direction(0));
 251
 252     # ifdef USE_OMP
 253     # pragma omp parallel num_threads(np.size())
 254     # endif
 255     {
 256         # ifdef USE_OMP
 257         const LongList<labelledPoint>& newPoints = np[omp_get_thread_num()];
 258         # else
 259         const LongList<labelledPoint>& newPoints = np[0];
 260         # endif
 261
 262         forAll(newPoints, i)
 263         {
 264             const labelledPoint& lp = newPoints[i];
 265             const label pointI = lp.pointLabel();
 266
 267             pts[pointI] = lp.coordinates();
 268             updateType[pointI] |= SMOOTH;
 269
 270             forAllRow(pointFacets, pointI, ptI)
 271             {
 272                 const labelledTri& tri = surf_[pointFacets(pointI, ptI)];
 273
 274                 if( pointType_[tri[2]] & FACECENTRE )
 275                     updateType[tri[2]] |= FACECENTRE;
 276             }
 277         }
 278     }
 279
 280     //- update coordinates of buffer layer points
 281     if( Pstream::parRun() )
 282     {
 283         const labelLongList& bufferLayerPoints = this->bufferLayerPoints();
 284         const VRWGraph& pProcs = this->pointAtProcs();
 285         const labelLongList& globalPointLabel = this->globalPointLabel();
 286         const Map<label>& globalToLocal = this->globalToLocalPointAddressing();
 287         const DynList<label>& neiProcs = this->neiProcs();
 288
 289         //- create the map
 290         std::map<label, LongList<labelledPoint> > exchangeData;
 291         forAll(neiProcs, i)
 292             exchangeData.insert
 293             (
 294                 std::make_pair(neiProcs[i], LongList<labelledPoint>())
 295             );
 296
 297         //- add points into the map
 298         forAll(bufferLayerPoints, pI)
 299         {
 300             const label pointI = bufferLayerPoints[pI];
 301
 302             if( !(updateType[pointI] & SMOOTH) )
 303                 continue;
 304
 305             forAllRow(pProcs, pointI, i)
 306             {
 307                 const label neiProc = pProcs(pointI, i);
 308
 309                 if( neiProc == Pstream::myProcNo() )
 310                     continue;
 311
 312                 exchangeData[neiProc].append
 313                 (
 314                     labelledPoint(globalPointLabel[pointI], pts[pointI])
 315                 );
 316             }
 317         }
 318
 319         LongList<labelledPoint> receivedData;
 320         help::exchangeMap(exchangeData, receivedData);
 321
 322         forAll(receivedData, i)
 323         {
 324             const labelledPoint& lp = receivedData[i];
 325
 326             const label triPointI = globalToLocal[lp.pointLabel()];
 327             pts[triPointI] = lp.coordinates();
 328
 329             forAllRow(pointFacets, triPointI, ptI)
 330             {
 331                 const labelledTri& tri = surf_[pointFacets(triPointI, ptI)];
 332
 333                 if( pointType_[tri[2]] & FACECENTRE )
 334                     updateType[tri[2]] |= FACECENTRE;
 335             }
 336         }
 337     }
 338
 339     //- update coordinates of FACECENTRE vertices
 340     # ifdef USE_OMP
 341     # pragma omp parallel for schedule(dynamic, 20)
 342     # endif
 343     forAll(updateType, pI)
 344     {
 345         if( updateType[pI] & FACECENTRE )
 346         {
 347             point centre(vector::zero);
 348             scalar faceArea(0.0);
 349             forAllRow(pointFacets, pI, ptI)
 350             {
 351                 const labelledTri& tri = surf_[pointFacets(pI, ptI)];
 352                 point c(vector::zero);
 353                 for(label i=0;i<3;++i)
 354                     c += pts[tri[i]];
 355                 c /= 3;
 356                 const scalar area = tri.mag(pts) + VSMALL;
 357
 358                 centre += c * area;
 359                 faceArea += area;
 360             }
 361
 362             pts[pI] = centre / faceArea;
 363         }
 364     }
 365
 366     if( Pstream::parRun() )
 367         updateBufferLayers();
 368 }
 369
 370 void partTriMesh::updateVertices()
 371 {
 372     const meshSurfaceEngine& mse = mPart_.surfaceEngine();
 373     const labelList& bPoints = mse.boundaryPoints();
 374     labelLongList movedPoints(bPoints.size());
 375     forAll(movedPoints, i)
 376         movedPoints[i] = i;
 377
 378     updateVertices(movedPoints);
 379 }
 380
 381 void partTriMesh::updateVertices(const labelLongList& movedPoints)
 382 {
 383     const meshSurfaceEngine& mse = mPart_.surfaceEngine();
 384     const pointFieldPMG& points = mse.points();
 385     const labelList& bPoints = mse.boundaryPoints();
 386
 387     triSurfModifier sMod(surf_);
 388     pointField& pts = sMod.pointsAccess();
 389     const VRWGraph& pointFacets = surf_.pointFacets();
 390
 391     List<direction> updateType(pts.size(), direction(0));
 392
 393     //- update coordinates of vertices which exist in the surface
 394     //- of the volume mesh
 395     # ifdef USE_OMP
 396     # pragma omp parallel for schedule(dynamic, 50)
 397     # endif
 398     forAll(movedPoints, i)
 399     {
 400         const label bpI = movedPoints[i];
 401         const label pointI = bPoints[bpI];
 402         const label triPointI = meshSurfacePointLabelInTriMesh_[bpI];
 403
 404         if( triPointI < 0 )
 405             continue;
 406
 407         pts[triPointI] = points[pointI];
 408         updateType[triPointI] |= SMOOTH;
 409
 410         forAllRow(pointFacets, triPointI, ptI)
 411         {
 412             const labelledTri& tri = surf_[pointFacets(triPointI, ptI)];
 413
 414             if( pointType_[tri[2]] & FACECENTRE )
 415                 updateType[tri[2]] |= FACECENTRE;
 416         }
 417     }
 418
 419     //- update coordinates of buffer layer points
 420     if( Pstream::parRun() )
 421     {
 422         const labelLongList& bufferLayerPoints = this->bufferLayerPoints();
 423         const VRWGraph& pProcs = this->pointAtProcs();
 424         const labelLongList& globalPointLabel = this->globalPointLabel();
 425         const Map<label>& globalToLocal = this->globalToLocalPointAddressing();
 426         const DynList<label>& neiProcs = this->neiProcs();
 427
 428         //- create the map
 429         std::map<label, LongList<labelledPoint> > exchangeData;
 430         forAll(neiProcs, i)
 431             exchangeData.insert
 432             (
 433                 std::make_pair(neiProcs[i], LongList<labelledPoint>())
 434             );
 435
 436         //- add points into the map
 437         forAll(bufferLayerPoints, pI)
 438         {
 439             const label pointI = bufferLayerPoints[pI];
 440
 441             if( !(updateType[pointI] & SMOOTH) )
 442                 continue;
 443
 444             forAllRow(pProcs, pointI, i)
 445             {
 446                 const label neiProc = pProcs(pointI, i);
 447
 448                 if( neiProc == Pstream::myProcNo() )
 449                     continue;
 450
 451                 exchangeData[neiProc].append
 452                 (
 453                     labelledPoint(globalPointLabel[pointI], pts[pointI])
 454                 );
 455             }
 456         }
 457
 458         LongList<labelledPoint> receivedData;
 459         help::exchangeMap(exchangeData, receivedData);
 460
 461         forAll(receivedData, i)
 462         {
 463             const labelledPoint& lp = receivedData[i];
 464
 465             const label triPointI = globalToLocal[lp.pointLabel()];
 466             pts[triPointI] = lp.coordinates();
 467
 468             forAllRow(pointFacets, triPointI, ptI)
 469             {
 470                 const labelledTri& tri = surf_[pointFacets(triPointI, ptI)];
 471
 472                 if( pointType_[tri[2]] & FACECENTRE )
 473                     updateType[tri[2]] |= FACECENTRE;
 474             }
 475         }
 476     }
 477
 478     //- update coordinates of FACECENTRE vertices
 479     # ifdef USE_OMP
 480     # pragma omp parallel for schedule(dynamic, 20)
 481     # endif
 482     forAll(updateType, pI)
 483     {
 484         if( updateType[pI] & FACECENTRE )
 485         {
 486             point centre(vector::zero);
 487             scalar faceArea(0.0);
 488             forAllRow(pointFacets, pI, ptI)
 489             {
 490                 const labelledTri& tri = surf_[pointFacets(pI, ptI)];
 491                 point c(vector::zero);
 492                 for(label i=0;i<3;++i)
 493                     c += pts[tri[i]];
 494                 c /= 3;
 495                 const scalar area = tri.mag(pts) + VSMALL;
 496
 497                 centre += c * area;
 498                 faceArea += area;
 499             }
 500
 501             pts[pI] = centre / faceArea;
 502         }
 503     }
 504
 505     if( Pstream::parRun() )
 506         updateBufferLayers();
 507 }
 508
 509 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 510
 511 } // End namespace Foam
 512
 513 // ************************************************************************* //