src/INTERP_KERNEL/InterpolationCurve.txx

   1 // Copyright (C) 2007-2020  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 __INTERPOLATIONCURVE_TXX__
  21 #define __INTERPOLATIONCURVE_TXX__
  22
  23 #include "InterpolationCurve.hxx"
  24 #include "InterpolationOptions.hxx"
  25 #include "CurveIntersectorP0P0.txx"
  26 #include "CurveIntersectorP1P0.txx"
  27 #include "CurveIntersectorP0P1.txx"
  28 #include "CurveIntersectorP1P1.txx"
  29 #include "CurveIntersectorP1P1PL.txx"
  30
  31 #include <time.h>
  32 #include <memory>
  33
  34 namespace INTERP_KERNEL
  35 {
  36   /**
  37    * \defgroup interpolationCurve InterpolationCurve
  38    *
  39    * \class InterpolationCurve
  40    * \brief Class used to compute the coefficients of the interpolation matrix between
  41    * two local meshes in two dimensions.
  42    */
  43   template<class RealCurve>
  44   InterpolationCurve<RealCurve>::InterpolationCurve()
  45   {
  46   }
  47
  48   template<class RealCurve>
  49   InterpolationCurve<RealCurve>::InterpolationCurve (const InterpolationOptions& io)
  50     :Interpolation< InterpolationCurve<RealCurve> >(io)
  51   {
  52   }
  53
  54   /** \brief Main function to interpolate 1D meshes.
  55       \details  The algorithm proceeds 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 IntersectorPlanar 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       * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  62       * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  63       * which have a non-zero intersection volume with the target element. The vector has indices running from
  64       * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  65       * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  66       *
  67
  68       * @param myMeshS  Planar source mesh
  69       * @Param myMeshT  Planar target mesh
  70       * @return vector containing for each element i of the source mesh, a map giving for each element j
  71       *         of the target mesh which i intersects, the area of the intersection
  72       *
  73       */
  74   template<class RealCurve>
  75   template<class MyMeshType, class MatrixType>
  76   typename MyMeshType::MyConnType InterpolationCurve<RealCurve>::interpolateMeshesInternal (const MyMeshType&  myMeshS,
  77                                                                                             const MyMeshType&  myMeshT,
  78                                                                                             MatrixType&        result,
  79                                                                                             const std::string& method,
  80                                                                                             std::function< void(const BBTree< MyMeshType::MY_SPACEDIM ,
  81                                                                                                                 typename MyMeshType::MyConnType>&, const double*,
  82                                                                                                                 std::vector<typename MyMeshType::MyConnType>&) > bbtreeMethod
  83                                                                                             )
  84   {
  85     static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
  86     typedef typename MyMeshType::MyConnType ConnType;
  87     static const NumberingPolicy numPol = MyMeshType::My_numPol;
  88
  89     long global_start = clock();
  90     std::size_t counter=0;
  91
  92     ConnType nbMailleS = myMeshS.getNumberOfElements();
  93     ConnType nbMailleT = myMeshT.getNumberOfElements();
  94
  95     std::unique_ptr< CurveIntersector<MyMeshType,MatrixType> > intersector;
  96     if(method=="P0P0")
  97       {
  98         switch (InterpolationOptions::getIntersectionType())
  99           {
 100             case Triangulation:
 101               {
 102                 intersector.reset( new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
 103                                                                                    InterpolationOptions::getPrecision(),
 104                                                                                    InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 105                                                                                    InterpolationOptions::getMedianPlane(),
 106                                                                                    InterpolationOptions::getPrintLevel()) );
 107                 break;
 108               }
 109             default:
 110               throw INTERP_KERNEL::Exception("For P0P0 in 1D or 2D curve only Triangulation supported for the moment !");
 111           }
 112       }
 113     else if(method=="P0P1")
 114       {
 115         switch (InterpolationOptions::getIntersectionType())
 116           {
 117             case Triangulation:
 118               {
 119                 intersector.reset( new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
 120                                                                                    InterpolationOptions::getPrecision(),
 121                                                                                    InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 122                                                                                    InterpolationOptions::getMedianPlane(),
 123                                                                                    InterpolationOptions::getPrintLevel()) );
 124                 break;
 125               }
 126             default:
 127               throw INTERP_KERNEL::Exception("For P0P1 in 1D or 2D curve only Triangulation supported for the moment !");
 128           }
 129       }
 130     else if(method=="P1P0")
 131       {
 132         switch (InterpolationOptions::getIntersectionType())
 133           {
 134             case Triangulation:
 135               {
 136                 intersector.reset( new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
 137                                                                                    InterpolationOptions::getPrecision(),
 138                                                                                    InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 139                                                                                    InterpolationOptions::getMedianPlane(),
 140                                                                                    InterpolationOptions::getPrintLevel()) );
 141                 break;
 142               }
 143           default:
 144             throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
 145           }
 146       }
 147     else if(method=="P1P1")
 148       {
 149         switch (InterpolationOptions::getIntersectionType())
 150           {
 151           case Triangulation:
 152             intersector.reset( new CurveIntersectorP1P1<MyMeshType,MatrixType>
 153                                (myMeshT, myMeshS,
 154                                 InterpolationOptions::getPrecision(),
 155                                 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 156                                 InterpolationOptions::getMedianPlane(),
 157                                 InterpolationOptions::getPrintLevel()) );
 158             break;
 159           case PointLocator:
 160             intersector.reset( new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
 161                                (myMeshT, myMeshS,
 162                                 InterpolationOptions::getPrecision(),
 163                                 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
 164                                 InterpolationOptions::getMedianPlane(),
 165                                 InterpolationOptions::getPrintLevel()) );
 166             break;
 167           default:
 168             throw INTERP_KERNEL::Exception("For P1P1 in 1D or 2D curve only Triangulation and PointLocator supported !");
 169           }
 170       }
 171     else
 172       throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
 173     /****************************************************************/
 174     /* Create a search tree based on the bounding boxes             */
 175     /* Instantiate the intersector and initialise the result vector */
 176     /****************************************************************/
 177
 178     long start_filtering=clock();
 179
 180     std::vector<double> bbox;
 181     intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
 182     intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustmentAbs());
 183     BBTree<SPACEDIM,ConnType> my_tree(&bbox[0], 0, 0, nbMailleS);//creating the search structure
 184
 185     long end_filtering = clock();
 186
 187     result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
 188
 189     /****************************************************/
 190     /* Loop on the target cells - core of the algorithm */
 191     /****************************************************/
 192     long start_intersection = clock();
 193     const ConnType *connIndxT = myMeshT.getConnectivityIndexPtr();
 194     for(ConnType iT=0; iT<nbMailleT; iT++)
 195       {
 196         ConnType nb_nodesT = connIndxT[iT+1] - connIndxT[iT];
 197         std::vector<ConnType> intersecting_elems;
 198         double bb[2*SPACEDIM];
 199         intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
 200         bbtreeMethod(my_tree,bb,intersecting_elems);
 201         intersector->intersectCells(iT,intersecting_elems,result);
 202         counter += intersecting_elems.size();
 203       }
 204
 205     if (InterpolationOptions::getPrintLevel() >= 1)
 206       {
 207         long end_intersection=clock();
 208         std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
 209         std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
 210         long global_end =clock();
 211         std::cout << "Number of computed intersections = " << counter << std::endl;
 212         std::cout << "Global time= " << global_end - global_start << std::endl;
 213       }
 214     return intersector->getNumberOfColsOfResMatrix();
 215   }
 216 }
 217
 218 #endif