src/INTERP_KERNEL/InterpolationCurve.txx

   1 // Copyright (C) 2007-2012  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 __INTERPOLATIONCURVE_TXX__
  20 #define __INTERPOLATIONCURVE_TXX__
  21
  22 #include "InterpolationCurve.hxx"
  23 #include "InterpolationOptions.hxx"
  24 #include "CurveIntersectorP0P0.txx"
  25 #include "CurveIntersectorP1P0.txx"
  26 #include "CurveIntersectorP0P1.txx"
  27 #include "CurveIntersectorP1P1.txx"
  28 #include "BBTree.txx"
  29
  30 #include <time.h>
  31
  32 namespace INTERP_KERNEL
  33 {
  34   /**
  35    * \defgroup interpolationCurve InterpolationCurve
  36    *
  37    * \class InterpolationCurve
  38    * \brief Class used to compute the coefficients of the interpolation matrix between
  39    * two local meshes in two dimensions.
  40    */
  41   template<class RealCurve>
  42   InterpolationCurve<RealCurve>::InterpolationCurve()
  43   {
  44   }
  45
  46   template<class RealCurve>
  47   InterpolationCurve<RealCurve>::InterpolationCurve (const InterpolationOptions& io)
  48     :Interpolation< InterpolationCurve<RealCurve> >(io)
  49   {
  50   }
  51
  52   /** \brief Main function to interpolate 1D meshes.
  53       \details  The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
  54       * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  55       * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
  56       * intersection matrix.
  57       *
  58       * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  59       * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  60       * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  61       * which have a non-zero intersection volume with the target element. The vector has indices running from
  62       * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  63       * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  64       *
  65
  66       * @param myMeshS  Planar source mesh
  67       * @Param myMeshT  Planar target mesh
  68       * @return vector containing for each element i of the source mesh, a map giving for each element j
  69       *         of the target mesh which i intersects, the area of the intersection
  70       *
  71       */
  72   template<class RealCurve>
  73   template<class MyMeshType, class MatrixType>
  74   int InterpolationCurve<RealCurve>::interpolateMeshes (const MyMeshType& myMeshS,
  75                                                         const MyMeshType& myMeshT,
  76                                                         MatrixType&       result,
  77                                                         const char *      method)
  78   {
  79     static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
  80     typedef typename MyMeshType::MyConnType ConnType;
  81     static const NumberingPolicy numPol = MyMeshType::My_numPol;
  82
  83     long global_start = clock();
  84     int counter=0;
  85
  86     long nbMailleS = myMeshS.getNumberOfElements();
  87     long nbMailleT = myMeshT.getNumberOfElements();
  88
  89     CurveIntersector<MyMeshType,MatrixType>* intersector=0;
  90     std::string meth(method);
  91     if(meth=="P0P0")
  92       {
  93         intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>
  94           (myMeshT, myMeshS,
  95            InterpolationOptions::getPrecision(),
  96            InterpolationOptions::getBoundingBoxAdjustmentAbs(),
  97            InterpolationOptions::getMedianPlane(),
  98            InterpolationOptions::getPrintLevel());
  99       }
 100     else if(meth=="P0P1")
 101       {
 102         intersector = new CurveIntersectorP0P1<MyMeshType,MatrixType>
 103           (myMeshT, myMeshS,
 104            InterpolationOptions::getPrecision(),
 105            InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 106            InterpolationOptions::getMedianPlane(),
 107            InterpolationOptions::getPrintLevel());
 108       }
 109     else if(meth=="P1P0")
 110       {
 111         intersector = new CurveIntersectorP1P0<MyMeshType,MatrixType>
 112           (myMeshT, myMeshS,
 113            InterpolationOptions::getPrecision(),
 114            InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 115            InterpolationOptions::getMedianPlane(),
 116            InterpolationOptions::getPrintLevel());
 117       }
 118     else if(meth=="P1P1")
 119       {
 120         intersector = new CurveIntersectorP1P1<MyMeshType,MatrixType>
 121           (myMeshT, myMeshS,
 122            InterpolationOptions::getPrecision(),
 123            InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 124            InterpolationOptions::getMedianPlane(),
 125            InterpolationOptions::getPrintLevel());
 126       }
 127     else
 128       throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
 129     /****************************************************************/
 130     /* Create a search tree based on the bounding boxes             */
 131     /* Instanciate the intersector and initialise the result vector */
 132     /****************************************************************/
 133
 134     long start_filtering=clock();
 135
 136     std::vector<double> bbox;
 137     intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
 138     intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustmentAbs());
 139     BBTree<SPACEDIM,ConnType> my_tree(&bbox[0], 0, 0, nbMailleS);//creating the search structure
 140
 141     long end_filtering = clock();
 142
 143     result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
 144
 145     /****************************************************/
 146     /* Loop on the target cells - core of the algorithm */
 147     /****************************************************/
 148     long start_intersection = clock();
 149     const ConnType *connIndxT = myMeshT.getConnectivityIndexPtr();
 150     for(int iT=0; iT<nbMailleT; iT++)
 151       {
 152         int nb_nodesT = connIndxT[iT+1] - connIndxT[iT];
 153         std::vector<int> intersecting_elems;
 154         double bb[2*SPACEDIM];
 155         intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
 156         my_tree.getIntersectingElems(bb, intersecting_elems);
 157         intersector->intersectCells(iT,intersecting_elems,result);
 158         counter += intersecting_elems.size();
 159       }
 160     int ret = intersector->getNumberOfColsOfResMatrix();
 161     delete intersector;
 162
 163     if (InterpolationOptions::getPrintLevel() >= 1)
 164       {
 165         long end_intersection=clock();
 166         std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
 167         std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
 168         long global_end =clock();
 169         std::cout << "Number of computed intersections = " << counter << std::endl;
 170         std::cout << "Global time= " << global_end - global_start << std::endl;
 171       }
 172     return ret;
 173   }
 174 }
 175
 176 #endif