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