src/INTERP_KERNEL/InterpolationPlanar.txx

   1 // Copyright (C) 2007-2013  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 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __INTERPOLATIONPLANAR_TXX__
  22 #define __INTERPOLATIONPLANAR_TXX__
  23
  24 #include "InterpolationPlanar.hxx"
  25 #include "Interpolation.txx"
  26 #include "InterpolationOptions.hxx"
  27 #include "PlanarIntersector.hxx"
  28 #include "PlanarIntersector.txx"
  29 #include "TriangulationIntersector.hxx"
  30 #include "TriangulationIntersector.txx"
  31 #include "ConvexIntersector.hxx"
  32 #include "ConvexIntersector.txx"
  33 #include "Geometric2DIntersector.hxx"
  34 #include "Geometric2DIntersector.txx"
  35 #include "PointLocator2DIntersector.hxx"
  36 #include "PointLocator2DIntersector.txx"
  37 #include "PlanarIntersectorP0P1PL.hxx"
  38 #include "PlanarIntersectorP0P1PL.txx"
  39 #include "PlanarIntersectorP1P0PL.hxx"
  40 #include "PlanarIntersectorP1P0PL.txx"
  41 #include "PlanarIntersectorP1P1PL.hxx"
  42 #include "PlanarIntersectorP1P1PL.txx"
  43 #include "VectorUtils.hxx"
  44 #include "BBTree.txx"
  45
  46 #include <limits>
  47 #include <time.h>
  48
  49 namespace INTERP_KERNEL
  50 {
  51   /**
  52    * \defgroup interpolationPlanar InterpolationPlanar
  53    *
  54    * \class InterpolationPlanar
  55    * \brief Class used to compute the coefficients of the interpolation matrix between
  56    * two local meshes in two dimensions. Meshes can contain mixed triangular and quadrangular elements.
  57    */
  58   template<class RealPlanar>
  59   InterpolationPlanar<RealPlanar>::InterpolationPlanar():_dim_caracteristic(1)
  60
  61   {
  62   }
  63
  64   template<class RealPlanar>
  65   InterpolationPlanar<RealPlanar>::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar<RealPlanar> >(io),_dim_caracteristic(1)
  66
  67   {
  68   }
  69
  70   /**
  71    *  \brief  Function used to set the options for the intersection calculation
  72    * \details The following options can be modified:
  73    *  -# Intersection_type: the type of algorithm to be used in the computation of the cell-cell intersections.
  74    *   - Values: Triangle, Convex.
  75    *   - Default: Triangle.
  76    *  -# Precision: Level of precision of the computations is precision times the characteristic size of the mesh.
  77    *   - Values: positive real number.
  78    *   - Default: 1.0E-12.
  79    *  -# PrintLevel: Level of verboseness during the computations.
  80    *   - Values: interger between 0 and 3.
  81    *   - Default: 0.
  82    */
  83   template<class RealPlanar>
  84   void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation)
  85   {
  86     InterpolationOptions::setPrecision(precision);
  87     InterpolationOptions::setPrintLevel(printLevel);
  88     InterpolationOptions::setIntersectionType(intersectionType);
  89     InterpolationOptions::setOrientation(orientation);
  90   }
  91
  92
  93   /** \brief Main function to interpolate triangular or quadrangular meshes.
  94       \details  The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
  95       * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  96       * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
  97       * intersection matrix.
  98       *
  99       * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
 100       * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
 101       * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
 102       * which have a non-zero intersection volume with the target element. The vector has indices running from
 103       * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
 104       * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
 105       *
 106
 107       * @param myMeshS  Planar source mesh
 108       * @Param myMeshT  Planar target mesh
 109       * @return            vector containing for each element i of the source mesh, a map giving for each element j
 110       *                    of the target mesh which i intersects, the area of the intersection
 111       *
 112       */
 113   template<class RealPlanar>
 114   template<class MyMeshType, class MatrixType>
 115   int InterpolationPlanar<RealPlanar>::interpolateMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, MatrixType& result, const std::string& method)
 116   {
 117     static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 118     typedef typename MyMeshType::MyConnType ConnType;
 119     static const NumberingPolicy numPol=MyMeshType::My_numPol;
 120
 121     long global_start =clock();
 122     int counter=0;
 123     /***********************************************************/
 124     /* Check both meshes are made of triangles and quadrangles */
 125     /***********************************************************/
 126
 127     long nbMailleS=myMeshS.getNumberOfElements();
 128     long nbMailleT=myMeshT.getNumberOfElements();
 129
 130     /**************************************************/
 131     /* Search the characteristic size of the meshes   */
 132     /**************************************************/
 133
 134     double BoxS[2*SPACEDIM]; myMeshS.getBoundingBox(BoxS);
 135     double BoxT[2*SPACEDIM]; myMeshT.getBoundingBox(BoxT);
 136     double diagonalS,dimCaracteristicS=std::numeric_limits<double>::max();
 137     if(nbMailleS!=0)
 138       {
 139         diagonalS=getDistanceBtw2Pts<SPACEDIM>(BoxS+SPACEDIM,BoxS);
 140         dimCaracteristicS=diagonalS/nbMailleS;
 141       }
 142     double diagonalT,dimCaracteristicT=std::numeric_limits<double>::max();
 143     if(nbMailleT!=0)
 144       {
 145         diagonalT=getDistanceBtw2Pts<SPACEDIM>(BoxT+SPACEDIM,BoxT);
 146         dimCaracteristicT=diagonalT/nbMailleT;
 147       }
 148
 149     _dim_caracteristic=std::min(dimCaracteristicS, dimCaracteristicT);
 150     if (InterpolationOptions::getPrintLevel()>=1)
 151       {
 152         std::cout << "  - Characteristic size of the source mesh : " << dimCaracteristicS << std::endl;
 153         std::cout << "  - Characteristic size of the target mesh: " << dimCaracteristicT << std::endl;
 154         std::cout << "InterpolationPlanar::computation of the intersections" << std::endl;
 155       }
 156
 157     PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
 158     std::string meth = InterpolationOptions::filterInterpolationMethod(method);
 159     if(meth=="P0P0")
 160       {
 161         switch (InterpolationOptions::getIntersectionType())
 162           {
 163           case Triangulation:
 164             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
 165                                                                                                   InterpolationOptions::getPrecision(),
 166                                                                                                   InterpolationOptions::getMaxDistance3DSurfIntersect(),
 167                                                                                                   InterpolationOptions::getMedianPlane(),
 168                                                                                                   InterpolationOptions::getOrientation(),
 169                                                                                                   InterpolationOptions::getPrintLevel());
 170             break;
 171           case Convex:
 172             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
 173                                                                                            InterpolationOptions::getPrecision(),
 174                                                                                            InterpolationOptions::getMaxDistance3DSurfIntersect(),
 175                                                                                            InterpolationOptions::getMedianPlane(),
 176                                                                                            InterpolationOptions::getDoRotate(),
 177                                                                                            InterpolationOptions::getOrientation(),
 178                                                                                            InterpolationOptions::getPrintLevel());
 179             break;
 180           case Geometric2D:
 181             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
 182                                                                                                 InterpolationOptions::getMaxDistance3DSurfIntersect(),
 183                                                                                                 InterpolationOptions::getMedianPlane(),
 184                                                                                                 InterpolationOptions::getPrecision(),
 185                                                                                                 InterpolationOptions::getOrientation());
 186             break;
 187           case PointLocator:
 188             intersector=new PointLocator2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
 189                                                                                                    InterpolationOptions::getMaxDistance3DSurfIntersect(),
 190                                                                                                    InterpolationOptions::getMedianPlane(),
 191                                                                                                    InterpolationOptions::getPrecision(),
 192                                                                                                    InterpolationOptions::getOrientation());
 193             break;
 194           default:
 195             throw INTERP_KERNEL::Exception("For P0P0 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator !");
 196           }
 197       }
 198     else if(meth=="P0P1")
 199       {
 200         switch (InterpolationOptions::getIntersectionType())
 201           {
 202           case Triangulation:
 203             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
 204                                                                                                   InterpolationOptions::getPrecision(),
 205                                                                                                   InterpolationOptions::getMaxDistance3DSurfIntersect(),
 206                                                                                                   InterpolationOptions::getMedianPlane(),
 207                                                                                                   InterpolationOptions::getOrientation(),
 208                                                                                                   InterpolationOptions::getPrintLevel());
 209             break;
 210           case Convex:
 211             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
 212                                                                                            InterpolationOptions::getPrecision(),
 213                                                                                            InterpolationOptions::getMaxDistance3DSurfIntersect(),
 214                                                                                            InterpolationOptions::getMedianPlane(),
 215                                                                                            InterpolationOptions::getDoRotate(),
 216                                                                                            InterpolationOptions::getOrientation(),
 217                                                                                            InterpolationOptions::getPrintLevel());
 218             break;
 219           case Geometric2D:
 220             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT, myMeshS, _dim_caracteristic,
 221                                                                                                 InterpolationOptions::getMaxDistance3DSurfIntersect(),
 222                                                                                                 InterpolationOptions::getMedianPlane(),
 223                                                                                                 InterpolationOptions::getPrecision(),
 224                                                                                                 InterpolationOptions::getOrientation());
 225             break;
 226           case PointLocator:
 227             intersector=new PlanarIntersectorP0P1PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
 228                                                                            InterpolationOptions::getMaxDistance3DSurfIntersect(),
 229                                                                            InterpolationOptions::getMedianPlane(),
 230                                                                            InterpolationOptions::getPrecision(),
 231                                                                            InterpolationOptions::getOrientation());
 232             break;
 233           case Barycentric:
 234             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1Bary>(myMeshT,myMeshS,_dim_caracteristic,
 235                                                                                                       InterpolationOptions::getPrecision(),
 236                                                                                                       InterpolationOptions::getMaxDistance3DSurfIntersect(),
 237                                                                                                       InterpolationOptions::getMedianPlane(),
 238                                                                                                       InterpolationOptions::getOrientation(),
 239                                                                                                       InterpolationOptions::getPrintLevel());
 240             break;
 241           case BarycentricGeo2D:
 242             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1Bary>(myMeshT, myMeshS, _dim_caracteristic,
 243                                                                                                     InterpolationOptions::getMaxDistance3DSurfIntersect(),
 244                                                                                                     InterpolationOptions::getMedianPlane(),
 245                                                                                                     InterpolationOptions::getPrecision(),
 246                                                                                                     InterpolationOptions::getOrientation());
 247             break;
 248           default:
 249             throw INTERP_KERNEL::Exception("For P0P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator, Barycentric, BarycentricGeo2D !");
 250           }
 251       }
 252     else if(meth=="P1P0")
 253       {
 254         switch (InterpolationOptions::getIntersectionType())
 255           {
 256           case Triangulation:
 257             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
 258                                                                                                   InterpolationOptions::getPrecision(),
 259                                                                                                   InterpolationOptions::getMaxDistance3DSurfIntersect(),
 260                                                                                                   InterpolationOptions::getMedianPlane(),
 261                                                                                                   InterpolationOptions::getOrientation(),
 262                                                                                                   InterpolationOptions::getPrintLevel());
 263             break;
 264           case Convex:
 265             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
 266                                                                                            InterpolationOptions::getPrecision(),
 267                                                                                            InterpolationOptions::getMaxDistance3DSurfIntersect(),
 268                                                                                            InterpolationOptions::getMedianPlane(),
 269                                                                                            InterpolationOptions::getDoRotate(),
 270                                                                                            InterpolationOptions::getOrientation(),
 271                                                                                            InterpolationOptions::getPrintLevel());
 272             break;
 273           case Geometric2D:
 274             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT, myMeshS, _dim_caracteristic,
 275                                                                                                 InterpolationOptions::getMaxDistance3DSurfIntersect(),
 276                                                                                                 InterpolationOptions::getMedianPlane(),
 277                                                                                                 InterpolationOptions::getPrecision(),
 278                                                                                                 InterpolationOptions::getOrientation());
 279             break;
 280           case PointLocator:
 281             intersector=new PlanarIntersectorP1P0PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
 282                                                                        InterpolationOptions::getMaxDistance3DSurfIntersect(),
 283                                                                        InterpolationOptions::getMedianPlane(),
 284                                                                        InterpolationOptions::getPrecision(),
 285                                                                        InterpolationOptions::getOrientation());
 286             break;
 287           case Barycentric:
 288              intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0Bary>(myMeshT,myMeshS,_dim_caracteristic,
 289                                                                                                        InterpolationOptions::getPrecision(),
 290                                                                                                        InterpolationOptions::getMaxDistance3DSurfIntersect(),
 291                                                                                                        InterpolationOptions::getMedianPlane(),
 292                                                                                                        InterpolationOptions::getOrientation(),
 293                                                                                                        InterpolationOptions::getPrintLevel());
 294              break;
 295           case BarycentricGeo2D:
 296             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0Bary>(myMeshT, myMeshS, _dim_caracteristic,
 297                                                                                                     InterpolationOptions::getMaxDistance3DSurfIntersect(),
 298                                                                                                     InterpolationOptions::getMedianPlane(),
 299                                                                                                     InterpolationOptions::getPrecision(),
 300                                                                                                     InterpolationOptions::getOrientation());
 301             break;
 302           }
 303       }
 304     else if(meth=="P1P1")
 305       {
 306         switch (InterpolationOptions::getIntersectionType())
 307           {
 308           case Triangulation:
 309             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
 310                                                                                                   InterpolationOptions::getPrecision(),
 311                                                                                                   InterpolationOptions::getMaxDistance3DSurfIntersect(),
 312                                                                                                   InterpolationOptions::getMedianPlane(),
 313                                                                                                   InterpolationOptions::getOrientation(),
 314                                                                                                   InterpolationOptions::getPrintLevel());
 315             break;
 316           case Convex:
 317             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
 318                                                                                            InterpolationOptions::getPrecision(),
 319                                                                                            InterpolationOptions::getMaxDistance3DSurfIntersect(),
 320                                                                                            InterpolationOptions::getMedianPlane(),
 321                                                                                            InterpolationOptions::getDoRotate(),
 322                                                                                            InterpolationOptions::getOrientation(),
 323                                                                                            InterpolationOptions::getPrintLevel());
 324             break;
 325           case Geometric2D:
 326             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT, myMeshS, _dim_caracteristic,
 327                                                                                                 InterpolationOptions::getMaxDistance3DSurfIntersect(),
 328                                                                                                 InterpolationOptions::getMedianPlane(),
 329                                                                                                 InterpolationOptions::getPrecision(),
 330                                                                                                 InterpolationOptions::getOrientation());
 331             break;
 332           case PointLocator:
 333             intersector=new PlanarIntersectorP1P1PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
 334                                                                            InterpolationOptions::getMaxDistance3DSurfIntersect(),
 335                                                                            InterpolationOptions::getMedianPlane(),
 336                                                                            InterpolationOptions::getPrecision(),
 337                                                                            InterpolationOptions::getOrientation());
 338             break;
 339           default:
 340             throw INTERP_KERNEL::Exception("For P1P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator !");
 341           }
 342       }
 343     else
 344       throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
 345     /****************************************************************/
 346     /* Create a search tree based on the bounding boxes             */
 347     /* Instanciate the intersector and initialise the result vector */
 348     /****************************************************************/
 349
 350     long start_filtering=clock();
 351
 352     std::vector<double> bbox;
 353     intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
 354     performAdjustmentOfBB(intersector,bbox);
 355     const double *bboxPtr=0;
 356     if(nbMailleS>0)
 357       bboxPtr=&bbox[0];
 358     BBTree<SPACEDIM,ConnType> my_tree(bboxPtr, 0, 0,nbMailleS);//creating the search structure
 359
 360     long end_filtering=clock();
 361
 362     result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
 363
 364     /****************************************************/
 365     /* Loop on the target cells - core of the algorithm */
 366     /****************************************************/
 367     long start_intersection=clock();
 368     long nbelem_type=myMeshT.getNumberOfElements();
 369     const ConnType *connIndxT=myMeshT.getConnectivityIndexPtr();
 370     for(int iT=0; iT<nbelem_type; iT++)
 371       {
 372         int nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
 373         std::vector<int> intersecting_elems;
 374         double bb[2*SPACEDIM];
 375         intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
 376         my_tree.getIntersectingElems(bb, intersecting_elems);
 377         intersector->intersectCells(iT,intersecting_elems,result);
 378         counter+=intersecting_elems.size();
 379         intersecting_elems.clear();
 380       }
 381     int ret=intersector->getNumberOfColsOfResMatrix();
 382     delete intersector;
 383
 384     if (InterpolationOptions::getPrintLevel() >=1)
 385       {
 386         long end_intersection=clock();
 387         std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
 388         std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
 389         long global_end =clock();
 390         std::cout << "Number of computed intersections = " << counter << std::endl;
 391         std::cout << "Global time= " << global_end - global_start << std::endl;
 392       }
 393     return ret;
 394   }
 395 }
 396
 397 #endif