src/INTERP_KERNEL/Interpolation3D.txx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __INTERPOLATION3D_TXX__
  22 #define __INTERPOLATION3D_TXX__
  23
  24 #include "Interpolation3D.hxx"
  25 #include "Interpolation.txx"
  26 #include "MeshElement.txx"
  27 #include "TransformedTriangle.hxx"
  28 #include "PolyhedronIntersectorP0P0.txx"
  29 #include "PointLocator3DIntersectorP0P0.txx"
  30 #include "PolyhedronIntersectorP0P1.txx"
  31 #include "PointLocator3DIntersectorP0P1.txx"
  32 #include "PolyhedronIntersectorP1P0.txx"
  33 #include "PolyhedronIntersectorP1P0Bary.txx"
  34 #include "PointLocator3DIntersectorP1P0.txx"
  35 #include "PolyhedronIntersectorP1P1.txx"
  36 #include "PointLocator3DIntersectorP1P1.txx"
  37 #include "Barycentric3DIntersectorP1P1.txx"
  38 #include "MappedBarycentric3DIntersectorP1P1.txx"
  39 #include "Log.hxx"
  40 // If defined, use recursion to traverse the binary search tree, else use the BBTree class
  41 //#define USE_RECURSIVE_BBOX_FILTER
  42
  43 #ifdef USE_RECURSIVE_BBOX_FILTER
  44 #include "MeshRegion.txx"
  45 #include "RegionNode.hxx"
  46 #include <stack>
  47
  48 #else // use BBTree class
  49
  50 #include "InterpolationHelper.txx"
  51
  52 #endif
  53
  54 #include <memory>
  55
  56 namespace INTERP_KERNEL
  57 {
  58   /**
  59    * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
  60    * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
  61    * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  62    * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the
  63    * intersection matrix.
  64    *
  65    * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  66    * It can also be an INTERP_KERNEL::Matrix object.
  67    * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  68    * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  69    * which have a non-zero intersection volume with the target element. The vector has indices running from
  70    * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  71    * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  72    *
  73
  74    * @param srcMesh     3-dimensional source mesh
  75    * @param targetMesh  3-dimesional target mesh, containing only tetraedra
  76    * @param result      matrix in which the result is stored
  77    *
  78    */
  79   template<class MyMeshType, class MatrixType>
  80   typename MyMeshType::MyConnType Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
  81   {
  82     using ConnType = typename MyMeshType::MyConnType;
  83     // create MeshElement objects corresponding to each element of the two meshes
  84     const ConnType numTargetElems = targetMesh.getNumberOfElements();
  85
  86     LOG(2, "Target mesh has " << numTargetElems << " elements ");
  87
  88     std::unique_ptr<Intersector3D<MyMeshType,MatrixType>> intersector;
  89     std::string methC = InterpolationOptions::filterInterpolationMethod(method);
  90     if(methC=="P0P0")
  91       {
  92         switch(InterpolationOptions::getIntersectionType())
  93           {
  94           case Triangulation:
  95             intersector.reset( new PolyhedronIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
  96             break;
  97           case PointLocator:
  98             intersector.reset( new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
  99             break;
 100           default:
 101             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangle or PointLocator.");
 102           }
 103       }
 104     else if(methC=="P0P1")
 105       {
 106         switch(InterpolationOptions::getIntersectionType())
 107           {
 108           case Triangulation:
 109             intersector.reset( new PolyhedronIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
 110             break;
 111           case PointLocator:
 112             intersector.reset( new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
 113             break;
 114           default:
 115             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P1 interp specified : must be Triangle or PointLocator.");
 116           }
 117       }
 118     else if(methC=="P1P0")
 119       {
 120         switch(InterpolationOptions::getIntersectionType())
 121           {
 122           case Triangulation:
 123             intersector.reset( new PolyhedronIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
 124             break;
 125           case PointLocator:
 126             intersector.reset( new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
 127             break;
 128           case Barycentric:
 129             intersector.reset( new PolyhedronIntersectorP1P0Bary<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
 130             break;
 131           default:
 132             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P0 interp specified : must be Triangle, PointLocator or Barycentric.");
 133           }
 134       }
 135     else if(methC=="P1P1")
 136       {
 137         switch(InterpolationOptions::getIntersectionType())
 138           {
 139           case Triangulation:
 140             intersector.reset( new PolyhedronIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy()) );
 141             break;
 142           case PointLocator:
 143             intersector.reset( new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
 144             break;
 145           case Barycentric:
 146             intersector.reset( new Barycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
 147             break;
 148           case MappedBarycentric:
 149             intersector.reset( new MappedBarycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
 150             break;
 151           default:
 152             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle, PointLocator, Barycentric or MappedBarycentric.");
 153           }
 154       }
 155     else
 156       throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
 157     // create empty maps for all source elements
 158     result.resize(intersector->getNumberOfRowsOfResMatrix());
 159
 160 #ifdef USE_RECURSIVE_BBOX_FILTER
 161
 162     /*
 163      * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh
 164      * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming
 165      * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the
 166      * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a
 167      * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements
 168      * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh.
 169      */
 170
 171     // create initial RegionNode and fill up its source region with all the source mesh elements and
 172     // its target region with all the target mesh elements whose bounding box
 173     // intersects that of the source region
 174
 175     RegionNode<ConnType>* firstNode = new RegionNode<ConnType>();
 176
 177     MeshRegion<ConnType>& srcRegion = firstNode->getSrcRegion();
 178
 179     for(ConnType i = 0 ; i < numSrcElems ; ++i)
 180       {
 181         srcRegion.addElement(srcElems[i], srcMesh);
 182       }
 183
 184     MeshRegion<ConnType>& targetRegion = firstNode->getTargetRegion();
 185
 186     for(ConnType i = 0 ; i < numTargetElems ; ++i)
 187       {
 188         if(!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]) ))
 189           {
 190             targetRegion.addElement(targetElems[i], targetMesh);
 191           }
 192       }
 193
 194     // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and
 195     // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the
 196     // elements of the target mesh whose bounding box intersects the corresponding part
 197     // Continue until the source region contains only one element, at which point the intersection volumes are
 198     // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero.
 199
 200     std::stack< RegionNode<ConnType>* > nodes;
 201     nodes.push(firstNode);
 202
 203     while(!nodes.empty())
 204       {
 205         RegionNode<ConnType>* currNode = nodes.top();
 206         nodes.pop();
 207         LOG(4, "Popping node ");
 208
 209         if(currNode->getTargetRegion().getNumberOfElements() == 1)
 210           {
 211             // calculate volumes
 212             LOG(4, " - One element");
 213
 214             MeshElement<ConnType>* targetElement = *(currNode->getTargetRegion().getBeginElements());
 215             std::vector<ConnType> intersectElems;
 216             for(typename std::vector< MeshElement<ConnType>* >::const_iterator iter = currNode->getSrcRegion().getBeginElements();iter != currNode->getSrcRegion().getEndElements();++iter)
 217               intersectElems.push_back((*iter)->getIndex());
 218             intersector->intersectCells(targetElement->getIndex(),intersectElems,result);
 219           }
 220         else // recursion
 221           {
 222
 223             LOG(4, " - Recursion");
 224
 225             RegionNode<ConnType>* leftNode = new RegionNode<ConnType>();
 226             RegionNode<ConnType>* rightNode = new RegionNode<ConnType>();
 227
 228             // split current source region
 229             //} decide on axis
 230             static BoundingBox::BoxCoord axis = BoundingBox::XMAX;
 231
 232             currNode->getTargetRegion().split(leftNode->getTargetRegion(), rightNode->getTargetRegion(), axis, targetMesh);
 233
 234             LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements()
 235                 << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements()
 236                 << " elements");
 237
 238             // ugly hack to avoid problem with enum which does not start at 0
 239             // I guess I ought to implement ++ for it instead ...
 240             // Anyway, it basically chooses the next axis, cyclically
 241             axis = (axis != BoundingBox::ZMAX) ? static_cast<BoundingBox::BoxCoord>(axis + 1) : BoundingBox::XMAX;
 242
 243             // add source elements of current node that overlap the target regions of the new nodes
 244             LOG(5, " -- Adding source elements");
 245             ConnType numLeftElements = 0;
 246             ConnType numRightElements = 0;
 247             for(typename std::vector<MeshElement<ConnType>*>::const_iterator iter = currNode->getSrcRegion().getBeginElements() ;
 248                 iter != currNode->getSrcRegion().getEndElements() ; ++iter)
 249               {
 250                 LOG(6, " --- New target node");
 251
 252                 if(!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter))
 253                   {
 254                     leftNode->getSrcRegion().addElement(*iter, srcMesh);
 255                     ++numLeftElements;
 256                   }
 257
 258                 if(!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter))
 259                   {
 260                     rightNode->getSrcRegion().addElement(*iter, srcMesh);
 261                     ++numRightElements;
 262                   }
 263
 264               }
 265
 266             LOG(5, "Left src region has " << numLeftElements << " elements and right src region has "
 267                 << numRightElements << " elements");
 268
 269             // push new nodes on stack
 270             if(numLeftElements != 0)
 271               {
 272                 nodes.push(leftNode);
 273               }
 274             else
 275               {
 276                 delete leftNode;
 277               }
 278
 279             if(numRightElements != 0)
 280               {
 281                 nodes.push(rightNode);
 282               }
 283             else
 284               {
 285                 delete rightNode;
 286               }
 287           }
 288
 289         // all nodes are deleted here
 290         delete currNode;
 291
 292         LOG(4, "Next iteration. Nodes left : " << nodes.size());
 293       }
 294
 295 #else // Use BBTree
 296     // create BBTree structure
 297     BBTreeStandAlone<3,ConnType> tree( BuildBBTree(srcMesh) );
 298
 299     // for each target element, get source elements with which to calculate intersection
 300     // - calculate intersection by calling intersectCells
 301     for(ConnType i = 0; i < numTargetElems; ++i)
 302       {
 303         MeshElement<ConnType> trgMeshElem(i, targetMesh);
 304
 305         const BoundingBox *box = trgMeshElem.getBoundingBox();
 306
 307         // get target bbox in right order
 308         double targetBox[6];
 309         box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox);
 310
 311         std::vector<ConnType> intersectElems;
 312
 313         tree.getIntersectingElems(targetBox, intersectElems);
 314
 315         if ( !intersectElems.empty() )
 316           intersector->intersectCells(i,intersectElems,result);
 317       }
 318
 319 #endif
 320     return intersector->getNumberOfColsOfResMatrix();
 321   }
 322 }
 323
 324 #endif