src/medtool/src/INTERP_KERNEL/Interpolation3D2D.txx

   1 // Copyright (C) 2007-2015  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 #ifndef __INTERPOLATION3D2D_TXX__
  20 #define __INTERPOLATION3D2D_TXX__
  21
  22 #include "Interpolation3D2D.hxx"
  23 #include "Interpolation.txx"
  24 #include "MeshElement.txx"
  25 #include "TransformedTriangle.hxx"
  26 #include "Polyhedron3D2DIntersectorP0P0.txx"
  27 #include "PointLocator3DIntersectorP0P0.txx"
  28 #include "PolyhedronIntersectorP0P1.txx"
  29 #include "PointLocator3DIntersectorP0P1.txx"
  30 #include "PolyhedronIntersectorP1P0.txx"
  31 #include "PolyhedronIntersectorP1P0Bary.txx"
  32 #include "PointLocator3DIntersectorP1P0.txx"
  33 #include "PolyhedronIntersectorP1P1.txx"
  34 #include "PointLocator3DIntersectorP1P1.txx"
  35 #include "Log.hxx"
  36
  37 #include "BBTree.txx"
  38
  39 namespace INTERP_KERNEL
  40 {
  41   /**
  42    * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
  43    * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
  44    * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  45    * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the
  46    * intersection matrix.
  47    *
  48    * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  49    * It can also be an INTERP_KERNEL::Matrix object.
  50    * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  51    * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  52    * which have a non-zero intersection volume with the target element. The vector has indices running from
  53    * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  54    * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  55    *
  56
  57    * @param srcMesh     3-dimensional source mesh
  58    * @param targetMesh  3-dimesional target mesh, containing only tetraedra
  59    * @param matrix      matrix in which the result is stored
  60    *
  61    */
  62   template<class MyMeshType, class MyMatrixType>
  63   int Interpolation3D2D::interpolateMeshes(const MyMeshType& srcMesh,
  64                                            const MyMeshType& targetMesh,
  65                                            MyMatrixType& matrix,
  66                                            const std::string& method)
  67   {
  68     typedef typename MyMeshType::MyConnType ConnType;
  69     // create MeshElement objects corresponding to each element of the two meshes
  70     const unsigned long numSrcElems = srcMesh.getNumberOfElements();
  71     const unsigned long numTargetElems = targetMesh.getNumberOfElements();
  72
  73     LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
  74
  75     std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
  76     std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
  77
  78     std::map<MeshElement<ConnType>*, int> indices;
  79     DuplicateFacesType intersectFaces;
  80
  81     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
  82       srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
  83
  84     for(unsigned long i = 0 ; i < numTargetElems ; ++i)
  85       targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
  86
  87     Intersector3D<MyMeshType,MyMatrixType>* intersector=0;
  88     std::string methC = InterpolationOptions::filterInterpolationMethod(method);
  89     const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel());
  90     if(methC=="P0P0")
  91       {
  92         switch(InterpolationOptions::getIntersectionType())
  93           {
  94           case Triangulation:
  95             intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
  96                                                                                    srcMesh,
  97                                                                                    dimCaracteristic,
  98                                                                                    getPrecision(),
  99                                                                                    intersectFaces,
 100                                                                                    getSplittingPolicy());
 101             break;
 102           case PointLocator:
 103             intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,srcMesh,getPrecision());
 104             break;
 105           default:
 106             throw INTERP_KERNEL::Exception("Invalid 3D to 2D intersection type for P0P0 interp specified : must be Triangulation or PointLocator.");
 107           }
 108       }
 109     else
 110       throw Exception("Invalid method choosed must be in \"P0P0\".");
 111     // create empty maps for all source elements
 112     matrix.resize(intersector->getNumberOfRowsOfResMatrix());
 113
 114     // create BBTree structure
 115     // - get bounding boxes
 116     double* bboxes = new double[6 * numSrcElems];
 117     int* srcElemIdx = new int[numSrcElems];
 118     for(unsigned long i = 0; i < numSrcElems ; ++i)
 119       {
 120         // get source bboxes in right order
 121         const BoundingBox* box = srcElems[i]->getBoundingBox();
 122         bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
 123         bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
 124         bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
 125         bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
 126         bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
 127         bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
 128
 129         // source indices have to begin with zero for BBox, I think
 130         srcElemIdx[i] = srcElems[i]->getIndex();
 131       }
 132
 133     BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems, 0.);
 134
 135     // for each target element, get source elements with which to calculate intersection
 136     // - calculate intersection by calling intersectCells
 137     for(unsigned long i = 0; i < numTargetElems; ++i)
 138       {
 139         const BoundingBox* box = targetElems[i]->getBoundingBox();
 140         const int targetIdx = targetElems[i]->getIndex();
 141
 142         // get target bbox in right order
 143         double targetBox[6];
 144         targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
 145         targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
 146         targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
 147         targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
 148         targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
 149         targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
 150
 151         std::vector<ConnType> intersectElems;
 152
 153         tree.getIntersectingElems(targetBox, intersectElems);
 154
 155         if ( !intersectElems.empty() )
 156             intersector->intersectCells(targetIdx, intersectElems, matrix);
 157
 158       }
 159
 160     delete [] bboxes;
 161     delete [] srcElemIdx;
 162
 163     DuplicateFacesType::iterator iter;
 164     for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter)
 165       {
 166         if (iter->second.size() > 1)
 167           {
 168             _duplicate_faces.insert(std::make_pair(iter->first, iter->second));
 169           }
 170       }
 171
 172     // free allocated memory
 173     int ret=intersector->getNumberOfColsOfResMatrix();
 174
 175     delete intersector;
 176
 177     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
 178       {
 179         delete srcElems[i];
 180       }
 181     for(unsigned long i = 0 ; i < numTargetElems ; ++i)
 182       {
 183         delete targetElems[i];
 184       }
 185     return ret;
 186
 187   }
 188 }
 189
 190 #endif