src/ParaMEDMEM/OverlapInterpolationMatrix.cxx

   1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D
   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 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "OverlapInterpolationMatrix.hxx"
  22 #include "ParaMESH.hxx"
  23 #include "ParaFIELD.hxx"
  24 #include "ProcessorGroup.hxx"
  25 #include "TranslationRotationMatrix.hxx"
  26 #include "Interpolation.hxx"
  27 #include "Interpolation1D.txx"
  28 #include "Interpolation2DCurve.hxx"
  29 #include "Interpolation2D.txx"
  30 #include "Interpolation3DSurf.hxx"
  31 #include "Interpolation3D.txx"
  32 #include "Interpolation2D3D.txx"
  33 #include "Interpolation2D1D.txx"
  34 #include "MEDCouplingUMesh.hxx"
  35 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  36 #include "InterpolationOptions.hxx"
  37 #include "NormalizedUnstructuredMesh.hxx"
  38 #include "ElementLocator.hxx"
  39 #include "InterpKernelAutoPtr.hxx"
  40
  41 #include <algorithm>
  42
  43 using namespace std;
  44
  45 namespace MEDCoupling
  46 {
  47   OverlapInterpolationMatrix::OverlapInterpolationMatrix(ParaFIELD *source_field,
  48                                                          ParaFIELD *target_field,
  49                                                          const ProcessorGroup& group,
  50                                                          const DECOptions& dec_options,
  51                                                          const INTERP_KERNEL::InterpolationOptions& i_opt,
  52                                                          const OverlapElementLocator & locator):
  53     INTERP_KERNEL::InterpolationOptions(i_opt),
  54     DECOptions(dec_options),
  55     _source_field(source_field),
  56     _target_field(target_field),
  57     _source_support(source_field->getSupport()->getCellMesh()),
  58     _target_support(target_field->getSupport()->getCellMesh()),
  59     _mapping(group, locator)
  60   {
  61   }
  62
  63   void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayIdType *ids)
  64   {
  65     _mapping.keepTracksOfSourceIds(procId,ids);
  66   }
  67
  68   void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayIdType *ids)
  69   {
  70     _mapping.keepTracksOfTargetIds(procId,ids);
  71   }
  72
  73   OverlapInterpolationMatrix::~OverlapInterpolationMatrix()
  74   {
  75   }
  76
  77   // TODO? Merge with MEDCouplingRemapper::prepareInterpKernelOnlyUU() ?
  78   /**!
  79    * Local run (on this proc) of the sequential interpolation algorithm.
  80    *
  81    * @param srcIds is null if the source mesh is on the local proc
  82    * @param trgIds is null if the source mesh is on the local proc
  83    *
  84    * One of the 2 is necessarily null (the two can be null together)
  85    */
  86   void OverlapInterpolationMatrix::computeLocalIntersection(const MEDCouplingPointSet *src, const DataArrayIdType *srcIds, const std::string& srcMeth, int srcProcId,
  87                                                    const MEDCouplingPointSet *trg, const DataArrayIdType *trgIds, const std::string& trgMeth, int trgProcId)
  88   {
  89     std::string interpMethod(srcMeth);
  90     interpMethod+=trgMeth;
  91     //creating the interpolator structure
  92     vector<SparseDoubleVec > sparse_matrix_part;
  93     mcIdType colSize=0;
  94     //computation of the intersection volumes between source and target elements
  95     const MEDCouplingUMesh *trgC=dynamic_cast<const MEDCouplingUMesh *>(trg);
  96     const MEDCouplingUMesh *srcC=dynamic_cast<const MEDCouplingUMesh *>(src);
  97     if ( src->getMeshDimension() == -1 )
  98       {
  99         if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==2)
 100           {
 101             MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC);
 102             INTERP_KERNEL::Interpolation2D interpolation(*this);
 103             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth);
 104           }
 105         else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3)
 106           {
 107             MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC);
 108             INTERP_KERNEL::Interpolation3D interpolation(*this);
 109             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth);
 110           }
 111         else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3)
 112           {
 113             MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC);
 114             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 115             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth);
 116           }
 117         else
 118           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 119       }
 120     else if ( trg->getMeshDimension() == -1 )
 121       {
 122         if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==2)
 123           {
 124             MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC);
 125             INTERP_KERNEL::Interpolation2D interpolation(*this);
 126             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth);
 127           }
 128         else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3)
 129           {
 130             MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC);
 131             INTERP_KERNEL::Interpolation3D interpolation(*this);
 132             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth);
 133           }
 134         else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3)
 135           {
 136             MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC);
 137             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 138             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth);
 139           }
 140         else
 141           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
 142       }
 143     else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 3
 144               && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
 145       {
 146         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
 147         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
 148
 149         INTERP_KERNEL::Interpolation2D3D interpolator (*this);
 150         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 151       }
 152     else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2
 153               && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
 154       {
 155         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
 156         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
 157
 158         INTERP_KERNEL::Interpolation2D3D interpolator (*this);
 159         vector<SparseDoubleVec > matrixTranspose;
 160         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,sparse_matrix_part,interpMethod);//not a bug target in source.
 161         TransposeMatrix(matrixTranspose,colSize,sparse_matrix_part);
 162         colSize=ToIdType(matrixTranspose.size());
 163       }
 164     else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2
 165               && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
 166       {
 167         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
 168         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
 169
 170         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
 171         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 172       }
 173     else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 1
 174               && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
 175       {
 176         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
 177         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
 178
 179         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
 180         vector<SparseDoubleVec > matrixTranspose;
 181         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,matrixTranspose,interpMethod);//not a bug target in source.
 182         TransposeMatrix(matrixTranspose,colSize,sparse_matrix_part);
 183         colSize=ToIdType(matrixTranspose.size());
 184       }
 185     else if (trg->getMeshDimension() != _source_support->getMeshDimension())
 186       {
 187         throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
 188       }
 189     else if( src->getMeshDimension() == 1
 190              && src->getSpaceDimension() == 1 )
 191       {
 192         MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(trgC);
 193         MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC);
 194
 195         INTERP_KERNEL::Interpolation1D interpolation(*this);
 196         colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 197       }
 198     else if( trg->getMeshDimension() == 1
 199              && trg->getSpaceDimension() == 2 )
 200       {
 201         MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(trgC);
 202         MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC);
 203
 204         INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
 205         colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 206       }
 207     else if ( trg->getMeshDimension() == 2
 208               && trg->getSpaceDimension() == 3 )
 209       {
 210         MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(trgC);
 211         MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC);
 212
 213         INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
 214         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 215       }
 216     else if ( trg->getMeshDimension() == 2
 217               && trg->getSpaceDimension() == 2)
 218       {
 219         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
 220         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
 221
 222         INTERP_KERNEL::Interpolation2D interpolator (*this);
 223         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 224       }
 225     else if ( trg->getMeshDimension() == 3
 226               && trg->getSpaceDimension() == 3 )
 227       {
 228         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
 229         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
 230
 231         INTERP_KERNEL::Interpolation3D interpolator (*this);
 232         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
 233       }
 234     else
 235       {
 236         throw INTERP_KERNEL::Exception("No interpolator exists for these mesh and space dimensions!");
 237       }
 238     /* Fill distributed matrix:
 239        In sparse_matrix_part rows refer to target, and columns (=first param of map in SparseDoubleVec)
 240        refer to source.
 241      */
 242     _mapping.addContributionST(sparse_matrix_part,srcIds,srcProcId,trgIds,trgProcId);
 243   }
 244
 245   /*!
 246    * 'procsToSendField' gives the list of procs field data has to be sent to.
 247    * See OverlapElementLocator::computeBoundingBoxesAndTodoList()
 248    */
 249   void OverlapInterpolationMatrix::prepare(const std::vector< int >& procsToSendField)
 250   {
 251     if(_source_support)
 252       _mapping.prepare(procsToSendField,_target_field->getField()->getNumberOfTuplesExpected());
 253     else
 254       _mapping.prepare(procsToSendField,0);
 255   }
 256
 257   void OverlapInterpolationMatrix::computeSurfacesAndDeno()
 258   {
 259     if(_target_field->getField()->getNature()==IntensiveMaximum)
 260       _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected());
 261     else
 262       throw INTERP_KERNEL::Exception("OverlapDEC: Policy not implemented yet: only IntensiveMaximum!");
 263 //      {
 264 //      if(_target_field->getField()->getNature()==IntensiveConservation)
 265 //        {
 266 //          MCAuto<MEDCouplingFieldDouble> f;
 267 //          int orient = getOrientation(); // From InterpolationOptions inheritance
 268 //          if(orient == 2)  // absolute areas
 269 //             f = _target_support->getMeasureField(true);
 270 //          else
 271 //            if(orient == 0)  // relative areas
 272 //              f = _target_support->getMeasureField(false);
 273 //            else
 274 //              throw INTERP_KERNEL::Exception("OverlapDEC: orientation policy not impl. yet!");
 275 //          _mapping.computeDenoRevIntegral(*f->getArray());
 276 //        }
 277 //      else
 278 //        throw INTERP_KERNEL::Exception("OverlapDEC: Policy not implemented yet: only IntensiveMaximum and IntensiveConservation defined!");
 279 //      }
 280   }
 281
 282   void OverlapInterpolationMatrix::multiply(double default_val)
 283   {
 284     _mapping.multiply(_source_field->getField(),_target_field->getField(), default_val);
 285   }
 286
 287   void OverlapInterpolationMatrix::transposeMultiply()
 288   {
 289     _mapping.transposeMultiply(_target_field->getField(),_source_field->getField());
 290   }
 291
 292   void OverlapInterpolationMatrix::TransposeMatrix(const std::vector<SparseDoubleVec >& matIn,
 293                                                    mcIdType nbColsMatIn, std::vector<SparseDoubleVec >& matOut)
 294   {
 295     matOut.resize(nbColsMatIn);
 296     int id=0;
 297     for(std::vector<SparseDoubleVec >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
 298       for(SparseDoubleVec::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 299         matOut[(*iter2).first][id]=(*iter2).second;
 300   }
 301 }