]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM/DisjointDEC.cxx
Salome HOME
Update copyrights 2014.
[tools/medcoupling.git] / src / ParaMEDMEM / DisjointDEC.cxx
1 // Copyright (C) 2007-2014  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 "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(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(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         /* Strange part of code localgroup not used !
279         ProcessorGroup* localgroup;
280         if (_source_group->containsMyRank())
281           localgroup=_source_group;
282         else
283         localgroup=_target_group;*/
284         delete _icoco_field;
285         
286         _icoco_field=new ICoCo::MEDField(*const_cast<ICoCo::TrioField* >(triofield));
287         attachLocalField(_icoco_field);
288         return;
289       }
290     throw INTERP_KERNEL::Exception("incompatible field type");
291   }
292   
293   /*!
294     Computes the field norm over its support 
295     on the source side and renormalizes the field on the target side
296     so that the norms match.
297
298     \f[
299     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
300     \f]
301
302     \f[
303     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
304     \f]
305     
306     \f[
307     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
308     \f]
309
310   */
311   void DisjointDEC::renormalizeTargetField(bool isWAbs)
312   {
313     if (_source_group->containsMyRank())
314       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
315         {
316           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
317           double source_norm = total_norm;
318           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
319
320         }
321     if (_target_group->containsMyRank())
322       {
323         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
324           {
325             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
326             double source_norm=total_norm;
327             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
328
329             if (fabs(total_norm)>1e-100)
330               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
331           }
332       }
333   }
334   /*! @} */
335
336   bool DisjointDEC::isInSourceSide() const
337   {
338     if(!_source_group)
339       return false;
340     return _source_group->containsMyRank();
341   }
342
343   bool DisjointDEC::isInTargetSide() const
344   {
345     if(!_target_group)
346       return false;
347     return _target_group->containsMyRank();
348   }
349
350   bool DisjointDEC::isInUnion() const
351   {
352     if(!_union_group)
353       return false;
354     return _union_group->containsMyRank();
355   }
356
357   void DisjointDEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
358   {
359     if (_local_field)
360       {
361         TypeOfField entity = _local_field->getField()->getTypeOfField();
362         if ( getMethod() == "P0" )
363           {
364             if ( entity != ON_CELLS )
365               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
366                                              " For P0 interpolation, field must be on MED_CELL's");
367           }
368         else if ( getMethod() == "P1" )
369           {
370             if ( entity != ON_NODES )
371               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
372                                              " For P1 interpolation, field must be on MED_NODE's");
373           }
374         else if ( getMethod() == "P1d" )
375           {
376             if ( entity != ON_CELLS )
377               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
378                                              " For P1d interpolation, field must be on MED_CELL's");
379             if ( _target_group->containsMyRank() )
380               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
381           }
382         else
383           {
384             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
385           }
386       }
387   }
388
389   /*!
390     If way==true, source procs call sendData() and target procs call recvData().
391     if way==false, it's the other way round.
392   */
393   void DisjointDEC::sendRecvData(bool way)
394   {
395     if(!isInUnion())
396       return;
397     if(isInSourceSide())
398       {
399         if(way)
400           sendData();
401         else
402           recvData();
403       }
404     else if(isInTargetSide())
405       {
406         if(way)
407           recvData();
408         else
409           sendData();
410       }
411   }
412 }