1 // Copyright (C) 2007-2021 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)
20 #ifndef __INTERPOLATIONCURVE_TXX__
21 #define __INTERPOLATIONCURVE_TXX__
23 #include "InterpolationCurve.hxx"
24 #include "InterpolationOptions.hxx"
25 #include "CurveIntersectorP0P0.txx"
26 #include "CurveIntersectorP1P0.txx"
27 #include "CurveIntersectorP0P1.txx"
28 #include "CurveIntersectorP1P1.txx"
29 #include "CurveIntersectorP1P1PL.txx"
34 namespace INTERP_KERNEL
37 template<class RealCurve>
38 InterpolationCurve<RealCurve>::InterpolationCurve()
42 template<class RealCurve>
43 InterpolationCurve<RealCurve>::InterpolationCurve (const InterpolationOptions& io)
44 :Interpolation< InterpolationCurve<RealCurve> >(io)
48 /** \brief Main function to interpolate 1D meshes.
49 \details The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
50 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
51 * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
52 * intersection matrix.
54 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
55 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
56 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
57 * which have a non-zero intersection volume with the target element. The vector has indices running from
58 * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
59 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
62 * @param myMeshS Planar source mesh
63 * @Param myMeshT Planar target mesh
64 * @return vector containing for each element i of the source mesh, a map giving for each element j
65 * of the target mesh which i intersects, the area of the intersection
68 template<class RealCurve>
69 template<class MyMeshType, class MatrixType>
70 typename MyMeshType::MyConnType InterpolationCurve<RealCurve>::interpolateMeshesInternal (const MyMeshType& myMeshS,
71 const MyMeshType& myMeshT,
73 const std::string& method,
74 std::function< void(const BBTree< MyMeshType::MY_SPACEDIM ,
75 typename MyMeshType::MyConnType>&, const double*,
76 std::vector<typename MyMeshType::MyConnType>&) > bbtreeMethod
79 static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
80 typedef typename MyMeshType::MyConnType ConnType;
81 static const NumberingPolicy numPol = MyMeshType::My_numPol;
83 long global_start = clock();
84 std::size_t counter=0;
86 ConnType nbMailleS = myMeshS.getNumberOfElements();
87 ConnType nbMailleT = myMeshT.getNumberOfElements();
89 std::unique_ptr< CurveIntersector<MyMeshType,MatrixType> > intersector;
92 switch (InterpolationOptions::getIntersectionType())
96 intersector.reset( new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
97 InterpolationOptions::getPrecision(),
98 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
99 InterpolationOptions::getMedianPlane(),
100 InterpolationOptions::getPrintLevel()) );
104 throw INTERP_KERNEL::Exception("For P0P0 in 1D or 2D curve only Triangulation supported for the moment !");
107 else if(method=="P0P1")
109 switch (InterpolationOptions::getIntersectionType())
113 intersector.reset( new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
114 InterpolationOptions::getPrecision(),
115 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
116 InterpolationOptions::getMedianPlane(),
117 InterpolationOptions::getPrintLevel()) );
121 throw INTERP_KERNEL::Exception("For P0P1 in 1D or 2D curve only Triangulation supported for the moment !");
124 else if(method=="P1P0")
126 switch (InterpolationOptions::getIntersectionType())
130 intersector.reset( new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
131 InterpolationOptions::getPrecision(),
132 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
133 InterpolationOptions::getMedianPlane(),
134 InterpolationOptions::getPrintLevel()) );
138 throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
141 else if(method=="P1P1")
143 switch (InterpolationOptions::getIntersectionType())
146 intersector.reset( new CurveIntersectorP1P1<MyMeshType,MatrixType>
148 InterpolationOptions::getPrecision(),
149 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
150 InterpolationOptions::getMedianPlane(),
151 InterpolationOptions::getPrintLevel()) );
154 intersector.reset( new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
156 InterpolationOptions::getPrecision(),
157 InterpolationOptions::getBoundingBoxAdjustmentAbs(),
158 InterpolationOptions::getMedianPlane(),
159 InterpolationOptions::getPrintLevel()) );
162 throw INTERP_KERNEL::Exception("For P1P1 in 1D or 2D curve only Triangulation and PointLocator supported !");
166 throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
167 /****************************************************************/
168 /* Create a search tree based on the bounding boxes */
169 /* Instantiate the intersector and initialise the result vector */
170 /****************************************************************/
172 long start_filtering=clock();
174 std::vector<double> bbox;
175 intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
176 intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustmentAbs());
177 BBTree<SPACEDIM,ConnType> my_tree(&bbox[0], 0, 0, nbMailleS);//creating the search structure
179 long end_filtering = clock();
181 result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
183 /****************************************************/
184 /* Loop on the target cells - core of the algorithm */
185 /****************************************************/
186 long start_intersection = clock();
187 const ConnType *connIndxT = myMeshT.getConnectivityIndexPtr();
188 for(ConnType iT=0; iT<nbMailleT; iT++)
190 ConnType nb_nodesT = connIndxT[iT+1] - connIndxT[iT];
191 std::vector<ConnType> intersecting_elems;
192 double bb[2*SPACEDIM];
193 intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
194 bbtreeMethod(my_tree,bb,intersecting_elems);
195 intersector->intersectCells(iT,intersecting_elems,result);
196 counter += intersecting_elems.size();
199 if (InterpolationOptions::getPrintLevel() >= 1)
201 long end_intersection=clock();
202 std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
203 std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
204 long global_end =clock();
205 std::cout << "Number of computed intersections = " << counter << std::endl;
206 std::cout << "Global time= " << global_end - global_start << std::endl;
208 return intersector->getNumberOfColsOfResMatrix();