Back to index

salome-med  6.5.0
MEDSPLITTER_MESHCollection.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 "MEDMEM_ConnectZone.hxx"
00021 #include "MEDMEM_DriversDef.hxx"
00022 #include "MEDMEM_Mesh.hxx"
00023 #include "MEDMEM_Meshing.hxx"
00024 #include "MEDMEM_GaussLocalization.hxx"
00025 #include "MEDMEM_Field.hxx"
00026 #include "MEDMEM_CellModel.hxx"
00027 #include "MEDMEM_Group.hxx"
00028 #include "MEDMEM_MeshFuse.hxx"
00029 
00030 #include "MEDMEM_Exception.hxx"
00031 
00032 #include "MEDSPLITTER_utils.hxx" 
00033 
00034 #include "MEDSPLITTER_Graph.hxx"
00035 
00036 #include "MEDSPLITTER_Topology.hxx"
00037 #include "MEDSPLITTER_ParallelTopology.hxx"
00038 #include "MEDSPLITTER_SequentialTopology.hxx"
00039 #include "MEDSPLITTER_ParaDomainSelector.hxx"
00040 #include "MEDSPLITTER_MeshSendReceive.hxx"
00041 #include "MEDSPLITTER_JointExchangeData.hxx"
00042 
00043 #include "MEDSPLITTER_MESHCollection.hxx"
00044 #include "MEDSPLITTER_MESHCollectionDriver.hxx"
00045 #include "MEDSPLITTER_MESHCollectionMedXMLDriver.hxx"
00046 #include "MEDSPLITTER_MESHCollectionMedAsciiDriver.hxx"
00047 
00048 #include "MEDSPLITTER_UserGraph.hxx"
00049 
00050 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
00051 #include "MEDSPLITTER_METISGraph.hxx"
00052 #endif
00053 #ifdef MED_ENABLE_SCOTCH
00054 #include "MEDSPLITTER_SCOTCHGraph.hxx"
00055 #endif
00056 
00057 #include "InterpKernelHashMap.hxx"
00058 
00059 #include <vector>
00060 #include <string>
00061 #include <set>
00062 
00063 #include <iostream>
00064 #include <fstream>
00065 
00066 using namespace MEDSPLITTER;
00067 
00068 using namespace INTERP_KERNEL;
00069 
00070 //template inclusion
00071 #include "MEDSPLITTER_MESHCollection.H"
00072 
00073 MESHCollection::MESHCollection()
00074   : _topology(0),
00075     _owns_topology(false),
00076     _driver(0),
00077     _domain_selector( 0 ),
00078     _i_non_empty_mesh(-1),
00079     _driver_type(MEDSPLITTER::MedXML),
00080     _subdomain_boundary_creates(false),
00081     _family_splitting(false),
00082     _create_empty_groups(false)
00083 {
00084 }
00085 
00086 namespace
00087 {
00088   //================================================================================
00092   //================================================================================
00093 
00094   MEDMEM::MESH* newMesh(const std::string& name, int dim, int space, MEDMEM::MESH* meshToDelete=0)
00095   {
00096     delete meshToDelete;
00097     MEDMEM::MESHING* mesh = new MEDMEM::MeshFuse;
00098     mesh->setName( name );
00099     //mesh->setMeshDimension ( dim );
00100     //mesh->setSpaceDimension( space );
00101     return mesh;
00102   }
00103 }
00104 
00116 MESHCollection::MESHCollection(const MESHCollection& initial_collection, Topology* topology, bool family_splitting, bool create_empty_groups)
00117   : _topology(topology),
00118     _owns_topology(false),
00119     _cell_graph(topology->getGraph()),
00120     _driver(0),
00121     _domain_selector( initial_collection._domain_selector ),
00122     _i_non_empty_mesh(-1),
00123     _name(initial_collection._name),
00124     _driver_type(MEDSPLITTER::MedXML),
00125     _subdomain_boundary_creates(false),
00126     _family_splitting(family_splitting),
00127     _create_empty_groups(create_empty_groups)
00128 {
00129   string mesh_name = initial_collection.getName();
00130   _mesh.resize(_topology->nbDomain());
00131 
00132   int space_dim = initial_collection.getSpaceDimension();
00133   int mesh_dim  = initial_collection.getMeshDimension();
00134   if ( mesh_dim < 1 )
00135     space_dim = mesh_dim = initial_collection._driver->readMeshDimension();
00136 
00137   for (int idomain=0; idomain < _topology->nbDomain(); idomain++)
00138   {
00139     //creating the new mesh
00140     _mesh[idomain]= newMesh( MEDMEM::STRING(mesh_name)<<"_"<<idomain+1, mesh_dim, space_dim );
00141 
00142     createNodalConnectivity(initial_collection,idomain, MED_EN::MED_CELL);
00143 
00144     if ( _mesh[idomain]->getNumberOfNodes() > 0 )
00145       _i_non_empty_mesh = idomain;
00146   }
00147 
00148   _topology->createFaceMapping(initial_collection, *this ); 
00149   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
00150   {
00151     switch (mesh_dim)
00152     {
00153     case 3:
00154       createNodalConnectivity(initial_collection,idomain, MED_EN::MED_FACE);
00155       break;
00156     case 2:
00157       createNodalConnectivity(initial_collection,idomain, MED_EN::MED_EDGE);
00158       break;
00159     default :
00160       if ( !isParallelMode() || _domain_selector->isMyDomain( idomain ))
00161         cerr<<"MEDSPLITTER : Mesh dimension must be 2 or 3"<<endl;
00162     }
00163   }
00164 
00165   castFamilies(initial_collection);
00166 
00167   // Exchange domain parts
00168 
00169   if ( isParallelMode() )
00170   {
00171     _domain_selector->setNbDomains( _topology->nbDomain() );
00172 
00173     vector< MeshSendReceive > mesh_sender( _topology->nbDomain() );
00174     list<int> domainsToClear; // sent domains
00175     bool isSent;
00176     // first, send domains
00177     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
00178     {
00179       // get node numbers global over all procs
00180       vector<int> glob_nodes_all_proc( _topology->getNodeNumber( idomain )); // to fill in
00181       vector<int> glob_cells_all_proc( _topology->getCellNumber( idomain ));
00182       vector<int> glob_faces_all_proc( _topology->getFaceNumber( idomain ));
00183       if ( !glob_cells_all_proc.empty() )
00184       {
00185         // get ids global on this proc
00186         _topology->getNodeList( idomain, & glob_nodes_all_proc[0] );
00187         _topology->getCellList( idomain, & glob_cells_all_proc[0] );
00188         _topology->getFaceList( idomain, & glob_faces_all_proc[0] );
00189         // convert cell ids to ids global over all procs
00190         int cell_shift = _domain_selector->getProcShift();
00191         for ( int i = 0; i < glob_cells_all_proc.size(); ++i )
00192           glob_cells_all_proc[i] += cell_shift;
00193       }
00194       if ( _domain_selector->isMyDomain( idomain ))
00195       {
00196         // prepare to receiving other parts of the domain
00197         ((MEDMEM::MeshFuse*) _mesh[idomain])->setNodeNumbers( glob_nodes_all_proc );
00198         _topology->getFusedCellNumbers( idomain ) = glob_cells_all_proc;
00199         _topology->getFusedFaceNumbers( idomain ) = glob_faces_all_proc;
00200       }
00201       else
00202       {
00203         // sending
00204         int target_proc = _domain_selector->getProccessorID( idomain );
00205         mesh_sender[ idomain ].send( target_proc, idomain, _mesh[idomain],
00206                                      glob_cells_all_proc,
00207                                      glob_faces_all_proc,
00208                                      glob_nodes_all_proc );
00209         if ( !glob_nodes_all_proc.empty() )
00210           domainsToClear.push_back( idomain );
00211       }
00212       // clear just sent domain meshes
00213       for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); )
00214       {
00215         if (( isSent = mesh_sender[ *dom ].isSent() ))
00216           _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim, _mesh[*dom]);
00217         dom = isSent ? domainsToClear.erase( dom ) : ++dom;
00218       }
00219     }
00220 
00221     // then, receive domains
00222     MeshSendReceive mesh_receiver;
00223     int this_proc = _domain_selector->rank();
00224     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
00225     {
00226       if ( _domain_selector->isMyDomain( idomain ))
00227       {
00228         for (int iproc = 0; iproc < _domain_selector->nbProcs(); ++iproc)
00229         {
00230           if ( iproc == this_proc ) continue;
00231           vector<int>  nodes_other_proc, cells_other_proc, faces_other_proc;
00232           mesh_receiver.recv( iproc, idomain, cells_other_proc, faces_other_proc,nodes_other_proc);
00233           if ( MEDMEM::MESH* received_mesh = mesh_receiver.getMesh() )
00234           {
00235             // unite meshes and global node numbers stored in MeshFuse
00236             MEDMEM::MeshFuse* fuse = (MEDMEM::MeshFuse*) _mesh[idomain];
00237             fuse->concatenate( received_mesh, nodes_other_proc );
00238             delete received_mesh;
00239 
00240             // unite global element numbers
00241             fuse->append( MED_EN::MED_CELL,
00242                           _topology->getFusedCellNumbers( idomain ), cells_other_proc );
00243 
00244             fuse->append( mesh_dim==3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE,
00245                           _topology->getFusedFaceNumbers( idomain ), faces_other_proc );
00246 
00247             if ( _mesh[idomain]->getNumberOfNodes() > 0 )
00248               _i_non_empty_mesh = idomain;
00249           }
00250         }
00251       }
00252       // clear just sent domain meshes
00253       for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); )
00254       {
00255         if (( isSent = mesh_sender[ *dom ].isSent() ))
00256           _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim,_mesh[*dom]);
00257         dom = isSent ? domainsToClear.erase( dom ) : ++dom;
00258       }
00259     }
00260     // clear sent domain meshes
00261     mesh_sender.clear();
00262     for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); ++dom)
00263       _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim,_mesh[*dom]);
00264 
00265     _topology->recreateMappingAfterFusion( getMesh() );
00266   }
00267   if ( _i_non_empty_mesh < 0 ) // non of domains resides on this proc,
00268     _i_non_empty_mesh = 0; // in this case we need only dimension that is set to all meshes
00269 
00270 }
00271 
00276 MESHCollection::MESHCollection(const string& filename)
00277   : _topology(0),
00278     _owns_topology(true),
00279     _driver(0),
00280     _domain_selector( 0 ),
00281     _i_non_empty_mesh(-1),
00282     _driver_type(MEDSPLITTER::Undefined),
00283     _subdomain_boundary_creates(false),
00284     _family_splitting(false),
00285     _create_empty_groups(false)
00286 {
00287   char filenamechar[256];
00288   strcpy(filenamechar,filename.c_str());
00289   try
00290   {
00291     _driver=new MESHCollectionMedXMLDriver(this);
00292     _driver->read (filenamechar);
00293     _driver_type = MedXML;
00294 
00295   }
00296   catch(MEDEXCEPTION&){
00297     delete _driver;
00298     try
00299     {
00300       _driver=new MESHCollectionMedAsciiDriver(this);
00301       _driver->read (filenamechar);
00302       _driver_type=MedAscii;
00303     }
00304     catch(MEDEXCEPTION&)
00305     {
00306       delete _driver;
00307       throw MEDEXCEPTION("file does not comply with any recognized format");
00308     }
00309   }
00310   for ( int idomain = 0; idomain < _mesh.size(); ++idomain )
00311     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
00312       _i_non_empty_mesh = idomain;
00313 }
00314 
00320 MESHCollection::MESHCollection(const string& filename, ParaDomainSelector& domainSelector)
00321   : _topology(0),
00322     _owns_topology(true),
00323     _driver(0),
00324     _domain_selector( domainSelector.nbProcs() > 1 ? & domainSelector : 0 ),
00325     _i_non_empty_mesh(-1),
00326     _driver_type(MEDSPLITTER::Undefined),
00327     _subdomain_boundary_creates(false),
00328     _family_splitting(false),
00329     _create_empty_groups(false)
00330 {
00331   try
00332   {
00333     _driver=new MESHCollectionMedXMLDriver(this);
00334     _driver->read ( (char*)filename.c_str(), _domain_selector );
00335     _driver_type = MedXML;
00336   }
00337   catch(MEDEXCEPTION&)
00338   {
00339     delete _driver;
00340     try
00341     {
00342       _driver=new MESHCollectionMedAsciiDriver(this);
00343       _driver->read ( (char*)filename.c_str(), _domain_selector );
00344       _driver_type=MedAscii;
00345     }
00346     catch(MEDEXCEPTION&)
00347     {
00348       delete _driver;
00349       throw MEDEXCEPTION("file does not comply with any recognized format");
00350     }
00351   }
00352   if ( isParallelMode() )
00353     // to know nb of cells on each proc to compute global cell ids from locally global
00354     _domain_selector->gatherNbOf( MED_EN::MED_CELL, getMesh() );
00355 
00356   // find non-empty domain mesh
00357   for ( int idomain = 0; idomain < _mesh.size(); ++idomain )
00358     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
00359       _i_non_empty_mesh = idomain;
00360 }
00361 
00367 MESHCollection::MESHCollection(const string& filename, const string& meshname)
00368   : _topology(0),
00369     _owns_topology(true),
00370     _driver(0),
00371     _domain_selector( 0 ),
00372     _i_non_empty_mesh(-1),
00373     _name(meshname),
00374     _driver_type(MEDSPLITTER::MedXML),
00375     _subdomain_boundary_creates(false),
00376     _family_splitting(false),
00377     _create_empty_groups(false)
00378 {
00379   char filenamechar[256];
00380   char meshnamechar[256];
00381   strcpy(filenamechar,filename.c_str());
00382   strcpy(meshnamechar,meshname.c_str());
00383   try // avoid memory leak in case of inexistent filename
00384   {
00385     retrieveDriver()->readSeq (filenamechar,meshnamechar);
00386   }
00387   catch ( MED_EXCEPTION& e )
00388   {
00389     if ( _driver ) delete _driver; _driver=0;
00390     throw e;
00391   }
00392   if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
00393     _i_non_empty_mesh = 0;
00394 }
00395 
00396 MESHCollection::~MESHCollection()
00397 {
00398   for (unsigned i=0; i<_mesh.size();i++)
00399     if (_mesh[i]!=0) {/*delete*/ _mesh[i]->removeReference(); }
00400   for (unsigned i=0; i<_connect_zones.size();i++)
00401     if (_connect_zones[i]!=0) {delete _connect_zones[i];}
00402   if (_driver !=0) {delete _driver; _driver=0;}
00403   if (_topology!=0 && _owns_topology) {delete _topology; _topology=0;}
00404 }
00405 
00418 void MESHCollection::getNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
00419                                          MED_EN::medGeometryElement type, int* type_connectivity) const
00420 {
00421   int *local=new int[nb_cells];
00422   int *ip=new int[nb_cells];
00423   switch (entity)
00424   {
00425   case MED_EN::MED_CELL:
00426     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
00427     break;
00428   case MED_EN::MED_FACE:
00429   case MED_EN::MED_EDGE:
00430     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
00431     break;
00432   }
00433 
00434 
00435   //  int nbnode_per_type=(int)type%100;
00436   //  vector<int> number_of_types_array(_topology->nbDomain(),0);
00437   //  for (int i=0; i<_topology->nbDomain(); i++)
00438   //    number_of_types_array[i]=_mesh[i]->getNumberOfTypes(entity);
00439 
00440   //defining a connectivity table for different domains 
00441   vector  <const int*> conn_ip(_topology->nbDomain());
00442   vector  <const int*> conn_index_ip(_topology->nbDomain());
00443 
00444 
00445   vector< map <MED_EN::medGeometryElement, int> > offset;
00446   //  offset.resize(_topology->nbDomain());
00447 
00448   for (int i=0; i<_topology->nbDomain();i++)
00449   {
00450     if ( !_mesh[i] ) continue;
00451     int nb_elem = _mesh[i]->getNumberOfElements(entity,type);
00452     if (nb_elem>0)
00453     {
00454       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,entity,MED_EN::MED_ALL_ELEMENTS);
00455       conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,entity);
00456       //       global_index= _mesh[i]->getGlobalNumberingIndex(entity);
00457     }                                             
00458     else
00459     {
00460       conn_ip[i]=0;
00461       conn_index_ip[i]=0;
00462     }     
00463     //      int number_of_types = number_of_types_array[i];
00464     //      const MEDMEM::CELLMODEL* types =  _mesh[ip[icell]]->getCellsTypes(entity);
00465     //      for (int itype=0; itype<number_of_types; itype++) 
00466     //        offset[i][types[itype].getType()]=global_index[itype]-1;
00467   }
00468 
00469   int* type_connectivity_ptr=type_connectivity;
00470   for (int icell=0; icell<nb_cells; icell++)
00471   {
00472     //      int type_offset = offset[ip[icell]][type];
00473     const int* conn=conn_ip[ip[icell]]; 
00474     const int* conn_index=conn_index_ip[ip[icell]];
00475     for (int inode=conn_index[local[icell]-1]; inode<conn_index[local[icell]]; inode++)
00476     {
00477       *type_connectivity_ptr=
00478         _topology->convertNodeToGlobal(ip[icell],conn[inode-1]);
00479       type_connectivity_ptr++;                       
00480     }
00481   }
00482 
00483   delete[]local;
00484   delete[]ip;
00485 }
00486 
00496 void MESHCollection::getPolygonNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
00497                                                 vector<int>& type_connectivity, vector<int>& connectivity_index) const
00498 {
00499 
00500   int *local=new int[nb_cells];
00501   int *ip=new int[nb_cells];
00502   switch (entity)
00503   {
00504   case MED_EN::MED_CELL:
00505     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
00506     break;
00507   case MED_EN::MED_FACE:
00508   case MED_EN::MED_EDGE:
00509     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
00510     break;
00511   }
00512 
00513 
00514   //defining a connectivity table for different domains 
00515   vector  <const int*> conn_ip(_topology->nbDomain());
00516   vector  <const int*> conn_index_ip(_topology->nbDomain());
00517   vector <const int* > conn_face_index(_topology->nbDomain());
00518   vector<int> nb_plain_elems(_topology->nbDomain());
00519 
00520   vector< map <MED_EN::medGeometryElement, int> > offset;
00521 
00522   for (int i=0; i<_topology->nbDomain();i++)
00523   {
00524     int nb_elem = _mesh[i]->getNumberOfElements(entity,MED_EN::MED_POLYGON);
00525     if (nb_elem>0)
00526     {
00527       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,entity,MED_EN::MED_ALL_ELEMENTS);
00528       conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,entity);
00529     }
00530     else
00531     {
00532       conn_ip[i]=0;
00533       conn_index_ip[i]=0;
00534     }     
00535   }
00536 
00537   connectivity_index.resize(nb_cells+1);
00538   connectivity_index[0]=1;
00539   for (int icell=0; icell<nb_cells; icell++)
00540   {
00541     //int nb_plain= nb_plain_elems[ip[icell]];
00542     const int* conn=conn_ip[ip[icell]]; 
00543     const int* conn_index=conn_index_ip[ip[icell]];
00544     for (int inode=conn_index[local[icell]-1/*-nb_plain*/]; inode<conn_index[local[icell]/*-nb_plain*/]; inode++)
00545     {
00546       type_connectivity.push_back(
00547                                   _topology->convertNodeToGlobal(ip[icell],conn[inode-1]));
00548     }
00549     connectivity_index[icell+1]=connectivity_index[icell]
00550       -conn_index[local[icell]-1/*-nb_plain*/]+conn_index[local[icell]/*-nb_plain*/];
00551   }
00552 
00553   delete[]local;
00554   delete[]ip;
00555 }
00556 
00557 
00567 void MESHCollection::getPolyhedraNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
00568                                                   vector<int>& type_connectivity, vector<int>& connectivity_index/*, vector<int>& face_connectivity_index*/) const
00569 {
00570 
00571   int *local=new int[nb_cells];
00572   int *ip=new int[nb_cells];
00573   switch (entity)
00574   {
00575   case MED_EN::MED_CELL:
00576     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
00577     break;
00578   case MED_EN::MED_FACE:
00579   case MED_EN::MED_EDGE:
00580     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
00581     break;
00582   }
00583 
00584 
00585   //defining a connectivity table for different domains 
00586   vector  <const int*> conn_ip(_topology->nbDomain());
00587   vector  <const int*> conn_index_ip(_topology->nbDomain());
00588   vector<int> nb_plain_elems(_topology->nbDomain());
00589 
00590   vector< map <MED_EN::medGeometryElement, int> > offset;
00591 
00592   for (int i=0; i<_topology->nbDomain();i++)
00593   {
00594     nb_plain_elems[i] = _mesh[i]->getNumberOfElements(entity, MED_EN::MED_ALL_ELEMENTS);
00595     int nb_elem = _mesh[i]->getNumberOfElements(entity,MED_EN::MED_POLYHEDRA);
00596     if (nb_elem>0)
00597     {
00598       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
00599       conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
00600     }
00601     else
00602     {
00603       conn_ip[i]=0;
00604       conn_index_ip[i]=0;
00605     }
00606   }
00607 
00608   connectivity_index.resize(nb_cells+1);
00609   connectivity_index[0]=1;
00610   for (int icell=0; icell<nb_cells; icell++)
00611   {
00612     const int* conn=conn_ip[ip[icell]]; 
00613     const int* conn_index=conn_index_ip[ip[icell]];
00614     connectivity_index[icell+1]=connectivity_index[icell]+
00615       conn_index[local[icell]]-conn_index[local[icell]-1];
00616 
00617     for (int inode=conn_index[local[icell]-1]; inode<conn_index[local[icell]]; inode++)
00618     {
00619       if ( conn[inode-1] == -1 )
00620         type_connectivity.push_back( -1 );
00621       else
00622         type_connectivity.push_back(_topology->convertNodeToGlobal(ip[icell],conn[inode-1]));
00623     }
00624 
00625   }
00626 
00627   delete[]local;
00628   delete[]ip;
00629 }
00630 
00641 void MESHCollection::write(const string& filename)
00642 {
00643   //building the connect zones necessary for writing joints
00644   cout<<"Building Connect Zones"<<endl;
00645   if (_topology->nbDomain()>1)
00646     buildConnectZones();
00647   cout <<"End of connect zones building"<<endl;
00648   //suppresses link with driver so that it can be changed for writing
00649   if (_driver!=0)delete _driver;
00650   _driver=0;
00651 
00652   char filenamechar[256];
00653   strcpy(filenamechar,filename.c_str());
00654   retrieveDriver()->write (filenamechar, _domain_selector);
00655 }
00656 
00659 MESHCollectionDriver* MESHCollection::retrieveDriver()
00660 {
00661   if (_driver==0)
00662   {
00663     switch(_driver_type)
00664     {
00665     case MedXML:
00666       _driver=new MESHCollectionMedXMLDriver(this);
00667       break;
00668     case MedAscii:
00669       _driver=new MESHCollectionMedAsciiDriver(this);
00670       break;
00671     default:
00672       throw MEDEXCEPTION("Unrecognized driver");
00673     }
00674   }
00675 
00676   return _driver;
00677 }
00678 
00679 
00683 MESHCollectionDriver* MESHCollection::getDriver() const
00684 {
00685   return _driver;
00686 }
00687 
00688 
00695 void MESHCollection::getTypeList(int* cell_list,int nb_cells,
00696                                  MED_EN::medEntityMesh entity,
00697                                  MED_EN::medGeometryElement* type_list) const 
00698 {
00699   MESSAGE_MED (" Beginning of getTypeList with entity "<<entity);
00700   int *local=new int[nb_cells];
00701   int *ip=new int[nb_cells];
00702   switch (entity)
00703   {
00704   case MED_EN::MED_CELL:
00705     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
00706     break;
00707   case MED_EN::MED_FACE:
00708   case MED_EN::MED_EDGE:
00709     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
00710     break;
00711   }
00712 
00713   for (int icell=0; icell<nb_cells; icell++)
00714   {
00715     type_list[icell]=_mesh[ip[icell]]->getElementType(entity,local[icell]);
00716   }
00717   delete[]local;
00718   delete[]ip;
00719   MESSAGE_MED("end of getTypeList");
00720 }
00721 
00722 
00723 
00736 void MESHCollection::getFaceConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
00737                                          MED_EN::medGeometryElement type, int* type_connectivity) const
00738 {
00739   int *local=new int[nb_cells];
00740   int *ip=new int[nb_cells];
00741   switch (entity)
00742   {
00743   case MED_EN::MED_CELL:
00744     _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
00745     break;
00746   case MED_EN::MED_FACE:
00747   case MED_EN::MED_EDGE:
00748     _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
00749     break;
00750   }
00751 
00752 
00753   int nbface_per_type;
00754   switch (type){
00755   case 308:
00756     nbface_per_type=6;
00757     break;
00758   case 304:
00759     nbface_per_type=4;
00760     break;
00761   case 306:
00762     nbface_per_type=5;
00763     break;
00764   }
00765 
00766   vector<int> number_of_types_array(_topology->nbDomain(),0);
00767   for (int i=0; i<_topology->nbDomain(); i++)
00768     number_of_types_array[i]=_mesh[i]->getNumberOfTypes(entity);
00769 
00770   //defining a connectivity table for different domains 
00771   vector  <const int*> conn_ip(_topology->nbDomain());
00772   for (int i=0; i<_topology->nbDomain();i++)
00773   {
00774     int nb_elem = _mesh[i]->getNumberOfElements(entity,type);
00775     if (nb_elem>0)
00776       conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_DESCENDING,entity,type);
00777     else
00778       conn_ip[i]=0;                     
00779   }
00780 
00781   for (int icell=0; icell<nb_cells; icell++)
00782   {
00783     int number_of_types = number_of_types_array[ip[icell]];
00784     const MEDMEM::CELLMODEL* types =  _mesh[ip[icell]]->getCellsTypes(entity);
00785     int type_offset=0;
00786     for (int itype=0; itype< number_of_types; itype++)
00787     {
00788       if (types[itype].getType() < type)
00789         type_offset += _mesh[ip[icell]]->getNumberOfElements(entity,types[itype].getType());
00790     }
00791     const int* conn=conn_ip[ip[icell]];           
00792     for (int iface=0; iface<nbface_per_type; iface++)
00793     {
00794       type_connectivity[icell*nbface_per_type+iface] = _topology->convertFaceToGlobal
00795         (ip[icell], abs(conn[(local[icell] - type_offset - 1) * nbface_per_type + iface]));           
00796     }
00797   }
00798 
00799   delete[]local;
00800   delete[]ip;
00801 }
00802 
00814 void MESHCollection::getCoordinates(int* node_list,int nb_nodes, double* coordinates) const
00815 {
00816   int* local=new int[nb_nodes];
00817   int* ip=new int[nb_nodes];
00818   int space_dimension= getSpaceDimension();
00819   _topology->convertGlobalNodeList(node_list,nb_nodes,local,ip);
00820   for (int i=0; i< nb_nodes; i++)
00821   {
00822     const double* coord=_mesh[ip[i]]->getCoordinates(MED_EN::MED_FULL_INTERLACE);
00823     for (int icoord=0; icoord<space_dimension; icoord++)
00824       coordinates[i*space_dimension+icoord]=coord[(local[i]-1)*space_dimension+icoord];
00825   }
00826   delete[]local;
00827   delete[] ip;
00828 }
00830 MED_EN::medEntityMesh MESHCollection::getSubEntity() const
00831 {
00832   return getMeshDimension()==3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE;
00833 }
00834 
00836 int MESHCollection::getSpaceDimension() const
00837 {
00838   return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getSpaceDimension();
00839 }
00841 int MESHCollection::getMeshDimension() const
00842 {
00843   return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
00844 }
00845 
00847 string MESHCollection::getSystem() const
00848 {
00849   return _i_non_empty_mesh < 0 ? "" : _mesh[_i_non_empty_mesh]->getCoordinatesSystem();
00850 }
00851 
00853 string MESHCollection::getMeshName() const
00854 {
00855   return _i_non_empty_mesh < 0 ? (_mesh[0] ? _mesh[0]->getName() : "") : _mesh[_i_non_empty_mesh]->getName();
00856 }
00857 
00858 vector<MEDMEM::MESH*>& MESHCollection::getMesh() 
00859 {
00860   return _mesh;
00861 }
00862 
00863 MEDMEM::MESH* MESHCollection::getMesh(int idomain) const
00864 {
00865   return _mesh[idomain];
00866 }
00867 
00868 vector<MEDMEM::CONNECTZONE*>& MESHCollection::getCZ()
00869 {
00870   return _connect_zones;
00871 }
00872 
00873 Topology* MESHCollection::getTopology() const 
00874 {
00875   return _topology;
00876 }
00877 
00878 void MESHCollection::setTopology(Topology* topo)
00879 {
00880   if (_topology!=0)
00881   {
00882     throw MED_EXCEPTION(STRING("Erreur : topology is already set"));
00883   }
00884   else
00885     _topology = topo;
00886 }
00887 
00888 void MESHCollection::setIndivisibleGroup(const string& name)
00889 {
00890   _indivisible_regions.push_back(name);
00891 
00892 }
00893 
00905 void MESHCollection::treatIndivisibleRegions(int* indivisible_tag)
00906 {
00907   //tag 0 is positioned on all the cells that are not affected by these tags
00908   for (int i=0; i<_topology->nbCells(); i++)
00909     indivisible_tag[i]=0;
00910 
00911   //treating cell groups    
00912   for (int idomain=0; idomain<_topology->nbDomain();idomain++)
00913     for (int igroup=0; igroup<_mesh[idomain]->getNumberOfGroups(MED_EN::MED_CELL); igroup++)
00914       for (unsigned i=0; i<_indivisible_regions.size(); i++)
00915       {
00916         const MEDMEM::GROUP* group = _mesh[idomain]->getGroup(MED_EN::MED_CELL,igroup+1);
00917         string groupname = group->getName();
00918         if (trim(groupname)==trim(_indivisible_regions[i]))
00919         {
00920           int nbcells=group->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
00921           const int* numbers=group->getNumber(MED_EN::MED_ALL_ELEMENTS);
00922           int* global=new int[nbcells];
00923           _topology->convertCellToGlobal(idomain,numbers,nbcells,global);
00924           for (int icell=0; icell<nbcells; icell++)
00925             indivisible_tag[global[icell]-1]=i+1;
00926           delete[] global;
00927         } 
00928       }
00929 }
00930 
00931 //================================================================================
00935 //================================================================================
00936 
00937 template<class TID2CONN>
00938 void MESHCollection::fillGlobalConnectivity(TID2CONN & node2cell, TID2CONN & cell2node )
00939 {
00940   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
00941   {
00942     if ( !_mesh[idomain] ) continue;
00943     //  MED_EN::medGeometryElement* type_array;
00944     int nb_cells= _topology->nbCells(idomain);
00945     int* cell_list = new int[nb_cells];
00946 
00947     //retrieving global id list
00948     _topology->getCellList(idomain, cell_list);
00949 
00950     int nb_plain_cells = _mesh[idomain]->getNumberOfElements(MED_EN::MED_CELL,
00951                                                              MED_EN::MED_ALL_ELEMENTS);
00952     if (nb_plain_cells >0)
00953     {          
00954       const int* conn_index = _mesh[idomain]->getConnectivityIndex(MED_EN::MED_NODAL,
00955                                                                    MED_EN::MED_CELL);
00956 
00957       const int* conn =   _mesh[idomain]->getConnectivity(MED_EN::MED_NODAL,
00958                                                           MED_EN::MED_CELL,
00959                                                           MED_EN::MED_ALL_ELEMENTS);                                                                    
00960       int nbnodes = conn_index[nb_plain_cells]-1;
00961       int* global_nodes =new int [nbnodes];
00962       _topology->convertNodeToGlobal(idomain, conn, nbnodes, global_nodes);
00963       for (int icell=0; icell< nb_plain_cells; icell++)
00964       {
00965         for (int inode=conn_index[icell]; inode < conn_index[icell+1]; inode++)
00966         {
00967           int node_global_id = global_nodes[inode-1];
00968           if ( node_global_id > 0 )
00969           {
00970             int cell_global_id = cell_list[icell];
00971             cell2node [cell_global_id].push_back(node_global_id);
00972             node2cell [node_global_id].push_back(cell_global_id);
00973           }
00974         }
00975       }
00976       delete[] global_nodes;
00977     }
00978 
00979     delete[] cell_list;
00980   }
00981 }
00982 
00990 void MESHCollection::buildCellGraph(MEDMEM::MEDSKYLINEARRAY* & array,int *& edgeweights )
00991 {
00992 
00993   int cell_number=1;
00994   int node_number=1;
00995   for (int i=0; i<_topology->nbDomain(); i++)
00996   {
00997     cell_number+=_topology->getCellNumber(i);
00998     node_number+=_topology->getNodeNumber(i);
00999   }
01000   //list of cells for a given node
01001   //vector< vector<int> > node2cell(node_number);
01002   map< int, vector<int> > node2cell;
01003 
01004   //list of nodes for a given cell
01005   //vector< vector <int> > cell2node(cell_number);
01006   map< int, vector <int> > cell2node;
01007 
01008   //  map<MED_EN::medGeometryElement,int*> type_cell_list;
01009 
01010   //tagging for the indivisible regions
01011   int* indivisible_tag=0;
01012   bool has_indivisible_regions=false;
01013   if (!_indivisible_regions.empty())
01014   {
01015     has_indivisible_regions=true;
01016     indivisible_tag=new int[_topology->nbCells()];
01017     treatIndivisibleRegions(indivisible_tag);
01018   }
01019 
01020   fillGlobalConnectivity(node2cell, cell2node );
01021 
01022   cout << "beginning of skyline creation"<<endl;
01023   //creating the MEDMEMSKYLINEARRAY containing the graph
01024 
01025   int* size = new int[_topology->nbCells()];
01026   int** temp=new int*[_topology->nbCells()];
01027   int** temp_edgeweight=0;
01028   if (has_indivisible_regions)
01029     temp_edgeweight=new int*[_topology->nbCells()];
01030 
01031   int cell_glob_shift = 0;
01032 
01033   // Get connection to cells on other procs
01034   multimap< int, int > loc2dist; // global cell ids on this proc -> other proc cells
01035   if ( isParallelMode() )
01036   {
01037     cell_glob_shift = _domain_selector->getProcShift();
01038 
01039     set<int> loc_domains; // domains on this proc
01040     for ( int idom = 0; idom < _mesh.size(); ++idom )
01041       if ( _mesh[ idom ] )
01042         loc_domains.insert( idom );
01043 
01044     for ( int idom = 0; idom < _mesh.size(); ++idom )
01045     {
01046       if ( !_mesh[idom] ) continue;
01047       vector<int> loc2glob_corr; // pairs of corresponding cells (loc_loc & glob_dist)
01048       retrieveDriver()->readLoc2GlobCellConnect(idom, loc_domains, _domain_selector, loc2glob_corr);
01049       //MEDMEM::STRING out;
01050       for ( int i = 0; i < loc2glob_corr.size(); i += 2 )
01051       {
01052         int glob_here  = _topology->convertCellToGlobal(idom,loc2glob_corr[i]); 
01053         int glob_there = loc2glob_corr[i+1]; 
01054         loc2dist.insert ( make_pair( glob_here, glob_there));
01055         //out << glob_here << "-" << glob_there << " ";
01056       }
01057       //cout << "\nRank " << _domain_selector->rank() << ": BndCZ: " << out << endl;
01058     }
01059   }
01060 
01061   //going across all cells
01062 
01063   map<int,int> cells_neighbours;
01064   for (int i=0; i< _topology->nbCells(); i++)
01065   {
01066 
01067 
01068     vector<int> cells(50);
01069 
01070     for (unsigned inode=0; inode< cell2node[i+1].size(); inode++)
01071     {
01072       int nodeid=cell2node[i+1][inode];
01073       for (unsigned icell=0; icell<node2cell[nodeid].size();icell++)
01074         cells_neighbours[node2cell[nodeid][icell]]++;
01075     }
01076     size[i]=0;
01077     int dimension = getMeshDimension();
01078     cells.clear();
01079 
01080     for (map<int,int>::const_iterator iter=cells_neighbours.begin(); iter != cells_neighbours.end(); iter++)  
01081     {
01082       if (iter->second >= dimension && iter->first != i+1) 
01083       {
01084         cells.push_back(iter->first + cell_glob_shift);
01085         //       cells[isize++]=iter->first;
01086       }
01087     }
01088     // add neighbour cells from distant domains
01089     multimap< int, int >::iterator loc_dist = loc2dist.find( i+1 );
01090     for (; loc_dist!=loc2dist.end() && loc_dist->first==( i+1 ); ++loc_dist )
01091       cells.push_back( loc_dist->second );
01092 
01093     size[i]=cells.size();
01094 
01095     temp[i]=new int[size[i]];
01096     if (has_indivisible_regions)
01097       temp_edgeweight[i]=new int[size[i]];
01098     //
01099     int itemp=0;
01100 
01101     for (vector<int>::const_iterator iter=cells.begin(); iter!=cells.end();iter++)
01102     {
01103       temp[i][itemp]=*iter;
01104       if (has_indivisible_regions)
01105       {
01106         int tag1 = indivisible_tag[(i+1)-1];
01107         int tag2 = indivisible_tag[*iter-1];
01108         if (tag1==tag2 && tag1!=0)
01109           temp_edgeweight[i][itemp]=_topology->nbCells()*100000;
01110         else
01111           temp_edgeweight[i][itemp]=1;
01112       } 
01113       itemp++;
01114     }
01115     cells_neighbours.clear();
01116   }
01117   cout <<"end of graph definition"<<endl;
01118   int* index=new int[_topology->nbCells()+1];
01119   index[0]=1;
01120   for (int i=0; i<_topology->nbCells(); i++)
01121     index[i+1]=index[i]+size[i];
01122 
01123   node2cell.clear();
01124   cell2node.clear();
01125   if (indivisible_tag!=0) delete [] indivisible_tag;
01126 
01127   //SKYLINEARRAY structure holding the cell graph
01128   array= new MEDMEM::MEDSKYLINEARRAY(_topology->nbCells(),index[_topology->nbCells()]-index[0]);
01129   array->setIndex(index);
01130 
01131   for (int i=0; i<_topology->nbCells(); i++)
01132   {
01133     array->setI(i+1,temp[i]);
01134     delete[]temp[i];
01135   }
01136 
01137   if (has_indivisible_regions)
01138   {
01139     edgeweights=new int[array->getLength()];
01140     for (int i=0; i<_topology->nbCells(); i++)
01141     {
01142       for (int j=index[i]; j<index[i+1];j++)
01143         edgeweights[j-1]=temp_edgeweight[i][j-index[i]];
01144       delete[] temp_edgeweight[i];  
01145     }
01146     delete[]temp_edgeweight;
01147   }
01148   delete[] index;
01149   delete[] temp;
01150   delete[] size;
01151 
01152   cout<< "end of graph creation"<<endl;
01153 }
01154 
01161 Topology* MESHCollection::createPartition(int nbdomain, 
01162                                           Graph::splitter_type split, 
01163                                           const string& options_string,
01164                                           int* user_edge_weights,
01165                                           int* user_vertices_weights)
01166 {
01167   if (nbdomain <1) throw MEDEXCEPTION("Number of subdomains must be >0");
01168   MEDMEM::MEDSKYLINEARRAY* array=0;
01169   int* edgeweights=0;
01170 
01171   MESSAGE_MED("Building cell graph");
01172   buildCellGraph(array,edgeweights);
01173 
01174   switch (split)
01175   {
01176   case Graph::METIS:
01177 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
01178     _cell_graph=boost::shared_ptr<Graph>(new METISGraph(array,edgeweights));
01179 #else
01180     throw MEDEXCEPTION("METIS Graph is not available. Check your products, please.");
01181 #endif
01182     break;
01183   case Graph::SCOTCH:
01184 #ifdef MED_ENABLE_SCOTCH
01185     _cell_graph=boost::shared_ptr<Graph>(new SCOTCHGraph(array,edgeweights));
01186 #else
01187     throw MEDEXCEPTION("SCOTCH Graph is not available. Check your products, please.");
01188 #endif
01189     break;
01190   }
01191 
01193   if (user_edge_weights!=0) 
01194     _cell_graph->setEdgesWeights(user_edge_weights);
01195   if (user_vertices_weights!=0)
01196     _cell_graph->setVerticesWeights(user_vertices_weights);
01197 
01198   MESSAGE_MED("Partitioning graph");
01199   _cell_graph->partGraph(nbdomain,options_string,_domain_selector);
01200 
01201   // DEBUG
01202 //   MEDMEM::STRING out("RESULT GRAPH #");
01203 //   out << (_domain_selector?_domain_selector->rank():0) << ": ";
01204 //   const int* part = _cell_graph->getPart();
01205 //   int n = _cell_graph->nbVertices();
01206 //   for ( int e=0; e < n; ++e )
01207 //     out << part[e] <<" ";
01208 //   cout << out << endl;
01209   
01210 
01211   MESSAGE_MED("Building new topology");
01212   //_cell_graph is a shared pointer 
01213   Topology* topology = new ParallelTopology (_cell_graph, nbdomain, getMeshDimension());
01214 
01215   //cleaning
01216   if (edgeweights!=0) delete[] edgeweights;
01217   //if (array!=0) delete array;
01218   MESSAGE_MED("End of partition creation");
01219   return topology;
01220 }
01221 
01228 Topology* MESHCollection::createPartition(const int* partition)
01229 {
01230   MEDMEM::MEDSKYLINEARRAY* array=0;
01231   int* edgeweights=0;
01232 
01233   buildCellGraph(array,edgeweights);
01234 
01235   set<int> domains;
01236   for (int i=0; i<_topology->nbCells(); i++)
01237   {
01238     domains.insert(partition[i]);
01239   }
01240   int nbdomain=domains.size();
01241 
01242   _cell_graph=boost::shared_ptr<Graph>(new UserGraph(array, partition, _topology->nbCells()));
01243 
01244   //_cell_graph is a shared pointer 
01245   Topology* topology = new ParallelTopology (_cell_graph, nbdomain, getMeshDimension());
01246 
01247   //if (array!=0) delete array;
01248   return topology;
01249 }
01250 
01251 
01261 // void MESHCollection::buildConnectZones(int idomain)
01262 // {
01263 //   // constructing node/node correspondencies
01264 //   vector<MEDMEM::MEDSKYLINEARRAY*> node_node_correspondency;
01265 //   node_node_correspondency.resize(_topology->nbDomain());
01266 
01267 //   cout << "Computing node/node corresp"<<endl;
01268 
01269 //   _topology->computeNodeNodeCorrespondencies(idomain, node_node_correspondency );
01270 
01271 //   for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
01272 //   {
01273 //     // on regarde si une correspondance noeud/noeud a été trouvée 
01274 //     // entre idomain et idistant
01275 //     // si oui, on crée une connectzone
01276 //     if (node_node_correspondency[idistant]!=0)
01277 //     {
01278 //       MEDMEM::CONNECTZONE* cz= new MEDMEM::CONNECTZONE();
01279 //       cz->setLocalMesh(_mesh[idomain]);
01280 //       cz->setDistantMesh(_mesh[idistant]);
01281 //       cz->setLocalDomainNumber(idomain);
01282 //       cz->setDistantDomainNumber(idistant);
01283 //       cz-> setName ("Connect zone defined by SPLITTER");
01284 //       cz->setNodeCorresp(node_node_correspondency[idistant]);
01285 //       _connect_zones.push_back(cz);  
01286 //     }
01287 //   }
01288 //   cout << "Computing node/node corresp"<<endl;
01289 
01290 //   vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency;
01291 //   cell_cell_correspondency.resize(_topology->nbDomain());
01292 //   _topology->computeCellCellCorrespondencies(idomain, cell_cell_correspondency, _cell_graph.get());
01293 
01294 //   for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
01295 //   {
01296 //     //the connect zone has been created by the node/node computation
01297 //     if (cell_cell_correspondency[idistant]!=0)
01298 //     {
01299 //       MEDMEM::CONNECTZONE* cz=0;
01300 //       for (int icz=0; icz<_connect_zones.size();icz++)
01301 //         if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
01302 //             _connect_zones[icz]->getDistantDomainNumber()==idistant)
01303 //           cz = _connect_zones[icz];
01304 //       if (cz!=0) 
01305 //         cz->setEntityCorresp(MED_EN::MED_CELL,MED_EN::MED_CELL, cell_cell_correspondency[idistant]);
01306 //       else 
01307 //         throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");   
01308 //       //delete cell_cell_correspondency[idistant];
01309 //     }
01310 
01311 //   }
01312 // }
01313 
01314 //================================================================================
01321 //================================================================================
01322 
01323 void MESHCollection::addJointGroup(const std::vector<int>& loc_face_ids, int idomain, int idistant)
01324 {
01325   MEDMEM::MESHING* meshing = dynamic_cast<MEDMEM::MESHING*> (_mesh[idomain]);
01326   MED_EN::medEntityMesh constituent_entity = getSubEntity();
01327 
01328   MEDMEM::STRING jointname("joint_");
01329   jointname<<idistant+1;
01330 
01331   MEDMEM::GROUP * tmp_grp = new GROUP, * joint_group = tmp_grp;
01332   // try to find already present group with such a name
01333   //  vector<MEDMEM::GROUP*> groups = meshing->getGroups( constituent_entity );
01334   //  for ( int g = 0; g < groups.size(); ++g )
01335   //    if ( groups[g]->getName() == jointname.str() )
01336   //    {
01337   //      joint_group = groups[g];
01338   //      break;
01339   //    }
01340   // assure uniqueness of group name
01341   bool unique = false;
01342   vector<MEDMEM::GROUP*> groups = meshing->getGroups( constituent_entity );
01343   do
01344   {
01345     unique = true;
01346     for ( int g = 0; unique && g < groups.size(); ++g )
01347       unique = ( groups[g]->getName() != jointname );
01348     if ( !unique )
01349       jointname << "_" << idomain+1;
01350   }
01351   while ( !unique );
01352   joint_group->setName(jointname);
01353   joint_group->setMesh(meshing);
01354   joint_group->setEntity(constituent_entity);
01355   map<MED_EN::medGeometryElement, vector<int> > joint_types;
01356 
01357   int nbfaces = loc_face_ids.size();
01358   for (int i=0; i<nbfaces; i++)
01359   {    
01360     MED_EN::medGeometryElement type = meshing->getElementType(constituent_entity,loc_face_ids[i]);
01361     joint_types[type].push_back(loc_face_ids[i]);
01362   }
01363   joint_group->setNumberOfGeometricType(joint_types.size());
01364   MED_EN::medGeometryElement* types=new MED_EN::medGeometryElement[joint_types.size()];
01365   int* nb_in_types=new int[joint_types.size()];
01366   int* group_index=new int[joint_types.size()+1];
01367 
01368   group_index[0]=1;
01369   int itype=0;
01370   int iface =0;
01371   int* group_value=new int[nbfaces];
01372   for (map<MED_EN::medGeometryElement, vector<int> >::const_iterator iterj=joint_types.begin();
01373        iterj != joint_types.end();
01374        iterj++)
01375   {
01376     nb_in_types[itype]=(iterj->second).size();
01377     types[itype]=iterj->first;
01378     itype++;
01379     group_index[itype]=group_index[itype-1]+(iterj->second).size();
01380     for (int i=0; i<  (iterj->second).size(); i++)
01381       group_value[iface++]=(iterj->second)[i];
01382   }
01383   joint_group->setGeometricType(types);
01384   joint_group->setNumberOfElements(nb_in_types);
01385   joint_group->setNumber(group_index, group_value, /*shallowCopy=*/true);
01386   delete[] types;
01387   delete[] nb_in_types;
01388 
01389   if ( joint_group == tmp_grp )
01390     meshing->addGroup(*tmp_grp);
01391   tmp_grp->removeReference();
01392 }
01393 
01398 void MESHCollection::buildConnectZones()
01399 {
01400   vector <map <MED_EN::medGeometryElement, vector<MEDSPLITTER_FaceModel*> > > face_map(_topology->nbDomain());
01401   map< pair<int,int>, MEDMEM::MEDSKYLINEARRAY*> cell_corresp_here;
01402 
01403   MED_EN::medEntityMesh constituent_entity = getSubEntity();
01404 
01405   if ( isParallelMode() )
01406   {
01407     buildConnectZonesBetweenProcs(face_map, cell_corresp_here);
01408   }
01409 
01410   cout << "Computing node/node corresp"<<endl;
01411 
01412   //Creating nodes
01413   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01414   {
01415 
01416     // constructing node/node correspondencies
01417     vector<MEDMEM::MEDSKYLINEARRAY*> node_node_correspondency(_topology->nbDomain());
01418     _topology->computeNodeNodeCorrespondencies(idomain, node_node_correspondency );
01419 
01420     for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
01421     {
01422       // on regarde si une correspondance noeud/noeud a été trouvée 
01423       // entre idomain et idistant
01424       // si oui, on crée une connectzone
01425       if (node_node_correspondency[idistant]!=0)
01426       {
01427         MEDMEM::CONNECTZONE* cz= new MEDMEM::CONNECTZONE();
01428         cz->setLocalMesh(_mesh[idomain]);
01429         cz->setDistantMesh(_mesh[idistant]);
01430         cz->setLocalDomainNumber(idomain);
01431         cz->setDistantDomainNumber(idistant);
01432         cz-> setName ("Connect zone defined by SPLITTER");
01433         cz->setNodeCorresp(node_node_correspondency[idistant]);
01434         _connect_zones.push_back(cz);  
01435       }
01436     }
01437   }
01438   cout << "Computing face corresp"<<endl;
01439 
01440   //creating faces if required 
01441   if (_subdomain_boundary_creates)
01442   {
01443     int global_face_id = _topology->getFaceNumber()+1;
01444     //int global_face_id = _topology->getMaxGlobalFace()+1;
01445 
01446     map <pair<int,int>, vector<int> > faces_in_joint;
01447 
01448     if ( !isParallelMode() )
01449       // taking faces that are already present in the mesh into account
01450       for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01451       {
01452         getFaces(idomain,face_map[idomain]); 
01453       }
01454 
01455     // creating faces that are located at the interface between
01456     // subdomains 
01457 
01458     vector <int> nb_added_groups( _topology->nbDomain(), 0 );
01459 
01460     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01461     {
01462       vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency( _topology->nbDomain() );
01463       if ( !isParallelMode() )
01464         _topology->computeCellCellCorrespondencies(idomain, cell_cell_correspondency, _cell_graph.get());
01465 
01466       for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
01467       {
01468         if (idistant <= idomain) continue;
01469 
01470         MEDMEM::MEDSKYLINEARRAY* cell_correspondency = 0;
01471         if ( isParallelMode() )
01472           cell_correspondency = cell_corresp_here[ make_pair (idomain,idistant)];
01473         else
01474           cell_correspondency = cell_cell_correspondency[idistant];
01475 
01476         //the connect zone has been created by the node/node computation
01477 
01478         if ( cell_correspondency )
01479         {
01480           int nbcells      = cell_correspondency->getNumberOf();
01481           const int* index = cell_correspondency->getIndex();
01482           const int* value = cell_correspondency->getValue();
01483           if ( isParallelMode() )
01484             global_face_id = _domain_selector->getFisrtGlobalIdOfSubentity( idomain, idistant );
01485 
01486           for (int ilocal=0; ilocal<nbcells; ilocal++)
01487           { 
01488             for (int icelldistant = index[ilocal]; icelldistant < index[ilocal+1]; icelldistant++)
01489             {
01490               int distant_id = value[icelldistant-1];
01491               MEDSPLITTER_FaceModel* face = getCommonFace(idomain,ilocal+1,idistant,distant_id,global_face_id);
01492               face_map[idomain][face->getType()].push_back(face);
01493               MEDSPLITTER_FaceModel* face2 = getCommonFace(idistant,distant_id,idomain, ilocal+1,global_face_id);
01494               face_map[idistant][face->getType()].push_back(face2);
01495               faces_in_joint[make_pair(idomain,idistant)].push_back(global_face_id);
01496               global_face_id++;
01497             } 
01498           }
01499         }
01500 
01501       }
01502       //cleaning up
01503       for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
01504         delete cell_cell_correspondency[idistant];
01505     }
01506 
01507 
01508     _topology->recreateFaceMapping(face_map);
01509 
01510     //transforming the face_map into a constituent entity connectivity
01511     for (int idomain=0; idomain< _topology->nbDomain();idomain++) 
01512     {
01513       int nbtypes = face_map[idomain].size();
01514       vector<medGeometryElement> types;
01515       vector <int> nb_elems;
01516       vector <int*> conn;
01517 
01518       MEDMEM::MESHING* meshing = dynamic_cast<MEDMEM::MESHING*> (_mesh[idomain]);
01519       if ( !meshing->getConnectivityptr() )
01520         continue; // no cells in idomain
01521 
01522       for (map <medGeometryElement, vector<MEDSPLITTER_FaceModel*> >::const_iterator iter= face_map[idomain].begin();
01523            iter != face_map[idomain].end(); iter ++)
01524       {
01525         types.push_back(iter->first);
01526         int nb_elem_in_type = (iter->second).size();
01527         nb_elems.push_back(nb_elem_in_type);
01528         int nb_node_per_type=(iter->first)%100;
01529         int* connectivity= new int [nb_node_per_type*nb_elem_in_type];
01530         for (int ielem=0; ielem<nb_elem_in_type; ielem++)
01531         {
01532           for (int inode=0;  inode<nb_node_per_type; inode++)
01533             connectivity[ielem*nb_node_per_type+inode]=(*(iter->second)[ielem])[inode];
01534         }
01535         conn.push_back(connectivity);
01536 
01537       }
01538       //setting the faces in the mesh
01539       meshing->setNumberOfTypes(nbtypes,constituent_entity);
01540       meshing->setTypes(&types[0],constituent_entity);
01541       meshing->setNumberOfElements(&nb_elems[0],constituent_entity);
01542 
01543       for (int itype=0; itype<nbtypes; itype++)
01544       {
01545         meshing->setConnectivity( constituent_entity, types[itype], conn[itype] );
01546         delete[]conn[itype];
01547       }
01548       for (int idistant =0; idistant<_topology->nbDomain(); idistant++)
01549       {
01550         map <pair<int,int>, vector<int> >::iterator iter;
01551         iter = faces_in_joint.find(make_pair(idomain,idistant));
01552         if (iter == faces_in_joint.end())
01553         {
01554           iter = faces_in_joint.find (make_pair(idistant,idomain));
01555           if (iter == faces_in_joint.end()) 
01556             continue;
01557         }
01558 
01559         int nbfaces = (iter->second).size();   
01560         vector<int> face_joint(nbfaces*2);
01561         MEDMEM::CONNECTZONE* cz=0;
01562         for (unsigned icz=0; icz<_connect_zones.size();icz++)
01563           if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
01564               _connect_zones[icz]->getDistantDomainNumber()==idistant)
01565             cz = _connect_zones[icz];
01566 
01567         int nbtotalfaces= _topology->getFaceNumber(idomain);
01568 
01569         //creating arrays for the MEDSKYLINEARRAY structure containing the joint
01570         int* index =new int[nbtotalfaces+1];
01571         for (int i=0; i<nbtotalfaces+1;i++)
01572           index[i]=0;
01573         int*value=new int[nbfaces];
01574 
01575         map<int,int> faces;
01576         vector<int> local_faces( nbfaces );
01577         for (int iface=0; iface<nbfaces; iface++)
01578         {
01579           int iglobal = (iter->second)[iface];
01580           int localid=_topology->convertGlobalFace(iglobal,idomain);
01581           int distantid=_topology->convertGlobalFace(iglobal,idistant);
01582           faces.insert(make_pair(localid,distantid));
01583           local_faces[iface]=localid;
01584         }
01585 
01586         int iloc=0;
01587         index[0]=1;
01588         for (map<int,int>::const_iterator iter=faces.begin(); 
01589              iter != faces.end();
01590              iter++)
01591         {
01592           index[iter->first]=1;
01593           value[iloc++]=iter->second;            
01594         }
01595 
01596         for (int i=0; i<nbtotalfaces;i++)
01597           index[i+1]+=index[i];
01598         bool shallowcopy=true;  
01599         MEDMEM::MEDSKYLINEARRAY* skarray=new MEDMEM::MEDSKYLINEARRAY(nbtotalfaces,nbfaces,index,value,shallowcopy);  
01600 
01601         if (cz!=0)  
01602           cz->setEntityCorresp(constituent_entity,constituent_entity,skarray);              
01603         else 
01604           throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");            
01605         // Creating a group of the faces constituting the joint
01606         addJointGroup( local_faces, idomain, idistant );
01607         nb_added_groups[ idomain ]++;
01608       }
01609     }
01610 
01611     if ( isParallelMode() )
01612     {
01613       // Now all faces have got local ids and we can receive local ids from other procs.
01614       // Set face/face data to zones with other procs and create a group
01615       for (int icz=0; icz<_connect_zones.size();icz++)
01616       {
01617         MEDMEM::CONNECTZONE* cz=_connect_zones[icz];
01618         if ( _domain_selector->isMyDomain( cz->getDistantDomainNumber()) ) continue;
01619         
01620         int glob_id = _domain_selector->getFisrtGlobalIdOfSubentity( cz->getLocalDomainNumber(),
01621                                                                      cz->getDistantDomainNumber());
01622         int nb_cz_faces = _domain_selector->getNbCellPairs( cz->getDistantDomainNumber(),
01623                                                             cz->getLocalDomainNumber());
01624         vector< int > loc_ids_here( nb_cz_faces );
01625         for ( int i = 0; i < nb_cz_faces; ++i )
01626           loc_ids_here[i] = _topology->convertGlobalFace(glob_id++,cz->getLocalDomainNumber());
01627 
01628         int* loc_ids_dist = _domain_selector->exchangeSubentityIds( cz->getLocalDomainNumber(),
01629                                                                     cz->getDistantDomainNumber(),
01630                                                                     loc_ids_here );
01631         int nb_faces_here= _topology->getFaceNumber(cz->getLocalDomainNumber());
01632         int* face_index = new int[ nb_faces_here+1 ];
01633         face_index[0]=1;
01634         for ( int loc_id = 0, i = 0; loc_id < nb_faces_here; ++loc_id)
01635         {
01636           face_index[ loc_id+1 ] = face_index[ loc_id ];
01637           if ( i < loc_ids_here.size() && loc_ids_here[i] == loc_id+1 )
01638           {
01639             face_index[ loc_id+1 ]++;
01640             i++;
01641           }
01642         }
01643         MEDMEM::MEDSKYLINEARRAY* skarray=
01644           new MEDMEM::MEDSKYLINEARRAY(nb_faces_here, nb_cz_faces, face_index, loc_ids_dist, true);
01645         cz->setEntityCorresp(constituent_entity,constituent_entity,skarray);
01646 
01647         addJointGroup( loc_ids_here, cz->getLocalDomainNumber(), cz->getDistantDomainNumber());
01648         nb_added_groups[ cz->getLocalDomainNumber() ]++;
01649       }
01650     }
01651 
01652     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01653     {
01654       // delete face_map
01655       for (map <medGeometryElement, vector<MEDSPLITTER_FaceModel*> >::const_iterator iter= face_map[idomain].begin();
01656            iter != face_map[idomain].end(); iter ++)
01657         for (unsigned i=0; i<(iter->second).size();i++)
01658           delete (iter->second)[i];
01659 
01660       if ( nb_added_groups[ idomain ] > 0 &&
01661            _mesh[idomain]->getNumberOfFamilies( constituent_entity ) > 0 )
01662         // needed because if there were face families before, driver won't
01663         // create families from just added groups (see MEDMEM_MedMeshDriver.cxx:3330),
01664         // actually it is a bug of driver - it must check presence of groups in families
01665         _mesh[idomain]->createFamilies(); 
01666     }
01667   }
01668 
01669   if ( isParallelMode() )
01670     // Excange info on types of constituent_entity needed while writing joints
01671     // to get ids local in geom type for distant procs
01672     _domain_selector->gatherEntityTypesInfo( _mesh, constituent_entity );
01673 
01674   cout << "Computing cell/cell corresp"<<endl;
01675 
01676   //Creating cell/cell correspondencies
01677   for (int idomain=0;idomain<_topology->nbDomain();idomain++)
01678   {
01679     vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency( _topology->nbDomain() );
01680     if ( !isParallelMode() )
01681       _topology->computeCellCellCorrespondencies(idomain,cell_cell_correspondency,_cell_graph.get());
01682 
01683     for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
01684     {
01685       MEDMEM::MEDSKYLINEARRAY* cell_correspondency = 0;
01686       if ( isParallelMode() )
01687         cell_correspondency = cell_corresp_here[ make_pair (idomain,idistant)];
01688       else
01689         cell_correspondency = cell_cell_correspondency[idistant];
01690 
01691       //the connect zone has been created by the node/node computation
01692       if ( cell_correspondency )
01693       {
01694         MEDMEM::CONNECTZONE* cz=0;
01695         for (unsigned icz=0; icz<_connect_zones.size();icz++)
01696           if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
01697               _connect_zones[icz]->getDistantDomainNumber()==idistant)
01698             cz = _connect_zones[icz];
01699         if (cz!=0)  
01700           cz->setEntityCorresp(MED_EN::MED_CELL,MED_EN::MED_CELL, cell_correspondency);
01701         else 
01702           throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");   
01703       }
01704     }
01705   }
01706 }
01707 
01712 void MESHCollection::buildConnectZonesBetweenProcs(TGeom2FacesByDomian & face_map,
01713                                                    map< pair<int,int>, MEDMEM::MEDSKYLINEARRAY*> & cell_cell_correspondency_here)
01714 {
01715   using namespace MED_EN;
01716 
01717   // graph over all procs
01718   auto_ptr<Graph> global_graph( _domain_selector->gatherGraph( _cell_graph.get() ));
01719 
01720   vector< vector< JointExchangeData > > joints_of_domain( _topology->nbDomain() );
01721 
01722   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01723   {
01724     if ( !_domain_selector->isMyDomain( idomain )) continue;
01725 
01726     vector< JointExchangeData > & joints = joints_of_domain[ idomain ];
01727     joints.resize( _topology->nbDomain() );
01728 
01729     // Find corresponding cells on other procs
01730 
01731     const int* gra_index = global_graph->getGraph()->getIndex();
01732     const int* gra_value = global_graph->getGraph()->getValue();
01733     const int* partition = global_graph->getPart();
01734     const int dj = gra_index[0];
01735 
01736     vector< int > glob_cells_here( _topology->getCellNumber( idomain ));
01737     _topology->getCellList( idomain, & glob_cells_here[0]);
01738     for ( int loc_here = 0; loc_here < glob_cells_here.size(); ++loc_here )
01739     {
01740       int glob_here = glob_cells_here[ loc_here ];
01741       for ( int j = gra_index[ glob_here-1 ]; j < gra_index[ glob_here ]; ++j )
01742       {
01743         int glob_neighbor = gra_value[ j-dj ];
01744         int neighbor_dom = partition[ glob_neighbor-1 ];
01745         if ( neighbor_dom == idomain ) continue;
01746 
01747         if ( _domain_selector->isMyDomain( neighbor_dom ))
01748         {
01749           joints[ neighbor_dom ].addCellCorrespondence
01750             (_mesh[idomain], neighbor_dom, idomain, glob_neighbor, glob_here, loc_here + 1,
01751              _topology->convertGlobalCell(glob_neighbor).second );
01752         }
01753         else
01754         {
01755           joints[ neighbor_dom ].addCellCorrespondence
01756             (_mesh[idomain], neighbor_dom, idomain, glob_neighbor, glob_here, loc_here + 1 );
01757         }
01758       }
01759     }
01760   }
01761   global_graph.reset(); // free memory
01762 
01763   // set joints in a queue to exchange
01764   typedef map< int, JointExchangeData* > TOrderedJoints;
01765   TOrderedJoints queue;
01766   for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01767   {
01768     if ( !_domain_selector->isMyDomain( idomain )) continue;
01769 
01770     vector< JointExchangeData > & joints = joints_of_domain[ idomain ];
01771     for (int idist=0; idist<_topology->nbDomain(); ++idist )
01772     {
01773       JointExchangeData& joint = joints[idist];
01774 
01775       int nb_cell_pairs = joint.nbCellPairs();
01776       if ( nb_cell_pairs == 0 )
01777         continue;
01778       else
01779         _domain_selector->setNbCellPairs( nb_cell_pairs, idist, idomain );
01780 
01781       joint.setMeshes( idist, _mesh[idist], idomain, _mesh[idomain] );
01782 
01783       if ( _domain_selector->isMyDomain( idist ))
01784       {
01785         // a joint on this proc
01786         cell_cell_correspondency_here[ make_pair( idomain, idist )] = joint.makeCellCorrespArray();
01787       }
01788       else
01789       {
01790         // a joint with distant proc
01791         joint.setConnectivity( & ((MEDMEM::MeshFuse*)_mesh[idomain])->getNodeNumbers()[0] );
01792         int order = _domain_selector->jointId( idomain, idist );
01793         queue[ order ] = & joint;
01794       }
01795     }
01796   }
01797   // gather info on cell geom types needed to exchange joints
01798   _domain_selector->gatherEntityTypesInfo( _mesh, MED_EN::MED_CELL );
01799 
01800   // gather info on nb of sub-entities to compute their global numbers for joints
01801   _domain_selector->gatherNbOf( getSubEntity(), _mesh );
01802   _domain_selector->gatherNbCellPairs();
01803   if ( _subdomain_boundary_creates )
01804   {
01805     // taking faces that are already present in the mesh into account
01806     for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
01807       if ( _domain_selector->isMyDomain( idomain ))
01808         getFaces(idomain,face_map[idomain]);
01809   }
01810   else
01811   {
01812     face_map.clear(); // mark for the joint not to create faces
01813   }
01814 
01815   // exchange joint data with other procs and make CONNECTZONEs
01816   TOrderedJoints::iterator ord_joint = queue.begin();
01817   for ( ; ord_joint != queue.end(); ++ord_joint )
01818   {
01819     JointExchangeData* joint = ord_joint->second;
01820 
01821     _domain_selector->exchangeJoint( joint );
01822     if ( _subdomain_boundary_creates )
01823     {
01824       int first_sub_id = _domain_selector->getFisrtGlobalIdOfSubentity( joint->localDomain(),
01825                                                                         joint->distantDomain() );
01826       joint->setFisrtGlobalIdOfSubentity( first_sub_id );
01827     }
01828     _connect_zones.push_back ( joint->makeConnectZone( face_map ));
01829   }
01830 }
01831 
01834 void MESHCollection::castFamilies(const MESHCollection& old_collection)
01835 {
01836   vector <list<int> > element_array  (_topology->nbDomain());
01837 
01838   //loop on old domains to create groups out of the existing families
01839   if (_family_splitting)
01840     for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
01841       old_collection.getMesh(idomain)->createGroups();
01842 
01843   //definition of the entities array which 
01844   //defines the entities over which the information is cast
01845   MED_EN::medEntityMesh entities[3];
01846   entities[0]=MED_EN::MED_NODE;
01847   entities[1]=getSubEntity();
01848   entities[2]=MED_EN::MED_CELL;
01849 
01850   for (int ientity=0; ientity<=2;ientity++)
01851   {
01852 
01853     //int nbgroups = old_collection.getMesh(0)->getNumberOfGroups(entities[ientity]);
01854 
01855     map <string, set<int> > group_map;
01856     for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
01857     {
01858       if ( !old_collection.getMesh(idomain) ) continue;
01859       for (int igroup=0; igroup<old_collection.getMesh(idomain)->getNumberOfGroups(entities[ientity]); igroup++)
01860       {
01861         //retrieves a group
01862         MEDMEM::GROUP* group = (old_collection.getMesh(idomain)->getGroups(entities[ientity]))[igroup];
01863         //increments the number of groups if it is a new group
01864         //if (group_map.find(group->getName())==group_map.end())
01865 
01866         group_map[group->getName()].insert(idomain);
01867         //   group_map.insert(make_pair(group->getName(), idomain);
01868 
01869       }   
01870     }
01871     int nbgroups=group_map.size();
01872     vector <int> igroupold(old_collection._topology->nbDomain(),0);
01873     map<string,set<int> >::const_iterator iter=group_map.begin();
01874 
01875     for (int igroup=0; igroup<nbgroups; igroup++)
01876     {
01877       vector <const MEDMEM::SUPPORT*> old_supports(old_collection._topology->nbDomain());
01878       string group_name = iter->first;
01879       iter++; 
01880 
01881       //parameters stored for passing group description
01882       // from the old meshes to the new ones
01883 
01884       for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
01885       {
01886         //                for (set<int>::iterator iter=group_map[group_name].begin(); iter!=group_map[group_name].end(); iter++)
01887         //                cout << *iter<<" ";
01888         //                cout <<endl;
01889         if (group_map[group_name].find(idomain)==group_map[group_name].end()) continue;
01890 
01891         //retrieves the group igroup on domain idomain
01892         MEDMEM::GROUP* group = (old_collection.getMesh(idomain)->getGroups(entities[ientity]))[igroupold[idomain]];
01893         old_supports[idomain] = static_cast<const MEDMEM::SUPPORT*> (group);
01894         igroupold[idomain]++;
01895       }
01896 
01897       vector <MEDMEM::GROUP*>new_groups(_topology->nbDomain());
01898       vector <MEDMEM::SUPPORT*> new_supports(_topology->nbDomain());
01899       for (int i=0; i<_topology->nbDomain(); i++)
01900       {
01901         new_groups[i]=new MEDMEM::GROUP();
01902         new_supports[i]=static_cast<MEDMEM::SUPPORT*>(new_groups[i]);
01903       }
01904       castSupport(old_collection,old_supports,new_supports);      
01905 
01906       //creating new groups from the previous list of elements
01907       for (int idomain=0; idomain <_topology->nbDomain(); idomain++)
01908       {
01909         MEDMEM::MESHING* mesh_builder=static_cast<MEDMEM::MESHING*> (_mesh[idomain]);
01910         if ( new_supports[idomain] )
01911           mesh_builder->addGroup(*new_groups[idomain]);
01912       }
01913       //groups are copied by the addGroup method,
01914       //so they can be safely deleted here
01915       for (int i=0; i<_topology->nbDomain(); i++)
01916       {
01917         if ( new_supports[i] ) new_groups[i]->removeReference();
01918       }
01919 
01920     }// on groups
01921   }//on entities
01922 }
01923 
01924 
01925 void MESHCollection::castSupport(const MESHCollection& old_collection, vector<const MEDMEM::SUPPORT*>& old_support, vector<MEDMEM::SUPPORT*>& new_support)
01926 {
01927 
01928   if (old_collection._topology->nbDomain() != (int)old_support.size())
01929   {
01930     throw MED_EXCEPTION(STRING("Error : wrong call to MESHCollection::castSupport"));
01931   }
01932   vector <list<int> > element_array  (_topology->nbDomain());
01933 
01934   //parameters stored for passing description
01935   // from the old meshes to the new ones
01936   string name;
01937   string description;
01938   MED_EN::medEntityMesh entity;
01939   vector <string> support_name(1);
01940   support_name[0]="support";
01941   for (int inew=0; inew< _topology->nbDomain(); inew++)
01942     element_array[inew].clear();
01943 
01944   for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
01945   {
01946     //retrieves the group igroup on domain idomain
01947     const MEDMEM::SUPPORT* support = old_support[idomain];
01948     if (old_support[idomain]==0) continue;
01949     name = support->getName();
01950     description=support->getDescription();
01951     int nbelem = support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
01952     if (nbelem==0 && !_create_empty_groups) continue;
01953 
01954     int* list_of_elems;
01955     if (support->isOnAllElements())
01956     {
01957       list_of_elems = new int[nbelem];
01958       for (int i=0; i<nbelem;i++)
01959         list_of_elems[i]=i+1;
01960     }
01961     else
01962       list_of_elems = const_cast<int*> (support->getNumber(MED_EN::MED_ALL_ELEMENTS));
01963 
01964     int* array=new int[nbelem];
01965     int* ip=0;
01966     int* local=0;
01967     int* full_array=0;
01968     entity = support->getEntity();
01969     int size;
01970 
01971     switch (entity)
01972     {
01973     case MED_EN::MED_CELL :
01974       ip=new int[nbelem];
01975       local= new int[nbelem];
01976       size=nbelem;
01977       old_collection.getTopology()->convertCellToGlobal(idomain,list_of_elems,nbelem,array);
01978       _topology->convertGlobalCellList(array,nbelem,local,ip);
01979       for (int i=0; i<nbelem; i++)
01980         //              cell_arrays[ip[i]][local[i]]=id;
01981       {
01982         //          cout <<"(glob,ip,iloc)/nbelem"<<array[i]<<" "<<ip[i]<<" "<<local[i]<<"/"<<nbelem<<endl;
01983         element_array[ip[i]].push_back(local[i]);
01984       }
01985       break;
01986     case MED_EN::MED_FACE :
01987     case MED_EN::MED_EDGE :
01988       old_collection.getTopology()->convertFaceToGlobal(idomain,list_of_elems,nbelem,array);
01989       _topology->convertGlobalFaceListWithTwins(array,nbelem,local,ip,full_array,size);
01990       for (int i=0; i<size; i++)
01991         element_array[ip[i]].push_back(local[i]);
01992       delete[] full_array;  
01993       break;
01994     case MED_EN::MED_NODE :
01995       old_collection.getTopology()->convertNodeToGlobal(idomain,list_of_elems,nbelem,array);
01996       _topology->convertGlobalNodeListWithTwins(array,nbelem,local,ip,full_array,size);
01997       for (int i=0; i<size; i++)
01998         element_array[ip[i]].push_back(local[i]);
01999       delete[] full_array;
02000       break;
02001 
02002     }
02003     delete[] ip;
02004     delete[] local;
02005     delete[] array;
02006 
02007     if (support->isOnAllElements()) delete[] list_of_elems;
02008   }
02009 
02010   //creating new groups from the previous list of elements
02011   for (int idomain=0; idomain <_topology->nbDomain(); idomain++)
02012   {
02013     if ( _mesh[idomain]->getNumberOfNodes() < 1 || 
02014          (element_array[idomain].empty() && !_create_empty_groups))
02015     {
02016       new_support[idomain]->removeReference();
02017       new_support[idomain]=0;
02018       continue;
02019     }
02020     MEDMEM::SUPPORT* support= new_support[idomain];
02021     support->setName(name);
02022     support->setMesh(_mesh[idomain]);
02023     support->setDescription(description);
02024     support->setEntity(entity);
02025 
02026     if ( !element_array[idomain].empty() ) /* if() was added for issue 0021576
02027                                               to prevent creation of faces */
02028       {
02029         element_array[idomain].sort();
02030         element_array[idomain].unique();
02031 
02032         if (entity != MED_EN::MED_NODE)
02033           support->fillFromElementList(element_array[idomain]);
02034         else
02035           support->fillFromNodeList(element_array[idomain]);
02036       }
02037   }
02038 }
02039 
02040 void MESHCollection::castField(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber)
02041 {
02042   int type=old_collection.getDriver()->getFieldType(fieldname);
02043   char field_char[80];
02044   strcpy(field_char,fieldname.c_str());
02045 
02046   if (type ==0)
02047     castFields<int>(old_collection, field_char, itnumber, ordernumber);
02048   else
02049     castFields<double>(old_collection, field_char, itnumber, ordernumber);
02050 }
02051 
02052 void MESHCollection::castAllFields(const MESHCollection& initial_collection)
02053 {
02054   vector <string> field_names;
02055   vector <int> iternumber;
02056   vector <int> ordernumber;
02057   vector <int> types;
02058   initial_collection.getDriver()->readFileStruct(field_names,iternumber,ordernumber,types);
02059 
02060   for (unsigned i=0; i<field_names.size(); i++)
02061   {
02062     char field_char[80];
02063     strcpy(field_char,field_names[i].c_str());
02064 
02065     // choosing whether the field is of int or double type
02066     if (types[i] ==0)
02067       castFields<int>(initial_collection, field_char, iternumber[i], ordernumber[i]);
02068     else
02069       castFields<double>(initial_collection, field_char, iternumber[i], ordernumber[i]);
02070   }
02071 }
02072 
02073 void MESHCollection::createNodalConnectivity(const MESHCollection& initial_collection,int idomain, MED_EN::medEntityMesh entity)
02074 {
02075   MESSAGE_MED ("beginning of createNodalConnectivity for entity "<<entity);
02076   int dimension=0;
02077   int nb_elems=0;
02078   MEDMEM::MESHING* mesh_builder = static_cast<MEDMEM::MESHING*>(_mesh[idomain]);
02079 
02080 
02081   //number of elements per type
02082   std::map<MED_EN::medGeometryElement,int> type_numbers;
02083 
02084   //creating arrays for storing global numbers and cell types
02085   switch (entity)
02086   {
02087   case MED_EN::MED_CELL:
02088     dimension=initial_collection.getMeshDimension();
02089     nb_elems=_topology->getCellNumber(idomain);
02090     break;
02091   case MED_EN::MED_EDGE:
02092   case MED_EN::MED_FACE:
02093     dimension=initial_collection.getMeshDimension()-1;
02094     nb_elems=_topology->getFaceNumber(idomain); 
02095     break;
02096   default:
02097     nb_elems=0;
02098     break;
02099   }
02100 
02101   if (nb_elems == 0) return;
02102   SCRUTE_MED(nb_elems);
02103 
02104 
02105   int *list= new int[nb_elems];
02106   MED_EN::medGeometryElement *cell_type_list= new MED_EN::medGeometryElement[nb_elems];
02107 
02108 
02109   //  cout << "Beginning of retrieval "<<endl;
02110   //retrieving global id list
02111   switch (entity)
02112   {
02113   case MED_EN::MED_CELL:
02114     _topology->getCellList(idomain,list);
02115     break;
02116   case MED_EN::MED_EDGE:
02117   case MED_EN::MED_FACE:
02118     _topology->getFaceList(idomain,list);
02119     break;
02120   default:
02121 
02122     break;
02123   }
02124 
02125   //retrieving cell_types
02126   initial_collection.getTypeList(list,nb_elems,entity,cell_type_list);
02127   //  cout <<"end of type retrieval"<<endl;
02128   //vector containing the number of cells per type
02129   type_numbers.clear();
02130   for (int icell=0; icell<nb_elems; icell++)
02131   {
02132     map<MED_EN::medGeometryElement,int>::iterator iter= type_numbers.find(cell_type_list[icell]);
02133     if (iter!=type_numbers.end())
02134       (iter->second)++;
02135     else  
02136       type_numbers[cell_type_list[icell]]=1;
02137 
02138   }
02139   //cout << "Nombre de tetras"<<type_numbers[304]<<endl;
02140   int nb_present_types=type_numbers.size();
02141 
02142   //setting the list of cells for each type
02143   map<MED_EN::medGeometryElement,int> index;
02144 
02145   map<MED_EN::medGeometryElement,int*> type_cell_list;
02146 
02147   MED_EN::MESH_ENTITIES::const_iterator currentEntity;
02148   std::map<MED_EN::medGeometryElement,int>::const_iterator iter;
02149   //currentEntity  = MED_EN::meshEntities.find(entity);
02150   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)  
02151   {
02152     MED_EN::medGeometryElement type = iter->first;
02153     if (!isDimensionOK(type,dimension)) continue;
02154     //if (iter->second==0) continue;
02155     index[type]=0;
02156     type_cell_list[type]=new int[type_numbers[type]];
02157     // cout << "type :"<<type<<" nb:"<<type_numbers[type]<<endl;
02158   }
02159 
02160   for (int icell=0; icell<nb_elems; icell++)
02161   {
02162     type_cell_list[cell_type_list[icell]][index[cell_type_list[icell]]++]=list[icell];
02163   }
02164 
02165   delete[]list;
02166   delete[]cell_type_list;
02167 
02168   //setting the list of present types
02169   int* present_type_numbers=new int[nb_present_types];
02170   MED_EN::medGeometryElement* type_array = new MED_EN::medGeometryElement[nb_present_types];
02171   MESSAGE_MED("Nb de types presents "<<nb_present_types);
02172   int itype=0;
02173   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)  
02174   {
02175     MED_EN::medGeometryElement type = iter->first;
02176     if (!isDimensionOK(type,dimension)) continue;
02177 
02178     type_array[itype]=type;
02179 
02180     present_type_numbers[itype]=type_numbers[type];
02181 
02182     MESSAGE_MED("Nombre d'elements de type "<<type<<" : "<<type_numbers[type]);
02183     itype++;
02184   }
02185 
02186   //retrieving connectivity in global numbering for each type
02187   map<MED_EN::medGeometryElement,int*> type_connectivity;
02188   vector<int> polygon_conn;
02189   vector<int> polygon_conn_index;
02190   vector<int> polyhedron_conn;
02191   vector<int> polyhedron_conn_index;
02192   vector<int> polyhedron_face_index;
02193 
02194   //Treating nodes
02195 
02196 
02197   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)  
02198   {
02199     MED_EN::medGeometryElement type = iter->first;
02200 
02201 
02202     if (!isDimensionOK(type,dimension)) continue;
02203     if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
02204     { 
02205       int nbnode_per_type = (int)type%100;
02206       type_connectivity[type]=new int[type_numbers[type]*nbnode_per_type];
02207       initial_collection.getNodeConnectivity(type_cell_list[type],type_numbers[type],entity,type,type_connectivity[type]);
02208     }
02209     else if (type == MED_EN::MED_POLYGON && dimension==2)
02210     {
02211       initial_collection.getPolygonNodeConnectivity(type_cell_list[type],type_numbers[type],entity,polygon_conn,polygon_conn_index);
02212     }
02213     else if (type == MED_EN::MED_POLYHEDRA && dimension==3)
02214     {
02215       initial_collection.getPolyhedraNodeConnectivity(type_cell_list[type],type_numbers[type],entity,polyhedron_conn,polyhedron_conn_index);
02216     }
02217     delete[] type_cell_list[type];
02218   } 
02219 
02220   //creating node mapping 
02222   if (entity==MED_EN::MED_CELL)
02223   {
02224     _topology->createNodeMapping(type_connectivity,type_numbers,polygon_conn,polygon_conn_index,
02225                                  polyhedron_conn,polyhedron_conn_index,polyhedron_face_index,idomain);
02226   } 
02227 
02228   //converting node global numberings to local numberings
02229   //for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
02230   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)  
02231   {
02232     MED_EN::medGeometryElement type = iter->first;
02233 
02234     if (!isDimensionOK(type, dimension)) continue;
02235     if (type_numbers[type]==0) continue;
02236     if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
02237     { 
02238       int nbnode_per_type = (int)type%100;
02239       _topology->convertToLocal2ndVersion(type_connectivity[type],type_numbers[type]*nbnode_per_type,idomain);
02240     }
02241     else if (type == MED_EN::MED_POLYGON && dimension==2)
02242     {
02243       int nbpoly = type_numbers[type]; 
02244       _topology->convertToLocal2ndVersion(&polygon_conn[0], polygon_conn_index[nbpoly]-1, idomain);  
02245     }
02246     else if (type == MED_EN::MED_POLYHEDRA && dimension==3)
02247     {
02248       int nbpoly = type_numbers[type]; 
02249       _topology->convertToLocal2ndVersion(&polyhedron_conn[0], polyhedron_face_index[polyhedron_conn_index[nbpoly]-1]-1, idomain);  
02250     }
02251 
02252   } 
02253 
02254 
02255   //writing coordinates
02256   if (entity==MED_EN::MED_CELL) 
02257   {
02258     //setting coordinates from initial_collection coordinates
02259     int nbnode=_topology->getNodeNumber(idomain);
02260     MESSAGE_MED("Number of nodes on domain "<< idomain <<" : "<<nbnode);
02261 
02262     double* coordinates=new double[initial_collection.getSpaceDimension()*nbnode];
02263     int* node_list=new int[nbnode];
02264     _topology->getNodeList(idomain,node_list);
02265     initial_collection.getCoordinates(node_list,nbnode,coordinates);
02266     delete[] node_list;
02267 
02268     // redundant specification of number of nodes is required!! MED imperfection, sorry...  
02269 
02270     //TODO : change MEDMEM so that it accepts a direct setting of coordinates
02271     // (in the present version, it is deep-copied)
02272     mesh_builder->setCoordinates(initial_collection.getSpaceDimension(),
02273                                  nbnode, coordinates, initial_collection.getSystem(),
02274                                  MED_EN::MED_FULL_INTERLACE);
02275     delete [] coordinates;
02276   }
02277 
02278   int nb_plain_types=0;
02279   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++) 
02280   { 
02281     MED_EN::medGeometryElement type = iter->first;
02282 
02283     if (!isDimensionOK(type, dimension)) continue;
02284     if (type_numbers[type]==0) continue;
02285     nb_plain_types++;
02286   }
02287   mesh_builder->setNumberOfTypes(nb_plain_types,entity);
02288   mesh_builder->setTypes(type_array,entity);
02289   mesh_builder->setNumberOfElements(present_type_numbers,entity);
02290 
02291   delete[]present_type_numbers;
02292   delete[]type_array;
02293   //setting node connectivities
02294   for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)  
02295   {
02296     MED_EN::medGeometryElement type = iter->first;
02297 
02298     if (!isDimensionOK(type,dimension)) continue;
02299     if (type_numbers[type]==0) continue;
02300 
02301     if (type != MED_EN::MED_POLYHEDRA && type != MED_EN::MED_POLYGON)
02302     {
02303       mesh_builder->setConnectivity(entity,type,type_connectivity[type]);
02304       delete[] type_connectivity[type];
02305     }
02306     else if (type == MED_EN::MED_POLYGON && dimension ==2)
02307     {
02308       mesh_builder->setConnectivity(entity,type,&polygon_conn[0],&polygon_conn_index[0]);
02309     }
02310     else if (type == MED_EN::MED_POLYHEDRA && dimension ==3)
02311     {
02312       mesh_builder->setConnectivity(entity,type,&polyhedron_conn[0],&polyhedron_conn_index[0]);
02313     }
02314   }
02315   MESSAGE_MED("end of createNodalConnectivity");
02316 }
02317 
02318 
02325 void MESHCollection::getFaces(int idomain, 
02326                               map<MED_EN::medGeometryElement, vector<MEDSPLITTER_FaceModel*> >& face_map)                     
02327 {
02328   MED_EN::medEntityMesh constituent_entity = getSubEntity();
02329   const medGeometryElement* types;
02330   try
02331   {
02332     types = _mesh[idomain]->getTypes(constituent_entity);
02333     if ( !types ) return;
02334   }
02335   catch(MEDEXCEPTION&){ return;}
02336 
02337   int nbtypes  = _mesh[idomain]->getNumberOfTypes(constituent_entity);
02338   const int* global_numbering= _mesh[idomain]->getGlobalNumberingIndex(constituent_entity);
02339   int* conn = const_cast<int*> (_mesh[idomain]->getConnectivity(MED_EN::MED_NODAL,constituent_entity, MED_EN::MED_ALL_ELEMENTS));
02340   for (int itype=0; itype<nbtypes; itype++)
02341   {
02342     for (int iface=global_numbering[itype]; iface<global_numbering[itype+1]; iface++)
02343     {
02344       MEDSPLITTER_FaceModel* face_model = new MEDSPLITTER_FaceModel();
02345       MED_EN::medGeometryElement type =  types[itype];
02346       face_model->setType(type);
02347       int nbnodes = type%100;
02348       face_model->setNbNodes(nbnodes);
02349       face_model->setGlobal(_topology->convertFaceToGlobal(idomain,iface));
02350       for (int i=0; i<nbnodes; i++)
02351       {
02352         (*face_model)[i]=*conn++;
02353       }
02354       face_map[type].push_back(face_model);
02355     }
02356   }
02357 }
02358 
02367 MEDSPLITTER_FaceModel* MESHCollection::getCommonFace(int ip1,int ilocal1,int ip2,int ilocal2,int face_index)
02368 {
02369   MED_EN::medGeometryElement type1 = _mesh[ip1]->getElementType(MED_EN::MED_CELL,ilocal1);
02370   MEDMEM::CELLMODEL celltype1 (type1);
02371 
02372   const int* conn_index1 =  _mesh[ip1]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
02373   const int* conn1 = _mesh[ip1]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
02374 
02375   // MED_EN::medGeometryElement type2 = _mesh[ip2]->getElementType(MED_EN::MED_CELL,ilocal2);
02376   //MEDMEM::CELLTYPE celltype2 (type2);
02377   const int* conn_index2 =  _mesh[ip2]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
02378   const int* conn2 = _mesh[ip2]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
02379 
02380   vector<int> nodes1, nodes1_local;
02381   vector<int> nodes2;
02382   for (int i=  conn_index1[ilocal1-1]; i<conn_index1[ilocal1]; i++)
02383   {
02384     nodes1.push_back(_topology->convertNodeToGlobal(ip1,*(conn1+i-1)));
02385     nodes1_local.push_back( conn1[i-1] );
02386   }
02387   for (int i=  conn_index2[ilocal2-1]; i<conn_index2[ilocal2]; i++)
02388     nodes2.push_back(_topology->convertNodeToGlobal(ip2,*(conn2+i-1)));
02389 
02390   return MEDSPLITTER_FaceModel::getCommonFace( &nodes1[0], &nodes1_local[0], celltype1,
02391                                                &nodes2[0], nodes2.size(),  face_index);
02392 }
02393 
02394 //================================================================================
02404 //================================================================================
02405 
02406 MEDSPLITTER_FaceModel*
02407 MEDSPLITTER_FaceModel::getCommonFace(const int*               nodes1,
02408                                      const int*               nodes1_local,
02409                                      const MEDMEM::CELLMODEL& celltype1,
02410                                      const int*               nodes2,
02411                                      int                      nb_nodes2,
02412                                      int                      global_id)
02413 {
02414   int nbfaces= celltype1.getNumberOfConstituents(1);
02415   int ** faces = celltype1.getConstituents(1);
02416   MED_EN::medGeometryElement* types = celltype1.getConstituentsType(1);
02417   int iface=0;
02418   int dimension=celltype1.getDimension();
02419 
02420   while (iface<nbfaces)
02421   {
02422     //SCRUTE_MED (iface);
02423     int nbnodes= types[iface]%100;
02424     const int* nodes = celltype1.getNodesConstituent(1,iface+1);
02425     int common_nodes=0;
02426     for (int i=0; i<nbnodes;i++)
02427     {
02428       for (int i2=0; i2<nb_nodes2; i2++)
02429       {
02430         if (nodes1[nodes[i]-1]==nodes2[i2]) common_nodes++;
02431       }
02432     }
02433     if (common_nodes>=dimension) break;
02434     iface++;
02435   }
02436 
02437   if (iface==nbfaces)
02438     throw MEDEXCEPTION("MEDSPLITTER::getCommonFace - No common face found !");
02439 
02440   MEDSPLITTER_FaceModel* face_model = new MEDSPLITTER_FaceModel;
02441   face_model->setType(types[iface]);
02442   int nbnodes = types[iface]%100;
02443   face_model->setNbNodes(nbnodes);
02444   face_model->setGlobal(global_id); 
02445   for (int i=0; i<nbnodes; i++)
02446     (*face_model)[i]=nodes1_local[faces[iface][i]-1];
02447 
02448   return face_model;
02449 }