src/ParaMEDMEM/OverlapInterpolationMatrix.cxx

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