Back to index

salome-med  6.5.0
DisjointDEC.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "DisjointDEC.hxx"
00021 #include "CommInterface.hxx"
00022 #include "Topology.hxx"
00023 #include "BlockTopology.hxx"
00024 #include "ComponentTopology.hxx"
00025 #include "ParaFIELD.hxx"
00026 #include "ParaMESH.hxx"
00027 #include "ICoCoField.hxx"
00028 #include "ICoCoMEDField.hxx"
00029 #include "ICoCoTrioField.hxx"
00030 #include "MPIProcessorGroup.hxx"
00031 
00032 #include <cmath>
00033 #include <iostream>
00034 
00072 namespace ParaMEDMEM
00073 {
00074 
00075 
00079   DisjointDEC::DisjointDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):_local_field(0), 
00080                                                                                        _source_group(&source_group),
00081                                                                                        _target_group(&target_group),
00082                                                                                        _owns_field(false),
00083                                                                                        _owns_groups(false),
00084                                                                                        _icoco_field(0)
00085   {
00086     _union_group = source_group.fuse(target_group);  
00087   }
00088   
00089   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)
00090   {
00091     copyInstance(s);
00092   }
00093 
00094   DisjointDEC & DisjointDEC::operator=(const DisjointDEC& s)
00095   {
00096     cleanInstance();
00097     copyInstance(s);
00098     return *this;
00099      
00100   }
00101 
00102   void DisjointDEC::copyInstance(const DisjointDEC& other)
00103   {
00104     DEC::copyFrom(other);
00105     if(other._target_group)
00106       {
00107         _target_group=other._target_group->deepCpy();
00108         _owns_groups=true;
00109       }
00110     if(other._source_group)
00111       {
00112         _source_group=other._source_group->deepCpy();
00113         _owns_groups=true;
00114       }
00115     if (_source_group && _target_group)
00116       _union_group = _source_group->fuse(*_target_group);
00117   }
00118 
00119   DisjointDEC::DisjointDEC(const std::set<int>& source_ids, const std::set<int>& target_ids, const MPI_Comm& world_comm):_local_field(0), 
00120                                                                                                                          _owns_field(false),
00121                                                                                                                          _owns_groups(true),
00122                                                                                                                          _icoco_field(0)
00123   {
00124     ParaMEDMEM::CommInterface comm;
00125     // Create the list of procs including source and target
00126     std::set<int> union_ids; // source and target ids in world_comm
00127     union_ids.insert(source_ids.begin(),source_ids.end());
00128     union_ids.insert(target_ids.begin(),target_ids.end());
00129     if(union_ids.size()!=(source_ids.size()+target_ids.size()))
00130       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 !");
00131     int* union_ranks_world=new int[union_ids.size()]; // ranks of sources and targets in world_comm
00132     std::copy(union_ids.begin(), union_ids.end(), union_ranks_world);
00133 
00134     // Create a communicator on these procs
00135     MPI_Group union_group,world_group;
00136     comm.commGroup(world_comm,&world_group);
00137     comm.groupIncl(world_group,union_ids.size(),union_ranks_world,&union_group);
00138     MPI_Comm union_comm;
00139     comm.commCreate(world_comm,union_group,&union_comm);
00140     delete[] union_ranks_world;
00141 
00142     if (union_comm==MPI_COMM_NULL)
00143       { // This process is not in union
00144         _source_group=0;
00145         _target_group=0;
00146         _union_group=0;
00147         return;
00148       }
00149 
00150     // Translate source_ids and target_ids from world_comm to union_comm
00151     int* source_ranks_world=new int[source_ids.size()]; // ranks of sources in world_comm
00152     std::copy(source_ids.begin(), source_ids.end(),source_ranks_world);
00153     int* source_ranks_union=new int[source_ids.size()]; // ranks of sources in union_comm
00154     int* target_ranks_world=new int[target_ids.size()]; // ranks of targets in world_comm
00155     std::copy(target_ids.begin(), target_ids.end(),target_ranks_world);
00156     int* target_ranks_union=new int[target_ids.size()]; // ranks of targets in union_comm
00157     MPI_Group_translate_ranks(world_group,source_ids.size(),source_ranks_world,union_group,source_ranks_union);
00158     MPI_Group_translate_ranks(world_group,target_ids.size(),target_ranks_world,union_group,target_ranks_union);
00159     std::set<int> source_ids_union;
00160     for (int i=0;i<(int)source_ids.size();i++)
00161       source_ids_union.insert(source_ranks_union[i]);
00162     std::set<int> target_ids_union;
00163     for (int i=0;i<(int)target_ids.size();i++)
00164       target_ids_union.insert(target_ranks_union[i]);
00165     delete [] source_ranks_world;
00166     delete [] source_ranks_union;
00167     delete [] target_ranks_world;
00168     delete [] target_ranks_union;
00169 
00170     // Create the MPIProcessorGroups
00171     _source_group = new MPIProcessorGroup(comm,source_ids_union,union_comm);
00172     _target_group = new MPIProcessorGroup(comm,target_ids_union,union_comm);
00173     _union_group = _source_group->fuse(*_target_group);
00174 
00175   }
00176 
00177   DisjointDEC::~DisjointDEC()
00178   {
00179     cleanInstance();
00180   }
00181 
00182   void DisjointDEC::cleanInstance()
00183   {
00184     if(_owns_field)
00185       {
00186         delete _local_field;
00187       }
00188     _local_field=0;
00189     _owns_field=false;
00190     if(_owns_groups)
00191       {
00192         delete _source_group;
00193         delete _target_group;
00194       }
00195     _owns_groups=false;
00196     _source_group=0;
00197     _target_group=0;
00198     delete _icoco_field;
00199     _icoco_field=0;
00200     delete _union_group;
00201     _union_group=0;
00202   }
00203 
00204   void DisjointDEC::setNature(NatureOfField nature)
00205   {
00206     if(_local_field)
00207       _local_field->getField()->setNature(nature);
00208   }
00209 
00215   void DisjointDEC::attachLocalField(const ParaFIELD* field, bool ownPt) 
00216   {
00217     if(!isInUnion())
00218       return ;
00219     if(_owns_field)
00220       delete _local_field;
00221     _local_field=field;
00222     _owns_field=ownPt;
00223     _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
00224     compareFieldAndMethod();
00225   }
00226 
00237   void DisjointDEC::attachLocalField(MEDCouplingFieldDouble* field) 
00238   {
00239     if(!isInUnion())
00240       return ;
00241     ProcessorGroup* local_group;
00242     if (_source_group->containsMyRank())
00243       local_group=_source_group;
00244     else if (_target_group->containsMyRank())
00245       local_group=_target_group;
00246     else
00247       throw INTERP_KERNEL::Exception("Invalid procgroup for field attachment to DEC");
00248     ParaMESH *paramesh=new ParaMESH((MEDCouplingPointSet *)field->getMesh(),*local_group,field->getMesh()->getName());
00249     ParaFIELD *tmp=new ParaFIELD(field, paramesh, *local_group);
00250     tmp->setOwnSupport(true);
00251     attachLocalField(tmp,true);
00252     //_comm_interface=&(local_group->getCommInterface());
00253   }
00254 
00265   void DisjointDEC::attachLocalField(const ICoCo::Field* field)
00266   {
00267     if(!isInUnion())
00268       return ;
00269     const ICoCo::MEDField* medfield=dynamic_cast<const ICoCo::MEDField*> (field);
00270     if(medfield !=0)
00271       {
00272         attachLocalField(medfield->getField());
00273         return;
00274       }
00275     const ICoCo::TrioField* triofield=dynamic_cast<const ICoCo::TrioField*> (field);
00276     if (triofield !=0)
00277       {
00278         ProcessorGroup* localgroup;
00279         if (_source_group->containsMyRank())
00280           localgroup=_source_group;
00281         else
00282           localgroup=_target_group;
00283         delete _icoco_field;
00284         
00285         _icoco_field=new ICoCo::MEDField(*const_cast<ICoCo::TrioField* >(triofield));
00286         attachLocalField(_icoco_field);
00287         return;
00288       }
00289     throw INTERP_KERNEL::Exception("incompatible field type");
00290   }
00291   
00310   void DisjointDEC::renormalizeTargetField(bool isWAbs)
00311   {
00312     if (_source_group->containsMyRank())
00313       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
00314         {
00315           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
00316           double source_norm = total_norm;
00317           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
00318 
00319         }
00320     if (_target_group->containsMyRank())
00321       {
00322         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
00323           {
00324             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
00325             double source_norm=total_norm;
00326             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
00327 
00328             if (fabs(total_norm)>1e-100)
00329               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
00330           }
00331       }
00332   }
00335   bool DisjointDEC::isInSourceSide() const
00336   {
00337     if(!_source_group)
00338       return false;
00339     return _source_group->containsMyRank();
00340   }
00341 
00342   bool DisjointDEC::isInTargetSide() const
00343   {
00344     if(!_target_group)
00345       return false;
00346     return _target_group->containsMyRank();
00347   }
00348 
00349   bool DisjointDEC::isInUnion() const
00350   {
00351     if(!_union_group)
00352       return false;
00353     return _union_group->containsMyRank();
00354   }
00355 
00356   void DisjointDEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
00357   {
00358     if (_local_field)
00359       {
00360         TypeOfField entity = _local_field->getField()->getTypeOfField();
00361         if ( getMethod() == "P0" )
00362           {
00363             if ( entity != ON_CELLS )
00364               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
00365                                              " For P0 interpolation, field must be on MED_CELL's");
00366           }
00367         else if ( getMethod() == "P1" )
00368           {
00369             if ( entity != ON_NODES )
00370               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
00371                                              " For P1 interpolation, field must be on MED_NODE's");
00372           }
00373         else if ( getMethod() == "P1d" )
00374           {
00375             if ( entity != ON_CELLS )
00376               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
00377                                              " For P1d interpolation, field must be on MED_CELL's");
00378             if ( _target_group->containsMyRank() )
00379               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
00380           }
00381         else
00382           {
00383             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
00384           }
00385       }
00386   }
00387 
00392   void DisjointDEC::sendRecvData(bool way)
00393   {
00394     if(!isInUnion())
00395       return;
00396     if(isInSourceSide())
00397       {
00398         if(way)
00399           sendData();
00400         else
00401           recvData();
00402       }
00403     else if(isInTargetSide())
00404       {
00405         if(way)
00406           recvData();
00407         else
00408           sendData();
00409       }
00410   }
00411 }