src/INTERP_KERNEL/Interpolation3D.txx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   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 "BBTree.txx"
  51
  52 #endif
  53
  54 namespace INTERP_KERNEL
  55 {
  56   /**
  57    * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
  58    * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
  59    * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  60    * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the
  61    * intersection matrix.
  62    *
  63    * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  64    * It can also be an INTERP_KERNEL::Matrix object.
  65    * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  66    * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  67    * which have a non-zero intersection volume with the target element. The vector has indices running from
  68    * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  69    * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  70    *
  71
  72    * @param srcMesh     3-dimensional source mesh
  73    * @param targetMesh  3-dimesional target mesh, containing only tetraedra
  74    * @param result      matrix in which the result is stored
  75    *
  76    */
  77   template<class MyMeshType, class MatrixType>
  78   typename MyMeshType::MyConnType Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
  79   {
  80     typedef typename MyMeshType::MyConnType ConnType;
  81     // create MeshElement objects corresponding to each element of the two meshes
  82     const ConnType numSrcElems = srcMesh.getNumberOfElements();
  83     const ConnType numTargetElems = targetMesh.getNumberOfElements();
  84
  85     LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
  86
  87     std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
  88     std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
  89
  90     std::map<MeshElement<ConnType>*, ConnType> indices;
  91
  92     for(ConnType i = 0 ; i < numSrcElems ; ++i)
  93       srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
  94
  95     for(ConnType i = 0 ; i < numTargetElems ; ++i)
  96       targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
  97
  98     Intersector3D<MyMeshType,MatrixType>* intersector=0;
  99     std::string methC = InterpolationOptions::filterInterpolationMethod(method);
 100     if(methC=="P0P0")
 101       {
 102         switch(InterpolationOptions::getIntersectionType())
 103           {
 104           case Triangulation:
 105             intersector=new PolyhedronIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
 106             break;
 107           case PointLocator:
 108             intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
 109             break;
 110           default:
 111             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangle or PointLocator.");
 112           }
 113       }
 114     else if(methC=="P0P1")
 115       {
 116         switch(InterpolationOptions::getIntersectionType())
 117           {
 118           case Triangulation:
 119             intersector=new PolyhedronIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
 120             break;
 121           case PointLocator:
 122             intersector=new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
 123             break;
 124           default:
 125             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P1 interp specified : must be Triangle or PointLocator.");
 126           }
 127       }
 128     else if(methC=="P1P0")
 129       {
 130         switch(InterpolationOptions::getIntersectionType())
 131           {
 132           case Triangulation:
 133             intersector=new PolyhedronIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
 134             break;
 135           case PointLocator:
 136             intersector=new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
 137             break;
 138           case Barycentric:
 139             intersector=new PolyhedronIntersectorP1P0Bary<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
 140             break;
 141           default:
 142             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P0 interp specified : must be Triangle, PointLocator or Barycentric.");
 143           }
 144       }
 145     else if(methC=="P1P1")
 146       {
 147         switch(InterpolationOptions::getIntersectionType())
 148           {
 149           case Triangulation:
 150             intersector=new PolyhedronIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
 151             break;
 152           case PointLocator:
 153             intersector=new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
 154             break;
 155           case Barycentric:
 156             intersector=new Barycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
 157             break;
 158           case MappedBarycentric:
 159             intersector=new MappedBarycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
 160             break;
 161           default:
 162             throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle, PointLocator, Barycentric or MappedBarycentric.");
 163           }
 164       }
 165     else
 166       throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
 167     // create empty maps for all source elements
 168     result.resize(intersector->getNumberOfRowsOfResMatrix());
 169
 170 #ifdef USE_RECURSIVE_BBOX_FILTER
 171
 172     /*
 173      * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh
 174      * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming
 175      * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the
 176      * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a
 177      * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements
 178      * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh.
 179      */
 180
 181     // create initial RegionNode and fill up its source region with all the source mesh elements and
 182     // its target region with all the target mesh elements whose bounding box
 183     // intersects that of the source region
 184
 185     RegionNode<ConnType>* firstNode = new RegionNode<ConnType>();
 186
 187     MeshRegion<ConnType>& srcRegion = firstNode->getSrcRegion();
 188
 189     for(ConnType i = 0 ; i < numSrcElems ; ++i)
 190       {
 191         srcRegion.addElement(srcElems[i], srcMesh);
 192       }
 193
 194     MeshRegion<ConnType>& targetRegion = firstNode->getTargetRegion();
 195
 196     for(ConnType i = 0 ; i < numTargetElems ; ++i)
 197       {
 198         if(!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]) ))
 199           {
 200             targetRegion.addElement(targetElems[i], targetMesh);
 201           }
 202       }
 203
 204     // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and
 205     // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the
 206     // elements of the target mesh whose bounding box intersects the corresponding part
 207     // Continue until the source region contains only one element, at which point the intersection volumes are
 208     // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero.
 209
 210     std::stack< RegionNode<ConnType>* > nodes;
 211     nodes.push(firstNode);
 212
 213     while(!nodes.empty())
 214       {
 215         RegionNode<ConnType>* currNode = nodes.top();
 216         nodes.pop();
 217         LOG(4, "Popping node ");
 218
 219         if(currNode->getTargetRegion().getNumberOfElements() == 1)
 220           {
 221             // calculate volumes
 222             LOG(4, " - One element");
 223
 224             MeshElement<ConnType>* targetElement = *(currNode->getTargetRegion().getBeginElements());
 225             std::vector<ConnType> intersectElems;
 226             for(typename std::vector< MeshElement<ConnType>* >::const_iterator iter = currNode->getSrcRegion().getBeginElements();iter != currNode->getSrcRegion().getEndElements();++iter)
 227               intersectElems.push_back((*iter)->getIndex());
 228             intersector->intersectCells(targetElement->getIndex(),intersectElems,result);
 229           }
 230         else // recursion
 231           {
 232
 233             LOG(4, " - Recursion");
 234
 235             RegionNode<ConnType>* leftNode = new RegionNode<ConnType>();
 236             RegionNode<ConnType>* rightNode = new RegionNode<ConnType>();
 237
 238             // split current source region
 239             //} decide on axis
 240             static BoundingBox::BoxCoord axis = BoundingBox::XMAX;
 241
 242             currNode->getTargetRegion().split(leftNode->getTargetRegion(), rightNode->getTargetRegion(), axis, targetMesh);
 243
 244             LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements()
 245                 << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements()
 246                 << " elements");
 247
 248             // ugly hack to avoid problem with enum which does not start at 0
 249             // I guess I ought to implement ++ for it instead ...
 250             // Anyway, it basically chooses the next axis, cyclically
 251             axis = (axis != BoundingBox::ZMAX) ? static_cast<BoundingBox::BoxCoord>(axis + 1) : BoundingBox::XMAX;
 252
 253             // add source elements of current node that overlap the target regions of the new nodes
 254             LOG(5, " -- Adding source elements");
 255             ConnType numLeftElements = 0;
 256             ConnType numRightElements = 0;
 257             for(typename std::vector<MeshElement<ConnType>*>::const_iterator iter = currNode->getSrcRegion().getBeginElements() ;
 258                 iter != currNode->getSrcRegion().getEndElements() ; ++iter)
 259               {
 260                 LOG(6, " --- New target node");
 261
 262                 if(!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter))
 263                   {
 264                     leftNode->getSrcRegion().addElement(*iter, srcMesh);
 265                     ++numLeftElements;
 266                   }
 267
 268                 if(!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter))
 269                   {
 270                     rightNode->getSrcRegion().addElement(*iter, srcMesh);
 271                     ++numRightElements;
 272                   }
 273
 274               }
 275
 276             LOG(5, "Left src region has " << numLeftElements << " elements and right src region has "
 277                 << numRightElements << " elements");
 278
 279             // push new nodes on stack
 280             if(numLeftElements != 0)
 281               {
 282                 nodes.push(leftNode);
 283               }
 284             else
 285               {
 286                 delete leftNode;
 287               }
 288
 289             if(numRightElements != 0)
 290               {
 291                 nodes.push(rightNode);
 292               }
 293             else
 294               {
 295                 delete rightNode;
 296               }
 297           }
 298
 299         // all nodes are deleted here
 300         delete currNode;
 301
 302         LOG(4, "Next iteration. Nodes left : " << nodes.size());
 303       }
 304
 305 #else // Use BBTree
 306
 307       // create BBTree structure
 308       // - get bounding boxes
 309     double* bboxes = new double[6 * numSrcElems];
 310     ConnType* srcElemIdx = new ConnType[numSrcElems];
 311     for(ConnType i = 0; i < numSrcElems ; ++i)
 312       {
 313         // get source bboxes in right order
 314         const BoundingBox* box = srcElems[i]->getBoundingBox();
 315         bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
 316         bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
 317         bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
 318         bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
 319         bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
 320         bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
 321
 322         // source indices have to begin with zero for BBox, I think
 323         srcElemIdx[i] = srcElems[i]->getIndex();
 324       }
 325
 326     BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems);
 327
 328     // for each target element, get source elements with which to calculate intersection
 329     // - calculate intersection by calling intersectCells
 330     for(ConnType i = 0; i < numTargetElems; ++i)
 331       {
 332         const BoundingBox* box = targetElems[i]->getBoundingBox();
 333         const ConnType targetIdx = targetElems[i]->getIndex();
 334
 335         // get target bbox in right order
 336         double targetBox[6];
 337         targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
 338         targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
 339         targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
 340         targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
 341         targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
 342         targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
 343
 344         std::vector<ConnType> intersectElems;
 345
 346         tree.getIntersectingElems(targetBox, intersectElems);
 347
 348         if ( !intersectElems.empty() )
 349           intersector->intersectCells(targetIdx,intersectElems,result);
 350       }
 351
 352     delete [] bboxes;
 353     delete [] srcElemIdx;
 354
 355 #endif
 356     // free allocated memory
 357     ConnType ret=intersector->getNumberOfColsOfResMatrix();
 358
 359     delete intersector;
 360
 361     for(ConnType i = 0 ; i < numSrcElems ; ++i)
 362       {
 363         delete srcElems[i];
 364       }
 365     for(ConnType i = 0 ; i < numTargetElems ; ++i)
 366       {
 367         delete targetElems[i];
 368       }
 369     return ret;
 370
 371   }
 372 }
 373
 374 #endif