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 : 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, int dim, class MatrixType, class ConnType>
78 void fillMatrix(const std::list< Interference > inter_of_axis[dim],
80 const ConnType src_nb_cells[dim],
81 const ConnType tgt_nb_cells[dim],
82 const ConnType src_i_cell = 0,
83 const ConnType tgt_i_cell = 0,
84 const ConnType src_prev_area = 1,
85 const ConnType 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 ConnType src_i = i->_src_index * src_prev_area + src_i_cell;
97 ConnType tgt_i = i->_tgt_index * tgt_prev_area + tgt_i_cell;
99 result[ tgt_i ].insert( std::make_pair( _TMI( src_i ), value ));
104 ConnType src_prev_area_next = src_prev_area * src_nb_cells[ axis ];
105 ConnType 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 ConnType src_i = i->_src_index * src_prev_area + src_i_cell;
111 ConnType tgt_i = i->_tgt_index * tgt_prev_area + tgt_i_cell;
113 // call for the next axis
114 fillMatrix<MyMeshType, 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 typename MyMeshType::MyConnType 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 typedef typename MyMeshType::MyConnType ConnType;
156 const ConnType ret = src_mesh.getNumberOfElements();
158 const double eps = getPrecision();
159 const int dim = MyMeshType::MY_MESHDIM;
160 //const NumberingPolicy numPol = MyMeshType::My_numPol;
162 const double* src_coords[ dim ];
163 const double* tgt_coords[ dim ];
164 ConnType src_nb_cells[ dim ];
165 ConnType tgt_nb_cells[ dim ];
166 for ( int j = 0; j < dim; ++j )
168 int axis = static_cast<int>( _TMI( j ));
169 src_coords[ j ] = src_mesh.getCoordsAlongAxis( axis );
170 tgt_coords[ j ] = tgt_mesh.getCoordsAlongAxis( axis );
171 src_nb_cells[ j ] = src_mesh.nbCellsAlongAxis( axis );
172 tgt_nb_cells[ j ] = tgt_mesh.nbCellsAlongAxis( axis );
175 // ============================================
176 // Calculate cell interferences along the axes
177 // ============================================
179 std::list < Interference > interferences[ dim ];
181 for ( int j = 0; j < dim; ++j ) // loop on axes of castesian space
183 std::list < Interference >& axis_interferences = interferences[j];
186 double x1t, x2t, x1s, x2s; // left and right ordinates of target and source cells
188 // look for the first interference
189 // --------------------------------
190 bool intersection = false;
191 while ( !intersection && it < tgt_nb_cells[j] && is < src_nb_cells[j] )
193 x1s = src_coords[ j ][ is ];
194 x2t = tgt_coords[ j ][ it+1 ];
197 it++; // source lays on the right of target
200 x1t = tgt_coords[ j ][ it ];
201 x2s = src_coords[ j ][ is+1 ];
204 is++; // source lays on the left of target
209 if ( !intersection ) return ret; // no intersections
211 // get all interferences
212 // ----------------------
213 while ( intersection )
215 x1s = src_coords[ j ][ is ];
216 x1t = tgt_coords[ j ][ it ];
217 x2t = tgt_coords[ j ][ it+1 ];
218 x2s = src_coords[ j ][ is+1 ];
220 double x1 = std::max( x1s ,x1t );
221 double x2 = std::min( x2s ,x2t );
222 axis_interferences.push_back( Interference( is, it, x2 - x1 ));
224 // to the next target and/or source cell
225 double diff2 = x2s - x2t;
227 intersection = ( ++it < tgt_nb_cells[j] );
229 intersection = ( ++is < src_nb_cells[j] && intersection);
240 fillMatrix<MyMeshType,3>( interferences, result, src_nb_cells,tgt_nb_cells );
244 fillMatrix<MyMeshType,2>( interferences, result, src_nb_cells,tgt_nb_cells );
248 fillMatrix<MyMeshType,1>( interferences, result, src_nb_cells,tgt_nb_cells );