1 // Copyright (C) 2007-2013 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.
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 "VectorUtils.hxx"
49 namespace INTERP_KERNEL
52 * \defgroup interpolationPlanar InterpolationPlanar
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.
58 template<class RealPlanar>
59 InterpolationPlanar<RealPlanar>::InterpolationPlanar():_dim_caracteristic(1)
64 template<class RealPlanar>
65 InterpolationPlanar<RealPlanar>::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar<RealPlanar> >(io),_dim_caracteristic(1)
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.
79 * -# PrintLevel: Level of verboseness during the computations.
80 * - Values: interger between 0 and 3.
83 template<class RealPlanar>
84 void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation)
86 InterpolationOptions::setPrecision(precision);
87 InterpolationOptions::setPrintLevel(printLevel);
88 InterpolationOptions::setIntersectionType(intersectionType);
89 InterpolationOptions::setOrientation(orientation);
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.
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].
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
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)
117 static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
118 typedef typename MyMeshType::MyConnType ConnType;
119 static const NumberingPolicy numPol=MyMeshType::My_numPol;
121 long global_start =clock();
123 /***********************************************************/
124 /* Check both meshes are made of triangles and quadrangles */
125 /***********************************************************/
127 long nbMailleS=myMeshS.getNumberOfElements();
128 long nbMailleT=myMeshT.getNumberOfElements();
130 /**************************************************/
131 /* Search the characteristic size of the meshes */
132 /**************************************************/
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();
139 diagonalS=getDistanceBtw2Pts<SPACEDIM>(BoxS+SPACEDIM,BoxS);
140 dimCaracteristicS=diagonalS/nbMailleS;
142 double diagonalT,dimCaracteristicT=std::numeric_limits<double>::max();
145 diagonalT=getDistanceBtw2Pts<SPACEDIM>(BoxT+SPACEDIM,BoxT);
146 dimCaracteristicT=diagonalT/nbMailleT;
149 _dim_caracteristic=std::min(dimCaracteristicS, dimCaracteristicT);
150 if (InterpolationOptions::getPrintLevel()>=1)
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;
157 PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
158 std::string meth = InterpolationOptions::filterInterpolationMethod(method);
161 switch (InterpolationOptions::getIntersectionType())
164 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
165 InterpolationOptions::getPrecision(),
166 InterpolationOptions::getMaxDistance3DSurfIntersect(),
167 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
168 InterpolationOptions::getMedianPlane(),
169 InterpolationOptions::getOrientation(),
170 InterpolationOptions::getPrintLevel());
173 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
174 InterpolationOptions::getPrecision(),
175 InterpolationOptions::getMaxDistance3DSurfIntersect(),
176 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
177 InterpolationOptions::getMedianPlane(),
178 InterpolationOptions::getDoRotate(),
179 InterpolationOptions::getOrientation(),
180 InterpolationOptions::getPrintLevel());
183 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
184 InterpolationOptions::getMaxDistance3DSurfIntersect(),
185 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
186 InterpolationOptions::getMedianPlane(),
187 InterpolationOptions::getPrecision(),
188 InterpolationOptions::getOrientation());
191 intersector=new PointLocator2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
192 InterpolationOptions::getMaxDistance3DSurfIntersect(),
193 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
194 InterpolationOptions::getMedianPlane(),
195 InterpolationOptions::getPrecision(),
196 InterpolationOptions::getOrientation());
199 throw INTERP_KERNEL::Exception("For P0P0 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator !");
202 else if(meth=="P0P1")
204 switch (InterpolationOptions::getIntersectionType())
207 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
208 InterpolationOptions::getPrecision(),
209 InterpolationOptions::getMaxDistance3DSurfIntersect(),
210 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
211 InterpolationOptions::getMedianPlane(),
212 InterpolationOptions::getOrientation(),
213 InterpolationOptions::getPrintLevel());
216 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
217 InterpolationOptions::getPrecision(),
218 InterpolationOptions::getMaxDistance3DSurfIntersect(),
219 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
220 InterpolationOptions::getMedianPlane(),
221 InterpolationOptions::getDoRotate(),
222 InterpolationOptions::getOrientation(),
223 InterpolationOptions::getPrintLevel());
226 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT, myMeshS, _dim_caracteristic,
227 InterpolationOptions::getMaxDistance3DSurfIntersect(),
228 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
229 InterpolationOptions::getMedianPlane(),
230 InterpolationOptions::getPrecision(),
231 InterpolationOptions::getOrientation());
234 intersector=new PlanarIntersectorP0P1PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
235 InterpolationOptions::getMaxDistance3DSurfIntersect(),
236 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
237 InterpolationOptions::getMedianPlane(),
238 InterpolationOptions::getPrecision(),
239 InterpolationOptions::getOrientation());
242 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1Bary>(myMeshT,myMeshS,_dim_caracteristic,
243 InterpolationOptions::getPrecision(),
244 InterpolationOptions::getMaxDistance3DSurfIntersect(),
245 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
246 InterpolationOptions::getMedianPlane(),
247 InterpolationOptions::getOrientation(),
248 InterpolationOptions::getPrintLevel());
250 case BarycentricGeo2D:
251 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1Bary>(myMeshT, myMeshS, _dim_caracteristic,
252 InterpolationOptions::getMaxDistance3DSurfIntersect(),
253 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
254 InterpolationOptions::getMedianPlane(),
255 InterpolationOptions::getPrecision(),
256 InterpolationOptions::getOrientation());
259 throw INTERP_KERNEL::Exception("For P0P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator, Barycentric, BarycentricGeo2D !");
262 else if(meth=="P1P0")
264 switch (InterpolationOptions::getIntersectionType())
267 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
268 InterpolationOptions::getPrecision(),
269 InterpolationOptions::getMaxDistance3DSurfIntersect(),
270 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
271 InterpolationOptions::getMedianPlane(),
272 InterpolationOptions::getOrientation(),
273 InterpolationOptions::getPrintLevel());
276 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
277 InterpolationOptions::getPrecision(),
278 InterpolationOptions::getMaxDistance3DSurfIntersect(),
279 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
280 InterpolationOptions::getMedianPlane(),
281 InterpolationOptions::getDoRotate(),
282 InterpolationOptions::getOrientation(),
283 InterpolationOptions::getPrintLevel());
286 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT, myMeshS, _dim_caracteristic,
287 InterpolationOptions::getMaxDistance3DSurfIntersect(),
288 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
289 InterpolationOptions::getMedianPlane(),
290 InterpolationOptions::getPrecision(),
291 InterpolationOptions::getOrientation());
294 intersector=new PlanarIntersectorP1P0PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
295 InterpolationOptions::getMaxDistance3DSurfIntersect(),
296 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
297 InterpolationOptions::getMedianPlane(),
298 InterpolationOptions::getPrecision(),
299 InterpolationOptions::getOrientation());
302 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0Bary>(myMeshT,myMeshS,_dim_caracteristic,
303 InterpolationOptions::getPrecision(),
304 InterpolationOptions::getMaxDistance3DSurfIntersect(),
305 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
306 InterpolationOptions::getMedianPlane(),
307 InterpolationOptions::getOrientation(),
308 InterpolationOptions::getPrintLevel());
310 case BarycentricGeo2D:
311 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0Bary>(myMeshT, myMeshS, _dim_caracteristic,
312 InterpolationOptions::getMaxDistance3DSurfIntersect(),
313 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
314 InterpolationOptions::getMedianPlane(),
315 InterpolationOptions::getPrecision(),
316 InterpolationOptions::getOrientation());
320 else if(meth=="P1P1")
322 switch (InterpolationOptions::getIntersectionType())
325 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
326 InterpolationOptions::getPrecision(),
327 InterpolationOptions::getMaxDistance3DSurfIntersect(),
328 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
329 InterpolationOptions::getMedianPlane(),
330 InterpolationOptions::getOrientation(),
331 InterpolationOptions::getPrintLevel());
334 intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
335 InterpolationOptions::getPrecision(),
336 InterpolationOptions::getMaxDistance3DSurfIntersect(),
337 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
338 InterpolationOptions::getMedianPlane(),
339 InterpolationOptions::getDoRotate(),
340 InterpolationOptions::getOrientation(),
341 InterpolationOptions::getPrintLevel());
344 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT, myMeshS, _dim_caracteristic,
345 InterpolationOptions::getMaxDistance3DSurfIntersect(),
346 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
347 InterpolationOptions::getMedianPlane(),
348 InterpolationOptions::getPrecision(),
349 InterpolationOptions::getOrientation());
352 intersector=new PlanarIntersectorP1P1PL<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
353 InterpolationOptions::getMaxDistance3DSurfIntersect(),
354 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
355 InterpolationOptions::getMedianPlane(),
356 InterpolationOptions::getPrecision(),
357 InterpolationOptions::getOrientation());
360 throw INTERP_KERNEL::Exception("For P1P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator !");
364 throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
365 /****************************************************************/
366 /* Create a search tree based on the bounding boxes */
367 /* Instanciate the intersector and initialise the result vector */
368 /****************************************************************/
370 long start_filtering=clock();
372 std::vector<double> bbox;
373 intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
374 performAdjustmentOfBB(intersector,bbox);
375 const double *bboxPtr=0;
378 BBTree<SPACEDIM,ConnType> my_tree(bboxPtr, 0, 0,nbMailleS);//creating the search structure
380 long end_filtering=clock();
382 result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
384 /****************************************************/
385 /* Loop on the target cells - core of the algorithm */
386 /****************************************************/
387 long start_intersection=clock();
388 long nbelem_type=myMeshT.getNumberOfElements();
389 const ConnType *connIndxT=myMeshT.getConnectivityIndexPtr();
390 for(int iT=0; iT<nbelem_type; iT++)
392 int nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
393 std::vector<int> intersecting_elems;
394 double bb[2*SPACEDIM];
395 intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
396 my_tree.getIntersectingElems(bb, intersecting_elems);
397 intersector->intersectCells(iT,intersecting_elems,result);
398 counter+=intersecting_elems.size();
399 intersecting_elems.clear();
401 int ret=intersector->getNumberOfColsOfResMatrix();
404 if (InterpolationOptions::getPrintLevel() >=1)
406 long end_intersection=clock();
407 std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
408 std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
409 long global_end =clock();
410 std::cout << "Number of computed intersections = " << counter << std::endl;
411 std::cout << "Global time= " << global_end - global_start << std::endl;