1 // Copyright (C) 2007-2008 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 #ifndef __INTERPOLATIONPLANAR_TXX__
20 #define __INTERPOLATIONPLANAR_TXX__
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"
37 namespace INTERP_KERNEL
40 template<class RealPlanar>
41 const double InterpolationPlanar<RealPlanar>::DEFAULT_PRECISION=1.e-12;
44 * \defgroup interpolationPlanar InterpolationPlanar
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.
50 template<class RealPlanar>
51 InterpolationPlanar<RealPlanar>::InterpolationPlanar():_dim_caracteristic(1)
56 template<class RealPlanar>
57 InterpolationPlanar<RealPlanar>::InterpolationPlanar(const InterpolationOptions& io):Interpolation< InterpolationPlanar<RealPlanar> >(io),_dim_caracteristic(1)
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.
72 * -# PrintLevel: Level of verboseness during the computations.
73 * - Values: interger between 0 and 3.
76 template<class RealPlanar>
77 void InterpolationPlanar<RealPlanar>::setOptions(double precision, int printLevel, IntersectionType intersectionType, int orientation)
79 InterpolationOptions::setPrecision(precision);
80 InterpolationOptions::setPrintLevel(printLevel);
81 InterpolationOptions::setIntersectionType(intersectionType);
82 InterpolationOptions::setOrientation(orientation);
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.
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].
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
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)
110 static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
111 typedef typename MyMeshType::MyConnType ConnType;
112 static const NumberingPolicy numPol=MyMeshType::My_numPol;
114 long global_start =clock();
116 /***********************************************************/
117 /* Check both meshes are made of triangles and quadrangles */
118 /***********************************************************/
120 long nbMailleS=myMeshS.getNumberOfElements();
121 long nbMailleT=myMeshT.getNumberOfElements();
123 /**************************************************/
124 /* Search the characteristic size of the meshes */
125 /**************************************************/
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;
134 _dim_caracteristic=std::min(DimCaracteristicS, DimCaracteristicT);
135 if (InterpolationOptions::getPrintLevel()>=1)
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;
142 PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
143 std::string meth(method);
146 switch (InterpolationOptions::getIntersectionType())
149 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
150 InterpolationOptions::getPrecision(),
151 InterpolationOptions::getMedianPlane(),
152 InterpolationOptions::getOrientation(),
153 InterpolationOptions::getPrintLevel());
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());
164 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
165 InterpolationOptions::getMedianPlane(),
166 InterpolationOptions::getPrecision(),
167 InterpolationOptions::getOrientation());
171 else if(meth=="P0P1")
173 switch (InterpolationOptions::getIntersectionType())
176 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
177 InterpolationOptions::getPrecision(),
178 InterpolationOptions::getMedianPlane(),
179 InterpolationOptions::getOrientation(),
180 InterpolationOptions::getPrintLevel());
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());
191 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT, myMeshS, _dim_caracteristic,
192 InterpolationOptions::getMedianPlane(),
193 InterpolationOptions::getPrecision(),
194 InterpolationOptions::getOrientation());
198 else if(meth=="P1P0")
200 switch (InterpolationOptions::getIntersectionType())
203 intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
204 InterpolationOptions::getPrecision(),
205 InterpolationOptions::getMedianPlane(),
206 InterpolationOptions::getOrientation(),
207 InterpolationOptions::getPrintLevel());
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());
218 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT, myMeshS, _dim_caracteristic,
219 InterpolationOptions::getMedianPlane(),
220 InterpolationOptions::getPrecision(),
221 InterpolationOptions::getOrientation());
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 /****************************************************************/
232 long start_filtering=clock();
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
239 long end_filtering=clock();
241 result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
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++)
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();
260 int ret=intersector->getNumberOfColsOfResMatrix();
263 if (InterpolationOptions::getPrintLevel() >=1)
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;