1 // Copyright (C) 2009-2022 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 src_mesh cartesian source mesh
88 * @param tgt_mesh 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 typename MyCMeshType::MyConnType InterpolationCU::interpolateMeshes(const MyCMeshType& src_mesh,
96 const MyUMeshType& tgt_mesh,
100 typedef typename MyCMeshType::MyConnType CConnType;
101 typedef typename MyUMeshType::MyConnType UConnType;
103 if ( std::string("P0P0") != method )
104 throw Exception("Only P0P0 method is implemented so far");
105 if ( MyCMeshType::MY_SPACEDIM != MyUMeshType::MY_SPACEDIM ||
106 MyCMeshType::MY_SPACEDIM != MyUMeshType::MY_MESHDIM )
107 throw Exception("InterpolationCU::interpolateMeshes(): dimension of meshes must be same");
109 const double eps = getPrecision();
110 const int dim = MyCMeshType::MY_SPACEDIM;
112 TargetIntersector<MyCMeshType, MatrixType>* intersector = 0;
115 case 1: intersector = new IntersectorCU1D<MyCMeshType, MyUMeshType, MatrixType>( src_mesh, tgt_mesh ); break;
116 case 2: intersector = new IntersectorCU2D<MyCMeshType, MyUMeshType, MatrixType>( src_mesh, tgt_mesh ); break;
117 case 3: intersector = new IntersectorCU3D<MyCMeshType, MyUMeshType, MatrixType>( src_mesh, tgt_mesh, getSplittingPolicy() ); break;
119 // create empty maps for all target elements
120 result.resize( intersector->getNumberOfRowsOfResMatrix() );
121 const CConnType ret = intersector->getNumberOfColsOfResMatrix();
123 const double* src_coords[ dim ];
124 CConnType src_nb_coords[ dim ];
125 std::map< double, CConnType> src_coord_to_index[ dim ];
126 for ( int j = 0; j < dim; ++j )
128 int axis = static_cast<int>( _TMIC( j ));
129 src_coords [j] = src_mesh.getCoordsAlongAxis( axis );
130 src_nb_coords[j] = static_cast<CConnType>(src_mesh.nbCellsAlongAxis( axis )) + 1;
131 for (CConnType i = 0; i < src_nb_coords[j]; ++i )
132 src_coord_to_index[j].insert( std::make_pair( src_coords[j][i], i ));
135 const UConnType tgtu_nb_cells = tgt_mesh.getNumberOfElements();
137 IntersectorCU<MyCMeshType, MyUMeshType, MatrixType> bbHelper(src_mesh, tgt_mesh);
140 // loop on unstructured tgt cells
142 for(UConnType iT=0; iT<tgtu_nb_cells; iT++)
144 result[ iT ].clear();
146 // get bounding box of target cell
147 bbHelper.getUElemBB( bb, _TMIU(iT));
149 bool doItersect = true;
150 for ( int j = 0; j < dim && doItersect; ++j )
152 bb[j*2] < src_coords[j][ src_nb_coords[j]-1 ] - eps &&
153 bb[j*2+1] > src_coords[j][0] + eps;
155 continue; // no intersection
157 // find structured src cells intersecting iT cell
158 std::vector< std::vector< CConnType > > structIndices(1);
159 typename std::map< double, CConnType>::iterator coo_ind;
160 for ( int j = 0; j < dim; ++j )
162 coo_ind = src_coord_to_index[j].lower_bound( bb[2*j+1] - eps );
163 if ( coo_ind == src_coord_to_index[j].end() )
165 CConnType max_i = coo_ind->second;
167 coo_ind = src_coord_to_index[j].upper_bound( bb[2*j ] + eps );
168 if ( coo_ind != src_coord_to_index[j].begin() )
170 CConnType min_i = coo_ind->second;
172 std::vector< std::vector< CConnType > > newStructIndices;
173 for ( unsigned int iInd = 0; iInd < structIndices.size(); ++iInd )
175 for ( CConnType i = min_i; i < max_i; ++i )
177 std::vector< CConnType > index = structIndices[iInd];
178 index.push_back( i );
179 newStructIndices.push_back( index );
182 structIndices.swap( newStructIndices );
185 // perform intersection
187 for ( unsigned int iInd = 0; iInd < structIndices.size(); ++iInd )
188 intersector->intersectCells( iT, structIndices[iInd], result );
194 //================================================================================
196 * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
197 * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
198 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
199 * volume of intersection is calculated for the remaining pairs, and entered into the
200 * intersection matrix.
202 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
203 * It can also be an INTERP_KERNEL::Matrix object.
204 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
205 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
206 * which have a non-zero intersection volume with the target element. The vector has indices running from
207 * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
208 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
211 * @param meshS 2-dimesional unstructured target mesh
212 * @param meshT 2-dimensional cartesian source mesh
213 * @param result matrix in which the result is stored
214 * @param method interpolation method
216 //================================================================================
218 template<class MyUMeshType, class MyCMeshType, class MatrixType>
219 typename MyUMeshType::MyConnType InterpolationCU::interpolateMeshesRev(const MyUMeshType& meshS, const MyCMeshType& meshT, MatrixType& result, const char *method)
221 typedef typename MyCMeshType::MyConnType CConnType;
222 typedef typename MyUMeshType::MyConnType UConnType;
224 MatrixType revResult;
225 CConnType sizeT = interpolateMeshes( meshT, meshS, revResult, method );
226 UConnType sizeS = static_cast<UConnType>(revResult.size());
227 result.resize( sizeT );
229 for ( CConnType iS = 0; iS < sizeS; ++iS )
231 typename MatrixType::value_type & row = revResult[iS];
232 typename MatrixType::value_type::iterator iT_surf = row.begin();
233 for ( ; iT_surf != row.end(); ++iT_surf )
234 result[ iT_surf->first ][ iS ] = iT_surf->second;