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