src/INTERP_KERNEL/InterpolationCC.txx

   1 // Copyright (C) 2009-2012  OPEN CASCADE
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : InterpolationCC.txx
  20 // Created   : Fri Aug 14 11:39:27 2009
  21 // Author    : Edward AGAPOV (eap)
  22 //
  23
  24 #include "InterpolationCC.hxx"
  25 #include "InterpolationUtils.hxx"
  26
  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))
  31
  32 namespace INTERP_KERNEL
  33 {
  34   //================================================================================
  35   /*!
  36    * \brief Constructor does nothing
  37    */
  38   //================================================================================
  39   InterpolationCC::InterpolationCC()
  40   {
  41   }
  42
  43   //================================================================================
  44   /*!
  45    * \brief An 1D intersection result
  46    */
  47   //================================================================================
  48
  49   struct Interference
  50   {
  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){}
  55   };
  56
  57   //================================================================================
  58   /*!
  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
  70    */
  71   //================================================================================
  72
  73   template <class MyMeshType, class MatrixType, int dim>
  74   void fillMatrix(const std::list< Interference >  inter_of_axis[dim],
  75                   MatrixType&                      result,
  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,
  82                   const int                        axis = 0,
  83                   const double                     prev_value = 1.0)
  84   {
  85     typedef std::list < Interference >::const_iterator TIntIterator;
  86
  87     if ( axis + 1 == dim )
  88     {
  89       for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
  90       {
  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;
  94
  95         result[ tgt_i ].insert( std::make_pair( _TMI( src_i ), value ));
  96       }
  97     }
  98     else
  99     {
 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 ];
 102
 103       for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
 104       {
 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;
 108
 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,
 113                                                 axis+1, value );
 114       }
 115     }
 116   }
 117
 118   //================================================================================
 119   /*!
 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
 126    *
 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]
 136    */
 137   //================================================================================
 138
 139   template<class MyMeshType, class MatrixType>
 140   int InterpolationCC::interpolateMeshes(const MyMeshType& src_mesh,
 141                                          const MyMeshType& tgt_mesh,
 142                                          MatrixType&       result,
 143                                          const char *      method)
 144   {
 145     if ( std::string("P0P0") != method )
 146       throw Exception("Only P0P0 method is implemented so far");
 147
 148     // create empty maps for all target elements
 149     result.resize( tgt_mesh.getNumberOfElements() );
 150
 151     const int ret = src_mesh.getNumberOfElements();
 152
 153     const double eps = getPrecision();
 154     const int dim = MyMeshType::MY_MESHDIM;
 155     //const NumberingPolicy numPol = MyMeshType::My_numPol;
 156
 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 )
 162     {
 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 ));
 167     }
 168
 169     // ============================================
 170     // Calculate cell interferences along the axes
 171     // ============================================
 172
 173     std::list < Interference > interferences[ dim ];
 174
 175     for ( int j = 0; j < dim; ++j ) // loop on axes of castesian space
 176     {
 177       std::list < Interference >& axis_interferences = interferences[j];
 178
 179       int it = 0, is = 0;
 180       double x1t, x2t, x1s, x2s; // left and right ordinates of target and source cells
 181
 182       // look for the first interference
 183       // --------------------------------
 184       bool intersection = false;
 185       while ( !intersection && it < tgt_nb_cells[j] && is < src_nb_cells[j] )
 186       {
 187         x1s = src_coords[ j ][ is ];
 188         x2t = tgt_coords[ j ][ it+1 ];
 189         if ( x2t < x1s+eps )
 190         {
 191           it++; // source lays on the right of target
 192           continue;
 193         }
 194         x1t = tgt_coords[ j ][ it ];
 195         x2s = src_coords[ j ][ is+1 ];
 196         if ( x2s < x1t+eps )
 197         {
 198           is++; // source lays on the left of target
 199           continue;
 200         }
 201         intersection = true;
 202       }
 203       if ( !intersection ) return ret; // no intersections
 204
 205       // get all interferences
 206       // ----------------------
 207       while ( intersection )
 208       {
 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 ];
 213
 214         double x1 = std::max( x1s ,x1t );
 215         double x2 = std::min( x2s ,x2t );
 216         axis_interferences.push_back( Interference( is, it, x2 - x1 ));
 217
 218         // to the next target and/or source cell
 219         double diff2 = x2s - x2t;
 220         if ( diff2 > -eps )
 221           intersection = ( ++it < tgt_nb_cells[j] );
 222         if ( diff2 < eps )
 223           intersection = ( ++is < src_nb_cells[j] && intersection);
 224       }
 225     }
 226
 227     // ================
 228     // Fill the matrix
 229     // ================
 230
 231     switch ( dim )
 232     {
 233     case 3:
 234       fillMatrix<MyMeshType,MatrixType,3>( interferences, result, src_nb_cells,tgt_nb_cells );
 235       break;
 236
 237     case 2:
 238       fillMatrix<MyMeshType,MatrixType,2>( interferences, result, src_nb_cells,tgt_nb_cells );
 239       break;
 240
 241     case 1:
 242       fillMatrix<MyMeshType,MatrixType,1>( interferences, result, src_nb_cells,tgt_nb_cells );
 243       break;
 244     }
 245
 246     return ret;
 247   }
 248 }