src/INTERP_KERNEL/Interpolation3D.txx

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