src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H

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