src/ParaMEDMEM/OverlapInterpolationMatrix.cxx

   1 // Copyright (C) 2007-2015  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 "Interpolation3D2D.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 ParaMEDMEM
  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     INTERP_KERNEL::InterpolationOptions(i_opt),
  53     DECOptions(dec_options),
  54     _source_field(source_field),
  55     _target_field(target_field),
  56     _source_support(source_field->getSupport()->getCellMesh()),
  57     _target_support(target_field->getSupport()->getCellMesh()),
  58     _mapping(group),
  59     _group(group)
  60   {
  61     int nbelems = source_field->getField()->getNumberOfTuples();
  62     _row_offsets.resize(nbelems+1);
  63     _coeffs.resize(nbelems);
  64     _target_volume.resize(nbelems);
  65   }
  66
  67   void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayInt *ids)
  68   {
  69     _mapping.keepTracksOfSourceIds(procId,ids);
  70   }
  71
  72   void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayInt *ids)
  73   {
  74     _mapping.keepTracksOfTargetIds(procId,ids);
  75   }
  76
  77   OverlapInterpolationMatrix::~OverlapInterpolationMatrix()
  78   {
  79   }
  80
  81   void OverlapInterpolationMatrix::addContribution(const MEDCouplingPointSet *src, const DataArrayInt *srcIds, const std::string& srcMeth, int srcProcId,
  82                                                    const MEDCouplingPointSet *trg, const DataArrayInt *trgIds, const std::string& trgMeth, int trgProcId)
  83   {
  84     std::string interpMethod(srcMeth);
  85     interpMethod+=trgMeth;
  86     //creating the interpolator structure
  87     vector<map<int,double> > surfaces;
  88     int colSize=0;
  89     //computation of the intersection volumes between source and target elements
  90     const MEDCouplingUMesh *trgC=dynamic_cast<const MEDCouplingUMesh *>(trg);
  91     const MEDCouplingUMesh *srcC=dynamic_cast<const MEDCouplingUMesh *>(src);
  92     if ( src->getMeshDimension() == -1 )
  93       {
  94         if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==2)
  95           {
  96             MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC);
  97             INTERP_KERNEL::Interpolation2D interpolation(*this);
  98             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth);
  99           }
 100         else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3)
 101           {
 102             MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC);
 103             INTERP_KERNEL::Interpolation3D interpolation(*this);
 104             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth);
 105           }
 106         else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3)
 107           {
 108             MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC);
 109             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 110             colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth);
 111           }
 112         else
 113           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 114       }
 115     else if ( trg->getMeshDimension() == -1 )
 116       {
 117         if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==2)
 118           {
 119             MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC);
 120             INTERP_KERNEL::Interpolation2D interpolation(*this);
 121             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth);
 122           }
 123         else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3)
 124           {
 125             MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC);
 126             INTERP_KERNEL::Interpolation3D interpolation(*this);
 127             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth);
 128           }
 129         else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3)
 130           {
 131             MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC);
 132             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 133             colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth);
 134           }
 135         else
 136           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
 137       }
 138     else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 3
 139               && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
 140       {
 141         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
 142         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
 143
 144         INTERP_KERNEL::Interpolation3D2D interpolator (*this);
 145         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 146         target_wrapper.releaseTempArrays();
 147         source_wrapper.releaseTempArrays();
 148       }
 149     else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2
 150               && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
 151       {
 152         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
 153         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
 154
 155         INTERP_KERNEL::Interpolation3D2D interpolator (*this);
 156         vector<map<int,double> > surfacesTranspose;
 157         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);//not a bug target in source.
 158         TransposeMatrix(surfacesTranspose,colSize,surfaces);
 159         colSize=surfacesTranspose.size();
 160         target_wrapper.releaseTempArrays();
 161         source_wrapper.releaseTempArrays();
 162       }
 163     else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2
 164               && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
 165       {
 166         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
 167         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
 168
 169         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
 170         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 171         target_wrapper.releaseTempArrays();
 172         source_wrapper.releaseTempArrays();
 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<map<int,double> > surfacesTranspose;
 182         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfacesTranspose,interpMethod);//not a bug target in source.
 183         TransposeMatrix(surfacesTranspose,colSize,surfaces);
 184         colSize=surfacesTranspose.size();
 185         target_wrapper.releaseTempArrays();
 186         source_wrapper.releaseTempArrays();
 187       }
 188     else if (trg->getMeshDimension() != _source_support->getMeshDimension())
 189       {
 190         throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
 191       }
 192     else if( src->getMeshDimension() == 1
 193              && src->getSpaceDimension() == 1 )
 194       {
 195         MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(trgC);
 196         MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC);
 197
 198         INTERP_KERNEL::Interpolation1D interpolation(*this);
 199         colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 200         target_wrapper.releaseTempArrays();
 201         source_wrapper.releaseTempArrays();
 202       }
 203     else if( trg->getMeshDimension() == 1
 204              && trg->getSpaceDimension() == 2 )
 205       {
 206         MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(trgC);
 207         MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC);
 208
 209         INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
 210         colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 211         target_wrapper.releaseTempArrays();
 212         source_wrapper.releaseTempArrays();
 213       }
 214     else if ( trg->getMeshDimension() == 2
 215               && trg->getSpaceDimension() == 3 )
 216       {
 217         MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(trgC);
 218         MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC);
 219
 220         INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
 221         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 222         target_wrapper.releaseTempArrays();
 223         source_wrapper.releaseTempArrays();
 224       }
 225     else if ( trg->getMeshDimension() == 2
 226               && trg->getSpaceDimension() == 2)
 227       {
 228         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
 229         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
 230
 231         INTERP_KERNEL::Interpolation2D interpolator (*this);
 232         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 233         target_wrapper.releaseTempArrays();
 234         source_wrapper.releaseTempArrays();
 235       }
 236     else if ( trg->getMeshDimension() == 3
 237               && trg->getSpaceDimension() == 3 )
 238       {
 239         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
 240         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
 241
 242         INTERP_KERNEL::Interpolation3D interpolator (*this);
 243         colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod);
 244         target_wrapper.releaseTempArrays();
 245         source_wrapper.releaseTempArrays();
 246       }
 247     else
 248       {
 249         throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
 250       }
 251     bool needSourceSurf=isSurfaceComputationNeeded(srcMeth);
 252     MEDCouplingFieldDouble *source_triangle_surf=0;
 253     if(needSourceSurf)
 254       source_triangle_surf=src->getMeasureField(getMeasureAbsStatus());
 255     //
 256     fillDistributedMatrix(surfaces,srcIds,srcProcId,trgIds,trgProcId);
 257     //
 258     if(needSourceSurf)
 259       source_triangle_surf->decrRef();
 260   }
 261
 262   /*!
 263    * \b res rows refers to target and column (first param of map) to source.
 264    */
 265   void OverlapInterpolationMatrix::fillDistributedMatrix(const std::vector< std::map<int,double> >& res,
 266                                                          const DataArrayInt *srcIds, int srcProc,
 267                                                          const DataArrayInt *trgIds, int trgProc)
 268   {
 269     _mapping.addContributionST(res,srcIds,srcProc,trgIds,trgProc);
 270   }
 271
 272   /*!
 273    * 'procsInInteraction' gives the global view of interaction between procs.
 274    * In 'procsInInteraction' for a proc with id i, is in interaction with procs listed in procsInInteraction[i]
 275    */
 276   void OverlapInterpolationMatrix::prepare(const std::vector< std::vector<int> >& procsInInteraction)
 277   {
 278     if(_source_support)
 279       _mapping.prepare(procsInInteraction,_target_field->getField()->getNumberOfTuplesExpected());
 280     else
 281       _mapping.prepare(procsInInteraction,0);
 282   }
 283
 284   void OverlapInterpolationMatrix::computeDeno()
 285   {
 286     if(_target_field->getField()->getNature()==ConservativeVolumic)
 287       _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected());
 288     else
 289       throw INTERP_KERNEL::Exception("Policy Not implemented yet : only ConservativeVolumic defined !");
 290   }
 291
 292   void OverlapInterpolationMatrix::multiply()
 293   {
 294     _mapping.multiply(_source_field->getField(),_target_field->getField());
 295   }
 296
 297   void OverlapInterpolationMatrix::transposeMultiply()
 298   {
 299     _mapping.transposeMultiply(_target_field->getField(),_source_field->getField());
 300   }
 301
 302   bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
 303   {
 304     return method=="P0";
 305   }
 306
 307   void OverlapInterpolationMatrix::TransposeMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
 308   {
 309     matOut.resize(nbColsMatIn);
 310     int id=0;
 311     for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
 312       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 313         matOut[(*iter2).first][id]=(*iter2).second;
 314   }
 315 }