1 // Copyright (C) 2009-2012 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.
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 : InterpolationCC.txx
20 // Created : Fri Aug 14 11:39:27 2009
21 // Author : Edward AGAPOV (eap)
24 #include "InterpolationCC.hxx"
25 #include "InterpolationUtils.hxx"
27 // convert index "From Mesh Index"
28 #define _FMI(i) OTT<typename MyMeshType::MyConnType,MyMeshType::My_numPol>::ind2C((i))
29 // convert index "To Mesh Index"
30 #define _TMI(i) OTT<typename MyMeshType::MyConnType,MyMeshType::My_numPol>::indFC((i))
32 namespace INTERP_KERNEL
34 //================================================================================
36 * \brief Constructor does nothing
38 //================================================================================
39 InterpolationCC::InterpolationCC()
43 //================================================================================
45 * \brief An 1D intersection result
47 //================================================================================
51 int _src_index; // source cell index along an axis
52 int _tgt_index; // target cell index along an axis
53 double _length; // interference length
54 Interference(int is = -1, int it = -1, double l = 0):_src_index(is),_tgt_index(it),_length(l){}
57 //================================================================================
59 * \brief Fills the matrix by precomputed cell interferences along axes
60 * \param inter_of_axis - cell/cell interferences along each axis
61 * \param result - matrix to fill in
62 * \param src_nb_cells[] - nb of cells along each of axes in the source mesh
63 * \param tgt_nb_cells[] - nb of cells along each of axes in the target mesh
64 * \param src_i_cell - source cell number accumulated by previous axes
65 * \param tgt_i_cell - target cell number accumulated by previous axes
66 * \param src_prev_area - factor by which this axis icreases cell number
67 * \param tgt_prev_area - factor by which this axis icreases cell number
68 * \param axis - the axis to treat
69 * \param prev_value - intersection size computed by previous axes
71 //================================================================================
73 template <class MyMeshType, class MatrixType, int dim>
74 void fillMatrix(const std::list< Interference > inter_of_axis[dim],
76 const int src_nb_cells[dim],
77 const int tgt_nb_cells[dim],
78 const int src_i_cell = 0,
79 const int tgt_i_cell = 0,
80 const int src_prev_area = 1,
81 const int tgt_prev_area = 1,
83 const double prev_value = 1.0)
85 typedef std::list < Interference >::const_iterator TIntIterator;
87 if ( axis + 1 == dim )
89 for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
91 double value = i->_length * prev_value;
92 int src_i = i->_src_index * src_prev_area + src_i_cell;
93 int tgt_i = i->_tgt_index * tgt_prev_area + tgt_i_cell;
95 result[ tgt_i ].insert( std::make_pair( _TMI( src_i ), value ));
100 int src_prev_area_next = src_prev_area * src_nb_cells[ axis ];
101 int tgt_prev_area_next = tgt_prev_area * tgt_nb_cells[ axis ];
103 for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
105 double value = i->_length * prev_value;
106 int src_i = i->_src_index * src_prev_area + src_i_cell;
107 int tgt_i = i->_tgt_index * tgt_prev_area + tgt_i_cell;
109 // call for the next axis
110 fillMatrix<MyMeshType, MatrixType, dim>(inter_of_axis, result,
111 src_nb_cells, tgt_nb_cells, src_i, tgt_i,
112 src_prev_area_next, tgt_prev_area_next,
118 //================================================================================
120 * \brief Calculates the matrix of volumes of intersection between the elements of
121 * src_mesh and the elements of targetMesh
122 * \param src_mesh - source mesh
123 * \param tgt_mesh - target mesh
124 * \param result - matrix in which the result is stored
125 * \param method - interpolation method, not used as only "P0P0" is implemented so far
127 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
128 * It can also be an INTERP_KERNEL::Matrix object.
129 * The length of the vector is equal to the number of target elements - for each target
130 * element there is a map, regardless of whether the element intersects any source
131 * elements or not. But in the maps there are only entries for those source elements
132 * which have a non-zero intersection volume with the target element. The vector has
133 * indices running from 0 to (nb target elements - 1), meaning that the map for target
134 * element i is stored at index i - 1. In the maps, however, the indexing is more natural:
135 * the intersection volume of the target element i with source element j is found at matrix[i-1][j]
137 //================================================================================
139 template<class MyMeshType, class MatrixType>
140 int InterpolationCC::interpolateMeshes(const MyMeshType& src_mesh,
141 const MyMeshType& tgt_mesh,
145 if ( std::string("P0P0") != method )
146 throw Exception("Only P0P0 method is implemented so far");
148 // create empty maps for all target elements
149 result.resize( tgt_mesh.getNumberOfElements() );
151 const int ret = src_mesh.getNumberOfElements();
153 const double eps = getPrecision();
154 const int dim = MyMeshType::MY_MESHDIM;
155 //const NumberingPolicy numPol = MyMeshType::My_numPol;
157 const double* src_coords[ dim ];
158 const double* tgt_coords[ dim ];
159 int src_nb_cells[ dim ];
160 int tgt_nb_cells[ dim ];
161 for ( int j = 0; j < dim; ++j )
163 src_coords[ j ] = src_mesh.getCoordsAlongAxis( _TMI( j ));
164 tgt_coords[ j ] = tgt_mesh.getCoordsAlongAxis( _TMI( j ));
165 src_nb_cells[ j ] = src_mesh.nbCellsAlongAxis( _TMI( j ));
166 tgt_nb_cells[ j ] = tgt_mesh.nbCellsAlongAxis( _TMI( j ));
169 // ============================================
170 // Calculate cell interferences along the axes
171 // ============================================
173 std::list < Interference > interferences[ dim ];
175 for ( int j = 0; j < dim; ++j ) // loop on axes of castesian space
177 std::list < Interference >& axis_interferences = interferences[j];
180 double x1t, x2t, x1s, x2s; // left and right ordinates of target and source cells
182 // look for the first interference
183 // --------------------------------
184 bool intersection = false;
185 while ( !intersection && it < tgt_nb_cells[j] && is < src_nb_cells[j] )
187 x1s = src_coords[ j ][ is ];
188 x2t = tgt_coords[ j ][ it+1 ];
191 it++; // source lays on the right of target
194 x1t = tgt_coords[ j ][ it ];
195 x2s = src_coords[ j ][ is+1 ];
198 is++; // source lays on the left of target
203 if ( !intersection ) return ret; // no intersections
205 // get all interferences
206 // ----------------------
207 while ( intersection )
209 x1s = src_coords[ j ][ is ];
210 x1t = tgt_coords[ j ][ it ];
211 x2t = tgt_coords[ j ][ it+1 ];
212 x2s = src_coords[ j ][ is+1 ];
214 double x1 = std::max( x1s ,x1t );
215 double x2 = std::min( x2s ,x2t );
216 axis_interferences.push_back( Interference( is, it, x2 - x1 ));
218 // to the next target and/or source cell
219 double diff2 = x2s - x2t;
221 intersection = ( ++it < tgt_nb_cells[j] );
223 intersection = ( ++is < src_nb_cells[j] && intersection);
234 fillMatrix<MyMeshType,MatrixType,3>( interferences, result, src_nb_cells,tgt_nb_cells );
238 fillMatrix<MyMeshType,MatrixType,2>( interferences, result, src_nb_cells,tgt_nb_cells );
242 fillMatrix<MyMeshType,MatrixType,1>( interferences, result, src_nb_cells,tgt_nb_cells );