src/INTERP_KERNEL/InterpolationCC.txx

   1 // Copyright (C) 2009-2014  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, or (at your option) any later version.
   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   InterpolationCC::InterpolationCC(const InterpolationOptions& io):Interpolation<InterpolationCC>(io)
  44   {
  45   }
  46
  47   //================================================================================
  48   /*!
  49    * \brief An 1D intersection result
  50    */
  51   //================================================================================
  52
  53   struct Interference
  54   {
  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){}
  59   };
  60
  61   //================================================================================
  62   /*!
  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
  74    */
  75   //================================================================================
  76
  77   template <class MyMeshType, class MatrixType, int dim>
  78   void fillMatrix(const std::list< Interference >  inter_of_axis[dim],
  79                   MatrixType&                      result,
  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,
  86                   const int                        axis = 0,
  87                   const double                     prev_value = 1.0)
  88   {
  89     typedef std::list < Interference >::const_iterator TIntIterator;
  90
  91     if ( axis + 1 == dim )
  92     {
  93       for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
  94       {
  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;
  98
  99         result[ tgt_i ].insert( std::make_pair( _TMI( src_i ), value ));
 100       }
 101     }
 102     else
 103     {
 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 ];
 106
 107       for ( TIntIterator i = inter_of_axis[axis].begin(); i != inter_of_axis[axis].end(); ++i )
 108       {
 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;
 112
 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,
 117                                                 axis+1, value );
 118       }
 119     }
 120   }
 121
 122   //================================================================================
 123   /*!
 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
 130    *
 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]
 140    */
 141   //================================================================================
 142
 143   template<class MyMeshType, class MatrixType>
 144   int InterpolationCC::interpolateMeshes(const MyMeshType& src_mesh,
 145                                          const MyMeshType& tgt_mesh,
 146                                          MatrixType&       result,
 147                                          const char *      method)
 148   {
 149     if ( std::string("P0P0") != method )
 150       throw Exception("Only P0P0 method is implemented so far");
 151
 152     // create empty maps for all target elements
 153     result.resize( tgt_mesh.getNumberOfElements() );
 154
 155     const int ret = src_mesh.getNumberOfElements();
 156
 157     const double eps = getPrecision();
 158     const int dim = MyMeshType::MY_MESHDIM;
 159     //const NumberingPolicy numPol = MyMeshType::My_numPol;
 160
 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 )
 166     {
 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 ));
 171     }
 172
 173     // ============================================
 174     // Calculate cell interferences along the axes
 175     // ============================================
 176
 177     std::list < Interference > interferences[ dim ];
 178
 179     for ( int j = 0; j < dim; ++j ) // loop on axes of castesian space
 180     {
 181       std::list < Interference >& axis_interferences = interferences[j];
 182
 183       int it = 0, is = 0;
 184       double x1t, x2t, x1s, x2s; // left and right ordinates of target and source cells
 185
 186       // look for the first interference
 187       // --------------------------------
 188       bool intersection = false;
 189       while ( !intersection && it < tgt_nb_cells[j] && is < src_nb_cells[j] )
 190       {
 191         x1s = src_coords[ j ][ is ];
 192         x2t = tgt_coords[ j ][ it+1 ];
 193         if ( x2t < x1s+eps )
 194         {
 195           it++; // source lays on the right of target
 196           continue;
 197         }
 198         x1t = tgt_coords[ j ][ it ];
 199         x2s = src_coords[ j ][ is+1 ];
 200         if ( x2s < x1t+eps )
 201         {
 202           is++; // source lays on the left of target
 203           continue;
 204         }
 205         intersection = true;
 206       }
 207       if ( !intersection ) return ret; // no intersections
 208
 209       // get all interferences
 210       // ----------------------
 211       while ( intersection )
 212       {
 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 ];
 217
 218         double x1 = std::max( x1s ,x1t );
 219         double x2 = std::min( x2s ,x2t );
 220         axis_interferences.push_back( Interference( is, it, x2 - x1 ));
 221
 222         // to the next target and/or source cell
 223         double diff2 = x2s - x2t;
 224         if ( diff2 > -eps )
 225           intersection = ( ++it < tgt_nb_cells[j] );
 226         if ( diff2 < eps )
 227           intersection = ( ++is < src_nb_cells[j] && intersection);
 228       }
 229     }
 230
 231     // ================
 232     // Fill the matrix
 233     // ================
 234
 235     switch ( dim )
 236     {
 237     case 3:
 238       fillMatrix<MyMeshType,MatrixType,3>( interferences, result, src_nb_cells,tgt_nb_cells );
 239       break;
 240
 241     case 2:
 242       fillMatrix<MyMeshType,MatrixType,2>( interferences, result, src_nb_cells,tgt_nb_cells );
 243       break;
 244
 245     case 1:
 246       fillMatrix<MyMeshType,MatrixType,1>( interferences, result, src_nb_cells,tgt_nb_cells );
 247       break;
 248     }
 249
 250     return ret;
 251   }
 252 }