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