]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM/DisjointDEC.cxx
Salome HOME
typo-fix by Kunda
[tools/medcoupling.git] / src / ParaMEDMEM / DisjointDEC.cxx
1 // Copyright (C) 2007-2016  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
20 #include "DisjointDEC.hxx"
21 #include "CommInterface.hxx"
22 #include "Topology.hxx"
23 #include "BlockTopology.hxx"
24 #include "ComponentTopology.hxx"
25 #include "ParaFIELD.hxx"
26 #include "ParaMESH.hxx"
27 #include "ICoCoField.hxx"
28 #include "ICoCoMEDField.hxx"
29 #include "MPIProcessorGroup.hxx"
30
31 #include <cmath>
32 #include <iostream>
33
34
35 namespace MEDCoupling
36 {
37
38   /*!
39    * \anchor DisjointDEC-det
40    * \class DisjointDEC
41    *
42    * \section DisjointDEC-over Overview
43    *
44    * Abstract interface class representing a link between two
45    * processor groups for exchanging mesh or field data. The two processor groups must
46    * have a void intersection (\ref MEDCoupling::OverlapDEC "OverlapDEC" is to be considered otherwise).
47    * The %DEC is initialized by attaching a field on the receiving or on the
48    * sending side.
49    *
50    * The data is sent or received through calls to the (abstract) methods recvData() and sendData().
51    *
52    * One can attach either a \c MEDCoupling::ParaFIELD, or a
53    * \c ICoCo::Field, or directly a \c MEDCoupling::MEDCouplingFieldDouble instance.
54    * See the various signatures of the method DisjointDEC::attachLocalField()
55    *
56    * The derivations of this class should be considered for practical instantiation:
57    * - \ref InterpKernelDEC-det "InterpKernelDEC"
58    * - \ref ExplicitCoincidentDEC-det "ExplicitCoincidentDEC"
59    * - \ref StructuredCoincidentDEC-det "StructuredCoincidentDEC"
60    *
61    * \section DisjointDEC-options DisjointDEC options
62    * The options supported by %DisjointDEC objects are the same that the ones supported for all
63    * DECs in general and are all inherited from the class \ref MEDCoupling::DECOptions "DECOptions"
64    *
65   */
66
67   DisjointDEC::DisjointDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
68       _local_field(0),
69       _source_group(&source_group),
70       _target_group(&target_group),
71       _comm_interface(0),
72       _owns_field(false),
73       _owns_groups(false),
74       _union_comm(MPI_COMM_NULL)
75   {
76     checkPartitionGroup();
77     _union_group = source_group.fuse(target_group);  
78   }
79   
80   DisjointDEC::DisjointDEC(const DisjointDEC& s):
81       DEC(s),
82       _local_field(0),
83       _union_group(0),
84       _source_group(0),
85       _target_group(0),
86       _comm_interface(0),
87       _owns_field(false),
88       _owns_groups(false),
89       _union_comm(MPI_COMM_NULL)
90   {
91     copyInstance(s);
92   }
93
94   DisjointDEC & DisjointDEC::operator=(const DisjointDEC& s)
95   {
96     cleanInstance();
97     copyInstance(s);
98     return *this;
99      
100   }
101
102   void DisjointDEC::copyInstance(const DisjointDEC& other)
103   {
104     DEC::copyFrom(other);
105     if (other._union_comm != MPI_COMM_NULL)
106       {
107         // Tricky: the DEC is responsible for the management of _union_comm. And this comm is referenced by
108         // the MPIProcGroups (source/targets). In the case where _union_comm is not NULL we must take care of rebuilding the
109         // MPIProcGroups with a communicator that will survive the destruction of 'other'.
110         _owns_groups = true;
111         MPI_Comm_dup(other._union_comm, &_union_comm);
112 //        std::cout << "DUP union comm - new is "<< _union_comm << "\n";
113         _target_group = new MPIProcessorGroup(*_comm_interface, other._target_group->getProcIDs(), _union_comm);
114         _source_group = new MPIProcessorGroup(*_comm_interface, other._source_group->getProcIDs(), _union_comm);
115       }
116     else{
117       if(other._target_group)
118         {
119           _target_group=other._target_group->deepCopy();
120           _owns_groups=true;
121         }
122       if(other._source_group)
123         {
124           _source_group=other._source_group->deepCopy();
125           _owns_groups=true;
126         }
127     }
128     if (_source_group && _target_group)
129       _union_group = _source_group->fuse(*_target_group);
130   }
131
132   DisjointDEC::DisjointDEC(const std::set<int>& source_ids,
133                            const std::set<int>& target_ids,
134                            const MPI_Comm& world_comm):
135      _local_field(0),
136      _owns_field(false),
137      _owns_groups(true),
138      _comm_interface(0),
139      _union_comm(MPI_COMM_NULL)
140   {
141     MEDCoupling::CommInterface comm;
142     // Create the list of procs including source and target
143     std::set<int> union_ids; // source and target ids in world_comm
144     union_ids.insert(source_ids.begin(),source_ids.end());
145     union_ids.insert(target_ids.begin(),target_ids.end());
146     if(union_ids.size()!=(source_ids.size()+target_ids.size()))
147       throw INTERP_KERNEL::Exception("DisjointDEC constructor : source_ids and target_ids overlap partially or fully. This type of DEC does not support it! OverlapDEC class could be the solution!");
148     int* union_ranks_world=new int[union_ids.size()]; // ranks of sources and targets in world_comm
149     std::copy(union_ids.begin(), union_ids.end(), union_ranks_world);
150
151     // Create a communicator on these procs
152     MPI_Group union_group,world_group;
153     comm.commGroup(world_comm,&world_group);
154     comm.groupIncl(world_group,union_ids.size(),union_ranks_world,&union_group);
155     comm.commCreate(world_comm,union_group,&_union_comm);
156     delete[] union_ranks_world;
157     if (_union_comm==MPI_COMM_NULL)
158       { // This process is not in union
159         _source_group=0;
160         _target_group=0;
161         _union_group=0;
162         comm.groupFree(&union_group);
163         comm.groupFree(&world_group);
164         return;
165       }
166
167     // Translate source_ids and target_ids from world_comm to union_comm
168     int* source_ranks_world=new int[source_ids.size()]; // ranks of sources in world_comm
169     std::copy(source_ids.begin(), source_ids.end(),source_ranks_world);
170     int* source_ranks_union=new int[source_ids.size()]; // ranks of sources in union_comm
171     int* target_ranks_world=new int[target_ids.size()]; // ranks of targets in world_comm
172     std::copy(target_ids.begin(), target_ids.end(),target_ranks_world);
173     int* target_ranks_union=new int[target_ids.size()]; // ranks of targets in union_comm
174     MPI_Group_translate_ranks(world_group,source_ids.size(),source_ranks_world,union_group,source_ranks_union);
175     MPI_Group_translate_ranks(world_group,target_ids.size(),target_ranks_world,union_group,target_ranks_union);
176     std::set<int> source_ids_union;
177     for (int i=0;i<(int)source_ids.size();i++)
178       source_ids_union.insert(source_ranks_union[i]);
179     std::set<int> target_ids_union;
180     for (int i=0;i<(int)target_ids.size();i++)
181       target_ids_union.insert(target_ranks_union[i]);
182     delete [] source_ranks_world;
183     delete [] source_ranks_union;
184     delete [] target_ranks_world;
185     delete [] target_ranks_union;
186
187     // Create the MPIProcessorGroups
188     _source_group = new MPIProcessorGroup(comm,source_ids_union,_union_comm);
189     _target_group = new MPIProcessorGroup(comm,target_ids_union,_union_comm);
190     _union_group = _source_group->fuse(*_target_group);
191     comm.groupFree(&union_group);
192     comm.groupFree(&world_group);
193   }
194
195   DisjointDEC::~DisjointDEC()
196   {
197     cleanInstance();
198   }
199
200   void DisjointDEC::cleanInstance()
201   {
202     if(_owns_field)
203       {
204         delete _local_field;
205       }
206     _local_field=0;
207     _owns_field=false;
208     if(_owns_groups)
209       {
210         delete _source_group;
211         delete _target_group;
212       }
213     _owns_groups=false;
214     _source_group=0;
215     _target_group=0;
216     delete _union_group;
217     _union_group=0;
218     if (_union_comm != MPI_COMM_NULL)
219       _comm_interface->commFree(&_union_comm);
220   }
221
222   /**
223    * Check that the sources and targets procs form a partition of the world communicator referenced in the groups.
224    * This world communicator is not necessarily MPI_WORLD_COMM, but it has to be covered completely for the DECs to work.
225    */
226   void DisjointDEC::checkPartitionGroup() const
227   {
228     int size = -1;
229     MPIProcessorGroup * tgt = static_cast<MPIProcessorGroup *>(_target_group);
230     MPIProcessorGroup * src = static_cast<MPIProcessorGroup *>(_source_group);
231     MPI_Comm comm_t = tgt->getWorldComm();
232     MPI_Comm comm_s = src->getWorldComm();
233     if (comm_t != comm_s)
234       throw INTERP_KERNEL::Exception("DisjointDEC constructor: Inconsistent world communicator when building DisjointDEC");
235     MPI_Comm_size(comm_t, &size);
236
237     std::set<int> union_ids; // source and target ids in world_comm
238     union_ids.insert(src->getProcIDs().begin(),src->getProcIDs().end());
239     union_ids.insert(tgt->getProcIDs().begin(),tgt->getProcIDs().end());
240     if(union_ids.size()!=size)
241       throw INTERP_KERNEL::Exception("DisjointDEC constructor: source_ids and target_ids do not form a partition of the communicator! Restrain the world communicator passed to MPIProcessorGroup ctor.");
242   }
243
244   void DisjointDEC::setNature(NatureOfField nature)
245   {
246     if(_local_field)
247       _local_field->getField()->setNature(nature);
248   }
249
250   /*! Attaches a local field to a DEC.
251     If the processor is on the receiving end of the DEC, the field
252     will be updated by a recvData() call.
253     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
254   */
255   void DisjointDEC::attachLocalField(const ParaFIELD *field, bool ownPt)
256   {
257     if(!isInUnion())
258       return ;
259     if(_owns_field)
260       delete _local_field;
261     _local_field=field;
262     _owns_field=ownPt;
263     _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
264     compareFieldAndMethod();
265   }
266
267   /*! Attaches a local field to a DEC. The method will test whether the processor
268     is on the source or the target side and will associate the mesh underlying the 
269     field to the local side.
270
271     If the processor is on the receiving end of the DEC, the field
272     will be updated by a recvData() call.
273     Reversely, if the processor is on the sending end, the field will be read, possibly transformed,
274     and sent appropriately to the other side.
275   */
276
277   void DisjointDEC::attachLocalField(MEDCouplingFieldDouble *field)
278   {
279     if(!isInUnion())
280       return ;
281     ProcessorGroup* local_group;
282     if (_source_group->containsMyRank())
283       local_group=_source_group;
284     else if (_target_group->containsMyRank())
285       local_group=_target_group;
286     else
287       throw INTERP_KERNEL::Exception("Invalid procgroup for field attachment to DEC");
288     ParaMESH *paramesh=new ParaMESH(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(field->getMesh())),*local_group,field->getMesh()->getName());
289     ParaFIELD *tmp=new ParaFIELD(field, paramesh, *local_group);
290     tmp->setOwnSupport(true);
291     attachLocalField(tmp,true);
292     //_comm_interface=&(local_group->getCommInterface());
293   }
294
295   /*! 
296     Attaches a local field to a DEC.
297     If the processor is on the receiving end of the DEC, the field
298     will be updated by a recvData() call.
299     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
300     The field type is a generic ICoCo Field, so that the DEC can couple a number of different fields :
301     - a ICoCo::MEDField, that is created from a MEDCoupling structure
302     
303   */
304   void DisjointDEC::attachLocalField(const ICoCo::MEDField *field)
305   {
306     if(!isInUnion())
307       return ;
308     if(!field)
309       throw INTERP_KERNEL::Exception("DisjointDEC::attachLocalField : ICoCo::MEDField pointer is NULL !");
310     attachLocalField(field->getField());
311   }
312   
313   /*!
314     Computes the field norm over its support 
315     on the source side and renormalizes the field on the target side
316     so that the norms match.
317
318     \f[
319     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
320     \f]
321
322     \f[
323     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
324     \f]
325     
326     \f[
327     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
328     \f]
329
330   */
331   void DisjointDEC::renormalizeTargetField(bool isWAbs)
332   {
333     if (_source_group->containsMyRank())
334       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
335         {
336           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
337           double source_norm = total_norm;
338           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
339
340         }
341     if (_target_group->containsMyRank())
342       {
343         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
344           {
345             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
346             double source_norm=total_norm;
347             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
348
349             if (fabs(total_norm)>1e-100)
350               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
351           }
352       }
353   }
354
355   bool DisjointDEC::isInSourceSide() const
356   {
357     if(!_source_group)
358       return false;
359     return _source_group->containsMyRank();
360   }
361
362   bool DisjointDEC::isInTargetSide() const
363   {
364     if(!_target_group)
365       return false;
366     return _target_group->containsMyRank();
367   }
368
369   bool DisjointDEC::isInUnion() const
370   {
371     if(!_union_group)
372       return false;
373     return _union_group->containsMyRank();
374   }
375
376   void DisjointDEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
377   {
378     if (_local_field)
379       {
380         TypeOfField entity = _local_field->getField()->getTypeOfField();
381         if ( getMethod() == "P0" )
382           {
383             if ( entity != ON_CELLS )
384               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
385                                              " For P0 interpolation, field must be on MED_CELL's");
386           }
387         else if ( getMethod() == "P1" )
388           {
389             if ( entity != ON_NODES )
390               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
391                                              " For P1 interpolation, field must be on MED_NODE's");
392           }
393         else if ( getMethod() == "P1d" )
394           {
395             if ( entity != ON_CELLS )
396               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
397                                              " For P1d interpolation, field must be on MED_CELL's");
398             if ( _target_group->containsMyRank() )
399               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
400           }
401         else
402           {
403             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
404           }
405       }
406   }
407
408   /*!
409     If way==true, source procs call sendData() and target procs call recvData().
410     if way==false, it's the other way round.
411   */
412   void DisjointDEC::sendRecvData(bool way)
413   {
414     if(!isInUnion())
415       return;
416     if(isInSourceSide())
417       {
418         if(way)
419           sendData();
420         else
421           recvData();
422       }
423     else if(isInTargetSide())
424       {
425         if(way)
426           recvData();
427         else
428           sendData();
429       }
430   }
431 }