src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H

   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 Class
  25     GGIInterpolation
  26
  27 Description
  28     Mass-conservative face interpolation of face data between two
  29     primitivePatches
  30
  31 Author
  32     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  33
  34 Contributor:
  35     Martin Beaudoin, Hydro-Quebec, (2008)
  36
  37 SourceFiles
  38     GGIInterpolation.C
  39     GGIInterpolate.C
  40     GGIInterpolationWeights.C
  41
  42 \*---------------------------------------------------------------------------*/
  43
  44 #ifndef GGIInterpolation_H
  45 #define GGIInterpolation_H
  46
  47 #include "className.H"
  48 #include "labelList.H"
  49 #include "scalarField.H"
  50 #include "pointField.H"
  51 #include "FieldFields.H"
  52 #include "faceList.H"
  53 #include "intersection.H"
  54 #include "point2D.H"
  55 #include "NamedEnum.H"
  56 #include "tolerancesSwitch.H"
  57 #include "optimisationSwitch.H"
  58
  59 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  60
  61 namespace Foam
  62 {
  63
  64 /*---------------------------------------------------------------------------*\
  65                    Class GGIInterpolationName Declaration
  66 \*---------------------------------------------------------------------------*/
  67
  68 class GGIInterpolationName
  69 {
  70 public:
  71
  72     // Public enumerations
  73
  74         //- Define a return type for isVertexInsidePolygon.
  75         //  The names are rather obvious
  76         // on vertex, or on edge or inside == inside
  77         enum insideOutside
  78         {
  79             ALL_OUTSIDE,
  80             ALL_INSIDE,
  81             PARTIALLY_OVERLAPPING
  82         };
  83
  84         //- Quick reject overlap search algorithm
  85         enum quickReject
  86         {
  87             THREE_D_DISTANCE,
  88             AABB,
  89             BB_OCTREE,
  90             N_SQUARED
  91         };
  92
  93
  94     // Static data
  95
  96         ClassName("GGIInterpolation");
  97
  98         //- Quick reject names
  99         static const NamedEnum<quickReject, 4> quickRejectNames_;
 100
 101
 102     // Constructors
 103
 104         //- Construct null
 105         GGIInterpolationName()
 106         {}
 107
 108 };
 109
 110
 111 /*---------------------------------------------------------------------------*\
 112                       Class GGIInterpolation Declaration
 113 \*---------------------------------------------------------------------------*/
 114
 115 template<class MasterPatch, class SlavePatch>
 116 class GGIInterpolation
 117 :
 118     public GGIInterpolationName
 119 {
 120     // Private data types
 121
 122         // Definition for a 3D orthoNormalBasis type
 123         typedef VectorSpace<Vector<vector>, vector, 3> orthoNormalBasis;
 124
 125
 126     // Private data
 127
 128         //- Reference to the master patch
 129         const MasterPatch& masterPatch_;
 130
 131         //- Reference to the slave patch
 132         const SlavePatch& slavePatch_;
 133
 134         //- Slave-to-master transformation tensor.  Transforms slave data to
 135         //  master plane.  Size equals number of slave faces; zero length
 136         //  indicates no transform.  Size 1 indicates constant transform
 137         //  HJ, 6/Jan/2009
 138         tensorField forwardT_;
 139
 140         //- Master-to-slave transformation tensor.  Transforms slave data to
 141         //  master plane.  Size equals number of master faces; zero length
 142         //  indicates no transform.  Size 1 indicates constant transform
 143         //  HJ, 6/Jan/2009
 144         tensorField reverseT_;
 145
 146         //- Slave-to-master separation vector. Translation of slave data to
 147         //  master plane.  Size equals number of slave faces; zero length
 148         //  indicates no translation.  MB, 28/Jan/2009
 149         vectorField forwardSep_;
 150
 151         //- Master non-overlap face tolerance factor
 152         const scalar masterNonOverlapFaceTol_;
 153
 154         //- Slave non-overlap face tolerance factor
 155         const scalar slaveNonOverlapFaceTol_;
 156
 157         //- Rescale the GGI weighting factors or not
 158         const bool rescaleGGIWeightingFactors_;
 159
 160         //- Quick reject search algorithm
 161         quickReject reject_;
 162
 163
 164     // Demand-driven data
 165
 166         // Master-to-slave interpolation
 167
 168             //- Master to slave addressing
 169             mutable labelListList* masterAddrPtr_;
 170
 171             //- Master to slave  weights
 172             mutable scalarListList* masterWeightsPtr_;
 173
 174
 175         // Slave-to-master interpolation
 176
 177             //- Slave to master addressing
 178             mutable labelListList* slaveAddrPtr_;
 179
 180             //- Slave to master weights
 181             mutable scalarListList* slaveWeightsPtr_;
 182
 183             //- List of uncovered master faces
 184             mutable labelList* uncoveredMasterAddrPtr_;
 185
 186             //- List of uncovered slave faces
 187             mutable labelList* uncoveredSlaveAddrPtr_;
 188
 189
 190     // Private static data
 191
 192         //- Facet area error tolerance
 193         static const debug::tolerancesSwitch areaErrorTol_;
 194
 195         //- Facet normal featureCos criteria
 196         static const debug::tolerancesSwitch featureCosTol_;
 197
 198         //- Facet bound box extension factor
 199         static const debug::tolerancesSwitch faceBoundBoxExtendSpanFraction_;
 200
 201        //- The next 3 attributes are parameters controlling the creation
 202        //  of an octree search engine for the GGI facets neighbourhood
 203        //  determination.
 204        //
 205        //  Tests using GGI patches of up to ~100K facets have shown that
 206        //  the following values gave the best compromise between the
 207        //  "octree creation time" vs "octree search time":
 208        //
 209        //     octreeSearchMinNLevel_     : 3
 210        //     octreeSearchMaxLeafRatio_  : 3
 211        //     octreeSearchMaxShapeRatio_ : 1
 212        //
 213        //  For GGI patches larger than ~100K facets, your mileage may vary.
 214        //  So these 3 control parameters are adjustable using the following
 215        //  global optimization switches:
 216        //
 217        //     GGIOctreeSearchMinNLevel
 218        //     GGIOctreeSearchMaxLeafRatio
 219        //     GGIOctreeSearchMaxShapeRatio
 220        //
 221
 222        //- Octree search: minNlevel parameter for octree constructor
 223        static const debug::optimisationSwitch octreeSearchMinNLevel_;
 224
 225        //- Octree search: maxLeafRatio parameter for octree constructor
 226        static const debug::optimisationSwitch octreeSearchMaxLeafRatio_;
 227
 228        //- Octree search: maxShapeRatio parameter for octree constructor
 229        static const debug::optimisationSwitch octreeSearchMaxShapeRatio_;
 230
 231
 232     // Private Member Functions
 233
 234         //- Disallow default bitwise copy construct
 235         GGIInterpolation(const GGIInterpolation&);
 236
 237         //- Disallow default bitwise assignment
 238         void operator=(const GGIInterpolation&);
 239
 240
 241     // Helper functions for demand-driven data
 242
 243         //- Evaluate faces neighborhood based of sphere defined by faces BB
 244         void findNeighbours3D(labelListList& result) const;
 245
 246         //- Evaluate faces neighborhood based of faces Axis Aligned BB
 247         void findNeighboursAABB(labelListList& result) const;
 248
 249         //- Evaluate faces neighborhood based of faces BB and an octree
 250         //  search engine
 251         void findNeighboursBBOctree(labelListList& result) const;
 252
 253         //- Projects a list of points onto a plane located at
 254         //  planeOrig, oriented along planeNormal
 255         tmp<pointField> projectPointsOnPlane
 256         (
 257             const pointField& lpoints,
 258             const vector& planeOrig,
 259             const vector& planeDirection,
 260             scalarField& distanceProjection
 261         ) const;
 262
 263         //- Compute an orthonormal basis. Useful for 2D projection
 264         orthoNormalBasis computeOrthonormalBasis
 265         (
 266             const vector& normalVectorCentre,
 267             const vector& normalVector,
 268             const pointField& pointsOnPlane
 269         ) const;
 270
 271         //- Compute the 2D projection of a 3D list of point into the
 272         //  orthonormal basis orthoBase
 273         List<point2D> projectPoints3Dto2D
 274         (
 275             const orthoNormalBasis& orthoBase,
 276             const vector& orthoBaseOffset,
 277             const pointField& pointsIn3D,
 278             scalarField& distanceProjection
 279         ) const;
 280
 281         //- Compute the 1D projection of a 2D list of points onto the
 282         //  projectionDir direction
 283         scalarField projectPoints2Dto1D
 284         (
 285             const vector2D& normalizedProjectionDir,
 286             const point2D& normalizedProjectionDirOffset,
 287             const List<point2D>& lPoints2D
 288         ) const;
 289
 290         //- Detect the intersection of 2 2D polygons using the method of
 291         //  Separating Axes Theorem
 292         bool detect2dPolygonsOverlap
 293         (
 294             const List<point2D>& poly1,
 295             const List<point2D>& poly2,
 296             const scalar& tolFactor,
 297             const bool firstCall = true // Optional parameter: first call
 298         ) const;
 299
 300         //- Compute the intersection area between 2 2D polygons
 301         scalar polygonIntersection
 302         (
 303             const List<point2D>& poly1,
 304             const List<point2D>& poly2
 305         ) const;
 306
 307         //- Detect is the points from a subject polygon are inside or
 308         //  outside a clipping polygon
 309         insideOutside isVertexInsidePolygon
 310         (
 311             const List<point2D>& clippingPolygon,
 312             const List<point2D>& subjectPolygon,
 313             List<bool>& subjectVertexInside
 314         ) const;
 315
 316         //- Compute the intersection between two 2D polygons using the
 317         //  Sutherland Hodgman algorithm
 318         List<point2D> clipPolygon2DSutherlandHodgman
 319         (
 320             const List<point2D>& clippingPolygon,
 321             const List<point2D>& subjectPolygon
 322         ) const;
 323
 324         //- Compute the intersection between two 2D polygons using the
 325         //  Greiner Hormann algorithm
 326         List<point2D> clipPolygon2DGreinerHormann
 327         (
 328             const List<point2D>& clippingPolygon,
 329             const List<point2D>& subjectPolygon,
 330             const List<bool>& subjectVertexInside
 331         ) const;
 332
 333         //- Compute the area of a 2D polygon
 334         scalar area2D
 335         (
 336             const List<point2D>& polygon
 337         ) const;
 338
 339
 340         //- Calculate addressing and weights
 341         void calcAddressing() const;
 342
 343         //- Rescale GGI weighting factors
 344         void rescaleWeightingFactors() const;
 345
 346         //- Find non-overlapping faces
 347         tmp<labelField> findNonOverlappingFaces
 348         (
 349             const scalarListList& patchWeights,
 350             const scalar& nonOverlapFaceTol
 351         ) const;
 352
 353         //- Clear all geometry and addressing
 354         void clearOut();
 355
 356
 357     // Interpolation data access
 358
 359         //- Interpolate given source and target, addressing and weights
 360         template<class Type>
 361         static void interpolate
 362         (
 363             const Field<Type>& ff,
 364             Field<Type>& result,
 365             const labelListList& addr,
 366             const scalarListList& weights
 367         );
 368
 369         //- Interpolate given source and target, addressing and weights
 370         //  for masked faces only
 371         template<class Type>
 372         static void maskedInterpolate
 373         (
 374             const Field<Type>& ff,
 375             Field<Type>& result,
 376             const labelList& mask,
 377             const labelListList& addr,
 378             const scalarListList& weights
 379         );
 380
 381         //- Bridge uncovered faces given addressing
 382         template<class Type>
 383         static void bridge
 384         (
 385             const Field<Type>& bridgeField,
 386             Field<Type>& ff,
 387             const labelList& addr
 388         );
 389
 390         //- Bridge uncovered faces given addressing
 391         //  for masked faces only
 392         template<class Type>
 393         static void maskedBridge
 394         (
 395             const Field<Type>& bridgeField,
 396             Field<Type>& ff,
 397             const labelList& mask,
 398             const labelList& uncoveredFaces
 399         );
 400
 401         //- Is a transform required?
 402         inline bool doTransform() const
 403         {
 404             return forwardT_.size() > 0;
 405         }
 406
 407         //- Is a separation required?
 408         inline bool doSeparation() const
 409         {
 410             return forwardSep_.size() > 0;
 411         }
 412
 413
 414 public:
 415
 416     // Constructors
 417
 418         //- Construct from components
 419         GGIInterpolation
 420         (
 421             const MasterPatch& masterPatch,
 422             const SlavePatch&  slavePatch,
 423             const tensorField& forwardT,
 424             const tensorField& reverseT,
 425             const vectorField& forwardSep,
 426             const scalar masterFaceNonOverlapFaceTol = 0,
 427             const scalar slaveFaceNonOverlapFaceTol = 0,
 428             const bool rescaleGGIWeightingFactors = true,
 429             const quickReject reject = AABB
 430         );
 431
 432
 433     // Destructor
 434
 435         ~GGIInterpolation();
 436
 437
 438     // Member Functions
 439
 440         // Access
 441
 442             //- Reference to the master patch
 443             const MasterPatch& masterPatch() const
 444             {
 445                 return masterPatch_;
 446             }
 447
 448             //- Reference to the slave patch
 449             const SlavePatch& slavePatch() const
 450             {
 451                 return slavePatch_;
 452             }
 453
 454             //- Return reference to master addressing
 455             const labelListList& masterAddr() const;
 456
 457             //- Return reference to master weights
 458             const scalarListList& masterWeights() const;
 459
 460             //- Return reference to slave addressing
 461             const labelListList& slaveAddr() const;
 462
 463             //- Return reference to slave weights
 464             const scalarListList& slaveWeights() const;
 465
 466             //- Return reference to the master list of non-overlap faces
 467             const labelList& uncoveredMasterFaces() const;
 468
 469             //- Return reference to the slave list of non-overlap faces
 470             const labelList& uncoveredSlaveFaces() const;
 471
 472
 473     // Interpolation functions
 474
 475         //- Interpolate from master to slave
 476         template<class Type>
 477         tmp<Field<Type> > masterToSlave(const Field<Type>& ff) const;
 478
 479         template<class Type>
 480         tmp<Field<Type> > masterToSlave
 481         (
 482             const tmp<Field<Type> >& tff
 483         ) const;
 484
 485         //- Interpolate from master to slave, only for marked master faces
 486         //  Addressing picks the faces from patch that are selected
 487         //  for collection into the result field
 488         template<class Type>
 489         void maskedMasterToSlave
 490         (
 491             const Field<Type>& ff,
 492             Field<Type>& result,
 493             const labelList& mask
 494         ) const;
 495
 496
 497         //- Interpolate from slave to master
 498         template<class Type>
 499         tmp<Field<Type> > slaveToMaster(const Field<Type>& ff) const;
 500
 501         template<class Type>
 502         tmp<Field<Type> > slaveToMaster
 503         (
 504             const tmp<Field<Type> >& tff
 505         ) const;
 506
 507         //- Interpolate from slave to master, only for marked slave faces
 508         //  Addressing picks the faces from patch that are selected
 509         //  for collection into the result field
 510         template<class Type>
 511         void maskedSlaveToMaster
 512         (
 513             const Field<Type>& ff,
 514             Field<Type>& result,
 515             const labelList& mask
 516         ) const;
 517
 518
 519         //- Bridge uncovered master patch field
 520         template<class Type>
 521         void bridgeMaster
 522         (
 523             const Field<Type>& bridgeField,
 524             Field<Type>& ff
 525         ) const;
 526
 527         //- Bridge uncovered master patch field, only for marked master faces
 528         template<class Type>
 529         void maskedBridgeMaster
 530         (
 531             const Field<Type>& bridgeField,
 532             Field<Type>& ff,
 533             const labelList& mask
 534         ) const;
 535
 536         //- Bridge uncovered slave patch field
 537         template<class Type>
 538         void bridgeSlave
 539         (
 540             const Field<Type>& bridgeField,
 541             Field<Type>& ff
 542         ) const;
 543
 544         //- Bridge uncovered slave patch field, only for marked slave faces
 545         template<class Type>
 546         void maskedBridgeSlave
 547         (
 548             const Field<Type>& bridgeField,
 549             Field<Type>& ff,
 550             const labelList& mask
 551         ) const;
 552
 553
 554     // Edit
 555
 556         //- Correct weighting factors for moving mesh.
 557         bool movePoints
 558         (
 559             const tensorField& forwardT,
 560             const tensorField& reverseT,
 561             const vectorField& forwardSep
 562         );
 563 };
 564
 565
 566 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 567
 568 } // End namespace Foam
 569
 570 #ifdef NoRepository
 571 #   include "GGIInterpolation.C"
 572 #endif
 573
 574 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 575
 576 #endif
 577
 578 // ************************************************************************* //