src/INTERP_KERNEL/Interpolation3D.txx

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