1 // Copyright (C) 2007-2015 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 __INTERPOLATION2D1D_TXX__
21 #define __INTERPOLATION2D1D_TXX__
23 #include "Interpolation2D1D.hxx"
25 namespace INTERP_KERNEL
28 /** \brief Main function to interpolate triangular or quadrangular meshes.
29 \details The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the
30 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
31 * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the
32 * intersection matrix.
34 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
35 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
36 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
37 * which have a non-zero intersection volume with the target element. The vector has indices running from
38 * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
39 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
42 * @param myMeshS Planar source mesh
43 * @Param myMeshT Planar target mesh
44 * @return vector containing for each element i of the source mesh, a map giving for each element j
45 * of the target mesh which i intersects, the area of the intersection
48 template<class MyMeshType, class MatrixType>
49 int Interpolation2D1D::interpolateMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, MatrixType& result, const std::string& method)
51 static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
52 typedef typename MyMeshType::MyConnType ConnType;
53 static const NumberingPolicy numPol=MyMeshType::My_numPol;
55 long global_start =clock();
57 /***********************************************************/
58 /* Check both meshes are made of triangles and quadrangles */
59 /***********************************************************/
61 long nbMailleS=myMeshS.getNumberOfElements();
63 /**************************************************/
64 /* Search the characteristic size of the meshes */
65 /**************************************************/
67 int printLevel = InterpolationOptions::getPrintLevel();
68 _dim_caracteristic = CalculateCharacteristicSizeOfMeshes(myMeshS, myMeshT, printLevel);
71 std::cout << "Interpolation2D1D::computation of the intersections" << std::endl;
74 PlanarIntersector<MyMeshType,MatrixType>* intersector=0;
75 std::string meth = InterpolationOptions::filterInterpolationMethod(method);
78 switch (InterpolationOptions::getIntersectionType())
81 intersector=new Geometric2DIntersector<MyMeshType,MatrixType,Planar2D1DIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
82 InterpolationOptions::getMaxDistance3DSurfIntersect(),
83 InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
84 InterpolationOptions::getMedianPlane(),
85 InterpolationOptions::getPrecision(),
86 InterpolationOptions::getOrientation());
89 throw INTERP_KERNEL::Exception("Invalid intersection type ! Must be : Geometric2D");
93 throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be : \"P0P0\"");
95 /****************************************************************/
96 /* Create a search tree based on the bounding boxes */
97 /* Instanciate the intersector and initialise the result vector */
98 /****************************************************************/
100 long start_filtering=clock();
102 std::vector<double> bbox;
103 intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
104 const double *bboxPtr=0;
107 BBTree<SPACEDIM,ConnType> my_tree(bboxPtr, 0, 0,nbMailleS, -getPrecision());//creating the search structure
109 long end_filtering=clock();
111 result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise.
113 /****************************************************/
114 /* Loop on the target cells - core of the algorithm */
115 /****************************************************/
116 long start_intersection=clock();
117 long nbelem_type=myMeshT.getNumberOfElements();
118 const ConnType *connIndxT=myMeshT.getConnectivityIndexPtr();
119 for(int iT=0; iT<nbelem_type; iT++)
121 int nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
122 std::vector<int> intersecting_elems;
123 double bb[2*SPACEDIM];
124 intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
125 my_tree.getIntersectingElems(bb, intersecting_elems);
126 intersector->intersectCells(iT,intersecting_elems,result);
127 counter+=intersecting_elems.size();
128 intersecting_elems.clear();
130 int ret=intersector->getNumberOfColsOfResMatrix();
132 const DuplicateFacesType& intersectFaces = *intersector->getIntersectFaces();
133 DuplicateFacesType::const_iterator iter;
134 for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter)
136 if (iter->second.size() > 1)
138 _duplicate_faces.insert(std::make_pair(iter->first, iter->second));
144 if (InterpolationOptions::getPrintLevel() >=1)
146 long end_intersection=clock();
147 std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl;
148 std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl;
149 long global_end =clock();
150 std::cout << "Number of computed intersections = " << counter << std::endl;
151 std::cout << "Global time= " << global_end - global_start << std::endl;