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 InterpolationCC::InterpolationCC(const InterpolationOptions& io):Interpolation<InterpolationCC>(io)
47 //================================================================================
49 * \brief An 1D intersection result
51 //================================================================================
55 int _src_index; // source cell index along an axis
56 int _tgt_index; // target cell index along an axis
57 double _length; // interference length
58 Interference(int is = -1, int it = -1, double l = 0):_src_index(is),_tgt_index(it),_length(l){}
61 //================================================================================
63 * \brief Fills the matrix by precomputed cell interferences along axes
64 * \param inter_of_axis - cell/cell interferences along each axis
65 * \param result - matrix to fill in
66 * \param src_nb_cells[] - nb of cells along each of axes in the source mesh
67 * \param tgt_nb_cells[] - nb of cells along each of axes in the target mesh
68 * \param src_i_cell - source cell number accumulated by previous axes
69 * \param tgt_i_cell - target cell number accumulated by previous axes
70 * \param src_prev_area - factor by which this axis icreases cell number
71 * \param tgt_prev_area - factor by which this axis icreases cell number
72 * \param axis - the axis to treat
73 * \param prev_value - intersection size computed by previous axes
75 //================================================================================
77 template <class MyMeshType, class MatrixType, int dim>
78 void fillMatrix(const std::list< Interference > inter_of_axis[dim],
80 const int src_nb_cells[dim],
81 const int tgt_nb_cells[dim],
82 const int src_i_cell = 0,
83 const int tgt_i_cell = 0,
84 const int src_prev_area = 1,
85 const int tgt_prev_area = 1,
87 const double prev_value = 1.0)
89 typedef std::list < Interference >::const_iterator TIntIterator;
91 if ( axis + 1 == dim )
93 for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
95 double value = i->_length * prev_value;
96 int src_i = i->_src_index * src_prev_area + src_i_cell;
97 int tgt_i = i->_tgt_index * tgt_prev_area + tgt_i_cell;
99 result[ tgt_i ].insert( std::make_pair( _TMI( src_i ), value ));
104 int src_prev_area_next = src_prev_area * src_nb_cells[ axis ];
105 int tgt_prev_area_next = tgt_prev_area * tgt_nb_cells[ axis ];
107 for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
109 double value = i->_length * prev_value;
110 int src_i = i->_src_index * src_prev_area + src_i_cell;
111 int tgt_i = i->_tgt_index * tgt_prev_area + tgt_i_cell;
113 // call for the next axis
114 fillMatrix<MyMeshType, MatrixType, dim>(inter_of_axis, result,
115 src_nb_cells, tgt_nb_cells, src_i, tgt_i,
116 src_prev_area_next, tgt_prev_area_next,
122 //================================================================================
124 * \brief Calculates the matrix of volumes of intersection between the elements of
125 * src_mesh and the elements of targetMesh
126 * \param src_mesh - source mesh
127 * \param tgt_mesh - target mesh
128 * \param result - matrix in which the result is stored
129 * \param method - interpolation method, not used as only "P0P0" is implemented so far
131 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
132 * It can also be an INTERP_KERNEL::Matrix object.
133 * The length of the vector is equal to the number of target elements - for each target
134 * element there is a map, regardless of whether the element intersects any source
135 * elements or not. But in the maps there are only entries for those source elements
136 * which have a non-zero intersection volume with the target element. The vector has
137 * indices running from 0 to (nb target elements - 1), meaning that the map for target
138 * element i is stored at index i - 1. In the maps, however, the indexing is more natural:
139 * the intersection volume of the target element i with source element j is found at matrix[i-1][j]
141 //================================================================================
143 template<class MyMeshType, class MatrixType>
144 int InterpolationCC::interpolateMeshes(const MyMeshType& src_mesh,
145 const MyMeshType& tgt_mesh,
149 if ( std::string("P0P0") != method )
150 throw Exception("Only P0P0 method is implemented so far");
152 // create empty maps for all target elements
153 result.resize( tgt_mesh.getNumberOfElements() );
155 const int ret = src_mesh.getNumberOfElements();
157 const double eps = getPrecision();
158 const int dim = MyMeshType::MY_MESHDIM;
159 //const NumberingPolicy numPol = MyMeshType::My_numPol;
161 const double* src_coords[ dim ];
162 const double* tgt_coords[ dim ];
163 int src_nb_cells[ dim ];
164 int tgt_nb_cells[ dim ];
165 for ( int j = 0; j < dim; ++j )
167 src_coords[ j ] = src_mesh.getCoordsAlongAxis( _TMI( j ));
168 tgt_coords[ j ] = tgt_mesh.getCoordsAlongAxis( _TMI( j ));
169 src_nb_cells[ j ] = src_mesh.nbCellsAlongAxis( _TMI( j ));
170 tgt_nb_cells[ j ] = tgt_mesh.nbCellsAlongAxis( _TMI( j ));
173 // ============================================
174 // Calculate cell interferences along the axes
175 // ============================================
177 std::list < Interference > interferences[ dim ];
179 for ( int j = 0; j < dim; ++j ) // loop on axes of castesian space
181 std::list < Interference >& axis_interferences = interferences[j];
184 double x1t, x2t, x1s, x2s; // left and right ordinates of target and source cells
186 // look for the first interference
187 // --------------------------------
188 bool intersection = false;
189 while ( !intersection && it < tgt_nb_cells[j] && is < src_nb_cells[j] )
191 x1s = src_coords[ j ][ is ];
192 x2t = tgt_coords[ j ][ it+1 ];
195 it++; // source lays on the right of target
198 x1t = tgt_coords[ j ][ it ];
199 x2s = src_coords[ j ][ is+1 ];
202 is++; // source lays on the left of target
207 if ( !intersection ) return ret; // no intersections
209 // get all interferences
210 // ----------------------
211 while ( intersection )
213 x1s = src_coords[ j ][ is ];
214 x1t = tgt_coords[ j ][ it ];
215 x2t = tgt_coords[ j ][ it+1 ];
216 x2s = src_coords[ j ][ is+1 ];
218 double x1 = std::max( x1s ,x1t );
219 double x2 = std::min( x2s ,x2t );
220 axis_interferences.push_back( Interference( is, it, x2 - x1 ));
222 // to the next target and/or source cell
223 double diff2 = x2s - x2t;
225 intersection = ( ++it < tgt_nb_cells[j] );
227 intersection = ( ++is < src_nb_cells[j] && intersection);
238 fillMatrix<MyMeshType,MatrixType,3>( interferences, result, src_nb_cells,tgt_nb_cells );
242 fillMatrix<MyMeshType,MatrixType,2>( interferences, result, src_nb_cells,tgt_nb_cells );
246 fillMatrix<MyMeshType,MatrixType,1>( interferences, result, src_nb_cells,tgt_nb_cells );