src/INTERP_KERNEL/InterpolationPlanar.txx

   1 //  Copyright (C) 2007-2008  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 "InterpolationOptions.hxx"
  24 #include "PlanarIntersector.hxx"
  25 #include "PlanarIntersector.txx"
  26 #include "TriangulationIntersector.hxx"
  27 #include "TriangulationIntersector.txx"
  28 #include "ConvexIntersector.hxx"
  29 #include "ConvexIntersector.txx"
  30 #include "Geometric2DIntersector.hxx"
  31 #include "Geometric2DIntersector.txx"
  32 #include "VectorUtils.hxx"
  33 #include "BBTree.txx"
  34
  35 #include <time.h>
  36
  37 namespace INTERP_KERNEL
  38 {
  39
  40   template<class RealPlanar>
  41   const double InterpolationPlanar<RealPlanar>::DEFAULT_PRECISION=1.e-12;
  42
  43   /**
  44    * \defgroup interpolationPlanar InterpolationPlanar
  45    *
  46    * \class InterpolationPlanar
  47    * \brief Class used to compute the coefficients of the interpolation matrix between
  48    * two local meshes in two dimensions. Meshes can contain mixed triangular and quadrangular elements.
  49    */
  50   template<class RealPlanar>
  51   InterpolationPlanar<RealPlanar>::InterpolationPlanar():_dim_caracteristic(1)
  52
  53   {
  54   }
  55
  56   template<class RealPlanar>
  57   InterpolationPlanar<RealPlanar>::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar<RealPlanar> >(io),_dim_caracteristic(1)
  58
  59   {
  60   }
  61
  62
  63   /**
  64    *  \brief  Function used to set the options for the intersection calculation
  65    * \details The following options can be modified:
  66    *  -# Intersection_type: the type of algorithm to be used in the computation of the cell-cell intersections.
  67    *   - Values: Triangle, Convex.
  68    *   - Default: Triangle.
  69    *  -# Precision: Level of precision of the computations is precision times the characteristic size of the mesh.
  70    *   - Values: positive real number.
  71    *   - Default: 1.0E-12.
  72    *  -# PrintLevel: Level of verboseness during the computations.
  73    *   - Values: interger between 0 and 3.
  74    *   - Default: 0.
  75    */
  76   template<class RealPlanar>
  77   void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation)
  78   {
  79     InterpolationOptions::setPrecision(precision);
  80     InterpolationOptions::setPrintLevel(printLevel);
  81     InterpolationOptions::setIntersectionType(intersectionType);
  82     InterpolationOptions::setOrientation(orientation);
  83   }
  84
  85
  86   /** \brief Main function to interpolate triangular or quadrangular meshes.
  87       \details  The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
  88       * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  89       * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
  90       * intersection matrix.
  91       *
  92       * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  93       * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  94       * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  95       * which have a non-zero intersection volume with the target element. The vector has indices running from
  96       * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  97       * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  98       *
  99
 100       * @param myMeshS  Planar source mesh
 101       * @Param myMeshT  Planar target mesh
 102       * @return            vector containing for each element i of the source mesh, a map giving for each element j
 103       *                    of the target mesh which i intersects, the area of the intersection
 104       *
 105       */
 106   template<class RealPlanar>
 107   template<class MyMeshType, class MatrixType>
 108   int InterpolationPlanar<RealPlanar>::interpolateMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, MatrixType& result, const char *method)
 109   {
 110     static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 111     typedef typename MyMeshType::MyConnType ConnType;
 112     static const NumberingPolicy numPol=MyMeshType::My_numPol;
 113
 114     long global_start =clock();
 115     int counter=0;
 116     /***********************************************************/
 117     /* Check both meshes are made of triangles and quadrangles */
 118     /***********************************************************/
 119
 120     long nbMailleS=myMeshS.getNumberOfElements();
 121     long nbMailleT=myMeshT.getNumberOfElements();
 122
 123     /**************************************************/
 124     /* Search the characteristic size of the meshes   */
 125     /**************************************************/
 126
 127     double BoxS[2*SPACEDIM]; myMeshS.getBoundingBox(BoxS);
 128     double BoxT[2*SPACEDIM]; myMeshT.getBoundingBox(BoxT);
 129     double diagonalS=getDistanceBtw2Pts<SPACEDIM>(BoxS+SPACEDIM,BoxS);
 130     double DimCaracteristicS=diagonalS/nbMailleS;
 131     double diagonalT=getDistanceBtw2Pts<SPACEDIM>(BoxT+SPACEDIM,BoxT);
 132     double DimCaracteristicT=diagonalT/nbMailleT;
 133
 134     _dim_caracteristic=std::min(DimCaracteristicS, DimCaracteristicT);
 135     if (InterpolationOptions::getPrintLevel()>=1)
 136       {
 137         std::cout << "  - Characteristic size of the source mesh : " << DimCaracteristicS << std::endl;
 138         std::cout << "  - Characteristic size of the target mesh: " << DimCaracteristicT << std::endl;
 139         std::cout << "InterpolationPlanar::computation of the intersections" << std::endl;
 140       }
 141
 142     PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
 143     std::string meth(method);
 144     if(meth=="P0P0")
 145       {
 146         switch (InterpolationOptions::getIntersectionType())
 147           {
 148           case Triangulation:
 149             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
 150                                                                                                   InterpolationOptions::getPrecision(),
 151                                                                                                   InterpolationOptions::getMedianPlane(),
 152                                                                                                   InterpolationOptions::getOrientation(),
 153                                                                                                   InterpolationOptions::getPrintLevel());
 154             break;
 155           case Convex:
 156             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
 157                                                                                            InterpolationOptions::getPrecision(),
 158                                                                                            InterpolationOptions::getDoRotate(),
 159                                                                                            InterpolationOptions::getMedianPlane(),
 160                                                                                            InterpolationOptions::getOrientation(),
 161                                                                                            InterpolationOptions::getPrintLevel());
 162             break;
 163           case Geometric2D:
 164             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
 165                                                                                                 InterpolationOptions::getMedianPlane(),
 166                                                                                                 InterpolationOptions::getPrecision(),
 167                                                                                                 InterpolationOptions::getOrientation());
 168             break;
 169           }
 170       }
 171     else if(meth=="P0P1")
 172       {
 173         switch (InterpolationOptions::getIntersectionType())
 174           {
 175           case Triangulation:
 176             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
 177                                                                                                   InterpolationOptions::getPrecision(),
 178                                                                                                   InterpolationOptions::getMedianPlane(),
 179                                                                                                   InterpolationOptions::getOrientation(),
 180                                                                                                   InterpolationOptions::getPrintLevel());
 181             break;
 182           case Convex:
 183             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
 184                                                                                            InterpolationOptions::getPrecision(),
 185                                                                                            InterpolationOptions::getDoRotate(),
 186                                                                                            InterpolationOptions::getMedianPlane(),
 187                                                                                            InterpolationOptions::getOrientation(),
 188                                                                                            InterpolationOptions::getPrintLevel());
 189             break;
 190           case Geometric2D:
 191             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT, myMeshS, _dim_caracteristic,
 192                                                                                                 InterpolationOptions::getMedianPlane(),
 193                                                                                                 InterpolationOptions::getPrecision(),
 194                                                                                                 InterpolationOptions::getOrientation());
 195             break;
 196           }
 197       }
 198     else if(meth=="P1P0")
 199       {
 200         switch (InterpolationOptions::getIntersectionType())
 201           {
 202           case Triangulation:
 203             intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
 204                                                                                                   InterpolationOptions::getPrecision(),
 205                                                                                                   InterpolationOptions::getMedianPlane(),
 206                                                                                                   InterpolationOptions::getOrientation(),
 207                                                                                                   InterpolationOptions::getPrintLevel());
 208             break;
 209           case Convex:
 210             intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
 211                                                                                            InterpolationOptions::getPrecision(),
 212                                                                                            InterpolationOptions::getDoRotate(),
 213                                                                                            InterpolationOptions::getMedianPlane(),
 214                                                                                            InterpolationOptions::getOrientation(),
 215                                                                                            InterpolationOptions::getPrintLevel());
 216             break;
 217           case Geometric2D:
 218             intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT, myMeshS, _dim_caracteristic,
 219                                                                                                 InterpolationOptions::getMedianPlane(),
 220                                                                                                 InterpolationOptions::getPrecision(),
 221                                                                                                 InterpolationOptions::getOrientation());
 222             break;
 223           }
 224       }
 225     else
 226       throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" or \"P1P0\"");
 227     /****************************************************************/
 228     /* Create a search tree based on the bounding boxes             */
 229     /* Instanciate the intersector and initialise the result vector */
 230     /****************************************************************/
 231
 232     long start_filtering=clock();
 233
 234     std::vector<double> bbox;
 235     intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
 236     performAdjustmentOfBB(intersector,bbox);
 237     BBTree<SPACEDIM,ConnType> my_tree(&bbox[0], 0, 0,nbMailleS);//creating the search structure
 238
 239     long end_filtering=clock();
 240
 241     result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
 242
 243     /****************************************************/
 244     /* Loop on the target cells - core of the algorithm */
 245     /****************************************************/
 246     long start_intersection=clock();
 247     long nbelem_type=myMeshT.getNumberOfElements();
 248     const ConnType *connIndxT=myMeshT.getConnectivityIndexPtr();
 249     for(int iT=0; iT<nbelem_type; iT++)
 250       {
 251         int nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
 252         std::vector<int> intersecting_elems;
 253         double bb[2*SPACEDIM];
 254         intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
 255         my_tree.getIntersectingElems(bb, intersecting_elems);
 256         intersector->intersectCells(iT,intersecting_elems,result);
 257         counter+=intersecting_elems.size();
 258         intersecting_elems.clear();
 259       }
 260     int ret=intersector->getNumberOfColsOfResMatrix();
 261     delete intersector;
 262
 263     if (InterpolationOptions::getPrintLevel() >=1)
 264       {
 265         long end_intersection=clock();
 266         std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
 267         std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
 268         long global_end =clock();
 269         std::cout << "Number of computed intersections = " << counter << std::endl;
 270         std::cout << "Global time= " << global_end - global_start << std::endl;
 271       }
 272     return ret;
 273   }
 274 }
 275
 276 #endif