src/INTERP_KERNEL/Interpolation2D1D.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 // Author : Anthony Geay (CEA/DEN)
  20 #ifndef __INTERPOLATION2D1D_TXX__
  21 #define __INTERPOLATION2D1D_TXX__
  22
  23 #include "Interpolation2D1D.hxx"
  24
  25 namespace INTERP_KERNEL
  26 {
  27
  28   /** \brief Main function to interpolate triangular or quadrangular meshes.
  29       \details  The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
  30       * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  31       * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
  32       * intersection matrix.
  33       *
  34       * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  35       * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  36       * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  37       * which have a non-zero intersection volume with the target element. The vector has indices running from
  38       * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  39       * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  40       *
  41
  42       * @param myMeshS  Planar source mesh
  43       * @Param myMeshT  Planar target mesh
  44       * @return            vector containing for each element i of the source mesh, a map giving for each element j
  45       *                    of the target mesh which i intersects, the area of the intersection
  46       *
  47       */
  48   template<class MyMeshType, class MatrixType>
  49   int Interpolation2D1D::interpolateMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, MatrixType& result, const std::string& method)
  50   {
  51     static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
  52     typedef typename MyMeshType::MyConnType ConnType;
  53     static const NumberingPolicy numPol=MyMeshType::My_numPol;
  54
  55     long global_start =clock();
  56     int counter=0;
  57     /***********************************************************/
  58     /* Check both meshes are made of triangles and quadrangles */
  59     /***********************************************************/
  60
  61     long nbMailleS=myMeshS.getNumberOfElements();
  62
  63     /**************************************************/
  64     /* Search the characteristic size of the meshes   */
  65     /**************************************************/
  66
  67     int printLevel = InterpolationOptions::getPrintLevel();
  68     _dim_caracteristic = CalculateCharacteristicSizeOfMeshes(myMeshS, myMeshT, printLevel);
  69     if (printLevel>=1)
  70       {
  71         std::cout << "Interpolation2D1D::computation of the intersections" << std::endl;
  72       }
  73
  74     PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
  75     std::string meth = InterpolationOptions::filterInterpolationMethod(method);
  76     if(meth=="P0P0")
  77       {
  78         switch (InterpolationOptions::getIntersectionType())
  79           {
  80           case Geometric2D:
  81             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,Planar2D1DIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
  82                                                                                                     InterpolationOptions::getMaxDistance3DSurfIntersect(),
  83                                                                                                     InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
  84                                                                                                     InterpolationOptions::getMedianPlane(),
  85                                                                                                     InterpolationOptions::getPrecision(),
  86                                                                                                     InterpolationOptions::getOrientation());
  87             break;
  88           default:
  89             throw INTERP_KERNEL::Exception("Invalid intersection type ! Must be : Geometric2D");
  90           }
  91       }
  92     else
  93       throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be : \"P0P0\"");
  94
  95     /****************************************************************/
  96     /* Create a search tree based on the bounding boxes             */
  97     /* Instanciate the intersector and initialise the result vector */
  98     /****************************************************************/
  99
 100     long start_filtering=clock();
 101
 102     std::vector<double> bbox;
 103     intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
 104     const double *bboxPtr=0;
 105     if(nbMailleS>0)
 106       bboxPtr=&bbox[0];
 107     BBTree<SPACEDIM,ConnType> my_tree(bboxPtr, 0, 0,nbMailleS, -getPrecision());//creating the search structure
 108
 109     long end_filtering=clock();
 110
 111     result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
 112
 113     /****************************************************/
 114     /* Loop on the target cells - core of the algorithm */
 115     /****************************************************/
 116     long start_intersection=clock();
 117     long nbelem_type=myMeshT.getNumberOfElements();
 118     const ConnType *connIndxT=myMeshT.getConnectivityIndexPtr();
 119     for(int iT=0; iT<nbelem_type; iT++)
 120       {
 121         int nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
 122         std::vector<int> intersecting_elems;
 123         double bb[2*SPACEDIM];
 124         intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
 125         my_tree.getIntersectingElems(bb, intersecting_elems);
 126         intersector->intersectCells(iT,intersecting_elems,result);
 127         counter+=intersecting_elems.size();
 128         intersecting_elems.clear();
 129       }
 130     int ret=intersector->getNumberOfColsOfResMatrix();
 131
 132     const DuplicateFacesType& intersectFaces = *intersector->getIntersectFaces();
 133     DuplicateFacesType::const_iterator iter;
 134     for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter)
 135       {
 136         if (iter->second.size() > 1)
 137           {
 138             _duplicate_faces.insert(std::make_pair(iter->first, iter->second));
 139           }
 140       }
 141
 142     delete intersector;
 143
 144     if (InterpolationOptions::getPrintLevel() >=1)
 145       {
 146         long end_intersection=clock();
 147         std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
 148         std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
 149         long global_end =clock();
 150         std::cout << "Number of computed intersections = " << counter << std::endl;
 151         std::cout << "Global time= " << global_end - global_start << std::endl;
 152       }
 153     return ret;
 154   }
 155
 156 }
 157 #endif