]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM/DisjointDEC.cxx
Salome HOME
add mpi part
[tools/medcoupling.git] / src / ParaMEDMEM / DisjointDEC.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
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 ParaMEDMEM
36 {
37
38   /*!
39    * \anchor DisjointDEC-det
40    * \class DisjointDEC
41    *
42    * Interface class for creation of a link between two
43    * processor groups for exhanging mesh or field data.
44    * The \c DEC is defined by attaching a field on the receiving or on the
45    * sending side.
46    * On top of attaching a \c ParaMEDMEM::ParaFIELD, it is possible to
47    * attach a ICoCo::Field. This class is an abstract class that enables
48    * coupling of codes that respect the ICoCo interface \ref icoco. It has two implementations:
49    * one for codes that express their fields as \ref fields "MEDCoupling fields" (ICoCo::MEDField).
50    *
51    * \section dec_options DEC Options
52    * Options supported by \c DEC objects are
53    *
54    * <TABLE BORDER=1 >
55    * <TR><TD>Option</TD><TD>Description</TD><TD>Default value</TD></TR>
56    * <TR><TD>ForcedRenormalization</TD><TD>After receiving data, the target field is renormalized so that L2-norms of the source and target fields match.</TD><TD> false </TD></TR>
57    *</TABLE>
58
59
60    The following code excerpt shows how to set options for an object that inherits from \c DEC :
61
62    \code
63    InterpKernelDEC dec(source_group,target_group);
64    dec.setOptions("ForcedRenormalization",true);
65    dec.attachLocalField(field);
66    dec.synchronize();
67    if (source_group.containsMyRank())
68      dec.sendData();
69    else
70      dec.recvData();
71    \endcode
72   */
73
74
75   DisjointDEC::DisjointDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
76       _local_field(0),
77       _source_group(&source_group),
78       _target_group(&target_group),
79       _comm_interface(0),
80       _owns_field(false),
81       _owns_groups(false),
82       _union_comm(MPI_COMM_NULL)
83   {
84     _union_group = source_group.fuse(target_group);  
85   }
86   
87   DisjointDEC::DisjointDEC(const DisjointDEC& s):
88       DEC(s),
89       _local_field(0),
90       _union_group(0),
91       _source_group(0),
92       _target_group(0),
93       _comm_interface(0),
94       _owns_field(false),
95       _owns_groups(false),
96       _union_comm(MPI_COMM_NULL)
97   {
98     copyInstance(s);
99   }
100
101   DisjointDEC & DisjointDEC::operator=(const DisjointDEC& s)
102   {
103     cleanInstance();
104     copyInstance(s);
105     return *this;
106      
107   }
108
109   void DisjointDEC::copyInstance(const DisjointDEC& other)
110   {
111     DEC::copyFrom(other);
112     if(other._target_group)
113       {
114         _target_group=other._target_group->deepCpy();
115         _owns_groups=true;
116       }
117     if(other._source_group)
118       {
119         _source_group=other._source_group->deepCpy();
120         _owns_groups=true;
121       }
122     if (_source_group && _target_group)
123       _union_group = _source_group->fuse(*_target_group);
124   }
125
126   DisjointDEC::DisjointDEC(const std::set<int>& source_ids,
127                            const std::set<int>& target_ids,
128                            const MPI_Comm& world_comm):
129      _local_field(0),
130      _owns_field(false),
131      _owns_groups(true),
132      _comm_interface(0),
133      _union_comm(MPI_COMM_NULL)
134   {
135     ParaMEDMEM::CommInterface comm;
136     // Create the list of procs including source and target
137     std::set<int> union_ids; // source and target ids in world_comm
138     union_ids.insert(source_ids.begin(),source_ids.end());
139     union_ids.insert(target_ids.begin(),target_ids.end());
140     if(union_ids.size()!=(source_ids.size()+target_ids.size()))
141       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 !");
142     int* union_ranks_world=new int[union_ids.size()]; // ranks of sources and targets in world_comm
143     std::copy(union_ids.begin(), union_ids.end(), union_ranks_world);
144
145     // Create a communicator on these procs
146     MPI_Group union_group,world_group;
147     comm.commGroup(world_comm,&world_group);
148     comm.groupIncl(world_group,union_ids.size(),union_ranks_world,&union_group);
149     comm.commCreate(world_comm,union_group,&_union_comm);
150     delete[] union_ranks_world;
151     if (_union_comm==MPI_COMM_NULL)
152       { // This process is not in union
153         _source_group=0;
154         _target_group=0;
155         _union_group=0;
156         comm.groupFree(&union_group);
157         comm.groupFree(&world_group);
158         return;
159       }
160
161     // Translate source_ids and target_ids from world_comm to union_comm
162     int* source_ranks_world=new int[source_ids.size()]; // ranks of sources in world_comm
163     std::copy(source_ids.begin(), source_ids.end(),source_ranks_world);
164     int* source_ranks_union=new int[source_ids.size()]; // ranks of sources in union_comm
165     int* target_ranks_world=new int[target_ids.size()]; // ranks of targets in world_comm
166     std::copy(target_ids.begin(), target_ids.end(),target_ranks_world);
167     int* target_ranks_union=new int[target_ids.size()]; // ranks of targets in union_comm
168     MPI_Group_translate_ranks(world_group,source_ids.size(),source_ranks_world,union_group,source_ranks_union);
169     MPI_Group_translate_ranks(world_group,target_ids.size(),target_ranks_world,union_group,target_ranks_union);
170     std::set<int> source_ids_union;
171     for (int i=0;i<(int)source_ids.size();i++)
172       source_ids_union.insert(source_ranks_union[i]);
173     std::set<int> target_ids_union;
174     for (int i=0;i<(int)target_ids.size();i++)
175       target_ids_union.insert(target_ranks_union[i]);
176     delete [] source_ranks_world;
177     delete [] source_ranks_union;
178     delete [] target_ranks_world;
179     delete [] target_ranks_union;
180
181     // Create the MPIProcessorGroups
182     _source_group = new MPIProcessorGroup(comm,source_ids_union,_union_comm);
183     _target_group = new MPIProcessorGroup(comm,target_ids_union,_union_comm);
184     _union_group = _source_group->fuse(*_target_group);
185     comm.groupFree(&union_group);
186     comm.groupFree(&world_group);
187   }
188
189   DisjointDEC::~DisjointDEC()
190   {
191     cleanInstance();
192   }
193
194   void DisjointDEC::cleanInstance()
195   {
196     if(_owns_field)
197       {
198         delete _local_field;
199       }
200     _local_field=0;
201     _owns_field=false;
202     if(_owns_groups)
203       {
204         delete _source_group;
205         delete _target_group;
206       }
207     _owns_groups=false;
208     _source_group=0;
209     _target_group=0;
210     delete _union_group;
211     _union_group=0;
212     if (_union_comm != MPI_COMM_NULL)
213       _comm_interface->commFree(&_union_comm);
214   }
215
216   void DisjointDEC::setNature(NatureOfField nature)
217   {
218     if(_local_field)
219       _local_field->getField()->setNature(nature);
220   }
221
222   /*! Attaches a local field to a DEC.
223     If the processor is on the receiving end of the DEC, the field
224     will be updated by a recvData() call.
225     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
226   */
227   void DisjointDEC::attachLocalField(const ParaFIELD *field, bool ownPt)
228   {
229     if(!isInUnion())
230       return ;
231     if(_owns_field)
232       delete _local_field;
233     _local_field=field;
234     _owns_field=ownPt;
235     _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
236     compareFieldAndMethod();
237   }
238
239   /*! Attaches a local field to a DEC. The method will test whether the processor
240     is on the source or the target side and will associate the mesh underlying the 
241     field to the local side.
242
243     If the processor is on the receiving end of the DEC, the field
244     will be updated by a recvData() call.
245     Reversely, if the processor is on the sending end, the field will be read, possibly transformed,
246     and sent appropriately to the other side.
247   */
248
249   void DisjointDEC::attachLocalField(MEDCouplingFieldDouble *field)
250   {
251     if(!isInUnion())
252       return ;
253     ProcessorGroup* local_group;
254     if (_source_group->containsMyRank())
255       local_group=_source_group;
256     else if (_target_group->containsMyRank())
257       local_group=_target_group;
258     else
259       throw INTERP_KERNEL::Exception("Invalid procgroup for field attachment to DEC");
260     ParaMESH *paramesh=new ParaMESH(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(field->getMesh())),*local_group,field->getMesh()->getName());
261     ParaFIELD *tmp=new ParaFIELD(field, paramesh, *local_group);
262     tmp->setOwnSupport(true);
263     attachLocalField(tmp,true);
264     //_comm_interface=&(local_group->getCommInterface());
265   }
266
267   /*! 
268     Attaches a local field to a DEC.
269     If the processor is on the receiving end of the DEC, the field
270     will be updated by a recvData() call.
271     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
272     The field type is a generic ICoCo Field, so that the DEC can couple a number of different fields :
273     - a ICoCo::MEDField, that is created from a MEDCoupling structure
274     
275   */
276   void DisjointDEC::attachLocalField(const ICoCo::MEDField *field)
277   {
278     if(!isInUnion())
279       return ;
280     if(!field)
281       throw INTERP_KERNEL::Exception("DisjointDEC::attachLocalField : ICoCo::MEDField pointer is NULL !");
282     attachLocalField(field->getField());
283   }
284   
285   /*!
286     Computes the field norm over its support 
287     on the source side and renormalizes the field on the target side
288     so that the norms match.
289
290     \f[
291     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
292     \f]
293
294     \f[
295     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
296     \f]
297     
298     \f[
299     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
300     \f]
301
302   */
303   void DisjointDEC::renormalizeTargetField(bool isWAbs)
304   {
305     if (_source_group->containsMyRank())
306       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
307         {
308           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
309           double source_norm = total_norm;
310           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
311
312         }
313     if (_target_group->containsMyRank())
314       {
315         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
316           {
317             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
318             double source_norm=total_norm;
319             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
320
321             if (fabs(total_norm)>1e-100)
322               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
323           }
324       }
325   }
326
327   bool DisjointDEC::isInSourceSide() const
328   {
329     if(!_source_group)
330       return false;
331     return _source_group->containsMyRank();
332   }
333
334   bool DisjointDEC::isInTargetSide() const
335   {
336     if(!_target_group)
337       return false;
338     return _target_group->containsMyRank();
339   }
340
341   bool DisjointDEC::isInUnion() const
342   {
343     if(!_union_group)
344       return false;
345     return _union_group->containsMyRank();
346   }
347
348   void DisjointDEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
349   {
350     if (_local_field)
351       {
352         TypeOfField entity = _local_field->getField()->getTypeOfField();
353         if ( getMethod() == "P0" )
354           {
355             if ( entity != ON_CELLS )
356               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
357                                              " For P0 interpolation, field must be on MED_CELL's");
358           }
359         else if ( getMethod() == "P1" )
360           {
361             if ( entity != ON_NODES )
362               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
363                                              " For P1 interpolation, field must be on MED_NODE's");
364           }
365         else if ( getMethod() == "P1d" )
366           {
367             if ( entity != ON_CELLS )
368               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
369                                              " For P1d interpolation, field must be on MED_CELL's");
370             if ( _target_group->containsMyRank() )
371               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
372           }
373         else
374           {
375             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
376           }
377       }
378   }
379
380   /*!
381     If way==true, source procs call sendData() and target procs call recvData().
382     if way==false, it's the other way round.
383   */
384   void DisjointDEC::sendRecvData(bool way)
385   {
386     if(!isInUnion())
387       return;
388     if(isInSourceSide())
389       {
390         if(way)
391           sendData();
392         else
393           recvData();
394       }
395     else if(isInTargetSide())
396       {
397         if(way)
398           recvData();
399         else
400           sendData();
401       }
402   }
403 }