1 // Copyright (C) 2009-2015 OPEN CASCADE
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 // File : InterpolationCU.txx
20 // Created : Mon Dec 14 17:30:25 2009
21 // Author : Edward AGAPOV (eap)
23 #ifndef __InterpolationCU_TXX__
24 #define __InterpolationCU_TXX__
26 #include "InterpolationCU.hxx"
28 #include "Interpolation.txx"
29 #include "IntersectorCU1D.txx"
30 #include "IntersectorCU2D.txx"
31 #include "IntersectorCU3D.txx"
35 // // convert index "From Mesh Index"
36 #define _FMIU(i) OTT<typename MyUMeshType::MyConnType,MyUMeshType::My_numPol>::ind2C((i))
37 #define _FMIC(i) OTT<typename MyCMeshType::MyConnType,MyCMeshType::My_numPol>::ind2C((i))
38 // convert index "To Mesh Index"
39 #define _TMIU(i) OTT<typename MyUMeshType::MyConnType,MyUMeshType::My_numPol>::indFC((i))
40 #define _TMIC(i) OTT<typename MyCMeshType::MyConnType,MyCMeshType::My_numPol>::indFC((i))
41 // convert coord "From Mesh Coord"
42 #define _FMCOO(i) OTT<typename MyUMeshType::MyConnType,MyUMeshType::My_numPol>::coo2C((i))
43 // convert connectivity "From Mesh Connectivity"
44 #define _FMCON(i) OTT<typename MyUMeshType::MyConnType,MyUMeshType::My_numPol>::conn2C((i))
46 namespace INTERP_KERNEL
49 * \defgroup InterpolationCU InterpolationCU
50 * \class InterpolationCU
51 * \brief Class used to calculate the volumes of intersection between the elements of a cartesian and an unstructured meshes.
54 //================================================================================
59 //================================================================================
61 InterpolationCU::InterpolationCU()
65 InterpolationCU::InterpolationCU(const InterpolationOptions & io)
66 :Interpolation<InterpolationCU>(io)
70 //================================================================================
72 * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
73 * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
74 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
75 * volume of intersection is calculated for the remaining pairs, and entered into the
76 * intersection matrix.
78 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
79 * It can also be an INTERP_KERNEL::Matrix object.
80 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
81 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
82 * which have a non-zero intersection volume with the target element. The vector has indices running from
83 * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
84 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
87 * @param srcMesh cartesian source mesh
88 * @param targetMesh unstructured target mesh
89 * @param result matrix in which the result is stored
90 * @param method interpolation method
92 //================================================================================
94 template<class MyCMeshType, class MyUMeshType, class MatrixType>
95 int InterpolationCU::interpolateMeshes(const MyCMeshType& src_mesh,
96 const MyUMeshType& tgt_mesh,
100 typedef typename MyCMeshType::MyConnType CConnType;
102 if ( std::string("P0P0") != method )
103 throw Exception("Only P0P0 method is implemented so far");
104 if ( MyCMeshType::MY_SPACEDIM != MyUMeshType::MY_SPACEDIM ||
105 MyCMeshType::MY_SPACEDIM != MyUMeshType::MY_MESHDIM )
106 throw Exception("InterpolationCU::interpolateMeshes(): dimension of meshes must be same");
108 const double eps = getPrecision();
109 const int dim = MyCMeshType::MY_SPACEDIM;
111 TargetIntersector<MyCMeshType, MatrixType>* intersector = 0;
114 case 1: intersector = new IntersectorCU1D<MyCMeshType, MyUMeshType, MatrixType>( src_mesh, tgt_mesh ); break;
115 case 2: intersector = new IntersectorCU2D<MyCMeshType, MyUMeshType, MatrixType>( src_mesh, tgt_mesh ); break;
116 case 3: intersector = new IntersectorCU3D<MyCMeshType, MyUMeshType, MatrixType>( src_mesh, tgt_mesh, getSplittingPolicy() ); break;
118 // create empty maps for all target elements
119 result.resize( intersector->getNumberOfRowsOfResMatrix() );
120 const int ret = intersector->getNumberOfColsOfResMatrix();
122 const double* src_coords[ dim ];
123 int src_nb_coords[ dim ];
124 std::map< double, int> src_coord_to_index[ dim ];
125 for ( int j = 0; j < dim; ++j )
127 src_coords [j] = src_mesh.getCoordsAlongAxis( _TMIC( j ));
128 src_nb_coords[j] = src_mesh.nbCellsAlongAxis ( _TMIC( j )) + 1;
129 for (int i = 0; i < src_nb_coords[j]; ++i )
130 src_coord_to_index[j].insert( std::make_pair( src_coords[j][i], i ));
133 const unsigned long tgtu_nb_cells = tgt_mesh.getNumberOfElements();
135 IntersectorCU<MyCMeshType, MyUMeshType, MatrixType> bbHelper(src_mesh, tgt_mesh);
138 // loop on unstructured tgt cells
140 for(unsigned int iT=0; iT<tgtu_nb_cells; iT++)
142 result[ iT ].clear();
144 // get bounding box of target cell
145 bbHelper.getUElemBB( bb, _TMIU(iT));
147 bool doItersect = true;
148 for ( int j = 0; j < dim && doItersect; ++j )
150 bb[j*2] < src_coords[j][ src_nb_coords[j]-1 ] - eps &&
151 bb[j*2+1] > src_coords[j][0] + eps;
153 continue; // no intersection
155 // find structured src cells intersecting iT cell
156 std::vector< std::vector< CConnType > > structIndices(1);
157 std::map< double, int>::iterator coo_ind;
158 for ( int j = 0; j < dim; ++j )
160 coo_ind = src_coord_to_index[j].lower_bound( bb[2*j+1] - eps );
161 if ( coo_ind == src_coord_to_index[j].end() )
163 int max_i = coo_ind->second;
165 coo_ind = src_coord_to_index[j].upper_bound( bb[2*j ] + eps );
166 if ( coo_ind != src_coord_to_index[j].begin() )
168 int min_i = coo_ind->second;
170 std::vector< std::vector< CConnType > > newStructIndices;
171 for ( unsigned int iInd = 0; iInd < structIndices.size(); ++iInd )
173 for ( int i = min_i; i < max_i; ++i )
175 std::vector< CConnType > index = structIndices[iInd];
176 index.push_back( i );
177 newStructIndices.push_back( index );
180 structIndices.swap( newStructIndices );
183 // perform intersection
185 for ( unsigned int iInd = 0; iInd < structIndices.size(); ++iInd )
186 intersector->intersectCells( iT, structIndices[iInd], result );
192 //================================================================================
194 * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
195 * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
196 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
197 * volume of intersection is calculated for the remaining pairs, and entered into the
198 * intersection matrix.
200 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
201 * It can also be an INTERP_KERNEL::Matrix object.
202 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
203 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
204 * which have a non-zero intersection volume with the target element. The vector has indices running from
205 * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
206 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
209 * @param srcMesh 2-dimesional unstructured target mesh
210 * @param targetMesh 2-dimensional cartesian source mesh
211 * @param result matrix in which the result is stored
212 * @param method interpolation method
214 //================================================================================
216 template<class MyUMeshType, class MyCMeshType, class MatrixType>
217 int InterpolationCU::interpolateMeshesRev(const MyUMeshType& meshS, const MyCMeshType& meshT, MatrixType& result, const char *method)
219 MatrixType revResult;
220 int sizeT = interpolateMeshes( meshT, meshS, revResult, method );
221 int sizeS = revResult.size();
222 result.resize( sizeT );
224 for ( int iS = 0; iS < sizeS; ++iS )
226 typename MatrixType::value_type & row = revResult[iS];
227 typename MatrixType::value_type::iterator iT_surf = row.begin();
228 for ( ; iT_surf != row.end(); ++iT_surf )
229 result[ iT_surf->first ][ iS ] = iT_surf->second;