1 // Copyright (C) 2007-2020 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #ifndef __INTERPOLATIONPLANAR_TXX__
22 #define __INTERPOLATIONPLANAR_TXX__
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 "MappedBarycentric2DIntersectorP1P1.hxx"
44 #include "MappedBarycentric2DIntersectorP1P1.txx"
45 #include "VectorUtils.hxx"
51 namespace INTERP_KERNEL
53 template<class RealPlanar>
54 InterpolationPlanar<RealPlanar>::InterpolationPlanar():_dim_caracteristic(1)
59 template<class RealPlanar>
60 InterpolationPlanar<RealPlanar>::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar<RealPlanar> >(io),_dim_caracteristic(1)
66 * \brief Function used to set the options for the intersection calculation
67 * \details The following options can be modified:
68 * -# Intersection_type: the type of algorithm to be used in the computation of the cell-cell intersections.
69 * - Values: Triangle, Convex.
70 * - Default: Triangle.
71 * -# Precision: Level of precision of the computations is precision times the characteristic size of the mesh.
72 * - Values: positive real number.
74 * -# PrintLevel: Level of verboseness during the computations.
75 * - Values: integer between 0 and 3.
78 template<class RealPlanar>
79 void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation)
81 InterpolationOptions::setPrecision(precision);
82 InterpolationOptions::setPrintLevel(printLevel);
83 InterpolationOptions::setIntersectionType(intersectionType);
84 InterpolationOptions::setOrientation(orientation);
88 /** \brief Main function to interpolate triangular or quadrangular meshes.
89 \details The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
90 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
91 * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
92 * intersection matrix.
94 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
95 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
96 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
97 * which have a non-zero intersection volume with the target element. The vector has indices running from
98 * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
99 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
102 * @param myMeshS Planar source mesh
103 * @Param myMeshT Planar target mesh
104 * @return vector containing for each element i of the source mesh, a map giving for each element j
105 * of the target mesh which i intersects, the area of the intersection
108 template<class RealPlanar>
109 template<class MyMeshType, class MatrixType>
110 typename MyMeshType::MyConnType InterpolationPlanar<RealPlanar>::interpolateMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, MatrixType& result, const std::string& method)
112 static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
113 typedef typename MyMeshType::MyConnType ConnType;
114 static const NumberingPolicy numPol=MyMeshType::My_numPol;
116 long global_start =clock();
117 std::size_t counter=0;
118 /***********************************************************/
119 /* Check both meshes are made of triangles and quadrangles */
120 /***********************************************************/
122 ConnType nbMailleS=myMeshS.getNumberOfElements();
123 ConnType nbMailleT=myMeshT.getNumberOfElements();
125 /**************************************************/
126 /* Search the characteristic size of the meshes */
127 /**************************************************/
129 double BoxS[2*SPACEDIM]; myMeshS.getBoundingBox(BoxS);
130 double BoxT[2*SPACEDIM]; myMeshT.getBoundingBox(BoxT);
131 double diagonalS,dimCaracteristicS=std::numeric_limits<double>::max();
134 diagonalS=getDistanceBtw2Pts<SPACEDIM>(BoxS+SPACEDIM,BoxS);
135 dimCaracteristicS=diagonalS/(double)(nbMailleS);
137 double diagonalT,dimCaracteristicT=std::numeric_limits<double>::max();
140 diagonalT=getDistanceBtw2Pts<SPACEDIM>(BoxT+SPACEDIM,BoxT);
141 dimCaracteristicT=diagonalT/(double)(nbMailleT);
144 _dim_caracteristic=std::min(dimCaracteristicS, dimCaracteristicT);
145 if (InterpolationOptions::getPrintLevel()>=1)
147 std::cout << " - Characteristic size of the source mesh : " << dimCaracteristicS << std::endl;
148 std::cout << " - Characteristic size of the target mesh: " << dimCaracteristicT << std::endl;
149 std::cout << "InterpolationPlanar::computation of the intersections" << std::endl;
152 PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
153 std::string meth = InterpolationOptions::filterInterpolationMethod(method);
156 switch (InterpolationOptions::getIntersectionType())
159 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
160 InterpolationOptions::getPrecision(),
161 InterpolationOptions::getMaxDistance3DSurfIntersect(),
162 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
163 InterpolationOptions::getMedianPlane(),
164 InterpolationOptions::getOrientation(),
165 InterpolationOptions::getPrintLevel());
168 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
169 InterpolationOptions::getPrecision(),
170 InterpolationOptions::getMaxDistance3DSurfIntersect(),
171 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
172 InterpolationOptions::getMedianPlane(),
173 InterpolationOptions::getDoRotate(),
174 InterpolationOptions::getOrientation(),
175 InterpolationOptions::getPrintLevel());
178 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
179 InterpolationOptions::getMaxDistance3DSurfIntersect(),
180 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
181 InterpolationOptions::getMedianPlane(),
182 InterpolationOptions::getPrecision(),
183 InterpolationOptions::getOrientation());
186 intersector=new PointLocator2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
187 InterpolationOptions::getMaxDistance3DSurfIntersect(),
188 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
189 InterpolationOptions::getMedianPlane(),
190 InterpolationOptions::getPrecision(),
191 InterpolationOptions::getOrientation());
194 throw INTERP_KERNEL::Exception("For P0P0 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator !");
197 else if(meth=="P0P1")
199 switch (InterpolationOptions::getIntersectionType())
202 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
203 InterpolationOptions::getPrecision(),
204 InterpolationOptions::getMaxDistance3DSurfIntersect(),
205 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
206 InterpolationOptions::getMedianPlane(),
207 InterpolationOptions::getOrientation(),
208 InterpolationOptions::getPrintLevel());
211 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
212 InterpolationOptions::getPrecision(),
213 InterpolationOptions::getMaxDistance3DSurfIntersect(),
214 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
215 InterpolationOptions::getMedianPlane(),
216 InterpolationOptions::getDoRotate(),
217 InterpolationOptions::getOrientation(),
218 InterpolationOptions::getPrintLevel());
221 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT, myMeshS, _dim_caracteristic,
222 InterpolationOptions::getMaxDistance3DSurfIntersect(),
223 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
224 InterpolationOptions::getMedianPlane(),
225 InterpolationOptions::getPrecision(),
226 InterpolationOptions::getOrientation());
229 intersector=new PlanarIntersectorP0P1PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
230 InterpolationOptions::getMaxDistance3DSurfIntersect(),
231 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
232 InterpolationOptions::getMedianPlane(),
233 InterpolationOptions::getPrecision(),
234 InterpolationOptions::getOrientation());
237 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1Bary>(myMeshT,myMeshS,_dim_caracteristic,
238 InterpolationOptions::getPrecision(),
239 InterpolationOptions::getMaxDistance3DSurfIntersect(),
240 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
241 InterpolationOptions::getMedianPlane(),
242 InterpolationOptions::getOrientation(),
243 InterpolationOptions::getPrintLevel());
245 case BarycentricGeo2D:
246 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1Bary>(myMeshT, myMeshS, _dim_caracteristic,
247 InterpolationOptions::getMaxDistance3DSurfIntersect(),
248 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
249 InterpolationOptions::getMedianPlane(),
250 InterpolationOptions::getPrecision(),
251 InterpolationOptions::getOrientation());
254 throw INTERP_KERNEL::Exception("For P0P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator, Barycentric, BarycentricGeo2D !");
257 else if(meth=="P1P0")
259 switch (InterpolationOptions::getIntersectionType())
262 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
263 InterpolationOptions::getPrecision(),
264 InterpolationOptions::getMaxDistance3DSurfIntersect(),
265 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
266 InterpolationOptions::getMedianPlane(),
267 InterpolationOptions::getOrientation(),
268 InterpolationOptions::getPrintLevel());
271 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
272 InterpolationOptions::getPrecision(),
273 InterpolationOptions::getMaxDistance3DSurfIntersect(),
274 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
275 InterpolationOptions::getMedianPlane(),
276 InterpolationOptions::getDoRotate(),
277 InterpolationOptions::getOrientation(),
278 InterpolationOptions::getPrintLevel());
281 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT, myMeshS, _dim_caracteristic,
282 InterpolationOptions::getMaxDistance3DSurfIntersect(),
283 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
284 InterpolationOptions::getMedianPlane(),
285 InterpolationOptions::getPrecision(),
286 InterpolationOptions::getOrientation());
289 intersector=new PlanarIntersectorP1P0PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
290 InterpolationOptions::getMaxDistance3DSurfIntersect(),
291 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
292 InterpolationOptions::getMedianPlane(),
293 InterpolationOptions::getPrecision(),
294 InterpolationOptions::getOrientation());
297 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0Bary>(myMeshT,myMeshS,_dim_caracteristic,
298 InterpolationOptions::getPrecision(),
299 InterpolationOptions::getMaxDistance3DSurfIntersect(),
300 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
301 InterpolationOptions::getMedianPlane(),
302 InterpolationOptions::getOrientation(),
303 InterpolationOptions::getPrintLevel());
305 case BarycentricGeo2D:
306 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0Bary>(myMeshT, myMeshS, _dim_caracteristic,
307 InterpolationOptions::getMaxDistance3DSurfIntersect(),
308 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
309 InterpolationOptions::getMedianPlane(),
310 InterpolationOptions::getPrecision(),
311 InterpolationOptions::getOrientation());
314 throw INTERP_KERNEL::Exception("For P1P0 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator, BarycentricGeo2D or Barycentric!");
317 else if(meth=="P1P1")
319 switch (InterpolationOptions::getIntersectionType())
322 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
323 InterpolationOptions::getPrecision(),
324 InterpolationOptions::getMaxDistance3DSurfIntersect(),
325 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
326 InterpolationOptions::getMedianPlane(),
327 InterpolationOptions::getOrientation(),
328 InterpolationOptions::getPrintLevel());
331 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
332 InterpolationOptions::getPrecision(),
333 InterpolationOptions::getMaxDistance3DSurfIntersect(),
334 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
335 InterpolationOptions::getMedianPlane(),
336 InterpolationOptions::getDoRotate(),
337 InterpolationOptions::getOrientation(),
338 InterpolationOptions::getPrintLevel());
341 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT, myMeshS, _dim_caracteristic,
342 InterpolationOptions::getMaxDistance3DSurfIntersect(),
343 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
344 InterpolationOptions::getMedianPlane(),
345 InterpolationOptions::getPrecision(),
346 InterpolationOptions::getOrientation());
349 intersector=new PlanarIntersectorP1P1PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
350 InterpolationOptions::getMaxDistance3DSurfIntersect(),
351 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
352 InterpolationOptions::getMedianPlane(),
353 InterpolationOptions::getPrecision(),
354 InterpolationOptions::getOrientation());
356 case MappedBarycentric:
357 intersector=new MappedBarycentric2DIntersectorP1P1<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
358 InterpolationOptions::getMaxDistance3DSurfIntersect(),
359 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
360 InterpolationOptions::getMedianPlane(),
361 InterpolationOptions::getPrecision(),
362 InterpolationOptions::getOrientation());
365 throw INTERP_KERNEL::Exception("For P1P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator, MappedBarycentric !");
369 throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
370 /****************************************************************/
371 /* Create a search tree based on the bounding boxes */
372 /* Instantiate the intersector and initialise the result vector */
373 /****************************************************************/
375 long start_filtering=clock();
377 std::vector<double> bbox;
378 intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
379 performAdjustmentOfBB(intersector,bbox);
380 const double *bboxPtr=0;
383 BBTree<SPACEDIM,ConnType> my_tree(bboxPtr, 0, 0,nbMailleS);//creating the search structure
385 long end_filtering=clock();
387 result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
389 /****************************************************/
390 /* Loop on the target cells - core of the algorithm */
391 /****************************************************/
392 long start_intersection=clock();
393 ConnType nbelem_type=myMeshT.getNumberOfElements();
394 const ConnType *connIndxT=myMeshT.getConnectivityIndexPtr();
395 for(ConnType iT=0; iT<nbelem_type; iT++)
397 ConnType nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
398 std::vector<ConnType> intersecting_elems;
399 double bb[2*SPACEDIM];
400 intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
401 my_tree.getIntersectingElems(bb, intersecting_elems);
402 intersector->intersectCells(iT,intersecting_elems,result);
403 counter+=intersecting_elems.size();
404 intersecting_elems.clear();
406 ConnType ret=intersector->getNumberOfColsOfResMatrix();
409 if (InterpolationOptions::getPrintLevel() >=1)
411 long end_intersection=clock();
412 std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
413 std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
414 long global_end =clock();
415 std::cout << "Number of computed intersections = " << counter << std::endl;
416 std::cout << "Global time= " << global_end - global_start << std::endl;