1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "MEDMEM_ConnectZone.hxx"
21 #include "MEDMEM_DriversDef.hxx"
22 #include "MEDMEM_Mesh.hxx"
23 #include "MEDMEM_Meshing.hxx"
24 #include "MEDMEM_GaussLocalization.hxx"
25 #include "MEDMEM_Field.hxx"
26 #include "MEDMEM_CellModel.hxx"
27 #include "MEDMEM_Group.hxx"
28 #include "MEDMEM_MeshFuse.hxx"
30 #include "MEDMEM_Exception.hxx"
32 #include "MEDSPLITTER_utils.hxx"
34 #include "MEDSPLITTER_Graph.hxx"
36 #include "MEDSPLITTER_Topology.hxx"
37 #include "MEDSPLITTER_ParallelTopology.hxx"
38 #include "MEDSPLITTER_SequentialTopology.hxx"
39 #include "MEDSPLITTER_ParaDomainSelector.hxx"
40 #include "MEDSPLITTER_MeshSendReceive.hxx"
41 #include "MEDSPLITTER_JointExchangeData.hxx"
43 #include "MEDSPLITTER_MESHCollection.hxx"
44 #include "MEDSPLITTER_MESHCollectionDriver.hxx"
45 #include "MEDSPLITTER_MESHCollectionMedXMLDriver.hxx"
46 #include "MEDSPLITTER_MESHCollectionMedAsciiDriver.hxx"
48 #include "MEDSPLITTER_UserGraph.hxx"
50 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
51 #include "MEDSPLITTER_METISGraph.hxx"
53 #ifdef MED_ENABLE_SCOTCH
54 #include "MEDSPLITTER_SCOTCHGraph.hxx"
57 #include "InterpKernelHashMap.hxx"
66 using namespace MEDSPLITTER;
68 using namespace INTERP_KERNEL;
71 #include "MEDSPLITTER_MESHCollection.H"
73 MESHCollection::MESHCollection()
75 _owns_topology(false),
77 _domain_selector( 0 ),
78 _i_non_empty_mesh(-1),
79 _driver_type(MEDSPLITTER::MedXML),
80 _subdomain_boundary_creates(false),
81 _family_splitting(false),
82 _create_empty_groups(false)
88 //================================================================================
90 * \brief Creates a new domain mesh
92 //================================================================================
94 MEDMEM::MESH* newMesh(const std::string& name, int dim, int space, MEDMEM::MESH* meshToDelete=0)
97 MEDMEM::MESHING* mesh = new MEDMEM::MeshFuse;
98 mesh->setName( name );
99 //mesh->setMeshDimension ( dim );
100 //mesh->setSpaceDimension( space );
105 /*!constructor creating a new mesh collection (mesh series + topology)
106 *from an old collection and a new topology
108 * On output, the constructor has built the meshes corresponding to the new mesh collection.
109 * The new topology has been updated so that face and node mappings are included.
110 * The families have been cast to their projections in the new topology.
112 * \param initial_collection collection from which the data (coordinates, connectivity) are taken
113 * \param topology topology containing the cell mappings
116 MESHCollection::MESHCollection(const MESHCollection& initial_collection, Topology* topology, bool family_splitting, bool create_empty_groups)
117 : _topology(topology),
118 _owns_topology(false),
119 _cell_graph(topology->getGraph()),
121 _domain_selector( initial_collection._domain_selector ),
122 _i_non_empty_mesh(-1),
123 _name(initial_collection._name),
124 _driver_type(MEDSPLITTER::MedXML),
125 _subdomain_boundary_creates(false),
126 _family_splitting(family_splitting),
127 _create_empty_groups(create_empty_groups)
129 string mesh_name = initial_collection.getName();
130 _mesh.resize(_topology->nbDomain());
132 int space_dim = initial_collection.getSpaceDimension();
133 int mesh_dim = initial_collection.getMeshDimension();
135 space_dim = mesh_dim = initial_collection._driver->readMeshDimension();
137 for (int idomain=0; idomain < _topology->nbDomain(); idomain++)
139 //creating the new mesh
140 _mesh[idomain]= newMesh( MEDMEM::STRING(mesh_name)<<"_"<<idomain+1, mesh_dim, space_dim );
142 createNodalConnectivity(initial_collection,idomain, MED_EN::MED_CELL);
144 if ( _mesh[idomain]->getNumberOfNodes() > 0 )
145 _i_non_empty_mesh = idomain;
148 _topology->createFaceMapping(initial_collection, *this );
149 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
154 createNodalConnectivity(initial_collection,idomain, MED_EN::MED_FACE);
157 createNodalConnectivity(initial_collection,idomain, MED_EN::MED_EDGE);
160 if ( !isParallelMode() || _domain_selector->isMyDomain( idomain ))
161 cerr<<"MEDSPLITTER : Mesh dimension must be 2 or 3"<<endl;
165 castFamilies(initial_collection);
167 // Exchange domain parts
169 if ( isParallelMode() )
171 _domain_selector->setNbDomains( _topology->nbDomain() );
173 vector< MeshSendReceive > mesh_sender( _topology->nbDomain() );
174 list<int> domainsToClear; // sent domains
176 // first, send domains
177 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
179 // get node numbers global over all procs
180 vector<int> glob_nodes_all_proc( _topology->getNodeNumber( idomain )); // to fill in
181 vector<int> glob_cells_all_proc( _topology->getCellNumber( idomain ));
182 vector<int> glob_faces_all_proc( _topology->getFaceNumber( idomain ));
183 if ( !glob_cells_all_proc.empty() )
185 // get ids global on this proc
186 _topology->getNodeList( idomain, & glob_nodes_all_proc[0] );
187 _topology->getCellList( idomain, & glob_cells_all_proc[0] );
188 _topology->getFaceList( idomain, & glob_faces_all_proc[0] );
189 // convert cell ids to ids global over all procs
190 int cell_shift = _domain_selector->getProcShift();
191 for ( int i = 0; i < glob_cells_all_proc.size(); ++i )
192 glob_cells_all_proc[i] += cell_shift;
194 if ( _domain_selector->isMyDomain( idomain ))
196 // prepare to receiving other parts of the domain
197 ((MEDMEM::MeshFuse*) _mesh[idomain])->setNodeNumbers( glob_nodes_all_proc );
198 _topology->getFusedCellNumbers( idomain ) = glob_cells_all_proc;
199 _topology->getFusedFaceNumbers( idomain ) = glob_faces_all_proc;
204 int target_proc = _domain_selector->getProccessorID( idomain );
205 mesh_sender[ idomain ].send( target_proc, idomain, _mesh[idomain],
208 glob_nodes_all_proc );
209 if ( !glob_nodes_all_proc.empty() )
210 domainsToClear.push_back( idomain );
212 // clear just sent domain meshes
213 for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); )
215 if (( isSent = mesh_sender[ *dom ].isSent() ))
216 _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim, _mesh[*dom]);
217 dom = isSent ? domainsToClear.erase( dom ) : ++dom;
221 // then, receive domains
222 MeshSendReceive mesh_receiver;
223 int this_proc = _domain_selector->rank();
224 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
226 if ( _domain_selector->isMyDomain( idomain ))
228 for (int iproc = 0; iproc < _domain_selector->nbProcs(); ++iproc)
230 if ( iproc == this_proc ) continue;
231 vector<int> nodes_other_proc, cells_other_proc, faces_other_proc;
232 mesh_receiver.recv( iproc, idomain, cells_other_proc, faces_other_proc,nodes_other_proc);
233 if ( MEDMEM::MESH* received_mesh = mesh_receiver.getMesh() )
235 // unite meshes and global node numbers stored in MeshFuse
236 MEDMEM::MeshFuse* fuse = (MEDMEM::MeshFuse*) _mesh[idomain];
237 fuse->concatenate( received_mesh, nodes_other_proc );
238 delete received_mesh;
240 // unite global element numbers
241 fuse->append( MED_EN::MED_CELL,
242 _topology->getFusedCellNumbers( idomain ), cells_other_proc );
244 fuse->append( mesh_dim==3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE,
245 _topology->getFusedFaceNumbers( idomain ), faces_other_proc );
247 if ( _mesh[idomain]->getNumberOfNodes() > 0 )
248 _i_non_empty_mesh = idomain;
252 // clear just sent domain meshes
253 for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); )
255 if (( isSent = mesh_sender[ *dom ].isSent() ))
256 _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim,_mesh[*dom]);
257 dom = isSent ? domainsToClear.erase( dom ) : ++dom;
260 // clear sent domain meshes
262 for ( list<int>::iterator dom = domainsToClear.begin(); dom != domainsToClear.end(); ++dom)
263 _mesh[*dom] = newMesh( _mesh[*dom]->getName(), mesh_dim, space_dim,_mesh[*dom]);
265 _topology->recreateMappingAfterFusion( getMesh() );
267 if ( _i_non_empty_mesh < 0 ) // non of domains resides on this proc,
268 _i_non_empty_mesh = 0; // in this case we need only dimension that is set to all meshes
272 /*! constructing the MESH collection from a distributed file
274 * \param filename name of the master file containing the list of all the MED files
276 MESHCollection::MESHCollection(const string& filename)
278 _owns_topology(true),
280 _domain_selector( 0 ),
281 _i_non_empty_mesh(-1),
282 _driver_type(MEDSPLITTER::Undefined),
283 _subdomain_boundary_creates(false),
284 _family_splitting(false),
285 _create_empty_groups(false)
287 char filenamechar[256];
288 strcpy(filenamechar,filename.c_str());
291 _driver=new MESHCollectionMedXMLDriver(this);
292 _driver->read (filenamechar);
293 _driver_type = MedXML;
296 catch(MEDEXCEPTION&){
300 _driver=new MESHCollectionMedAsciiDriver(this);
301 _driver->read (filenamechar);
302 _driver_type=MedAscii;
307 throw MEDEXCEPTION("file does not comply with any recognized format");
310 for ( int idomain = 0; idomain < _mesh.size(); ++idomain )
311 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
312 _i_non_empty_mesh = idomain;
315 /*! Constructing the MESH collection from selected domains of a distributed file
317 * \param filename - name of the master file containing the list of all the MED files
318 * \param domainSelector - selector of domains to load
320 MESHCollection::MESHCollection(const string& filename, ParaDomainSelector& domainSelector)
322 _owns_topology(true),
324 _domain_selector( domainSelector.nbProcs() > 1 ? & domainSelector : 0 ),
325 _i_non_empty_mesh(-1),
326 _driver_type(MEDSPLITTER::Undefined),
327 _subdomain_boundary_creates(false),
328 _family_splitting(false),
329 _create_empty_groups(false)
333 _driver=new MESHCollectionMedXMLDriver(this);
334 _driver->read ( (char*)filename.c_str(), _domain_selector );
335 _driver_type = MedXML;
342 _driver=new MESHCollectionMedAsciiDriver(this);
343 _driver->read ( (char*)filename.c_str(), _domain_selector );
344 _driver_type=MedAscii;
349 throw MEDEXCEPTION("file does not comply with any recognized format");
352 if ( isParallelMode() )
353 // to know nb of cells on each proc to compute global cell ids from locally global
354 _domain_selector->gatherNbOf( MED_EN::MED_CELL, getMesh() );
356 // find non-empty domain mesh
357 for ( int idomain = 0; idomain < _mesh.size(); ++idomain )
358 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
359 _i_non_empty_mesh = idomain;
362 /*! constructing the MESH collection from a sequential MED-file
364 * \param filename MED file
365 * \param meshname name of the mesh that is to be read
367 MESHCollection::MESHCollection(const string& filename, const string& meshname)
369 _owns_topology(true),
371 _domain_selector( 0 ),
372 _i_non_empty_mesh(-1),
374 _driver_type(MEDSPLITTER::MedXML),
375 _subdomain_boundary_creates(false),
376 _family_splitting(false),
377 _create_empty_groups(false)
379 char filenamechar[256];
380 char meshnamechar[256];
381 strcpy(filenamechar,filename.c_str());
382 strcpy(meshnamechar,meshname.c_str());
383 try // avoid memory leak in case of inexistent filename
385 retrieveDriver()->readSeq (filenamechar,meshnamechar);
387 catch ( MED_EXCEPTION& e )
389 if ( _driver ) delete _driver; _driver=0;
392 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
393 _i_non_empty_mesh = 0;
396 MESHCollection::~MESHCollection()
398 for (unsigned i=0; i<_mesh.size();i++)
399 if (_mesh[i]!=0) {/*delete*/ _mesh[i]->removeReference(); }
400 for (unsigned i=0; i<_connect_zones.size();i++)
401 if (_connect_zones[i]!=0) {delete _connect_zones[i];}
402 if (_driver !=0) {delete _driver; _driver=0;}
403 if (_topology!=0 && _owns_topology) {delete _topology; _topology=0;}
406 /*!gets the connectivity for a certain type
408 * The output array type_connectivity should have been allocated
409 * at dimension nbnode_per_type* nb_cells before the call
411 * \param cell_list list of elements (global cell numbers) for which the connectivity is required
412 * \param nb_cells number of elements
413 * \param entity type of entity for which the nodal connectivity is required
414 * \param type type of the elements for which the connectivity is required
415 * \param type_connectivity on output contains the connectivity of all the elements of the list
418 void MESHCollection::getNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
419 MED_EN::medGeometryElement type, int* type_connectivity) const
421 int *local=new int[nb_cells];
422 int *ip=new int[nb_cells];
425 case MED_EN::MED_CELL:
426 _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
428 case MED_EN::MED_FACE:
429 case MED_EN::MED_EDGE:
430 _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
435 // int nbnode_per_type=(int)type%100;
436 // vector<int> number_of_types_array(_topology->nbDomain(),0);
437 // for (int i=0; i<_topology->nbDomain(); i++)
438 // number_of_types_array[i]=_mesh[i]->getNumberOfTypes(entity);
440 //defining a connectivity table for different domains
441 vector <const int*> conn_ip(_topology->nbDomain());
442 vector <const int*> conn_index_ip(_topology->nbDomain());
445 vector< map <MED_EN::medGeometryElement, int> > offset;
446 // offset.resize(_topology->nbDomain());
448 for (int i=0; i<_topology->nbDomain();i++)
450 if ( !_mesh[i] ) continue;
451 int nb_elem = _mesh[i]->getNumberOfElements(entity,type);
454 conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,entity,MED_EN::MED_ALL_ELEMENTS);
455 conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,entity);
456 // global_index= _mesh[i]->getGlobalNumberingIndex(entity);
463 // int number_of_types = number_of_types_array[i];
464 // const MEDMEM::CELLMODEL* types = _mesh[ip[icell]]->getCellsTypes(entity);
465 // for (int itype=0; itype<number_of_types; itype++)
466 // offset[i][types[itype].getType()]=global_index[itype]-1;
469 int* type_connectivity_ptr=type_connectivity;
470 for (int icell=0; icell<nb_cells; icell++)
472 // int type_offset = offset[ip[icell]][type];
473 const int* conn=conn_ip[ip[icell]];
474 const int* conn_index=conn_index_ip[ip[icell]];
475 for (int inode=conn_index[local[icell]-1]; inode<conn_index[local[icell]]; inode++)
477 *type_connectivity_ptr=
478 _topology->convertNodeToGlobal(ip[icell],conn[inode-1]);
479 type_connectivity_ptr++;
487 /*!gets the connectivity for MED_POLYGON type
489 * \param cell_list list of elements (global cell numbers) for which the connectivity is required
490 * \param nb_cells number of elements
491 * \param entity type of entity for which the nodal connectivity is required
492 * \param type_connectivity on output contains the connectivity of all the elements of the list
493 * \param connectivity_index on output contains the connectivity index for all polygons
496 void MESHCollection::getPolygonNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
497 vector<int>& type_connectivity, vector<int>& connectivity_index) const
500 int *local=new int[nb_cells];
501 int *ip=new int[nb_cells];
504 case MED_EN::MED_CELL:
505 _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
507 case MED_EN::MED_FACE:
508 case MED_EN::MED_EDGE:
509 _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
514 //defining a connectivity table for different domains
515 vector <const int*> conn_ip(_topology->nbDomain());
516 vector <const int*> conn_index_ip(_topology->nbDomain());
517 vector <const int* > conn_face_index(_topology->nbDomain());
518 vector<int> nb_plain_elems(_topology->nbDomain());
520 vector< map <MED_EN::medGeometryElement, int> > offset;
522 for (int i=0; i<_topology->nbDomain();i++)
524 int nb_elem = _mesh[i]->getNumberOfElements(entity,MED_EN::MED_POLYGON);
527 conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,entity,MED_EN::MED_ALL_ELEMENTS);
528 conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,entity);
537 connectivity_index.resize(nb_cells+1);
538 connectivity_index[0]=1;
539 for (int icell=0; icell<nb_cells; icell++)
541 //int nb_plain= nb_plain_elems[ip[icell]];
542 const int* conn=conn_ip[ip[icell]];
543 const int* conn_index=conn_index_ip[ip[icell]];
544 for (int inode=conn_index[local[icell]-1/*-nb_plain*/]; inode<conn_index[local[icell]/*-nb_plain*/]; inode++)
546 type_connectivity.push_back(
547 _topology->convertNodeToGlobal(ip[icell],conn[inode-1]));
549 connectivity_index[icell+1]=connectivity_index[icell]
550 -conn_index[local[icell]-1/*-nb_plain*/]+conn_index[local[icell]/*-nb_plain*/];
558 /*!gets the connectivity for MED_POLYHEDRA type
560 * \param cell_list list of elements (global cell numbers) for which the connectivity is required
561 * \param nb_cells number of elements
562 * \param entity type of entity for which the nodal connectivity is required
563 * \param type_connectivity on output contains the connectivity of all the elements of the list
564 * \param connectivity_index on output contains the connectivity index for all polygons
567 void MESHCollection::getPolyhedraNodeConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
568 vector<int>& type_connectivity, vector<int>& connectivity_index/*, vector<int>& face_connectivity_index*/) const
571 int *local=new int[nb_cells];
572 int *ip=new int[nb_cells];
575 case MED_EN::MED_CELL:
576 _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
578 case MED_EN::MED_FACE:
579 case MED_EN::MED_EDGE:
580 _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
585 //defining a connectivity table for different domains
586 vector <const int*> conn_ip(_topology->nbDomain());
587 vector <const int*> conn_index_ip(_topology->nbDomain());
588 vector<int> nb_plain_elems(_topology->nbDomain());
590 vector< map <MED_EN::medGeometryElement, int> > offset;
592 for (int i=0; i<_topology->nbDomain();i++)
594 nb_plain_elems[i] = _mesh[i]->getNumberOfElements(entity, MED_EN::MED_ALL_ELEMENTS);
595 int nb_elem = _mesh[i]->getNumberOfElements(entity,MED_EN::MED_POLYHEDRA);
598 conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
599 conn_index_ip[i] = _mesh[i]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
608 connectivity_index.resize(nb_cells+1);
609 connectivity_index[0]=1;
610 for (int icell=0; icell<nb_cells; icell++)
612 const int* conn=conn_ip[ip[icell]];
613 const int* conn_index=conn_index_ip[ip[icell]];
614 connectivity_index[icell+1]=connectivity_index[icell]+
615 conn_index[local[icell]]-conn_index[local[icell]-1];
617 for (int inode=conn_index[local[icell]-1]; inode<conn_index[local[icell]]; inode++)
619 if ( conn[inode-1] == -1 )
620 type_connectivity.push_back( -1 );
622 type_connectivity.push_back(_topology->convertNodeToGlobal(ip[icell],conn[inode-1]));
631 /*! constructing the MESH collection from a file
633 * The method creates as many MED-files as there are domains in the
634 * collection. It also creates a master file that lists all the MED files.
635 * The MED files created in ths manner contain joints that describe the
636 * connectivity between subdomains.
638 * \param filename name of the master file that will contain the list of the MED files
641 void MESHCollection::write(const string& filename)
643 //building the connect zones necessary for writing joints
644 cout<<"Building Connect Zones"<<endl;
645 if (_topology->nbDomain()>1)
647 cout <<"End of connect zones building"<<endl;
648 //suppresses link with driver so that it can be changed for writing
649 if (_driver!=0)delete _driver;
652 char filenamechar[256];
653 strcpy(filenamechar,filename.c_str());
654 retrieveDriver()->write (filenamechar, _domain_selector);
657 /*! creates or gets the link to the collection driver
659 MESHCollectionDriver* MESHCollection::retrieveDriver()
666 _driver=new MESHCollectionMedXMLDriver(this);
669 _driver=new MESHCollectionMedAsciiDriver(this);
672 throw MEDEXCEPTION("Unrecognized driver");
680 /*! gets an existing driver
683 MESHCollectionDriver* MESHCollection::getDriver() const
689 /*! gets the list of types for global numbers cell_list
691 * \param cell_list list of global numbers
692 * \param entity entity type
693 * \param type_list on output, list of types for the cells given in cell_list
695 void MESHCollection::getTypeList(int* cell_list,int nb_cells,
696 MED_EN::medEntityMesh entity,
697 MED_EN::medGeometryElement* type_list) const
699 MESSAGE_MED (" Beginning of getTypeList with entity "<<entity);
700 int *local=new int[nb_cells];
701 int *ip=new int[nb_cells];
704 case MED_EN::MED_CELL:
705 _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
707 case MED_EN::MED_FACE:
708 case MED_EN::MED_EDGE:
709 _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
713 for (int icell=0; icell<nb_cells; icell++)
715 type_list[icell]=_mesh[ip[icell]]->getElementType(entity,local[icell]);
719 MESSAGE_MED("end of getTypeList");
724 /*!gets the descending connectivity for a certain type
726 * The output array type_connectivity should have been allocated
727 * at dimension nbnode_per_type* nb_cells before the call
729 * \param cell_list list of elements (global cell numbers) for which the connectivity is required
730 * \param nb_cells number of elements
731 * \param entity type of entity for which the nodal connectivity is required
732 * \param type type of the elements for which the connectivity is required
733 * \param type_connectivity on output contains the connectivity of all the elements of the list
736 void MESHCollection::getFaceConnectivity(const int* cell_list,int nb_cells,MED_EN::medEntityMesh entity,
737 MED_EN::medGeometryElement type, int* type_connectivity) const
739 int *local=new int[nb_cells];
740 int *ip=new int[nb_cells];
743 case MED_EN::MED_CELL:
744 _topology->convertGlobalCellList(cell_list,nb_cells,local,ip);
746 case MED_EN::MED_FACE:
747 case MED_EN::MED_EDGE:
748 _topology->convertGlobalFaceList(cell_list,nb_cells,local,ip);
766 vector<int> number_of_types_array(_topology->nbDomain(),0);
767 for (int i=0; i<_topology->nbDomain(); i++)
768 number_of_types_array[i]=_mesh[i]->getNumberOfTypes(entity);
770 //defining a connectivity table for different domains
771 vector <const int*> conn_ip(_topology->nbDomain());
772 for (int i=0; i<_topology->nbDomain();i++)
774 int nb_elem = _mesh[i]->getNumberOfElements(entity,type);
776 conn_ip[i]=_mesh[i]->getConnectivity(MED_EN::MED_DESCENDING,entity,type);
781 for (int icell=0; icell<nb_cells; icell++)
783 int number_of_types = number_of_types_array[ip[icell]];
784 const MEDMEM::CELLMODEL* types = _mesh[ip[icell]]->getCellsTypes(entity);
786 for (int itype=0; itype< number_of_types; itype++)
788 if (types[itype].getType() < type)
789 type_offset += _mesh[ip[icell]]->getNumberOfElements(entity,types[itype].getType());
791 const int* conn=conn_ip[ip[icell]];
792 for (int iface=0; iface<nbface_per_type; iface++)
794 type_connectivity[icell*nbface_per_type+iface] = _topology->convertFaceToGlobal
795 (ip[icell], abs(conn[(local[icell] - type_offset - 1) * nbface_per_type + iface]));
803 /*! gets the list of coordinates for a given list of global node numbers
805 * The vector containing the coordinates on output should
806 * have been allocated at a dimension _space_dimension * nb_nodes
809 * \param node_list list of global node numbers
810 * \param nb_nodes number of nodes in the list
811 * \param coordinates on output, contains the coordinates
814 void MESHCollection::getCoordinates(int* node_list,int nb_nodes, double* coordinates) const
816 int* local=new int[nb_nodes];
817 int* ip=new int[nb_nodes];
818 int space_dimension= getSpaceDimension();
819 _topology->convertGlobalNodeList(node_list,nb_nodes,local,ip);
820 for (int i=0; i< nb_nodes; i++)
822 const double* coord=_mesh[ip[i]]->getCoordinates(MED_EN::MED_FULL_INTERLACE);
823 for (int icoord=0; icoord<space_dimension; icoord++)
824 coordinates[i*space_dimension+icoord]=coord[(local[i]-1)*space_dimension+icoord];
829 /*! returns constituent entity*/
830 MED_EN::medEntityMesh MESHCollection::getSubEntity() const
832 return getMeshDimension()==3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE;
835 /*! retrieves the space dimension*/
836 int MESHCollection::getSpaceDimension() const
838 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getSpaceDimension();
840 /*! retrieves the mesh dimension*/
841 int MESHCollection::getMeshDimension() const
843 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
846 /*! retrieves the type of coordinates system*/
847 string MESHCollection::getSystem() const
849 return _i_non_empty_mesh < 0 ? "" : _mesh[_i_non_empty_mesh]->getCoordinatesSystem();
852 /*!retrieves the name of the mesh*/
853 string MESHCollection::getMeshName() const
855 return _i_non_empty_mesh < 0 ? (_mesh[0] ? _mesh[0]->getName() : "") : _mesh[_i_non_empty_mesh]->getName();
858 vector<MEDMEM::MESH*>& MESHCollection::getMesh()
863 MEDMEM::MESH* MESHCollection::getMesh(int idomain) const
865 return _mesh[idomain];
868 vector<MEDMEM::CONNECTZONE*>& MESHCollection::getCZ()
870 return _connect_zones;
873 Topology* MESHCollection::getTopology() const
878 void MESHCollection::setTopology(Topology* topo)
882 throw MED_EXCEPTION(STRING("Erreur : topology is already set"));
888 void MESHCollection::setIndivisibleGroup(const string& name)
890 _indivisible_regions.push_back(name);
894 /*! Browses the domains and the regions that have
895 * been marked as indivisible in order to create a vector
896 * the dimlension of which is the total number of cells, and
897 * that contains 0 if the cell belongs to no indivisible group
898 * and that contains an integer corresponding to the group otherwise.
900 * \param indivisible_tag on input is an int* allocated as int[nbcells]
901 * on output contains the tags
905 void MESHCollection::treatIndivisibleRegions(int* indivisible_tag)
907 //tag 0 is positioned on all the cells that are not affected by these tags
908 for (int i=0; i<_topology->nbCells(); i++)
909 indivisible_tag[i]=0;
911 //treating cell groups
912 for (int idomain=0; idomain<_topology->nbDomain();idomain++)
913 for (int igroup=0; igroup<_mesh[idomain]->getNumberOfGroups(MED_EN::MED_CELL); igroup++)
914 for (unsigned i=0; i<_indivisible_regions.size(); i++)
916 const MEDMEM::GROUP* group = _mesh[idomain]->getGroup(MED_EN::MED_CELL,igroup+1);
917 string groupname = group->getName();
918 if (trim(groupname)==trim(_indivisible_regions[i]))
920 int nbcells=group->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
921 const int* numbers=group->getNumber(MED_EN::MED_ALL_ELEMENTS);
922 int* global=new int[nbcells];
923 _topology->convertCellToGlobal(idomain,numbers,nbcells,global);
924 for (int icell=0; icell<nbcells; icell++)
925 indivisible_tag[global[icell]-1]=i+1;
931 //================================================================================
933 * \brief Create cell->node and node->cell connectivities for domains
935 //================================================================================
937 template<class TID2CONN>
938 void MESHCollection::fillGlobalConnectivity(TID2CONN & node2cell, TID2CONN & cell2node )
940 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
942 if ( !_mesh[idomain] ) continue;
943 // MED_EN::medGeometryElement* type_array;
944 int nb_cells= _topology->nbCells(idomain);
945 int* cell_list = new int[nb_cells];
947 //retrieving global id list
948 _topology->getCellList(idomain, cell_list);
950 int nb_plain_cells = _mesh[idomain]->getNumberOfElements(MED_EN::MED_CELL,
951 MED_EN::MED_ALL_ELEMENTS);
952 if (nb_plain_cells >0)
954 const int* conn_index = _mesh[idomain]->getConnectivityIndex(MED_EN::MED_NODAL,
957 const int* conn = _mesh[idomain]->getConnectivity(MED_EN::MED_NODAL,
959 MED_EN::MED_ALL_ELEMENTS);
960 int nbnodes = conn_index[nb_plain_cells]-1;
961 int* global_nodes =new int [nbnodes];
962 _topology->convertNodeToGlobal(idomain, conn, nbnodes, global_nodes);
963 for (int icell=0; icell< nb_plain_cells; icell++)
965 for (int inode=conn_index[icell]; inode < conn_index[icell+1]; inode++)
967 int node_global_id = global_nodes[inode-1];
968 if ( node_global_id > 0 )
970 int cell_global_id = cell_list[icell];
971 cell2node [cell_global_id].push_back(node_global_id);
972 node2cell [node_global_id].push_back(cell_global_id);
976 delete[] global_nodes;
983 /*! Method creating the cell graph
985 * \param array returns the pointer to the structure that contains the graph
986 * \param edgeweight returns the pointer to the table that contains the edgeweights
987 * (only used if indivisible regions are required)
990 void MESHCollection::buildCellGraph(MEDMEM::MEDSKYLINEARRAY* & array,int *& edgeweights )
995 for (int i=0; i<_topology->nbDomain(); i++)
997 cell_number+=_topology->getCellNumber(i);
998 node_number+=_topology->getNodeNumber(i);
1000 //list of cells for a given node
1001 //vector< vector<int> > node2cell(node_number);
1002 map< int, vector<int> > node2cell;
1004 //list of nodes for a given cell
1005 //vector< vector <int> > cell2node(cell_number);
1006 map< int, vector <int> > cell2node;
1008 // map<MED_EN::medGeometryElement,int*> type_cell_list;
1010 //tagging for the indivisible regions
1011 int* indivisible_tag=0;
1012 bool has_indivisible_regions=false;
1013 if (!_indivisible_regions.empty())
1015 has_indivisible_regions=true;
1016 indivisible_tag=new int[_topology->nbCells()];
1017 treatIndivisibleRegions(indivisible_tag);
1020 fillGlobalConnectivity(node2cell, cell2node );
1022 cout << "beginning of skyline creation"<<endl;
1023 //creating the MEDMEMSKYLINEARRAY containing the graph
1025 int* size = new int[_topology->nbCells()];
1026 int** temp=new int*[_topology->nbCells()];
1027 int** temp_edgeweight=0;
1028 if (has_indivisible_regions)
1029 temp_edgeweight=new int*[_topology->nbCells()];
1031 int cell_glob_shift = 0;
1033 // Get connection to cells on other procs
1034 multimap< int, int > loc2dist; // global cell ids on this proc -> other proc cells
1035 if ( isParallelMode() )
1037 cell_glob_shift = _domain_selector->getProcShift();
1039 set<int> loc_domains; // domains on this proc
1040 for ( int idom = 0; idom < _mesh.size(); ++idom )
1041 if ( _mesh[ idom ] )
1042 loc_domains.insert( idom );
1044 for ( int idom = 0; idom < _mesh.size(); ++idom )
1046 if ( !_mesh[idom] ) continue;
1047 vector<int> loc2glob_corr; // pairs of corresponding cells (loc_loc & glob_dist)
1048 retrieveDriver()->readLoc2GlobCellConnect(idom, loc_domains, _domain_selector, loc2glob_corr);
1049 //MEDMEM::STRING out;
1050 for ( int i = 0; i < loc2glob_corr.size(); i += 2 )
1052 int glob_here = _topology->convertCellToGlobal(idom,loc2glob_corr[i]);
1053 int glob_there = loc2glob_corr[i+1];
1054 loc2dist.insert ( make_pair( glob_here, glob_there));
1055 //out << glob_here << "-" << glob_there << " ";
1057 //cout << "\nRank " << _domain_selector->rank() << ": BndCZ: " << out << endl;
1061 //going across all cells
1063 map<int,int> cells_neighbours;
1064 for (int i=0; i< _topology->nbCells(); i++)
1068 vector<int> cells(50);
1070 for (unsigned inode=0; inode< cell2node[i+1].size(); inode++)
1072 int nodeid=cell2node[i+1][inode];
1073 for (unsigned icell=0; icell<node2cell[nodeid].size();icell++)
1074 cells_neighbours[node2cell[nodeid][icell]]++;
1077 int dimension = getMeshDimension();
1080 for (map<int,int>::const_iterator iter=cells_neighbours.begin(); iter != cells_neighbours.end(); iter++)
1082 if (iter->second >= dimension && iter->first != i+1)
1084 cells.push_back(iter->first + cell_glob_shift);
1085 // cells[isize++]=iter->first;
1088 // add neighbour cells from distant domains
1089 multimap< int, int >::iterator loc_dist = loc2dist.find( i+1 );
1090 for (; loc_dist!=loc2dist.end() && loc_dist->first==( i+1 ); ++loc_dist )
1091 cells.push_back( loc_dist->second );
1093 size[i]=cells.size();
1095 temp[i]=new int[size[i]];
1096 if (has_indivisible_regions)
1097 temp_edgeweight[i]=new int[size[i]];
1101 for (vector<int>::const_iterator iter=cells.begin(); iter!=cells.end();iter++)
1103 temp[i][itemp]=*iter;
1104 if (has_indivisible_regions)
1106 int tag1 = indivisible_tag[(i+1)-1];
1107 int tag2 = indivisible_tag[*iter-1];
1108 if (tag1==tag2 && tag1!=0)
1109 temp_edgeweight[i][itemp]=_topology->nbCells()*100000;
1111 temp_edgeweight[i][itemp]=1;
1115 cells_neighbours.clear();
1117 cout <<"end of graph definition"<<endl;
1118 int* index=new int[_topology->nbCells()+1];
1120 for (int i=0; i<_topology->nbCells(); i++)
1121 index[i+1]=index[i]+size[i];
1125 if (indivisible_tag!=0) delete [] indivisible_tag;
1127 //SKYLINEARRAY structure holding the cell graph
1128 array= new MEDMEM::MEDSKYLINEARRAY(_topology->nbCells(),index[_topology->nbCells()]-index[0]);
1129 array->setIndex(index);
1131 for (int i=0; i<_topology->nbCells(); i++)
1133 array->setI(i+1,temp[i]);
1137 if (has_indivisible_regions)
1139 edgeweights=new int[array->getLength()];
1140 for (int i=0; i<_topology->nbCells(); i++)
1142 for (int j=index[i]; j<index[i+1];j++)
1143 edgeweights[j-1]=temp_edgeweight[i][j-index[i]];
1144 delete[] temp_edgeweight[i];
1146 delete[]temp_edgeweight;
1152 cout<< "end of graph creation"<<endl;
1155 /*! Creates the partition corresponding to the cell graph and the partition number
1157 * \param nbdomain number of subdomains for the newly created graph
1159 * returns a topology based on the new graph
1161 Topology* MESHCollection::createPartition(int nbdomain,
1162 Graph::splitter_type split,
1163 const string& options_string,
1164 int* user_edge_weights,
1165 int* user_vertices_weights)
1167 if (nbdomain <1) throw MEDEXCEPTION("Number of subdomains must be >0");
1168 MEDMEM::MEDSKYLINEARRAY* array=0;
1171 MESSAGE_MED("Building cell graph");
1172 buildCellGraph(array,edgeweights);
1177 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
1178 _cell_graph=boost::shared_ptr<Graph>(new METISGraph(array,edgeweights));
1180 throw MEDEXCEPTION("METIS Graph is not available. Check your products, please.");
1184 #ifdef MED_ENABLE_SCOTCH
1185 _cell_graph=boost::shared_ptr<Graph>(new SCOTCHGraph(array,edgeweights));
1187 throw MEDEXCEPTION("SCOTCH Graph is not available. Check your products, please.");
1192 //!user-defined weights
1193 if (user_edge_weights!=0)
1194 _cell_graph->setEdgesWeights(user_edge_weights);
1195 if (user_vertices_weights!=0)
1196 _cell_graph->setVerticesWeights(user_vertices_weights);
1198 MESSAGE_MED("Partitioning graph");
1199 _cell_graph->partGraph(nbdomain,options_string,_domain_selector);
1202 // MEDMEM::STRING out("RESULT GRAPH #");
1203 // out << (_domain_selector?_domain_selector->rank():0) << ": ";
1204 // const int* part = _cell_graph->getPart();
1205 // int n = _cell_graph->nbVertices();
1206 // for ( int e=0; e < n; ++e )
1207 // out << part[e] <<" ";
1208 // cout << out << endl;
1211 MESSAGE_MED("Building new topology");
1212 //_cell_graph is a shared pointer
1213 Topology* topology = new ParallelTopology (_cell_graph, nbdomain, getMeshDimension());
1216 if (edgeweights!=0) delete[] edgeweights;
1217 //if (array!=0) delete array;
1218 MESSAGE_MED("End of partition creation");
1222 /*! Creates a topology for a partition specified by the user
1224 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
1226 * returns a topology based on the new partition
1228 Topology* MESHCollection::createPartition(const int* partition)
1230 MEDMEM::MEDSKYLINEARRAY* array=0;
1233 buildCellGraph(array,edgeweights);
1236 for (int i=0; i<_topology->nbCells(); i++)
1238 domains.insert(partition[i]);
1240 int nbdomain=domains.size();
1242 _cell_graph=boost::shared_ptr<Graph>(new UserGraph(array, partition, _topology->nbCells()));
1244 //_cell_graph is a shared pointer
1245 Topology* topology = new ParallelTopology (_cell_graph, nbdomain, getMeshDimension());
1247 //if (array!=0) delete array;
1252 /*! building Connect Zones for storing the informations
1253 * of the connectivity
1255 * The connect zones are created for every domain that has common nodes with
1258 * \param idomain domain number for which the connect zones are created
1261 // void MESHCollection::buildConnectZones(int idomain)
1263 // // constructing node/node correspondencies
1264 // vector<MEDMEM::MEDSKYLINEARRAY*> node_node_correspondency;
1265 // node_node_correspondency.resize(_topology->nbDomain());
1267 // cout << "Computing node/node corresp"<<endl;
1269 // _topology->computeNodeNodeCorrespondencies(idomain, node_node_correspondency );
1271 // for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1273 // // on regarde si une correspondance noeud/noeud a été trouvée
1274 // // entre idomain et idistant
1275 // // si oui, on crée une connectzone
1276 // if (node_node_correspondency[idistant]!=0)
1278 // MEDMEM::CONNECTZONE* cz= new MEDMEM::CONNECTZONE();
1279 // cz->setLocalMesh(_mesh[idomain]);
1280 // cz->setDistantMesh(_mesh[idistant]);
1281 // cz->setLocalDomainNumber(idomain);
1282 // cz->setDistantDomainNumber(idistant);
1283 // cz-> setName ("Connect zone defined by SPLITTER");
1284 // cz->setNodeCorresp(node_node_correspondency[idistant]);
1285 // _connect_zones.push_back(cz);
1288 // cout << "Computing node/node corresp"<<endl;
1290 // vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency;
1291 // cell_cell_correspondency.resize(_topology->nbDomain());
1292 // _topology->computeCellCellCorrespondencies(idomain, cell_cell_correspondency, _cell_graph.get());
1294 // for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1296 // //the connect zone has been created by the node/node computation
1297 // if (cell_cell_correspondency[idistant]!=0)
1299 // MEDMEM::CONNECTZONE* cz=0;
1300 // for (int icz=0; icz<_connect_zones.size();icz++)
1301 // if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
1302 // _connect_zones[icz]->getDistantDomainNumber()==idistant)
1303 // cz = _connect_zones[icz];
1305 // cz->setEntityCorresp(MED_EN::MED_CELL,MED_EN::MED_CELL, cell_cell_correspondency[idistant]);
1307 // throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");
1308 // //delete cell_cell_correspondency[idistant];
1314 //================================================================================
1316 * \brief Adds a group of joint faces
1317 * \param loc_face_ids - local numbers of faces
1318 * \param idomian - domain index where faces are local
1319 * \param idistant - the other domain index
1321 //================================================================================
1323 void MESHCollection::addJointGroup(const std::vector<int>& loc_face_ids, int idomain, int idistant)
1325 MEDMEM::MESHING* meshing = dynamic_cast<MEDMEM::MESHING*> (_mesh[idomain]);
1326 MED_EN::medEntityMesh constituent_entity = getSubEntity();
1328 MEDMEM::STRING jointname("joint_");
1329 jointname<<idistant+1;
1331 MEDMEM::GROUP * tmp_grp = new GROUP, * joint_group = tmp_grp;
1332 // try to find already present group with such a name
1333 // vector<MEDMEM::GROUP*> groups = meshing->getGroups( constituent_entity );
1334 // for ( int g = 0; g < groups.size(); ++g )
1335 // if ( groups[g]->getName() == jointname.str() )
1337 // joint_group = groups[g];
1340 // assure uniqueness of group name
1341 bool unique = false;
1342 vector<MEDMEM::GROUP*> groups = meshing->getGroups( constituent_entity );
1346 for ( int g = 0; unique && g < groups.size(); ++g )
1347 unique = ( groups[g]->getName() != jointname );
1349 jointname << "_" << idomain+1;
1352 joint_group->setName(jointname);
1353 joint_group->setMesh(meshing);
1354 joint_group->setEntity(constituent_entity);
1355 map<MED_EN::medGeometryElement, vector<int> > joint_types;
1357 int nbfaces = loc_face_ids.size();
1358 for (int i=0; i<nbfaces; i++)
1360 MED_EN::medGeometryElement type = meshing->getElementType(constituent_entity,loc_face_ids[i]);
1361 joint_types[type].push_back(loc_face_ids[i]);
1363 joint_group->setNumberOfGeometricType(joint_types.size());
1364 MED_EN::medGeometryElement* types=new MED_EN::medGeometryElement[joint_types.size()];
1365 int* nb_in_types=new int[joint_types.size()];
1366 int* group_index=new int[joint_types.size()+1];
1371 int* group_value=new int[nbfaces];
1372 for (map<MED_EN::medGeometryElement, vector<int> >::const_iterator iterj=joint_types.begin();
1373 iterj != joint_types.end();
1376 nb_in_types[itype]=(iterj->second).size();
1377 types[itype]=iterj->first;
1379 group_index[itype]=group_index[itype-1]+(iterj->second).size();
1380 for (int i=0; i< (iterj->second).size(); i++)
1381 group_value[iface++]=(iterj->second)[i];
1383 joint_group->setGeometricType(types);
1384 joint_group->setNumberOfElements(nb_in_types);
1385 joint_group->setNumber(group_index, group_value, /*shallowCopy=*/true);
1387 delete[] nb_in_types;
1389 if ( joint_group == tmp_grp )
1390 meshing->addGroup(*tmp_grp);
1391 tmp_grp->removeReference();
1394 /*! building Connect Zones for storing the informations
1395 * of the connectivity
1398 void MESHCollection::buildConnectZones()
1400 vector <map <MED_EN::medGeometryElement, vector<MEDSPLITTER_FaceModel*> > > face_map(_topology->nbDomain());
1401 map< pair<int,int>, MEDMEM::MEDSKYLINEARRAY*> cell_corresp_here;
1403 MED_EN::medEntityMesh constituent_entity = getSubEntity();
1405 if ( isParallelMode() )
1407 buildConnectZonesBetweenProcs(face_map, cell_corresp_here);
1410 cout << "Computing node/node corresp"<<endl;
1413 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1416 // constructing node/node correspondencies
1417 vector<MEDMEM::MEDSKYLINEARRAY*> node_node_correspondency(_topology->nbDomain());
1418 _topology->computeNodeNodeCorrespondencies(idomain, node_node_correspondency );
1420 for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1422 // on regarde si une correspondance noeud/noeud a été trouvée
1423 // entre idomain et idistant
1424 // si oui, on crée une connectzone
1425 if (node_node_correspondency[idistant]!=0)
1427 MEDMEM::CONNECTZONE* cz= new MEDMEM::CONNECTZONE();
1428 cz->setLocalMesh(_mesh[idomain]);
1429 cz->setDistantMesh(_mesh[idistant]);
1430 cz->setLocalDomainNumber(idomain);
1431 cz->setDistantDomainNumber(idistant);
1432 cz-> setName ("Connect zone defined by SPLITTER");
1433 cz->setNodeCorresp(node_node_correspondency[idistant]);
1434 _connect_zones.push_back(cz);
1438 cout << "Computing face corresp"<<endl;
1440 //creating faces if required
1441 if (_subdomain_boundary_creates)
1443 int global_face_id = _topology->getFaceNumber()+1;
1444 //int global_face_id = _topology->getMaxGlobalFace()+1;
1446 map <pair<int,int>, vector<int> > faces_in_joint;
1448 if ( !isParallelMode() )
1449 // taking faces that are already present in the mesh into account
1450 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1452 getFaces(idomain,face_map[idomain]);
1455 // creating faces that are located at the interface between
1458 vector <int> nb_added_groups( _topology->nbDomain(), 0 );
1460 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1462 vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency( _topology->nbDomain() );
1463 if ( !isParallelMode() )
1464 _topology->computeCellCellCorrespondencies(idomain, cell_cell_correspondency, _cell_graph.get());
1466 for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1468 if (idistant <= idomain) continue;
1470 MEDMEM::MEDSKYLINEARRAY* cell_correspondency = 0;
1471 if ( isParallelMode() )
1472 cell_correspondency = cell_corresp_here[ make_pair (idomain,idistant)];
1474 cell_correspondency = cell_cell_correspondency[idistant];
1476 //the connect zone has been created by the node/node computation
1478 if ( cell_correspondency )
1480 int nbcells = cell_correspondency->getNumberOf();
1481 const int* index = cell_correspondency->getIndex();
1482 const int* value = cell_correspondency->getValue();
1483 if ( isParallelMode() )
1484 global_face_id = _domain_selector->getFisrtGlobalIdOfSubentity( idomain, idistant );
1486 for (int ilocal=0; ilocal<nbcells; ilocal++)
1488 for (int icelldistant = index[ilocal]; icelldistant < index[ilocal+1]; icelldistant++)
1490 int distant_id = value[icelldistant-1];
1491 MEDSPLITTER_FaceModel* face = getCommonFace(idomain,ilocal+1,idistant,distant_id,global_face_id);
1492 face_map[idomain][face->getType()].push_back(face);
1493 MEDSPLITTER_FaceModel* face2 = getCommonFace(idistant,distant_id,idomain, ilocal+1,global_face_id);
1494 face_map[idistant][face->getType()].push_back(face2);
1495 faces_in_joint[make_pair(idomain,idistant)].push_back(global_face_id);
1503 for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1504 delete cell_cell_correspondency[idistant];
1508 _topology->recreateFaceMapping(face_map);
1510 //transforming the face_map into a constituent entity connectivity
1511 for (int idomain=0; idomain< _topology->nbDomain();idomain++)
1513 int nbtypes = face_map[idomain].size();
1514 vector<medGeometryElement> types;
1515 vector <int> nb_elems;
1518 MEDMEM::MESHING* meshing = dynamic_cast<MEDMEM::MESHING*> (_mesh[idomain]);
1519 if ( !meshing->getConnectivityptr() )
1520 continue; // no cells in idomain
1522 for (map <medGeometryElement, vector<MEDSPLITTER_FaceModel*> >::const_iterator iter= face_map[idomain].begin();
1523 iter != face_map[idomain].end(); iter ++)
1525 types.push_back(iter->first);
1526 int nb_elem_in_type = (iter->second).size();
1527 nb_elems.push_back(nb_elem_in_type);
1528 int nb_node_per_type=(iter->first)%100;
1529 int* connectivity= new int [nb_node_per_type*nb_elem_in_type];
1530 for (int ielem=0; ielem<nb_elem_in_type; ielem++)
1532 for (int inode=0; inode<nb_node_per_type; inode++)
1533 connectivity[ielem*nb_node_per_type+inode]=(*(iter->second)[ielem])[inode];
1535 conn.push_back(connectivity);
1538 //setting the faces in the mesh
1539 meshing->setNumberOfTypes(nbtypes,constituent_entity);
1540 meshing->setTypes(&types[0],constituent_entity);
1541 meshing->setNumberOfElements(&nb_elems[0],constituent_entity);
1543 for (int itype=0; itype<nbtypes; itype++)
1545 meshing->setConnectivity( constituent_entity, types[itype], conn[itype] );
1546 delete[]conn[itype];
1548 for (int idistant =0; idistant<_topology->nbDomain(); idistant++)
1550 map <pair<int,int>, vector<int> >::iterator iter;
1551 iter = faces_in_joint.find(make_pair(idomain,idistant));
1552 if (iter == faces_in_joint.end())
1554 iter = faces_in_joint.find (make_pair(idistant,idomain));
1555 if (iter == faces_in_joint.end())
1559 int nbfaces = (iter->second).size();
1560 vector<int> face_joint(nbfaces*2);
1561 MEDMEM::CONNECTZONE* cz=0;
1562 for (unsigned icz=0; icz<_connect_zones.size();icz++)
1563 if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
1564 _connect_zones[icz]->getDistantDomainNumber()==idistant)
1565 cz = _connect_zones[icz];
1567 int nbtotalfaces= _topology->getFaceNumber(idomain);
1569 //creating arrays for the MEDSKYLINEARRAY structure containing the joint
1570 int* index =new int[nbtotalfaces+1];
1571 for (int i=0; i<nbtotalfaces+1;i++)
1573 int*value=new int[nbfaces];
1576 vector<int> local_faces( nbfaces );
1577 for (int iface=0; iface<nbfaces; iface++)
1579 int iglobal = (iter->second)[iface];
1580 int localid=_topology->convertGlobalFace(iglobal,idomain);
1581 int distantid=_topology->convertGlobalFace(iglobal,idistant);
1582 faces.insert(make_pair(localid,distantid));
1583 local_faces[iface]=localid;
1588 for (map<int,int>::const_iterator iter=faces.begin();
1589 iter != faces.end();
1592 index[iter->first]=1;
1593 value[iloc++]=iter->second;
1596 for (int i=0; i<nbtotalfaces;i++)
1597 index[i+1]+=index[i];
1598 bool shallowcopy=true;
1599 MEDMEM::MEDSKYLINEARRAY* skarray=new MEDMEM::MEDSKYLINEARRAY(nbtotalfaces,nbfaces,index,value,shallowcopy);
1602 cz->setEntityCorresp(constituent_entity,constituent_entity,skarray);
1604 throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");
1605 // Creating a group of the faces constituting the joint
1606 addJointGroup( local_faces, idomain, idistant );
1607 nb_added_groups[ idomain ]++;
1611 if ( isParallelMode() )
1613 // Now all faces have got local ids and we can receive local ids from other procs.
1614 // Set face/face data to zones with other procs and create a group
1615 for (int icz=0; icz<_connect_zones.size();icz++)
1617 MEDMEM::CONNECTZONE* cz=_connect_zones[icz];
1618 if ( _domain_selector->isMyDomain( cz->getDistantDomainNumber()) ) continue;
1620 int glob_id = _domain_selector->getFisrtGlobalIdOfSubentity( cz->getLocalDomainNumber(),
1621 cz->getDistantDomainNumber());
1622 int nb_cz_faces = _domain_selector->getNbCellPairs( cz->getDistantDomainNumber(),
1623 cz->getLocalDomainNumber());
1624 vector< int > loc_ids_here( nb_cz_faces );
1625 for ( int i = 0; i < nb_cz_faces; ++i )
1626 loc_ids_here[i] = _topology->convertGlobalFace(glob_id++,cz->getLocalDomainNumber());
1628 int* loc_ids_dist = _domain_selector->exchangeSubentityIds( cz->getLocalDomainNumber(),
1629 cz->getDistantDomainNumber(),
1631 int nb_faces_here= _topology->getFaceNumber(cz->getLocalDomainNumber());
1632 int* face_index = new int[ nb_faces_here+1 ];
1634 for ( int loc_id = 0, i = 0; loc_id < nb_faces_here; ++loc_id)
1636 face_index[ loc_id+1 ] = face_index[ loc_id ];
1637 if ( i < loc_ids_here.size() && loc_ids_here[i] == loc_id+1 )
1639 face_index[ loc_id+1 ]++;
1643 MEDMEM::MEDSKYLINEARRAY* skarray=
1644 new MEDMEM::MEDSKYLINEARRAY(nb_faces_here, nb_cz_faces, face_index, loc_ids_dist, true);
1645 cz->setEntityCorresp(constituent_entity,constituent_entity,skarray);
1647 addJointGroup( loc_ids_here, cz->getLocalDomainNumber(), cz->getDistantDomainNumber());
1648 nb_added_groups[ cz->getLocalDomainNumber() ]++;
1652 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1655 for (map <medGeometryElement, vector<MEDSPLITTER_FaceModel*> >::const_iterator iter= face_map[idomain].begin();
1656 iter != face_map[idomain].end(); iter ++)
1657 for (unsigned i=0; i<(iter->second).size();i++)
1658 delete (iter->second)[i];
1660 if ( nb_added_groups[ idomain ] > 0 &&
1661 _mesh[idomain]->getNumberOfFamilies( constituent_entity ) > 0 )
1662 // needed because if there were face families before, driver won't
1663 // create families from just added groups (see MEDMEM_MedMeshDriver.cxx:3330),
1664 // actually it is a bug of driver - it must check presence of groups in families
1665 _mesh[idomain]->createFamilies();
1669 if ( isParallelMode() )
1670 // Excange info on types of constituent_entity needed while writing joints
1671 // to get ids local in geom type for distant procs
1672 _domain_selector->gatherEntityTypesInfo( _mesh, constituent_entity );
1674 cout << "Computing cell/cell corresp"<<endl;
1676 //Creating cell/cell correspondencies
1677 for (int idomain=0;idomain<_topology->nbDomain();idomain++)
1679 vector<MEDMEM::MEDSKYLINEARRAY*> cell_cell_correspondency( _topology->nbDomain() );
1680 if ( !isParallelMode() )
1681 _topology->computeCellCellCorrespondencies(idomain,cell_cell_correspondency,_cell_graph.get());
1683 for (int idistant=0; idistant< _topology->nbDomain(); idistant++)
1685 MEDMEM::MEDSKYLINEARRAY* cell_correspondency = 0;
1686 if ( isParallelMode() )
1687 cell_correspondency = cell_corresp_here[ make_pair (idomain,idistant)];
1689 cell_correspondency = cell_cell_correspondency[idistant];
1691 //the connect zone has been created by the node/node computation
1692 if ( cell_correspondency )
1694 MEDMEM::CONNECTZONE* cz=0;
1695 for (unsigned icz=0; icz<_connect_zones.size();icz++)
1696 if (_connect_zones[icz]->getLocalDomainNumber()==idomain &&
1697 _connect_zones[icz]->getDistantDomainNumber()==idistant)
1698 cz = _connect_zones[icz];
1700 cz->setEntityCorresp(MED_EN::MED_CELL,MED_EN::MED_CELL, cell_correspondency);
1702 throw MEDEXCEPTION("MESHCollection::buildConnectZones() -A connect zone should exist");
1708 /*! building Connect Zones for storing the informations
1709 * of the connectivity in the parallel mode
1712 void MESHCollection::buildConnectZonesBetweenProcs(TGeom2FacesByDomian & face_map,
1713 map< pair<int,int>, MEDMEM::MEDSKYLINEARRAY*> & cell_cell_correspondency_here)
1715 using namespace MED_EN;
1717 // graph over all procs
1718 auto_ptr<Graph> global_graph( _domain_selector->gatherGraph( _cell_graph.get() ));
1720 vector< vector< JointExchangeData > > joints_of_domain( _topology->nbDomain() );
1722 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1724 if ( !_domain_selector->isMyDomain( idomain )) continue;
1726 vector< JointExchangeData > & joints = joints_of_domain[ idomain ];
1727 joints.resize( _topology->nbDomain() );
1729 // Find corresponding cells on other procs
1731 const int* gra_index = global_graph->getGraph()->getIndex();
1732 const int* gra_value = global_graph->getGraph()->getValue();
1733 const int* partition = global_graph->getPart();
1734 const int dj = gra_index[0];
1736 vector< int > glob_cells_here( _topology->getCellNumber( idomain ));
1737 _topology->getCellList( idomain, & glob_cells_here[0]);
1738 for ( int loc_here = 0; loc_here < glob_cells_here.size(); ++loc_here )
1740 int glob_here = glob_cells_here[ loc_here ];
1741 for ( int j = gra_index[ glob_here-1 ]; j < gra_index[ glob_here ]; ++j )
1743 int glob_neighbor = gra_value[ j-dj ];
1744 int neighbor_dom = partition[ glob_neighbor-1 ];
1745 if ( neighbor_dom == idomain ) continue;
1747 if ( _domain_selector->isMyDomain( neighbor_dom ))
1749 joints[ neighbor_dom ].addCellCorrespondence
1750 (_mesh[idomain], neighbor_dom, idomain, glob_neighbor, glob_here, loc_here + 1,
1751 _topology->convertGlobalCell(glob_neighbor).second );
1755 joints[ neighbor_dom ].addCellCorrespondence
1756 (_mesh[idomain], neighbor_dom, idomain, glob_neighbor, glob_here, loc_here + 1 );
1761 global_graph.reset(); // free memory
1763 // set joints in a queue to exchange
1764 typedef map< int, JointExchangeData* > TOrderedJoints;
1765 TOrderedJoints queue;
1766 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1768 if ( !_domain_selector->isMyDomain( idomain )) continue;
1770 vector< JointExchangeData > & joints = joints_of_domain[ idomain ];
1771 for (int idist=0; idist<_topology->nbDomain(); ++idist )
1773 JointExchangeData& joint = joints[idist];
1775 int nb_cell_pairs = joint.nbCellPairs();
1776 if ( nb_cell_pairs == 0 )
1779 _domain_selector->setNbCellPairs( nb_cell_pairs, idist, idomain );
1781 joint.setMeshes( idist, _mesh[idist], idomain, _mesh[idomain] );
1783 if ( _domain_selector->isMyDomain( idist ))
1785 // a joint on this proc
1786 cell_cell_correspondency_here[ make_pair( idomain, idist )] = joint.makeCellCorrespArray();
1790 // a joint with distant proc
1791 joint.setConnectivity( & ((MEDMEM::MeshFuse*)_mesh[idomain])->getNodeNumbers()[0] );
1792 int order = _domain_selector->jointId( idomain, idist );
1793 queue[ order ] = & joint;
1797 // gather info on cell geom types needed to exchange joints
1798 _domain_selector->gatherEntityTypesInfo( _mesh, MED_EN::MED_CELL );
1800 // gather info on nb of sub-entities to compute their global numbers for joints
1801 _domain_selector->gatherNbOf( getSubEntity(), _mesh );
1802 _domain_selector->gatherNbCellPairs();
1803 if ( _subdomain_boundary_creates )
1805 // taking faces that are already present in the mesh into account
1806 for (int idomain=0; idomain<_topology->nbDomain(); idomain++)
1807 if ( _domain_selector->isMyDomain( idomain ))
1808 getFaces(idomain,face_map[idomain]);
1812 face_map.clear(); // mark for the joint not to create faces
1815 // exchange joint data with other procs and make CONNECTZONEs
1816 TOrderedJoints::iterator ord_joint = queue.begin();
1817 for ( ; ord_joint != queue.end(); ++ord_joint )
1819 JointExchangeData* joint = ord_joint->second;
1821 _domain_selector->exchangeJoint( joint );
1822 if ( _subdomain_boundary_creates )
1824 int first_sub_id = _domain_selector->getFisrtGlobalIdOfSubentity( joint->localDomain(),
1825 joint->distantDomain() );
1826 joint->setFisrtGlobalIdOfSubentity( first_sub_id );
1828 _connect_zones.push_back ( joint->makeConnectZone( face_map ));
1832 /*! projects old collection families on new collection families
1834 void MESHCollection::castFamilies(const MESHCollection& old_collection)
1836 vector <list<int> > element_array (_topology->nbDomain());
1838 //loop on old domains to create groups out of the existing families
1839 if (_family_splitting)
1840 for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1841 old_collection.getMesh(idomain)->createGroups();
1843 //definition of the entities array which
1844 //defines the entities over which the information is cast
1845 MED_EN::medEntityMesh entities[3];
1846 entities[0]=MED_EN::MED_NODE;
1847 entities[1]=getSubEntity();
1848 entities[2]=MED_EN::MED_CELL;
1850 for (int ientity=0; ientity<=2;ientity++)
1853 //int nbgroups = old_collection.getMesh(0)->getNumberOfGroups(entities[ientity]);
1855 map <string, set<int> > group_map;
1856 for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1858 if ( !old_collection.getMesh(idomain) ) continue;
1859 for (int igroup=0; igroup<old_collection.getMesh(idomain)->getNumberOfGroups(entities[ientity]); igroup++)
1862 MEDMEM::GROUP* group = (old_collection.getMesh(idomain)->getGroups(entities[ientity]))[igroup];
1863 //increments the number of groups if it is a new group
1864 //if (group_map.find(group->getName())==group_map.end())
1866 group_map[group->getName()].insert(idomain);
1867 // group_map.insert(make_pair(group->getName(), idomain);
1871 int nbgroups=group_map.size();
1872 vector <int> igroupold(old_collection._topology->nbDomain(),0);
1873 map<string,set<int> >::const_iterator iter=group_map.begin();
1875 for (int igroup=0; igroup<nbgroups; igroup++)
1877 vector <const MEDMEM::SUPPORT*> old_supports(old_collection._topology->nbDomain());
1878 string group_name = iter->first;
1881 //parameters stored for passing group description
1882 // from the old meshes to the new ones
1884 for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1886 // for (set<int>::iterator iter=group_map[group_name].begin(); iter!=group_map[group_name].end(); iter++)
1887 // cout << *iter<<" ";
1889 if (group_map[group_name].find(idomain)==group_map[group_name].end()) continue;
1891 //retrieves the group igroup on domain idomain
1892 MEDMEM::GROUP* group = (old_collection.getMesh(idomain)->getGroups(entities[ientity]))[igroupold[idomain]];
1893 old_supports[idomain] = static_cast<const MEDMEM::SUPPORT*> (group);
1894 igroupold[idomain]++;
1897 vector <MEDMEM::GROUP*>new_groups(_topology->nbDomain());
1898 vector <MEDMEM::SUPPORT*> new_supports(_topology->nbDomain());
1899 for (int i=0; i<_topology->nbDomain(); i++)
1901 new_groups[i]=new MEDMEM::GROUP();
1902 new_supports[i]=static_cast<MEDMEM::SUPPORT*>(new_groups[i]);
1904 castSupport(old_collection,old_supports,new_supports);
1906 //creating new groups from the previous list of elements
1907 for (int idomain=0; idomain <_topology->nbDomain(); idomain++)
1909 MEDMEM::MESHING* mesh_builder=static_cast<MEDMEM::MESHING*> (_mesh[idomain]);
1910 if ( new_supports[idomain] )
1911 mesh_builder->addGroup(*new_groups[idomain]);
1913 //groups are copied by the addGroup method,
1914 //so they can be safely deleted here
1915 for (int i=0; i<_topology->nbDomain(); i++)
1917 if ( new_supports[i] ) new_groups[i]->removeReference();
1925 void MESHCollection::castSupport(const MESHCollection& old_collection, vector<const MEDMEM::SUPPORT*>& old_support, vector<MEDMEM::SUPPORT*>& new_support)
1928 if (old_collection._topology->nbDomain() != (int)old_support.size())
1930 throw MED_EXCEPTION(STRING("Error : wrong call to MESHCollection::castSupport"));
1932 vector <list<int> > element_array (_topology->nbDomain());
1934 //parameters stored for passing description
1935 // from the old meshes to the new ones
1938 MED_EN::medEntityMesh entity;
1939 vector <string> support_name(1);
1940 support_name[0]="support";
1941 for (int inew=0; inew< _topology->nbDomain(); inew++)
1942 element_array[inew].clear();
1944 for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
1946 //retrieves the group igroup on domain idomain
1947 const MEDMEM::SUPPORT* support = old_support[idomain];
1948 if (old_support[idomain]==0) continue;
1949 name = support->getName();
1950 description=support->getDescription();
1951 int nbelem = support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1952 if (nbelem==0 && !_create_empty_groups) continue;
1955 if (support->isOnAllElements())
1957 list_of_elems = new int[nbelem];
1958 for (int i=0; i<nbelem;i++)
1959 list_of_elems[i]=i+1;
1962 list_of_elems = const_cast<int*> (support->getNumber(MED_EN::MED_ALL_ELEMENTS));
1964 int* array=new int[nbelem];
1968 entity = support->getEntity();
1973 case MED_EN::MED_CELL :
1975 local= new int[nbelem];
1977 old_collection.getTopology()->convertCellToGlobal(idomain,list_of_elems,nbelem,array);
1978 _topology->convertGlobalCellList(array,nbelem,local,ip);
1979 for (int i=0; i<nbelem; i++)
1980 // cell_arrays[ip[i]][local[i]]=id;
1982 // cout <<"(glob,ip,iloc)/nbelem"<<array[i]<<" "<<ip[i]<<" "<<local[i]<<"/"<<nbelem<<endl;
1983 element_array[ip[i]].push_back(local[i]);
1986 case MED_EN::MED_FACE :
1987 case MED_EN::MED_EDGE :
1988 old_collection.getTopology()->convertFaceToGlobal(idomain,list_of_elems,nbelem,array);
1989 _topology->convertGlobalFaceListWithTwins(array,nbelem,local,ip,full_array,size);
1990 for (int i=0; i<size; i++)
1991 element_array[ip[i]].push_back(local[i]);
1992 delete[] full_array;
1994 case MED_EN::MED_NODE :
1995 old_collection.getTopology()->convertNodeToGlobal(idomain,list_of_elems,nbelem,array);
1996 _topology->convertGlobalNodeListWithTwins(array,nbelem,local,ip,full_array,size);
1997 for (int i=0; i<size; i++)
1998 element_array[ip[i]].push_back(local[i]);
1999 delete[] full_array;
2007 if (support->isOnAllElements()) delete[] list_of_elems;
2010 //creating new groups from the previous list of elements
2011 for (int idomain=0; idomain <_topology->nbDomain(); idomain++)
2013 if ( _mesh[idomain]->getNumberOfNodes() < 1 ||
2014 (element_array[idomain].empty() && !_create_empty_groups))
2016 new_support[idomain]->removeReference();
2017 new_support[idomain]=0;
2020 MEDMEM::SUPPORT* support= new_support[idomain];
2021 support->setName(name);
2022 support->setMesh(_mesh[idomain]);
2023 support->setDescription(description);
2024 support->setEntity(entity);
2026 if ( !element_array[idomain].empty() ) /* if() was added for issue 0021576
2027 to prevent creation of faces */
2029 element_array[idomain].sort();
2030 element_array[idomain].unique();
2032 if (entity != MED_EN::MED_NODE)
2033 support->fillFromElementList(element_array[idomain]);
2035 support->fillFromNodeList(element_array[idomain]);
2040 void MESHCollection::castField(const MESHCollection& old_collection, const string& fieldname, int itnumber, int ordernumber)
2042 int type=old_collection.getDriver()->getFieldType(fieldname);
2043 char field_char[80];
2044 strcpy(field_char,fieldname.c_str());
2047 castFields<int>(old_collection, field_char, itnumber, ordernumber);
2049 castFields<double>(old_collection, field_char, itnumber, ordernumber);
2052 void MESHCollection::castAllFields(const MESHCollection& initial_collection)
2054 vector <string> field_names;
2055 vector <int> iternumber;
2056 vector <int> ordernumber;
2058 initial_collection.getDriver()->readFileStruct(field_names,iternumber,ordernumber,types);
2060 for (unsigned i=0; i<field_names.size(); i++)
2062 char field_char[80];
2063 strcpy(field_char,field_names[i].c_str());
2065 // choosing whether the field is of int or double type
2067 castFields<int>(initial_collection, field_char, iternumber[i], ordernumber[i]);
2069 castFields<double>(initial_collection, field_char, iternumber[i], ordernumber[i]);
2073 void MESHCollection::createNodalConnectivity(const MESHCollection& initial_collection,int idomain, MED_EN::medEntityMesh entity)
2075 MESSAGE_MED ("beginning of createNodalConnectivity for entity "<<entity);
2078 MEDMEM::MESHING* mesh_builder = static_cast<MEDMEM::MESHING*>(_mesh[idomain]);
2081 //number of elements per type
2082 std::map<MED_EN::medGeometryElement,int> type_numbers;
2084 //creating arrays for storing global numbers and cell types
2087 case MED_EN::MED_CELL:
2088 dimension=initial_collection.getMeshDimension();
2089 nb_elems=_topology->getCellNumber(idomain);
2091 case MED_EN::MED_EDGE:
2092 case MED_EN::MED_FACE:
2093 dimension=initial_collection.getMeshDimension()-1;
2094 nb_elems=_topology->getFaceNumber(idomain);
2101 if (nb_elems == 0) return;
2102 SCRUTE_MED(nb_elems);
2105 int *list= new int[nb_elems];
2106 MED_EN::medGeometryElement *cell_type_list= new MED_EN::medGeometryElement[nb_elems];
2109 // cout << "Beginning of retrieval "<<endl;
2110 //retrieving global id list
2113 case MED_EN::MED_CELL:
2114 _topology->getCellList(idomain,list);
2116 case MED_EN::MED_EDGE:
2117 case MED_EN::MED_FACE:
2118 _topology->getFaceList(idomain,list);
2125 //retrieving cell_types
2126 initial_collection.getTypeList(list,nb_elems,entity,cell_type_list);
2127 // cout <<"end of type retrieval"<<endl;
2128 //vector containing the number of cells per type
2129 type_numbers.clear();
2130 for (int icell=0; icell<nb_elems; icell++)
2132 map<MED_EN::medGeometryElement,int>::iterator iter= type_numbers.find(cell_type_list[icell]);
2133 if (iter!=type_numbers.end())
2136 type_numbers[cell_type_list[icell]]=1;
2139 //cout << "Nombre de tetras"<<type_numbers[304]<<endl;
2140 int nb_present_types=type_numbers.size();
2142 //setting the list of cells for each type
2143 map<MED_EN::medGeometryElement,int> index;
2145 map<MED_EN::medGeometryElement,int*> type_cell_list;
2147 MED_EN::MESH_ENTITIES::const_iterator currentEntity;
2148 std::map<MED_EN::medGeometryElement,int>::const_iterator iter;
2149 //currentEntity = MED_EN::meshEntities.find(entity);
2150 for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2152 MED_EN::medGeometryElement type = iter->first;
2153 if (!isDimensionOK(type,dimension)) continue;
2154 //if (iter->second==0) continue;
2156 type_cell_list[type]=new int[type_numbers[type]];
2157 // cout << "type :"<<type<<" nb:"<<type_numbers[type]<<endl;
2160 for (int icell=0; icell<nb_elems; icell++)
2162 type_cell_list[cell_type_list[icell]][index[cell_type_list[icell]]++]=list[icell];
2166 delete[]cell_type_list;
2168 //setting the list of present types
2169 int* present_type_numbers=new int[nb_present_types];
2170 MED_EN::medGeometryElement* type_array = new MED_EN::medGeometryElement[nb_present_types];
2171 MESSAGE_MED("Nb de types presents "<<nb_present_types);
2173 for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2175 MED_EN::medGeometryElement type = iter->first;
2176 if (!isDimensionOK(type,dimension)) continue;
2178 type_array[itype]=type;
2180 present_type_numbers[itype]=type_numbers[type];
2182 MESSAGE_MED("Nombre d'elements de type "<<type<<" : "<<type_numbers[type]);
2186 //retrieving connectivity in global numbering for each type
2187 map<MED_EN::medGeometryElement,int*> type_connectivity;
2188 vector<int> polygon_conn;
2189 vector<int> polygon_conn_index;
2190 vector<int> polyhedron_conn;
2191 vector<int> polyhedron_conn_index;
2192 vector<int> polyhedron_face_index;
2197 for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2199 MED_EN::medGeometryElement type = iter->first;
2202 if (!isDimensionOK(type,dimension)) continue;
2203 if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
2205 int nbnode_per_type = (int)type%100;
2206 type_connectivity[type]=new int[type_numbers[type]*nbnode_per_type];
2207 initial_collection.getNodeConnectivity(type_cell_list[type],type_numbers[type],entity,type,type_connectivity[type]);
2209 else if (type == MED_EN::MED_POLYGON && dimension==2)
2211 initial_collection.getPolygonNodeConnectivity(type_cell_list[type],type_numbers[type],entity,polygon_conn,polygon_conn_index);
2213 else if (type == MED_EN::MED_POLYHEDRA && dimension==3)
2215 initial_collection.getPolyhedraNodeConnectivity(type_cell_list[type],type_numbers[type],entity,polyhedron_conn,polyhedron_conn_index);
2217 delete[] type_cell_list[type];
2220 //creating node mapping
2221 //!TODO : compute the total number of nodes
2222 if (entity==MED_EN::MED_CELL)
2224 _topology->createNodeMapping(type_connectivity,type_numbers,polygon_conn,polygon_conn_index,
2225 polyhedron_conn,polyhedron_conn_index,polyhedron_face_index,idomain);
2228 //converting node global numberings to local numberings
2229 //for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
2230 for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2232 MED_EN::medGeometryElement type = iter->first;
2234 if (!isDimensionOK(type, dimension)) continue;
2235 if (type_numbers[type]==0) continue;
2236 if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
2238 int nbnode_per_type = (int)type%100;
2239 _topology->convertToLocal2ndVersion(type_connectivity[type],type_numbers[type]*nbnode_per_type,idomain);
2241 else if (type == MED_EN::MED_POLYGON && dimension==2)
2243 int nbpoly = type_numbers[type];
2244 _topology->convertToLocal2ndVersion(&polygon_conn[0], polygon_conn_index[nbpoly]-1, idomain);
2246 else if (type == MED_EN::MED_POLYHEDRA && dimension==3)
2248 int nbpoly = type_numbers[type];
2249 _topology->convertToLocal2ndVersion(&polyhedron_conn[0], polyhedron_face_index[polyhedron_conn_index[nbpoly]-1]-1, idomain);
2255 //writing coordinates
2256 if (entity==MED_EN::MED_CELL)
2258 //setting coordinates from initial_collection coordinates
2259 int nbnode=_topology->getNodeNumber(idomain);
2260 MESSAGE_MED("Number of nodes on domain "<< idomain <<" : "<<nbnode);
2262 double* coordinates=new double[initial_collection.getSpaceDimension()*nbnode];
2263 int* node_list=new int[nbnode];
2264 _topology->getNodeList(idomain,node_list);
2265 initial_collection.getCoordinates(node_list,nbnode,coordinates);
2268 // redundant specification of number of nodes is required!! MED imperfection, sorry...
2270 //TODO : change MEDMEM so that it accepts a direct setting of coordinates
2271 // (in the present version, it is deep-copied)
2272 mesh_builder->setCoordinates(initial_collection.getSpaceDimension(),
2273 nbnode, coordinates, initial_collection.getSystem(),
2274 MED_EN::MED_FULL_INTERLACE);
2275 delete [] coordinates;
2278 int nb_plain_types=0;
2279 for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2281 MED_EN::medGeometryElement type = iter->first;
2283 if (!isDimensionOK(type, dimension)) continue;
2284 if (type_numbers[type]==0) continue;
2287 mesh_builder->setNumberOfTypes(nb_plain_types,entity);
2288 mesh_builder->setTypes(type_array,entity);
2289 mesh_builder->setNumberOfElements(present_type_numbers,entity);
2291 delete[]present_type_numbers;
2293 //setting node connectivities
2294 for (iter = type_numbers.begin();iter != type_numbers.end(); iter++)
2296 MED_EN::medGeometryElement type = iter->first;
2298 if (!isDimensionOK(type,dimension)) continue;
2299 if (type_numbers[type]==0) continue;
2301 if (type != MED_EN::MED_POLYHEDRA && type != MED_EN::MED_POLYGON)
2303 mesh_builder->setConnectivity(entity,type,type_connectivity[type]);
2304 delete[] type_connectivity[type];
2306 else if (type == MED_EN::MED_POLYGON && dimension ==2)
2308 mesh_builder->setConnectivity(entity,type,&polygon_conn[0],&polygon_conn_index[0]);
2310 else if (type == MED_EN::MED_POLYHEDRA && dimension ==3)
2312 mesh_builder->setConnectivity(entity,type,&polyhedron_conn[0],&polyhedron_conn_index[0]);
2315 MESSAGE_MED("end of createNodalConnectivity");
2319 /*! retrieves the faces that are present in a mesh and stores them in a
2320 * dynamic structure made of a map of MEDSPLITTER_FaceModel
2322 * \param idomain domain id on which the faces are collected
2323 * \param face_map container storing the faces
2325 void MESHCollection::getFaces(int idomain,
2326 map<MED_EN::medGeometryElement, vector<MEDSPLITTER_FaceModel*> >& face_map)
2328 MED_EN::medEntityMesh constituent_entity = getSubEntity();
2329 const medGeometryElement* types;
2332 types = _mesh[idomain]->getTypes(constituent_entity);
2333 if ( !types ) return;
2335 catch(MEDEXCEPTION&){ return;}
2337 int nbtypes = _mesh[idomain]->getNumberOfTypes(constituent_entity);
2338 const int* global_numbering= _mesh[idomain]->getGlobalNumberingIndex(constituent_entity);
2339 int* conn = const_cast<int*> (_mesh[idomain]->getConnectivity(MED_EN::MED_NODAL,constituent_entity, MED_EN::MED_ALL_ELEMENTS));
2340 for (int itype=0; itype<nbtypes; itype++)
2342 for (int iface=global_numbering[itype]; iface<global_numbering[itype+1]; iface++)
2344 MEDSPLITTER_FaceModel* face_model = new MEDSPLITTER_FaceModel();
2345 MED_EN::medGeometryElement type = types[itype];
2346 face_model->setType(type);
2347 int nbnodes = type%100;
2348 face_model->setNbNodes(nbnodes);
2349 face_model->setGlobal(_topology->convertFaceToGlobal(idomain,iface));
2350 for (int i=0; i<nbnodes; i++)
2352 (*face_model)[i]=*conn++;
2354 face_map[type].push_back(face_model);
2359 /*! retrieves the face that is common to two cells located on two different processors
2361 * \param ip1 domain id for cell 1
2362 * \param ilocal1 cell id for cell 1
2363 * \param ip2 domain id for cell 2
2364 * \param ilocal2 cell id for cell 2
2365 * \param face_index global index for the newly created face
2367 MEDSPLITTER_FaceModel* MESHCollection::getCommonFace(int ip1,int ilocal1,int ip2,int ilocal2,int face_index)
2369 MED_EN::medGeometryElement type1 = _mesh[ip1]->getElementType(MED_EN::MED_CELL,ilocal1);
2370 MEDMEM::CELLMODEL celltype1 (type1);
2372 const int* conn_index1 = _mesh[ip1]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
2373 const int* conn1 = _mesh[ip1]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
2375 // MED_EN::medGeometryElement type2 = _mesh[ip2]->getElementType(MED_EN::MED_CELL,ilocal2);
2376 //MEDMEM::CELLTYPE celltype2 (type2);
2377 const int* conn_index2 = _mesh[ip2]->getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
2378 const int* conn2 = _mesh[ip2]->getConnectivity(MED_EN::MED_NODAL,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
2380 vector<int> nodes1, nodes1_local;
2382 for (int i= conn_index1[ilocal1-1]; i<conn_index1[ilocal1]; i++)
2384 nodes1.push_back(_topology->convertNodeToGlobal(ip1,*(conn1+i-1)));
2385 nodes1_local.push_back( conn1[i-1] );
2387 for (int i= conn_index2[ilocal2-1]; i<conn_index2[ilocal2]; i++)
2388 nodes2.push_back(_topology->convertNodeToGlobal(ip2,*(conn2+i-1)));
2390 return MEDSPLITTER_FaceModel::getCommonFace( &nodes1[0], &nodes1_local[0], celltype1,
2391 &nodes2[0], nodes2.size(), face_index);
2394 //================================================================================
2396 * \brief Makes a face common for two given cells
2397 * \param nodes1 - globl nodes of the first cell
2398 * \param nodes1_local - local nodes of the first cell
2399 * \param celltype1 - cell model of the first cell
2400 * \param nodes2 - globl nodes of the second cell
2401 * \param nb_nodes2 - nb of nodes of the second cell
2402 * \param global_id - id of the new common face
2404 //================================================================================
2406 MEDSPLITTER_FaceModel*
2407 MEDSPLITTER_FaceModel::getCommonFace(const int* nodes1,
2408 const int* nodes1_local,
2409 const MEDMEM::CELLMODEL& celltype1,
2414 int nbfaces= celltype1.getNumberOfConstituents(1);
2415 int ** faces = celltype1.getConstituents(1);
2416 MED_EN::medGeometryElement* types = celltype1.getConstituentsType(1);
2418 int dimension=celltype1.getDimension();
2420 while (iface<nbfaces)
2422 //SCRUTE_MED (iface);
2423 int nbnodes= types[iface]%100;
2424 const int* nodes = celltype1.getNodesConstituent(1,iface+1);
2426 for (int i=0; i<nbnodes;i++)
2428 for (int i2=0; i2<nb_nodes2; i2++)
2430 if (nodes1[nodes[i]-1]==nodes2[i2]) common_nodes++;
2433 if (common_nodes>=dimension) break;
2438 throw MEDEXCEPTION("MEDSPLITTER::getCommonFace - No common face found !");
2440 MEDSPLITTER_FaceModel* face_model = new MEDSPLITTER_FaceModel;
2441 face_model->setType(types[iface]);
2442 int nbnodes = types[iface]%100;
2443 face_model->setNbNodes(nbnodes);
2444 face_model->setGlobal(global_id);
2445 for (int i=0; i<nbnodes; i++)
2446 (*face_model)[i]=nodes1_local[faces[iface][i]-1];