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