src/INTERP_KERNEL/InterpolationCC.txx

   1 // Copyright (C) 2009-2020  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, int dim, class MatrixType, class ConnType>
  78   void fillMatrix(const std::list< Interference >  inter_of_axis[dim],
  79                   MatrixType&                      result,
  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,
  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         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;
  98
  99         result[ tgt_i ].insert( std::make_pair( _TMI( src_i ), value ));
 100       }
 101     }
 102     else
 103     {
 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 ];
 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         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;
 112
 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,
 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   typename MyMeshType::MyConnType 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     typedef typename MyMeshType::MyConnType ConnType;
 156     const ConnType ret = src_mesh.getNumberOfElements();
 157
 158     const double eps = getPrecision();
 159     const int dim = MyMeshType::MY_MESHDIM;
 160     //const NumberingPolicy numPol = MyMeshType::My_numPol;
 161
 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 )
 167     {
 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 );
 173     }
 174
 175     // ============================================
 176     // Calculate cell interferences along the axes
 177     // ============================================
 178
 179     std::list < Interference > interferences[ dim ];
 180
 181     for ( int j = 0; j < dim; ++j ) // loop on axes of castesian space
 182     {
 183       std::list < Interference >& axis_interferences = interferences[j];
 184
 185       int it = 0, is = 0;
 186       double x1t, x2t, x1s, x2s; // left and right ordinates of target and source cells
 187
 188       // look for the first interference
 189       // --------------------------------
 190       bool intersection = false;
 191       while ( !intersection && it < tgt_nb_cells[j] && is < src_nb_cells[j] )
 192       {
 193         x1s = src_coords[ j ][ is ];
 194         x2t = tgt_coords[ j ][ it+1 ];
 195         if ( x2t < x1s+eps )
 196         {
 197           it++; // source lays on the right of target
 198           continue;
 199         }
 200         x1t = tgt_coords[ j ][ it ];
 201         x2s = src_coords[ j ][ is+1 ];
 202         if ( x2s < x1t+eps )
 203         {
 204           is++; // source lays on the left of target
 205           continue;
 206         }
 207         intersection = true;
 208       }
 209       if ( !intersection ) return ret; // no intersections
 210
 211       // get all interferences
 212       // ----------------------
 213       while ( intersection )
 214       {
 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 ];
 219
 220         double x1 = std::max( x1s ,x1t );
 221         double x2 = std::min( x2s ,x2t );
 222         axis_interferences.push_back( Interference( is, it, x2 - x1 ));
 223
 224         // to the next target and/or source cell
 225         double diff2 = x2s - x2t;
 226         if ( diff2 > -eps )
 227           intersection = ( ++it < tgt_nb_cells[j] );
 228         if ( diff2 < eps )
 229           intersection = ( ++is < src_nb_cells[j] && intersection);
 230       }
 231     }
 232
 233     // ================
 234     // Fill the matrix
 235     // ================
 236
 237     switch ( dim )
 238     {
 239     case 3:
 240       fillMatrix<MyMeshType,3>( interferences, result, src_nb_cells,tgt_nb_cells );
 241       break;
 242
 243     case 2:
 244       fillMatrix<MyMeshType,2>( interferences, result, src_nb_cells,tgt_nb_cells );
 245       break;
 246
 247     case 1:
 248       fillMatrix<MyMeshType,1>( interferences, result, src_nb_cells,tgt_nb_cells );
 249       break;
 250     }
 251
 252     return ret;
 253   }
 254 }