src/ParaMEDMEM/OverlapInterpolationMatrix.cxx

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