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