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