Back to index

salome-med  6.5.0
NonCoincidentDEC.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 <mpi.h>
00021 #include "CommInterface.hxx"
00022 #include "Topology.hxx"
00023 #include "BlockTopology.hxx"
00024 #include "ComponentTopology.hxx"
00025 #include "ParaFIELD.hxx"
00026 #include "MPIProcessorGroup.hxx"
00027 #include "DEC.hxx"
00028 #include "NonCoincidentDEC.hxx"
00029 
00030 extern "C" {
00031 #include <fvm_parall.h>
00032 #include <fvm_nodal.h>
00033 #include <fvm_nodal_append.h>
00034 #include <fvm_locator.h>
00035 }
00036 
00037 namespace ParaMEDMEM
00038 {
00039 
00098   fvm_nodal_t*  medmemMeshToFVMMesh(const MEDMEM::MESH* mesh)
00099   {
00100     // create an FVM structure from the paramesh structure
00101     fvm_nodal_t * fvm_nodal = fvm_nodal_create(mesh->getName().c_str(),mesh->getMeshDimension());
00102       
00103     //loop on cell types
00104     int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL);
00105     const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL);
00106     for (int itype=0; itype<nbtypes; itype++)
00107       {
00108         fvm_element_t fvm_type;
00109         switch (types[itype]) 
00110           {
00111           case MED_EN::MED_TRIA3 :
00112             fvm_type=FVM_FACE_TRIA;
00113             break;
00114           case MED_EN::MED_QUAD4 :
00115             fvm_type=FVM_FACE_QUAD;
00116             break;
00117           case MED_EN::MED_TETRA4 :
00118             fvm_type=FVM_CELL_TETRA;
00119             break;
00120           case MED_EN::MED_HEXA8 :
00121             fvm_type=FVM_CELL_HEXA;
00122             break;
00123           default:
00124             throw MEDEXCEPTION(" MED type  conversion to fvm is not handled yet.");
00125             break;
00126 
00127           }
00128 
00129         fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
00130         fvm_lnum_t* conn = new fvm_lnum_t[nbelems*(types[itype]%100)];
00131         const int* mesh_conn =mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype]);
00132         for (int i=0; i<nbelems*(types[itype]%100); i++)
00133           conn[i]=mesh_conn[i]; 
00134         //swapping trias
00135         if (types[itype]==MED_EN::MED_TRIA3)
00136           {
00137             for (int i=0; i<nbelems;i++)
00138               {
00139                 int tmp=conn[3*i];
00140                 conn[3*i]=mesh_conn[3*i+1];
00141                 conn[3*i+1]=tmp;
00142               }
00143           }
00144         //swapping tetras
00145         if (types[itype]==MED_EN::MED_TETRA4)
00146           {
00147             for (int i=0; i<nbelems;i++)
00148               {
00149                 int tmp=conn[4*i];
00150                 conn[4*i]=mesh_conn[4*i+1];
00151                 conn[4*i+1]=tmp;
00152               }
00153           }
00154         fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0);
00155          
00156         int nbnodes= mesh->getNumberOfNodes();
00157         int spacedim=mesh->getSpaceDimension();
00158         fvm_coord_t* coords = new fvm_coord_t[nbnodes*spacedim];
00159         const double* mesh_coords=mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
00160         for (int i=0; i<nbnodes*spacedim; i++)
00161           coords[i]=mesh_coords[i];                  
00162         fvm_nodal_transfer_vertices(fvm_nodal,coords);
00163       }
00164     return fvm_nodal;
00165   }
00166   
00167   fvm_nodal_t*  medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support)
00168   {
00169 
00170     // create an FVM structure from the paramesh structure
00171     fvm_nodal_t * fvm_nodal = fvm_nodal_create(support->getName().c_str(),1);
00172       
00173     const MEDMEM::MESH* mesh= support->getMesh();
00174       
00175     //loop on cell types
00176     MED_EN::medEntityMesh entity = support->getEntity();
00177       
00178     int nbtypes = support->getNumberOfTypes();
00179     const MED_EN::medGeometryElement* types = support->getTypes();
00180     int ioffset=0;
00181     const int* type_offset = support->getNumberIndex();
00182       
00183     //browsing through all types
00184     for (int itype=0; itype<nbtypes; itype++)
00185       {
00186         fvm_element_t fvm_type;
00187         switch (types[itype]) 
00188           {
00189           case MED_EN::MED_TRIA3 :
00190             fvm_type=FVM_FACE_TRIA;
00191             break;
00192           case MED_EN::MED_QUAD4 :
00193             fvm_type=FVM_FACE_QUAD;
00194             break;
00195           case MED_EN::MED_TETRA4 :
00196             fvm_type=FVM_CELL_TETRA;
00197             break;
00198           case MED_EN::MED_HEXA8 :
00199             fvm_type=FVM_CELL_HEXA;
00200             break;
00201           default:
00202             throw MEDEXCEPTION(" MED type  conversion to fvm is not handled yet.");
00203             break;
00204 
00205           }
00206         fvm_lnum_t nbelems = support->getNumberOfElements(types[itype]);
00207          
00208         //for a partial support, defining the element numbers that are taken into
00209         //account in the support
00210         fvm_lnum_t* elem_numbers=0;
00211         if (!support->isOnAllElements())
00212           {
00213             elem_numbers = const_cast<fvm_lnum_t*> (support->getNumber(types[itype]));
00214            
00215             //creating work arrays to store list of elems for partial suports
00216             if (itype>0)
00217               {
00218                 fvm_lnum_t* temp = new int[nbelems];
00219                 for (int i=0; i< nbelems; i++)
00220                   temp[i] = elem_numbers [i]-ioffset;
00221                 ioffset+=type_offset[itype];
00222                 elem_numbers = temp;
00223               }
00224           }
00225         //retrieving original mesh connectivity
00226         fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype]));
00227        
00228         // adding the elements to the FVM structure 
00229         fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers);
00230          
00231         //cleaning work arrays (for partial supports)
00232         if (!support->isOnAllElements() && itype>0)
00233           delete[] elem_numbers;
00234       
00235       }
00236     return fvm_nodal;
00237   }
00238   
00239   NonCoincidentDEC::NonCoincidentDEC()
00240   {  
00241   }
00242 
00256   NonCoincidentDEC::NonCoincidentDEC(ProcessorGroup& source_group,
00257                                      ProcessorGroup& target_group)
00258     :DEC(source_group, target_group)
00259   {}
00260                                    
00261   NonCoincidentDEC::~NonCoincidentDEC()
00262   {
00263   }
00264 
00276   void NonCoincidentDEC::synchronize()
00277   {
00278   
00279     //initializing FVM parallel environment
00280     const MPI_Comm* comm=dynamic_cast<const MPIProcessorGroup*> (_union_group)->getComm();
00281     fvm_parall_set_mpi_comm(*const_cast<MPI_Comm*> (comm));
00282   
00283   
00284     //setting up the communication DEC on both sides
00285   
00286     if (_source_group->containsMyRank())
00287       {
00288         MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
00289         fvm_nodal_t* source_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh);
00290       
00291         int target_size = _target_group->size()  ;
00292         int start_rank=  _source_group->size();
00293         const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
00294       
00295         _locator =  fvm_locator_create(1e-6,
00296                                        *comm,
00297                                        target_size,
00298                                        start_rank);
00299       
00300         fvm_locator_set_nodal(_locator,
00301                               source_nodal,
00302                               mesh->getSpaceDimension(),
00303                               0,
00304                               NULL,
00305                               0);
00306 
00307       
00308         _nb_distant_points = fvm_locator_get_n_dist_points(_locator);
00309         _distant_coords = fvm_locator_get_dist_coords(_locator);
00310         _distant_locations = fvm_locator_get_dist_locations(_locator);
00311            
00312       }
00313     if (_target_group->containsMyRank())
00314       {
00315         MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
00316       
00317         fvm_nodal_t* target_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh);
00318         int source_size = _source_group->size();
00319         int start_rank=  0 ;
00320         const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
00321       
00322         _locator = fvm_locator_create(1e-6,
00323                                       *comm,
00324                                       source_size,
00325                                       start_rank);
00326         int nbcells = mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
00327         const MEDMEM::SUPPORT* support=_local_field->getField()->getSupport();
00328         MEDMEM::FIELD<double>* barycenter_coords = mesh->getBarycenter(support);
00329         const double* coords = barycenter_coords->getValue();
00330         fvm_locator_set_nodal(_locator,
00331                               target_nodal,
00332                               mesh->getSpaceDimension(),
00333                               nbcells,
00334                               NULL,
00335                               coords);  
00336         delete barycenter_coords;
00337       }
00338   }
00339 
00340 
00346   void NonCoincidentDEC::recvData()
00347   {
00348     int nbelems = _local_field->getField()->getSupport()->getMesh()->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
00349     int nbcomp =  _local_field->getField()->getNumberOfComponents();
00350     double* values = new double [nbelems*nbcomp];
00351     fvm_locator_exchange_point_var(_locator,
00352                                    0,
00353                                    values,
00354                                    0,
00355                                    sizeof(double),
00356                                    nbcomp,
00357                                    0);
00358     _local_field->getField()->setValue(values);
00359     if (_forced_renormalization_flag)
00360       renormalizeTargetField();
00361     delete[]values;
00362   }
00363 
00369   void NonCoincidentDEC::sendData()
00370   {
00371     const double* values=_local_field->getField()->getValue();
00372     int nbcomp = _local_field->getField()->getNumberOfComponents();
00373     double* distant_values = new double [_nb_distant_points*nbcomp];
00374 
00375     //cheap interpolation :  the value of the cell is transfered to the point
00376     for (int i=0; i<_nb_distant_points; i++)
00377       for (int j=0; j <nbcomp; j++)
00378         distant_values[i*nbcomp+j]=values[(_distant_locations[i]-1)*nbcomp+j];
00379   
00380     fvm_locator_exchange_point_var(_locator,
00381                                    distant_values,
00382                                    0,
00383                                    0,
00384                                    sizeof(double),
00385                                    nbcomp,
00386                                    0);
00387 
00388     delete [] distant_values;
00389     if (_forced_renormalization_flag)
00390       renormalizeTargetField();
00391 
00392   }
00396 }