MED33 porting

[tools/medcoupling.git] / src / ParaMEDMEM / OverlapDEC.cxx

diff --git a/src/ParaMEDMEM/OverlapDEC.cxx b/src/ParaMEDMEM/OverlapDEC.cxx
index c799c23aaeff315b73718fea1fc82766e7873943..75cbee3085ccb68d9d6859854bb5a43ec90407d1 100644 (file)
--- a/src/ParaMEDMEM/OverlapDEC.cxx
+++ b/src/ParaMEDMEM/OverlapDEC.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D

 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public


@@ -20,36 +20,42 @@
 
 #include "OverlapDEC.hxx"
 #include "CommInterface.hxx"
 
 #include "OverlapDEC.hxx"
 #include "CommInterface.hxx"

+#include "ParaMESH.hxx"

 #include "ParaFIELD.hxx"
 #include "MPIProcessorGroup.hxx"
 #include "OverlapElementLocator.hxx"
 #include "OverlapInterpolationMatrix.hxx"
 #include "ParaFIELD.hxx"
 #include "MPIProcessorGroup.hxx"
 #include "OverlapElementLocator.hxx"
 #include "OverlapInterpolationMatrix.hxx"

+#include "ICoCoMEDField.hxx"

 
 

-namespace ParaMEDMEM
+namespace MEDCoupling

 {
 /*!
     \anchor OverlapDEC-det
     \class OverlapDEC
 
 {
 /*!
     \anchor OverlapDEC-det
     \class OverlapDEC
 

+    \section OverlapDEC-over Overview
+

     The \c OverlapDEC enables the \ref InterpKerRemapGlobal "conservative remapping" of fields between
     two parallel codes. This remapping is based on the computation of intersection volumes on
     The \c OverlapDEC enables the \ref InterpKerRemapGlobal "conservative remapping" of fields between
     two parallel codes. This remapping is based on the computation of intersection volumes on

-    a \b same \b processor \b group. On this processor group are defined two field-templates called A
+    a \b single \b processor \b group. On this processor group are defined two field-templates called A

     and B. The computation is possible for 3D meshes, 2D meshes, 3D-surface meshes, 1D meshes and
     2D-curve meshes. Dimensions must be similar for the distribution templates A and B.
     and B. The computation is possible for 3D meshes, 2D meshes, 3D-surface meshes, 1D meshes and
     2D-curve meshes. Dimensions must be similar for the distribution templates A and B.

-    The main difference with \ref InterpKernelDEC-det is that this \ref para-dec "DEC" manages 2 field templates
-    on each processor of the processor group (A and B) called source and target.
-    Furthermore all processors in processor group cooperates in global interpolation matrix
-    computation. In this respect \ref InterpKernelDEC is a specialization of \c OverlapDEC.

 
 

-    \section ParaMEDMEMOverlapDECAlgorithmDescription Algorithm Description
+    The main difference with \ref InterpKernelDEC-det "InterpKernelDEC" is that this
+    \ref para-dec "DEC" works with a *single* processor group, in which processors will share the work.
+    Consequently each processor manages two \ref MEDCouplingFieldTemplatesPage "field templates" (A and B)
+    called source and target.
+    Furthermore all processors in the processor group cooperate in the global interpolation matrix
+    computation. In this respect \c InterpKernelDEC is a specialization of \c OverlapDEC.
+
+    \section ParaMEDMEMOverlapDECAlgorithmDescription Algorithm description

 
     Let's consider the following use case that is ran in ParaMEDMEMTest_OverlapDEC.cxx to describes
     the different steps of the computation. The processor group contains 3 processors.
     \anchor ParaMEDMEMOverlapDECImgTest1
 
     Let's consider the following use case that is ran in ParaMEDMEMTest_OverlapDEC.cxx to describes
     the different steps of the computation. The processor group contains 3 processors.
     \anchor ParaMEDMEMOverlapDECImgTest1

-    \image html OverlapDEC1.png "Example showing the use case in order to explain the different steps."
+    \image html OverlapDEC1.png "Example split of the source and target mesh among the 3 procs"

 
 

-    \subsection ParaMEDMEMOverlapDECAlgoStep1 Step 1 : Bounding box exchange and global interaction
-    between procs computation.
+    \subsection ParaMEDMEMOverlapDECAlgoStep1 Step 1 : Bounding box exchange and global interaction between procs computation.

 
     In order to reduce as much as possible the amount of communications between distant processors,
     every processor computes a bounding box for A and B. Then a AllToAll communication is performed
 
     In order to reduce as much as possible the amount of communications between distant processors,
     every processor computes a bounding box for A and B. Then a AllToAll communication is performed


@@ -66,7 +72,7 @@ namespace ParaMEDMEM
     Here the pair (0,2) does not appear because the bounding box of fieldtemplateA of proc#2 does
     not intersect that of fieldtemplate B on proc#0.
 
     Here the pair (0,2) does not appear because the bounding box of fieldtemplateA of proc#2 does
     not intersect that of fieldtemplate B on proc#0.
 

-    Stage performed by ParaMEDMEM::OverlapElementLocator::computeBoundingBoxes.
+    Stage performed by MEDCoupling::OverlapElementLocator::computeBoundingBoxes.

 
     \subsection ParaMEDMEMOverlapDECAlgoStep2 Step 2 : Computation of local TODO list
 
 
     \subsection ParaMEDMEMOverlapDECAlgoStep2 Step 2 : Computation of local TODO list
 


@@ -151,7 +157,7 @@ namespace ParaMEDMEM
     keep track of the ids sent to proc \#m for te matrix-vector computation.
     This is incarnated by OverlapMapping::keepTracksOfSourceIds in proc k.
 
     keep track of the ids sent to proc \#m for te matrix-vector computation.
     This is incarnated by OverlapMapping::keepTracksOfSourceIds in proc k.
 

-    This step is performed in ParaMEDMEM::OverlapElementLocator::exchangeMeshes method.
+    This step is performed in MEDCoupling::OverlapElementLocator::exchangeMeshes method.

 
     \subsection ParaMEDMEMOverlapDECAlgoStep4 Step 4 : Computation of the interpolation matrix
 
 
     \subsection ParaMEDMEMOverlapDECAlgoStep4 Step 4 : Computation of the interpolation matrix
 


@@ -161,7 +167,7 @@ namespace ParaMEDMEM
     the \b local TODO list per proc is expected to
     be as well balanced as possible.
 
     the \b local TODO list per proc is expected to
     be as well balanced as possible.
 

-    The interpolation is performed as \ref ParaMEDMEM::MEDCouplingRemapper "Remapper" does.
+    The interpolation is performed as the \ref MEDCoupling::MEDCouplingRemapper "remapper" does.

 
     This operation is performed by OverlapInterpolationMatrix::addContribution method.
 
 
     This operation is performed by OverlapInterpolationMatrix::addContribution method.
 


@@ -179,7 +185,7 @@ namespace ParaMEDMEM
     is equal to k.
 
     After this step, the matrix repartition is the following after a call to
     is equal to k.
 
     After this step, the matrix repartition is the following after a call to

-    ParaMEDMEM::OverlapMapping::prepare :
+    MEDCoupling::OverlapMapping::prepare :

 
     - proc\#0 : (0,0),(1,0),(2,0)
     - proc\#1 : (0,1),(2,1)
 
     - proc\#0 : (0,0),(1,0),(2,0)
     - proc\#1 : (0,1),(2,1)


@@ -189,35 +195,39 @@ namespace ParaMEDMEM
     "prepare". This is an example of item 0 in \ref ParaMEDMEMOverlapDECAlgoStep2 "Step2".
     Tuple (0,1) computed on proc 1 is stored in proc 1 too. This is an example of item 1 in \ref ParaMEDMEMOverlapDECAlgoStep2 "Step2".
 
     "prepare". This is an example of item 0 in \ref ParaMEDMEMOverlapDECAlgoStep2 "Step2".
     Tuple (0,1) computed on proc 1 is stored in proc 1 too. This is an example of item 1 in \ref ParaMEDMEMOverlapDECAlgoStep2 "Step2".
 

-    In the end ParaMEDMEM::OverlapMapping::_proc_ids_to_send_vector_st will contain :
+    In the end MEDCoupling::OverlapMapping::_proc_ids_to_send_vector_st will contain :

 
     - Proc\#0 : 0,1
     - Proc\#1 : 0,2
     - Proc\#2 : 0,1,2
 
 
     - Proc\#0 : 0,1
     - Proc\#1 : 0,2
     - Proc\#2 : 0,1,2
 

-    In the end ParaMEDMEM::OverlapMapping::_proc_ids_to_recv_vector_st will contain :
+    In the end MEDCoupling::OverlapMapping::_proc_ids_to_recv_vector_st will contain :

 
     - Proc\#0 : 0,1,2
     - Proc\#1 : 0,2
     - Proc\#2 : 1,2
 
 
     - Proc\#0 : 0,1,2
     - Proc\#1 : 0,2
     - Proc\#2 : 1,2
 

-    The method in charge to perform this is : ParaMEDMEM::OverlapMapping::prepare.
+    The method in charge to perform this is : MEDCoupling::OverlapMapping::prepare.

 */
 */

-  OverlapDEC::OverlapDEC(const std::set<int>& procIds, const MPI_Comm& world_comm):_own_group(true),_interpolation_matrix(0),
-                                                                                   _source_field(0),_own_source_field(false),
-                                                                                   _target_field(0),_own_target_field(false)
+  OverlapDEC::OverlapDEC(const std::set<int>& procIds, const MPI_Comm& world_comm):
+      _load_balancing_algo(1),
+      _own_group(true),_interpolation_matrix(0), _locator(0),
+      _source_field(0),_own_source_field(false),
+      _target_field(0),_own_target_field(false),
+      _default_field_value(0.0),
+      _comm(MPI_COMM_NULL)

   {
   {

-    ParaMEDMEM::CommInterface comm;
+    MEDCoupling::CommInterface comm;

     int *ranks_world=new int[procIds.size()]; // ranks of sources and targets in world_comm
     std::copy(procIds.begin(),procIds.end(),ranks_world);
     MPI_Group group,world_group;
     comm.commGroup(world_comm,&world_group);
     comm.groupIncl(world_group,procIds.size(),ranks_world,&group);
     delete [] ranks_world;
     int *ranks_world=new int[procIds.size()]; // ranks of sources and targets in world_comm
     std::copy(procIds.begin(),procIds.end(),ranks_world);
     MPI_Group group,world_group;
     comm.commGroup(world_comm,&world_group);
     comm.groupIncl(world_group,procIds.size(),ranks_world,&group);
     delete [] ranks_world;

-    MPI_Comm theComm;
-    comm.commCreate(world_comm,group,&theComm);
+    comm.commCreate(world_comm,group,&_comm);

     comm.groupFree(&group);
     comm.groupFree(&group);

-    if(theComm==MPI_COMM_NULL)
+    comm.groupFree(&world_group);
+    if(_comm==MPI_COMM_NULL)

       {
         _group=0;
         return ;
       {
         _group=0;
         return ;


@@ -225,7 +235,7 @@ namespace ParaMEDMEM
     std::set<int> idsUnion;
     for(std::size_t i=0;i<procIds.size();i++)
       idsUnion.insert(i);
     std::set<int> idsUnion;
     for(std::size_t i=0;i<procIds.size();i++)
       idsUnion.insert(i);

-    _group=new MPIProcessorGroup(comm,idsUnion,theComm);
+    _group=new MPIProcessorGroup(comm,idsUnion,_comm);

   }
 
   OverlapDEC::~OverlapDEC()
   }
 
   OverlapDEC::~OverlapDEC()


@@ -237,6 +247,12 @@ namespace ParaMEDMEM
     if(_own_target_field)
       delete _target_field;
     delete _interpolation_matrix;
     if(_own_target_field)
       delete _target_field;
     delete _interpolation_matrix;

+    delete _locator;
+    if (_comm != MPI_COMM_NULL)
+      {
+        MEDCoupling::CommInterface comm;
+        comm.commFree(&_comm);
+      }

   }
 
   void OverlapDEC::sendRecvData(bool way)
   }
 
   void OverlapDEC::sendRecvData(bool way)


@@ -249,7 +265,7 @@ namespace ParaMEDMEM
 
   void OverlapDEC::sendData()
   {
 
   void OverlapDEC::sendData()
   {

-    _interpolation_matrix->multiply();
+    _interpolation_matrix->multiply(_default_field_value);

   }
 
   void OverlapDEC::recvData()
   }
 
   void OverlapDEC::recvData()


@@ -262,24 +278,31 @@ namespace ParaMEDMEM
   {
     if(!isInGroup())
       return ;
   {
     if(!isInGroup())
       return ;

+    // Check number of components of field on both side (for now allowing void field/mesh on one proc is not allowed)
+    if (!_source_field || !_source_field->getField())
+      throw INTERP_KERNEL::Exception("OverlapDEC::synchronize(): currently, having a void source field on a proc is not allowed!");
+    if (!_target_field || !_target_field->getField())
+      throw INTERP_KERNEL::Exception("OverlapDEC::synchronize(): currently, having a void target field on a proc is not allowed!");
+    if (_target_field->getField()->getNumberOfComponents() != _source_field->getField()->getNumberOfComponents())
+      throw INTERP_KERNEL::Exception("OverlapDEC::synchronize(): source and target field have different number of components!");

     delete _interpolation_matrix;
     delete _interpolation_matrix;

-    _interpolation_matrix=new OverlapInterpolationMatrix(_source_field,_target_field,*_group,*this,*this);
-    OverlapElementLocator locator(_source_field,_target_field,*_group);
-    locator.copyOptions(*this);
-    locator.exchangeMeshes(*_interpolation_matrix);
-    std::vector< std::pair<int,int> > jobs=locator.getToDoList();
-    std::string srcMeth=locator.getSourceMethod();
-    std::string trgMeth=locator.getTargetMethod();
+    _locator = new OverlapElementLocator(_source_field,_target_field,*_group, getBoundingBoxAdjustmentAbs(), _load_balancing_algo);
+    _interpolation_matrix=new OverlapInterpolationMatrix(_source_field,_target_field,*_group,*this,*this, *_locator);
+    _locator->copyOptions(*this);
+    _locator->exchangeMeshes(*_interpolation_matrix);
+    std::vector< std::pair<int,int> > jobs=_locator->getToDoList();
+    std::string srcMeth=_locator->getSourceMethod();
+    std::string trgMeth=_locator->getTargetMethod();

     for(std::vector< std::pair<int,int> >::const_iterator it=jobs.begin();it!=jobs.end();it++)
       {
     for(std::vector< std::pair<int,int> >::const_iterator it=jobs.begin();it!=jobs.end();it++)
       {

-        const MEDCouplingPointSet *src=locator.getSourceMesh((*it).first);
-        const DataArrayInt *srcIds=locator.getSourceIds((*it).first);
-        const MEDCouplingPointSet *trg=locator.getTargetMesh((*it).second);
-        const DataArrayInt *trgIds=locator.getTargetIds((*it).second);
-        _interpolation_matrix->addContribution(src,srcIds,srcMeth,(*it).first,trg,trgIds,trgMeth,(*it).second);
+        const MEDCouplingPointSet *src=_locator->getSourceMesh((*it).first);
+        const DataArrayInt *srcIds=_locator->getSourceIds((*it).first);
+        const MEDCouplingPointSet *trg=_locator->getTargetMesh((*it).second);
+        const DataArrayInt *trgIds=_locator->getTargetIds((*it).second);
+        _interpolation_matrix->computeLocalIntersection(src,srcIds,srcMeth,(*it).first,trg,trgIds,trgMeth,(*it).second);

       }
       }

-    _interpolation_matrix->prepare(locator.getProcsInInteraction());
-    _interpolation_matrix->computeDeno();
+    _interpolation_matrix->prepare(_locator->getProcsToSendFieldData());
+    _interpolation_matrix->computeSurfacesAndDeno();

   }
 
   void OverlapDEC::attachSourceLocalField(ParaFIELD *field, bool ownPt)
   }
 
   void OverlapDEC::attachSourceLocalField(ParaFIELD *field, bool ownPt)


@@ -302,10 +325,49 @@ namespace ParaMEDMEM
     _own_target_field=ownPt;
   }
 
     _own_target_field=ownPt;
   }
 

+  void OverlapDEC::attachSourceLocalField(MEDCouplingFieldDouble *field)
+  {
+    if(!isInGroup())
+      return ;
+
+    ParaMESH *paramesh = new ParaMESH(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(field->getMesh())),
+                                      *_group,field->getMesh()->getName());
+    ParaFIELD *tmpField=new ParaFIELD(field, paramesh, *_group);
+    tmpField->setOwnSupport(true);
+    attachSourceLocalField(tmpField,true);
+  }
+
+  void OverlapDEC::attachTargetLocalField(MEDCouplingFieldDouble *field)
+  {
+    if(!isInGroup())
+      return ;
+
+    ParaMESH *paramesh = new ParaMESH(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(field->getMesh())),
+                                      *_group,field->getMesh()->getName());
+    ParaFIELD *tmpField=new ParaFIELD(field, paramesh, *_group);
+    tmpField->setOwnSupport(true);
+    attachTargetLocalField(tmpField,true);
+  }
+
+  void OverlapDEC::attachSourceLocalField(ICoCo::MEDField *field)
+  {
+    attachSourceLocalField(field->getField());
+  }
+
+  void OverlapDEC::attachTargetLocalField(ICoCo::MEDField *field)
+  {
+    attachTargetLocalField(field->getField());
+  }
+

   bool OverlapDEC::isInGroup() const
   {
     if(!_group)
       return false;
     return _group->containsMyRank();
   }
   bool OverlapDEC::isInGroup() const
   {
     if(!_group)
       return false;
     return _group->containsMyRank();
   }

+
+  void OverlapDEC::debugPrintWorkSharing(std::ostream & ostr) const
+  {
+    _locator->debugPrintWorkSharing(ostr);
+  }

 }
 }