src/mesh/cfMesh/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperCornersAndEdges.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 "meshOctree.H"
  29 #include "triSurf.H"
  30 #include "triSurfacePartitioner.H"
  31 #include "meshSurfaceMapper.H"
  32 #include "meshSurfaceEngine.H"
  33 #include "meshSurfaceEngineModifier.H"
  34 #include "meshSurfacePartitioner.H"
  35 #include "labelledScalar.H"
  36
  37 #include "helperFunctions.H"
  38
  39 # ifdef USE_OMP
  40 #include <omp.h>
  41 # endif
  42
  43 //#define DEBUGMapping
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 namespace Foam
  48 {
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52 void meshSurfaceMapper::findMappingDistance
  53 (
  54     const labelLongList& nodesToMap,
  55     scalarList& mappingDistance
  56 ) const
  57 {
  58     const vectorField& faceCentres = surfaceEngine_.faceCentres();
  59     const VRWGraph& pFaces = surfaceEngine_.pointFaces();
  60     const labelList& bPoints = surfaceEngine_.boundaryPoints();
  61     const pointFieldPMG& points = surfaceEngine_.points();
  62
  63     //- generate search distance for corner nodes
  64     mappingDistance.setSize(nodesToMap.size());
  65     # ifdef USE_OMP
  66     # pragma omp parallel for schedule(dynamic, 50)
  67     # endif
  68     forAll(nodesToMap, i)
  69     {
  70         const label bpI = nodesToMap[i];
  71
  72         mappingDistance[i] = 0.0;
  73
  74         const point& p = points[bPoints[bpI]];
  75         forAllRow(pFaces, bpI, pfI)
  76         {
  77             const scalar d = magSqr(faceCentres[pFaces(bpI, pfI)] - p);
  78             mappingDistance[i] = Foam::max(mappingDistance[i], d);
  79         }
  80
  81         //- safety factor
  82         mappingDistance[i] *= 4.0;
  83     }
  84
  85     if( Pstream::parRun() )
  86     {
  87         //- make sure that corner nodesd at parallel boundaries
  88         //- have the same range in which they accept the corners
  89         const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
  90         const labelList& globalPointLabel =
  91             surfaceEngine_.globalBoundaryPointLabel();
  92
  93         //- create the map for exchanging data
  94         std::map<label, DynList<labelledScalar> > exchangeData;
  95         const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
  96         forAll(neiProcs, i)
  97             exchangeData.insert
  98             (
  99                 std::make_pair(neiProcs[i], DynList<labelledScalar>())
 100             );
 101
 102         Map<label> globalToLocal;
 103
 104         forAll(nodesToMap, nI)
 105         {
 106             const label bpI = nodesToMap[nI];
 107
 108             if( bpAtProcs.sizeOfRow(bpI) != 0 )
 109                 globalToLocal.insert(globalPointLabel[bpI], nI);
 110
 111             forAllRow(bpAtProcs, bpI, i)
 112             {
 113                 const label neiProc = bpAtProcs(bpI, i);
 114                 if( neiProc == Pstream::myProcNo() )
 115                     continue;
 116
 117                 exchangeData[neiProc].append
 118                 (
 119                     labelledScalar(globalPointLabel[bpI], mappingDistance[nI])
 120                 );
 121             }
 122         }
 123
 124         //- exchange data between processors
 125         LongList<labelledScalar> receivedData;
 126         help::exchangeMap(exchangeData, receivedData);
 127
 128         //- select the maximum mapping distance for processor points
 129         forAll(receivedData, i)
 130         {
 131             const labelledScalar& ls = receivedData[i];
 132
 133             const label nI = globalToLocal[ls.scalarLabel()];
 134
 135             //- choose the maximum value for the mapping distance
 136             mappingDistance[nI] = Foam::max(mappingDistance[nI], ls.value());
 137         }
 138     }
 139 }
 140
 141 scalar meshSurfaceMapper::faceMetricInPatch
 142 (
 143     const label bfI,
 144     const label patchI
 145 ) const
 146 {
 147     const face& bf = surfaceEngine_.boundaryFaces()[bfI];
 148
 149     const pointFieldPMG& points = surfaceEngine_.points();
 150
 151     const point centre = bf.centre(points);
 152     const vector area = bf.normal(points);
 153
 154     point projCentre;
 155     scalar dSq;
 156     label nt;
 157
 158     meshOctree_.findNearestSurfacePointInRegion
 159     (
 160         projCentre,
 161         dSq,
 162         nt,
 163         patchI,
 164         centre
 165     );
 166
 167     DynList<point> projPoints(bf.size());
 168     forAll(bf, pI)
 169     {
 170         meshOctree_.findNearestSurfacePointInRegion
 171         (
 172             projPoints[pI],
 173             dSq,
 174             nt,
 175             patchI,
 176             points[bf[pI]]
 177         );
 178     }
 179
 180     vector projArea(vector::zero);
 181     forAll(bf, pI)
 182     {
 183         projArea +=
 184             triPointRef
 185             (
 186                 projPoints[pI],
 187                 projPoints[bf.fcIndex(pI)],
 188                 projCentre
 189             ).normal();
 190     }
 191
 192     return magSqr(centre - projCentre) + mag(mag(projArea) - mag(area));
 193 }
 194
 195 void meshSurfaceMapper::mapCorners(const labelLongList& nodesToMap)
 196 {
 197     const triSurfacePartitioner& sPartitioner = surfacePartitioner();
 198     const labelList& surfCorners = sPartitioner.corners();
 199     const List<DynList<label> >& cornerPatches = sPartitioner.cornerPatches();
 200
 201     const meshSurfacePartitioner& mPart = meshPartitioner();
 202     const labelHashSet& corners = mPart.corners();
 203     const VRWGraph& pPatches = mPart.pointPatches();
 204
 205     const pointFieldPMG& points = surfaceEngine_.points();
 206     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 207
 208     const triSurf& surf = meshOctree_.surface();
 209     const pointField& sPoints = surf.points();
 210
 211     //std::map<label, scalar> mappingDistance;
 212     scalarList mappingDistance;
 213     findMappingDistance(nodesToMap, mappingDistance);
 214
 215     //- for every corner in the mesh surface find the nearest corner in the
 216     //- triSurface
 217     meshSurfaceEngineModifier sMod(surfaceEngine_);
 218
 219     # ifdef USE_OMP
 220     # pragma omp parallel for schedule(dynamic, 50)
 221     # endif
 222     forAll(nodesToMap, cornerI)
 223     {
 224         const label bpI = nodesToMap[cornerI];
 225         if( !corners.found(bpI) )
 226             FatalErrorIn
 227             (
 228                 "meshSurfaceMapper::mapCorners(const labelLongList&)"
 229             ) << "Trying to map a point that is not a corner"
 230                 << abort(FatalError);
 231
 232         const point& p = points[bPoints[bpI]];
 233         const scalar maxDist = mappingDistance[cornerI];
 234
 235         //- find the nearest position to the given point patches
 236         const DynList<label> patches = pPatches[bpI];
 237
 238         point mapPointApprox(p);
 239         scalar distSqApprox;
 240
 241         label iter(0);
 242         while( iter++ < 20 )
 243         {
 244             point newP(vector::zero);
 245             forAll(patches, patchI)
 246             {
 247                 point np;
 248                 label nt;
 249                 meshOctree_.findNearestSurfacePointInRegion
 250                 (
 251                     np,
 252                     distSqApprox,
 253                     nt,
 254                     patches[patchI],
 255                     mapPointApprox
 256                 );
 257
 258                 newP += np;
 259             }
 260
 261             newP /= patches.size();
 262
 263             if( magSqr(newP - mapPointApprox) < 1e-8 * maxDist )
 264                 break;
 265
 266             mapPointApprox = newP;
 267         }
 268         distSqApprox = magSqr(mapPointApprox - p);
 269
 270         //- find the nearest triSurface corner for the given corner
 271         scalar distSq(mappingDistance[cornerI]);
 272         point mapPoint(p);
 273         forAll(surfCorners, scI)
 274         {
 275             const label cornerID = surfCorners[scI];
 276             const point& sp = sPoints[cornerID];
 277
 278             if( Foam::magSqr(sp - p) < distSq )
 279             {
 280                 bool store(true);
 281                 const DynList<label>& cPatches = cornerPatches[scI];
 282                 forAll(patches, i)
 283                 {
 284                     if( !cPatches.contains(patches[i]) )
 285                     {
 286                         store = false;
 287                         break;
 288                     }
 289                 }
 290
 291                 if( store )
 292                 {
 293                     mapPoint = sp;
 294                     distSq = Foam::magSqr(sp - p);
 295                 }
 296             }
 297         }
 298
 299         if( distSq > 1.2 * distSqApprox )
 300         {
 301             mapPoint = mapPointApprox;
 302         }
 303
 304         //- move the point to the nearest corner
 305         sMod.moveBoundaryVertexNoUpdate(bpI, mapPoint);
 306     }
 307
 308     sMod.updateGeometry(nodesToMap);
 309 }
 310
 311 void meshSurfaceMapper::mapEdgeNodes(const labelLongList& nodesToMap)
 312 {
 313     const pointFieldPMG& points = surfaceEngine_.points();
 314     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 315
 316     const meshSurfacePartitioner& mPart = meshPartitioner();
 317     const VRWGraph& pPatches = mPart.pointPatches();
 318
 319     //const triSurf& surf = meshOctree_.surface();
 320     //const pointField& sPoints = surf.points();
 321
 322     //- find mapping distance for selected vertices
 323     scalarList mappingDistance;
 324     findMappingDistance(nodesToMap, mappingDistance);
 325
 326     const VRWGraph* bpAtProcsPtr(NULL);
 327     if( Pstream::parRun() )
 328         bpAtProcsPtr = &surfaceEngine_.bpAtProcs();
 329
 330     LongList<parMapperHelper> parallelBndNodes;
 331
 332     meshSurfaceEngineModifier sMod(surfaceEngine_);
 333
 334     //- map point to the nearest vertex on the surface mesh
 335     # ifdef USE_OMP
 336     # pragma omp parallel for schedule(dynamic, 50)
 337     # endif
 338     forAll(nodesToMap, i)
 339     {
 340         const label bpI = nodesToMap[i];
 341         const point& p = points[bPoints[bpI]];
 342
 343         //- find patches at this edge point
 344         const DynList<label> patches = pPatches[bpI];
 345
 346         const scalar maxDist = mappingDistance[i];
 347
 348         //- find approximate position of the vertex on the edge
 349         point mapPointApprox(p);
 350         scalar distSqApprox;
 351         label iter(0);
 352         while( iter++ < 20 )
 353         {
 354             point newP(vector::zero);
 355
 356             forAll(patches, patchI)
 357             {
 358                 point np;
 359                 label nt;
 360                 meshOctree_.findNearestSurfacePointInRegion
 361                 (
 362                     np,
 363                     distSqApprox,
 364                     nt,
 365                     patches[patchI],
 366                     mapPointApprox
 367                 );
 368
 369                 newP += np;
 370             }
 371
 372             newP /= patches.size();
 373
 374             if( magSqr(newP - mapPointApprox) < 1e-8 * maxDist )
 375                 break;
 376
 377             mapPointApprox = newP;
 378         }
 379         distSqApprox = magSqr(mapPointApprox - p);
 380
 381         //- find the nearest vertex on the triSurface feature edge
 382         point mapPoint(mapPointApprox);
 383         scalar distSq(distSqApprox);
 384         label nse;
 385         meshOctree_.findNearestEdgePoint(mapPoint, distSq, nse, p, patches);
 386
 387         //- use the vertex with the smallest mapping distance
 388         if( distSq > 1.2 * distSqApprox )
 389         {
 390             mapPoint = mapPointApprox;
 391             distSq = distSqApprox;
 392         }
 393
 394         //- check if the mapping distance is within the given tolerances
 395         if( distSq > maxDist )
 396         {
 397             //- this indicates possible problems
 398             //- reduce the mapping distance
 399             const scalar f = Foam::sqrt(maxDist / distSq);
 400             distSq = mappingDistance[i];
 401             mapPoint = f * (mapPoint - p) + p;
 402         }
 403
 404         //- move the point to the nearest edge vertex
 405         sMod.moveBoundaryVertexNoUpdate(bpI, mapPoint);
 406
 407         if( bpAtProcsPtr && bpAtProcsPtr->sizeOfRow(bpI) )
 408         {
 409             # ifdef USE_OMP
 410             # pragma omp critical
 411             # endif
 412             {
 413                 parallelBndNodes.append
 414                 (
 415                     parMapperHelper
 416                     (
 417                         mapPoint,
 418                         distSq,
 419                         bpI,
 420                         -1
 421                     )
 422                 );
 423             }
 424         }
 425     }
 426
 427     sMod.updateGeometry(nodesToMap);
 428
 429     mapToSmallestDistance(parallelBndNodes);
 430 }
 431
 432 void meshSurfaceMapper::mapCornersAndEdges()
 433 {
 434     const meshSurfacePartitioner& mPart = meshPartitioner();
 435     const labelHashSet& cornerPoints = mPart.corners();
 436     labelLongList selectedPoints;
 437     forAllConstIter(labelHashSet, cornerPoints, it)
 438         selectedPoints.append(it.key());
 439
 440     mapCorners(selectedPoints);
 441
 442     selectedPoints.clear();
 443     const labelHashSet& edgePoints = mPart.edgePoints();
 444     forAllConstIter(labelHashSet, edgePoints, it)
 445         selectedPoints.append(it.key());
 446
 447     mapEdgeNodes(selectedPoints);
 448 }
 449
 450 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 451
 452 } // End namespace Foam
 453
 454 // ************************************************************************* //