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
23 #include "InterpKernelHashMap.hxx"
25 #include "MEDMEM_CellModel.hxx"
26 #include "MEDMEM_ConnectZone.hxx"
27 #include "MEDMEM_DriversDef.hxx"
28 #include "MEDMEM_Exception.hxx"
29 #include "MEDMEM_Mesh.hxx"
30 #include "MEDMEM_MeshFuse.hxx"
31 #include "MEDMEM_SkyLineArray.hxx"
32 #include "MEDMEM_Utilities.hxx"
34 #include "MEDSPLITTER_MESHCollection.hxx"
35 #include "MEDSPLITTER_Topology.hxx"
36 #include "MEDSPLITTER_Graph.hxx"
37 #include "MEDSPLITTER_ParallelTopology.hxx"
39 using namespace INTERP_KERNEL;
41 using namespace MEDSPLITTER;
44 ParallelTopology::ParallelTopology():m_nb_domain(0),m_mesh_dimension(0)
47 //!constructing topology according to mesh collection
48 ParallelTopology::ParallelTopology(const vector<MEDMEM::MESH*>& meshes,
49 const vector<MEDMEM::CONNECTZONE*>& cz,
50 vector<int*>& cellglobal,
51 vector<int*>& nodeglobal,
52 vector<int*>& faceglobal):m_nb_domain(meshes.size())/*,m_mesh_dimension(meshes[0]->getMeshDimension())*/
56 int index_node_global=0;
57 int index_face_global=0;
59 m_nb_cells.resize(m_nb_domain);
60 m_nb_nodes.resize(m_nb_domain);
61 m_nb_faces.resize(m_nb_domain);
63 m_loc_to_glob.resize(m_nb_domain);
64 m_node_loc_to_glob.resize(m_nb_domain);
65 m_face_loc_to_glob.resize(m_nb_domain);
67 MED_EN::medEntityMesh constituent_entity;
69 bool parallel_mode = false;
70 for (int idomain=0; !parallel_mode && idomain<m_nb_domain; idomain++)
71 parallel_mode = (!meshes[idomain]);
73 for (int idomain=0; idomain<m_nb_domain; idomain++)
75 if ( !meshes[idomain] )
77 m_mesh_dimension = meshes[idomain]->getMeshDimension();
78 constituent_entity = (m_mesh_dimension == 3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE );
81 m_nb_cells[idomain]=meshes[idomain]->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
82 // cout << "Nb cells (domain "<<idomain<<") = "<<m_nb_cells[idomain];
83 m_loc_to_glob[idomain].resize(m_nb_cells[idomain]);
85 if (cellglobal[idomain]==0 || parallel_mode)
87 MESSAGE_MED("Creating global numbering");
88 //creating global numbering from scratch
89 for (int i=0; i<m_nb_cells[idomain]; i++)
92 m_glob_to_loc[index_global]=make_pair(idomain,i+1);
93 //m_loc_to_glob[make_pair(idomain,i+1)]=index_global;
94 m_loc_to_glob[idomain][i]=index_global;
95 // cout<<"glob:"<<index_global<<" --> ("<<idomain<<","<<i+1<<")"<<endl;
98 //using global numbering coming from a previous numbering
101 MESSAGE_MED("Using former global numbering");
102 for (int i=0; i<m_nb_cells[idomain]; i++)
104 int global=cellglobal[idomain][i];
105 m_glob_to_loc[global]=make_pair(idomain,i+1);
106 //m_loc_to_glob[make_pair(idomain,i+1)]=global;
107 m_loc_to_glob[idomain][i]=global;
109 // cout<<"glob:"<<global<<" --> ("<<idomain<<","<<i+1<<")"<<endl;
116 m_nb_total_cells=index_global;
117 m_nb_cells[0]=index_global;
118 m_node_loc_to_glob[idomain].resize(meshes[idomain]->getNumberOfNodes());
119 for (int i=0; i<meshes[idomain]->getNumberOfNodes(); i++)
121 m_node_glob_to_loc.insert(make_pair(i+1,make_pair(0,i+1)));
122 //m_node_loc_to_glob.insert(make_pair(make_pair(0,i+1), i+1));
123 m_node_loc_to_glob[0][i]=i+1;
125 m_nb_total_nodes=meshes[idomain]->getNumberOfNodes();
126 m_nb_nodes[0]=m_nb_total_nodes;
128 // meshes[idomain]->getConnectivity( MED_EN::MED_DESCENDING, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
129 int nbfaces=meshes[idomain]->getNumberOfElements(constituent_entity,MED_EN::MED_ALL_ELEMENTS);
130 m_face_loc_to_glob[idomain].resize(nbfaces);
131 for (int i=0; i<nbfaces; i++)
133 m_face_glob_to_loc.insert(make_pair(i+1,make_pair(0,i+1)));
134 //m_face_loc_to_glob.insert(make_pair(make_pair(0,i+1), i+1));
135 m_face_loc_to_glob[0][i]=i+1;
137 m_nb_total_faces=nbfaces;
138 m_nb_faces[0]=nbfaces;
139 MESSAGE_MED ("nb total cells "<< m_nb_total_cells);
140 MESSAGE_MED("nb total nodes "<< m_nb_total_nodes);
141 MESSAGE_MED("nb total faces "<< m_nb_total_faces);
146 m_nb_nodes[idomain]=meshes[idomain]->getNumberOfNodes();
147 INTERP_KERNEL::HashMap <int,pair<int,int> > local2distant;
148 m_node_loc_to_glob[idomain].resize(m_nb_nodes[idomain]);
149 for (unsigned icz=0; icz<cz.size(); icz++)
151 if (cz[icz]->getLocalDomainNumber() == idomain &&
152 cz[icz]->getLocalDomainNumber()>cz[icz]->getDistantDomainNumber())
154 int nb_node= cz[icz]->getNodeNumber();
155 const int* node_corresp=cz[icz]->getNodeCorrespValue();
156 int distant_ip = cz[icz]->getDistantDomainNumber();
157 for (int i=0; i< nb_node; i++)
159 int local= node_corresp[i*2];
160 int distant = node_corresp[i*2+1];
161 local2distant.insert(make_pair(local, make_pair(distant_ip,distant)));
165 // setting mappings for all nodes
166 if (nodeglobal[idomain]==0)
168 for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
170 if (local2distant.find(inode+1)==local2distant.end())
173 m_node_glob_to_loc.insert(make_pair(index_node_global,make_pair(idomain,inode+1)));
174 //m_node_loc_to_glob[make_pair(idomain,inode+1)]=index_node_global;
175 m_node_loc_to_glob[idomain][inode]=index_node_global;
179 int ip = (local2distant.find(inode+1)->second).first;
180 int distant = (local2distant.find(inode+1)->second).second;
181 //int global_number=m_loc_to_glob[make_pair(ip,distant)];
182 int global_number=m_loc_to_glob[ip][distant-1];
183 m_node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode+1)));
184 //m_node_loc_to_glob[make_pair(idomain,inode+1)]=global_number;
185 m_node_loc_to_glob[idomain][inode]=global_number;
189 //using former node numbering
191 {// cout << "("<<idomain<<","<<i+1<<")->"<<i+1<<endl;
192 for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
194 int global_number=nodeglobal[idomain][inode];
195 // cout << "global_number "<<global_number<<endl;
196 m_node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode+1)));
197 //m_node_loc_to_glob[make_pair(idomain,inode+1)]=global_number;
198 m_node_loc_to_glob[idomain][inode]=global_number;
203 //creating dimension d-1 component mappings
205 // meshes[idomain]->getConnectivity( MED_EN::MED_DESCENDING, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
206 m_nb_faces[idomain]=meshes[idomain]->getNumberOfElements(constituent_entity,MED_EN::MED_ALL_ELEMENTS);
207 MESSAGE_MED ("Nb faces domain " << idomain<<m_nb_faces[idomain]);
208 m_face_loc_to_glob[idomain].resize(m_nb_faces[idomain]);
209 local2distant.clear();
210 for (unsigned icz=0; icz<cz.size(); icz++)
212 if (cz[icz]->getLocalDomainNumber() == idomain &&
213 cz[icz]->getLocalDomainNumber()>cz[icz]->getDistantDomainNumber())
215 int nb_face= cz[icz]->getFaceNumber();
216 const int* face_corresp=cz[icz]->getFaceCorrespValue();
217 int distant_ip = cz[icz]->getDistantDomainNumber();
218 for (int i=0; i< nb_face; i++)
220 int local= face_corresp[i*2];
221 int distant = face_corresp[i*2+1];
222 local2distant.insert(make_pair(local, make_pair(distant_ip,distant)));
226 // setting mappings for all faces
227 if (faceglobal[idomain]==0)
229 for (int iface=0; iface<m_nb_faces[idomain]; iface++)
231 if (local2distant.find(iface+1)==local2distant.end())
234 m_face_glob_to_loc.insert(make_pair(index_face_global,make_pair(idomain,iface+1)));
235 //m_face_loc_to_glob[make_pair(idomain,iface+1)]=index_face_global;
236 m_face_loc_to_glob[idomain][iface]=index_face_global;
240 int ip = (local2distant.find(iface+1)->second).first;
241 int distant = (local2distant.find(iface+1)->second).second;
242 //int global_number=m_loc_to_glob[make_pair(ip,distant)];
243 int global_number=m_loc_to_glob[ip][distant-1];
244 m_face_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,iface+1)));
245 //m_face_loc_to_glob[make_pair(idomain,iface+1)]=global_number;
246 m_face_loc_to_glob[idomain][iface]=global_number;
250 //using former face numbering
253 for (int iface=0; iface<m_nb_faces[idomain]; iface++)
255 int global_number=faceglobal[idomain][iface];
256 //SCRUTE_MED(global_number);
257 m_face_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,iface+1)));
258 //m_face_loc_to_glob[make_pair(idomain,iface+1)]=global_number;
259 m_face_loc_to_glob[idomain][iface]=global_number;
264 m_nb_total_cells=index_global;
265 m_nb_total_nodes=index_node_global;
266 m_nb_total_faces=index_face_global;
267 SCRUTE_MED(m_nb_total_cells);
268 SCRUTE_MED(m_nb_total_faces);
269 SCRUTE_MED(m_nb_total_nodes);
274 //!constructing ParallelTopology from an old topology and a graph
275 ParallelTopology::ParallelTopology(boost::shared_ptr<Graph> graph, int nb_domain, int mesh_dimension):
276 m_nb_domain(nb_domain),
277 m_mesh_dimension(mesh_dimension),
278 m_nb_cells(graph->nbVertices()),
281 m_nb_cells.resize(m_nb_domain);
282 m_nb_nodes.resize(m_nb_domain);
283 m_nb_faces.resize(m_nb_domain);
285 m_loc_to_glob.resize(m_nb_domain);
286 m_node_loc_to_glob.resize(m_nb_domain);
287 m_face_loc_to_glob.resize(m_nb_domain);
289 // used in parallel mode only
290 m_cell_loc_to_glob_fuse.resize(m_nb_domain);
291 m_face_loc_to_glob_fuse.resize(m_nb_domain);
293 for (int i=0; i<m_nb_domain; i++)
296 const int* part = graph-> getPart();
297 m_nb_total_cells= graph->nbVertices();
299 for (int icell=0; icell<m_nb_total_cells; icell++)
301 int idomain = part[icell];
302 m_nb_cells[idomain]++;
303 //m_loc_to_glob[make_pair(idomain,m_nb_cells[idomain])]=icell+1;
304 m_loc_to_glob[idomain].push_back(icell+1);
305 m_glob_to_loc[icell+1]=make_pair(idomain,m_nb_cells[idomain]);
308 for (int idomain=0; idomain<m_nb_domain; idomain++)
309 MESSAGE_MED("Nombre de cellules dans le domaine "<< idomain <<" : "<<m_nb_cells[idomain]);
311 SCRUTE_MED(m_nb_total_cells);
315 ParallelTopology::~ParallelTopology()
320 /*!Converts a list of global node numbers
321 * to a distributed array with local cell numbers.
323 * If a node in the list is represented on several domains,
324 * only the first value is returned
326 void ParallelTopology::convertGlobalNodeList(const int* node_list, int nbnode, int* local, int* ip)
328 if (m_node_glob_to_loc.empty())
329 throw MEDMEM::MEDEXCEPTION("convertGlobalNodeList - Node mapping has not yet been built");
330 for (int i=0; i< nbnode; i++)
332 pair<int,int> local_node = m_node_glob_to_loc.find(node_list[i])->second;
333 ip[i]=local_node.first;
334 local[i]=local_node.second;
338 /*!Converts a list of global node numbers on domain ip
339 * to a distributed array with local cell numbers.
341 * If a node in the list is represented on several domains,
342 * only the value with domain ip is returned
345 void ParallelTopology::convertGlobalNodeList(const int* node_list, int nbnode, int* local, int ip)
347 if (m_node_glob_to_loc.empty())
348 throw MEDMEM::MEDEXCEPTION("convertGlobalNodeList - Node mapping has not yet been built");
350 for (int i=0; i< nbnode; i++)
352 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
353 pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(node_list[i]);
354 for (mmiter it=range.first; it !=range.second; it++)
356 int ipfound=(it->second).first;
358 local[i]=(it->second).second;
363 /*!Converts a list of global node numbers
364 * to a distributed array with local cell numbers.
366 * If a node in the list is represented on several domains,
367 * all the values are put in the array
369 void ParallelTopology::convertGlobalNodeListWithTwins(const int* node_list, int nbnode, int*& local, int*& ip,int*& full_array, int& size)
371 if (m_node_glob_to_loc.empty())
372 throw MEDMEM::MEDEXCEPTION("convertGlobalNodeList - Node mapping has not yet been built");
375 for (int i=0; i< nbnode; i++)
377 int count= m_node_glob_to_loc.count(node_list[i]);
379 // cout << "noeud " << node_list[i]<< " doublon d'ordre " << count<<endl;
385 full_array=new int[size];
386 for (int i=0; i< nbnode; i++)
388 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
389 pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(node_list[i]);
390 for (mmiter it=range.first; it !=range.second; it++)
392 ip[index]=(it->second).first;
393 local[index]=(it->second).second;
394 full_array [index]=node_list[i];
401 /*!Converts a list of global face numbers
402 * to a distributed array with local face numbers.
404 * If a face in the list is represented on several domains,
405 * all the values are put in the array
407 void ParallelTopology::convertGlobalFaceListWithTwins(const int* face_list, int nbface, int*& local, int*& ip, int*& full_array,int& size)
410 for (int i=0; i< nbface; i++)
412 //int count = m_face_glob_to_loc.count(face_list[i]);
413 //if (count >1) MESSAGE_MED("face en doublon "<<face_list[i]);
414 size+= m_face_glob_to_loc.count(face_list[i]);
419 full_array=new int[size];
420 for (int i=0; i< nbface; i++)
422 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
423 pair<mmiter,mmiter> range=m_face_glob_to_loc.equal_range(face_list[i]);
424 for (mmiter it=range.first; it !=range.second; it++)
426 ip[index]=(it->second).first;
427 local[index]=(it->second).second;
428 full_array[index]=face_list[i];
435 //!converts a list of global cell numbers
436 //!to a distributed array with local cell numbers
437 void ParallelTopology::convertGlobalCellList(const int* cell_list, int nbcell, int* local, int* ip)
439 for (int i=0; i< nbcell; i++)
441 INTERP_KERNEL::HashMap<int, pair<int,int> >::const_iterator iter = m_glob_to_loc.find(cell_list[i]);
442 ip[i]=(iter->second).first;
443 local[i]=(iter->second).second;
447 /*!Converts a list of global face numbers
448 * to a distributed array with local face numbers
450 void ParallelTopology::convertGlobalFaceList(const int* face_list, int nbface, int* local, int* ip)
452 for (int i=0; i< nbface; i++)
454 INTERP_KERNEL::HashMap<int, pair<int,int> >::const_iterator iter = m_face_glob_to_loc.find(face_list[i]);
455 if (iter == m_face_glob_to_loc.end())
457 throw MED_EXCEPTION("convertGlobalFaceList - Face not found");
459 ip[i]=(iter->second).first;
460 local[i]=(iter->second).second;
461 // cout << " in convertGlobalFAceList face global "<<face_list[i]<<" -> ("<<ip[i]<<","<<local[i]<<")"<<endl;
465 /*!Converts a list of global node numbers on domain ip
466 * to a distributed array with local cell numbers.
468 * If a node in the list is represented on several domains,
469 * only the value with domain ip is returned
472 void ParallelTopology::convertGlobalFaceList(const int* face_list, int nbface, int* local, int ip)
474 for (int i=0; i< nbface; i++)
476 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
477 pair<mmiter,mmiter> range=m_face_glob_to_loc.equal_range(face_list[i]);
478 for (mmiter it=range.first; it !=range.second; it++)
480 int ipfound=(it->second).first;
482 local[i]=(it->second).second;
489 ////creating node mapping
490 void ParallelTopology::createNodeMapping(map<MED_EN::medGeometryElement,int*>& type_connectivity,
491 map<MED_EN::medGeometryElement,int>& present_type_numbers,
492 vector<int>& polygon_conn,
493 vector<int>& polygon_conn_index,
494 vector<int>& polyhedron_conn,
495 vector<int>& polyhedron_conn_index,
496 vector<int>& polyhedron_face_index,
499 set<int> local_numbers;
502 map<MED_EN::medGeometryElement,int>::const_iterator iter;
503 for (iter = present_type_numbers.begin(); iter!=present_type_numbers.end();iter++)
505 int type=iter->first;
506 int nodes_per_type= type%100;
508 if (!((type/100==m_mesh_dimension)
509 ||(type==MED_EN::MED_POLYGON && m_mesh_dimension==2)
510 ||(type==MED_EN::MED_POLYHEDRA && m_mesh_dimension==3))) continue;
512 if (type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
514 for (int icell=0; icell<present_type_numbers[type]; icell++)
516 for (int inode=0; inode<nodes_per_type; inode++)
518 int global=type_connectivity[type][icell*nodes_per_type+inode];
519 if(local_numbers.find(global)==local_numbers.end())
522 local_numbers.insert(global);
523 //m_node_loc_to_glob.insert(make_pair(make_pair(idomain,local_index),global));
524 m_node_loc_to_glob[idomain].push_back(global);
525 m_node_glob_to_loc.insert(make_pair(global,make_pair(idomain,local_index)));
526 // cout << "node : global ="<<global<<" local =("<<idomain<<","<<local_index<<")"<<endl;
532 else if (type== MED_EN::MED_POLYGON)
534 for ( unsigned i = 0; i < polygon_conn.size(); ++i )
536 int global=polygon_conn[i];
537 if(local_numbers.find(global)==local_numbers.end())
540 local_numbers.insert(global);
541 m_node_loc_to_glob[idomain].push_back(global);
542 m_node_glob_to_loc.insert(make_pair(global,make_pair(idomain,local_index)));
546 else if (type==MED_EN::MED_POLYHEDRA)
548 for ( unsigned i = 0; i < polyhedron_conn.size(); ++i )
550 int global=polyhedron_conn[i];
551 if(local_numbers.find(global)==local_numbers.end())
554 local_numbers.insert(global);
555 m_node_loc_to_glob[idomain].push_back(global);
556 m_node_glob_to_loc.insert(make_pair(global,make_pair(idomain,local_index)));
562 m_nb_nodes[idomain]=local_index;
565 //================================================================================
567 * \brief Return true if the domain mesh contains a cell based on given global nodes
569 //================================================================================
571 bool ParallelTopology::hasCellWithNodes( const MESHCollection& new_collection,
573 const set<int>& globNodes)
575 // convert global nodes to local in the given domain
577 set<int>::const_iterator n = globNodes.begin();
578 for ( ; n != globNodes.end(); ++n )
579 nodes.insert( convertGlobalNode( *n, domain ));
581 const MED_EN::medConnectivity connType = MED_EN::MED_NODAL;
582 const MED_EN::medEntityMesh entity = MED_EN::MED_CELL;
584 // loop on all types of cells
585 const MEDMEM::MESH* mesh = new_collection.getMesh( domain );
586 int nbTypes = mesh->getNumberOfTypes( entity );
587 const MED_EN::medGeometryElement * types = mesh->getTypes( entity );
588 for ( int t = 0; t < nbTypes; ++t )
591 if ( !mesh->existConnectivity( connType, entity ))
593 int nbCell = mesh->getNumberOfElements( entity, types[t] );
594 const int *conn, *index;
595 conn = mesh->getConnectivity(connType, entity, types[t]);
596 index = mesh->getConnectivityIndex(connType, entity);
597 // find a cell containing the first of given nodes,
598 // then check if the found cell contains all the given nodes
599 const int firstNode = *nodes.begin();
600 for ( int i = 0; i < nbCell; ++i )
602 for ( int j = index[i]-1; j < index[i+1]-1; ++j )
603 if ( conn[j] == firstNode )
606 for ( j = index[i]-1; j < index[i+1]-1; ++j )
607 nbSame += nodes.count( conn[j] );
608 if ( nbSame == nodes.size() )
617 ////creating face mapping
618 void ParallelTopology::createFaceMapping(const MESHCollection& initial_collection,
619 const MESHCollection& new_collection)
620 // map<MED_EN::medGeometryElement,int*>& type_list,
621 // map<MED_EN::medGeometryElement,int>& present_type_numbers,
625 // containers for the new topology
626 vector<int> new_counts(m_nb_domain,0);
627 vector<int> domain_counts(m_nb_domain,0);
628 const Topology* old_topology=initial_collection.getTopology();
629 int nb_domain_old=old_topology->nbDomain();
630 int global_index=old_topology->getFaceNumber();
631 //cout << "nb faces = " << global_index << endl;
632 set <pair<int, pair<int,int> > > global_treated;
634 //definition of the d-1 constituent for the considered mesh dimension
635 MED_EN::medEntityMesh constituent_entity;
636 switch (m_mesh_dimension)
639 constituent_entity= MED_EN::MED_FACE;
642 constituent_entity = MED_EN::MED_EDGE;
646 for (int iold=0; iold<nb_domain_old;iold++)
648 if ( !initial_collection.getMesh(iold) ) continue;
649 int nbtotalface = initial_collection.getMesh(iold)->getNumberOfElements(constituent_entity,MED_EN::MED_ALL_ELEMENTS);
650 SCRUTE_MED(nbtotalface);
651 const int* face_conn = 0;
652 const int* face_offset = 0;
655 face_conn = initial_collection.getMesh(iold)->getConnectivity(
656 MED_EN::MED_NODAL,constituent_entity,MED_EN::MED_ALL_ELEMENTS);
657 face_offset = initial_collection.getMesh(iold)->getConnectivityIndex(MED_EN::MED_NODAL,constituent_entity);
659 for (int iface=0;iface<nbtotalface; iface++)
661 int global_face_number = old_topology->convertFaceToGlobal(iold,iface+1);
663 // int inode = face_offset[iface];
664 for (int i=0; i<m_nb_domain; i++) domain_counts[i]=0;
668 nbnodes=face_offset[iface+1]-face_offset[iface];
669 for (int inode= face_offset[iface];inode < face_offset[iface+1]; inode++)
671 int node=face_conn[inode-1];
673 int global = old_topology->convertNodeToGlobal(iold,node);
674 // cout << "global node "<<global<<"ip "<<iold<< "noeud"<<node<<endl;
675 nodes.insert(global);
676 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
677 pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(global);
680 for (mmiter it=range.first; it !=range.second; it++)
682 ip=(it->second).first;
688 set<int>::const_iterator iter_node = nodes.begin();
690 numbers[2]=0; // for segments
691 for (int i=0; i<nbnodes; i++)
693 numbers[i]=*iter_node;
696 set <pair<int, pair<int,int> > > ::iterator iter_triplets;
697 pair<int, pair<int,int> > triplet = make_pair(numbers[0],make_pair(numbers[1],numbers[2]));
698 iter_triplets=global_treated.find(triplet);
699 if (iter_triplets==global_treated.end())
701 global_treated.insert(triplet);
702 // int nbnodes=face_offset[iface+1]-face_offset[iface];
703 if (global_face_number == -1)
706 global_face_number=global_index;
709 // SCRUTE_MED(nbnodes);
711 for (int inew=0;inew<m_nb_domain;inew++)
713 // SCRUTE_MED(domain_counts[inew]);
714 if(domain_counts[inew]==nbnodes)
716 if ( !hasCellWithNodes( new_collection, inew, nodes ))
717 continue; // 0020861: EDF 1387 MED: Result of medsplitter gives standalone triangles
720 m_face_glob_to_loc.insert(make_pair(global_face_number,make_pair(inew,new_counts[inew])));
721 //m_face_loc_to_glob.insert(make_pair(make_pair(inew,new_counts[inew]),global_face_number));
722 m_face_loc_to_glob[inew].push_back(global_face_number);
725 //cout << "global_face_number = " << global_face_number << endl;
730 for (int inew=0;inew<m_nb_domain;inew++)
732 m_nb_faces[inew]=new_counts[inew];
733 MESSAGE_MED(" Nb faces ["<<inew<<"]="<<m_nb_faces[inew]);
735 MESSAGE_MED(" total number of faces"<<getFaceNumber());
738 ////creating node mapping
739 void ParallelTopology::createFaceMapping2ndversion(const MESHCollection& initial_collection)
742 // containers for the new topology
743 vector<int> new_counts(m_nb_domain,0);
744 vector<int> domain_counts(m_nb_domain,0);
745 const Topology* old_topology=initial_collection.getTopology();
746 int nb_domain_old=old_topology->nbDomain();
747 //int global_index=old_topology->getFaceNumber();
748 // set <pair<int, pair<int,int> > > global_treated;
750 //definition of the d-1 constituent for the considered mesh dimension
751 MED_EN::medEntityMesh constituent_entity;
752 switch (m_mesh_dimension)
755 constituent_entity= MED_EN::MED_FACE;
758 constituent_entity = MED_EN::MED_EDGE;
762 for (int iold=0; iold<nb_domain_old;iold++)
764 int nbcell=old_topology->getCellNumber(iold);
766 const int* face_conn = initial_collection.getMesh(iold)->
767 getConnectivity(MED_EN::MED_DESCENDING,MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
768 const int* face_offset = initial_collection.getMesh(iold)->
769 getConnectivityIndex(MED_EN::MED_DESCENDING,MED_EN::MED_CELL);
770 MESSAGE_MED("end of connectivity calculation");
771 set<int> global_treated;
772 for (int icell=0; icell<nbcell; icell++)
774 int global_cell_number=old_topology->convertCellToGlobal(iold,icell+1);
775 int inew = getCellDomainNumber(global_cell_number);
777 for (int iface = face_offset[icell]; iface < face_offset[icell+1]; iface++)
779 int global_face_number=old_topology->convertFaceToGlobal(iold,abs(face_conn[iface-1]));
780 if (global_treated.find (global_face_number)==global_treated.end())
783 m_face_glob_to_loc.insert(make_pair(global_face_number,make_pair(inew,new_counts[inew])));
785 //m_face_loc_to_glob.insert(make_pair(make_pair(inew,new_counts[inew]),global_face_number));
786 m_face_loc_to_glob[inew].push_back(global_face_number);
787 global_treated.insert(global_face_number);
795 for (int inew=0;inew<m_nb_domain;inew++)
797 m_nb_faces[inew]=new_counts[inew];
798 // cout << " Nb faces ["<<inew<<"]="<<m_nb_faces[inew]<<endl;
800 MESSAGE_MED(" total number of faces"<<getFaceNumber());
804 //replacing a table of global numbering with a table with local numberings
805 // type_connectivity contains global connectivity for each type in input
806 // type_connectivity contains local connectivity for each type in output
807 void ParallelTopology::convertToLocal(map<MED_EN::medGeometryElement,int*>& type_connectivity,
808 map<MED_EN::medGeometryElement,int>& present_type_numbers,
810 MED_EN::medEntityMesh entity)
815 case MED_EN::MED_CELL:
816 dimension=m_mesh_dimension;
818 case MED_EN::MED_FACE:
821 case MED_EN::MED_EDGE:
826 MED_EN::MESH_ENTITIES::const_iterator currentEntity;
827 list<MED_EN::medGeometryElement>::const_iterator iter;
828 currentEntity = MED_EN::meshEntities.find(MED_EN::MED_CELL);
830 for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
832 MED_EN::medGeometryElement type = (*iter);
833 if (type/100 != dimension) continue;
834 for (int inode=0; inode<present_type_numbers[type]*(type%100); inode++)
836 // cout <<" inode :"<<inode<< " global = "<<type_connectivity[type][inode];
837 int global = type_connectivity[type][inode];
838 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
839 pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(global);
840 for (mmiter it=range.first; it !=range.second; it++)
842 if ((it->second).first==idomain)
843 type_connectivity[type][inode]=(it->second).second;
845 // cout << " local = " <<type_connectivity[type][inode]<<endl;
851 //replacing a table of global numbering with a table with local numberings
852 // type_connectivity contains global connectivity for each type in input
853 // type_connectivity contains local connectivity for each type in output
854 void ParallelTopology::convertToLocal2ndVersion(int* nodes, int nbnodes, int idomain)
856 for (int inode=0; inode<nbnodes; inode++)
858 // cout <<" inode :"<<inode<< " global = "<<type_connectivity[type][inode];
859 int global = nodes[inode];
860 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::iterator mmiter;
861 pair<mmiter,mmiter> range=m_node_glob_to_loc.equal_range(global);
862 for (mmiter it=range.first; it !=range.second; it++)
864 if ((it->second).first==idomain)
865 nodes[inode]=(it->second).second;
871 //! computing arrays with node/node correspondencies
872 void ParallelTopology::computeNodeNodeCorrespondencies(int idomain, vector<MEDMEM::MEDSKYLINEARRAY*>& corr) const
874 vector <int> sizes (m_nb_domain,0);
875 vector <int*> values (m_nb_domain);
876 for (int i=0; i<m_nb_domain; i++)
881 for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
883 //int global = (m_node_loc_to_glob.find(make_pair(idomain,inode+1)))->second;
884 int global = m_node_loc_to_glob[idomain][inode];
885 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::const_iterator mm;
886 pair<mm,mm> range = m_node_glob_to_loc.equal_range(global);
887 for (mm it=range.first; it !=range.second; it++)
889 int id = (it->second).first;
897 for (int ip=0; ip<m_nb_domain; ip++)
900 values[ip]=new int[2*sizes[ip]];
905 for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
907 //int global = (m_node_loc_to_glob.find(make_pair(idomain,inode+1)))->second;
908 int global = m_node_loc_to_glob[idomain][inode];
909 typedef INTERP_KERNEL::HashMultiMap<int,pair<int,int> >::const_iterator mm;
910 pair<mm,mm> range = m_node_glob_to_loc.equal_range(global);
911 for (mm it=range.first; it !=range.second; it++)
913 int id = (it->second).first;
916 values[id][sizes[id]*2]=inode+1;
917 values[id][sizes[id]*2+1]=(it->second).second;
923 for (int i=0; i< m_nb_domain; i++)
927 int* index = new int[sizes[i]+1];
928 for (int j=0; j<sizes[i]+1; j++)
930 corr[i]=new MEDMEM::MEDSKYLINEARRAY(sizes[i],2*sizes[i],index,values[i]);
937 //================================================================================
939 * \brief computing arrays with cell/cell correspondencies
940 * \param idomain - domain for which to compute correspondencies
941 * \param corr - out correspondencies
942 * \param graph - graph containing new cell distribution among domains
943 * \param id_shift - shitf to get cell id global on this processor from id global
946 //================================================================================
948 void ParallelTopology::computeCellCellCorrespondencies(int idomain, vector<MEDMEM::MEDSKYLINEARRAY*>& corr, const Graph* graph) const
950 vector <int> sizes (m_nb_domain,0);
951 vector <int*> values (m_nb_domain);
952 for (int i=0; i<m_nb_domain; i++)
957 vector <INTERP_KERNEL::HashMultiMap<int,int> > cell_corresp;
958 //TODO : remplacer int* par une map <int,int>
959 // vector<int* > number_of_connections(m_nb_domain);
960 // vector<map<int,int> > number_of_connections;
961 vector<map<int,int> > number_of_connections;
962 cell_corresp.resize(m_nb_domain);
963 number_of_connections.resize(m_nb_domain);
964 // for (int i=0; i<m_nb_domain; i++)
966 // // cell_corresp[i]=new multimap<int,int>;
967 // if (m_nb_cells[i] >0)
969 // number_of_connections[i]=new int[m_nb_cells[idomain]];
970 // for (int j=0; j<m_nb_cells[idomain]; j++)
971 // number_of_connections[i][j]=0;
975 const MEDMEM::MEDSKYLINEARRAY* skylinegraph = graph->getGraph();
977 const int* index=skylinegraph->getIndex();
978 const int* value=skylinegraph->getValue();
980 TGlob2DomainLoc::const_iterator gl_do_lo_end = m_glob_to_loc.end();
982 for (int icell=0; icell<m_nb_cells[idomain]; icell++)
984 int global= m_loc_to_glob[idomain][icell];
985 for (int ii=index[global-1]-1; ii<index[global]-1; ii++)
987 int distant_global=value[ii];
989 const pair<int,int>& local = m_glob_to_loc.find(distant_global)->second;
991 if (local.first != idomain)
993 cell_corresp[local.first].insert(make_pair(icell+1,local.second));
994 // number_of_connections[local.first][icell]++;
995 if (number_of_connections[local.first].find(icell)==number_of_connections[local.first].end())
996 number_of_connections[local.first].insert(make_pair(icell,1));
998 number_of_connections[local.first][icell]++;
1004 for (int inew=0; inew<m_nb_domain; inew++)
1006 if (inew==idomain || number_of_connections[inew].empty()) continue;
1008 int* new_index=new int[m_nb_cells[idomain]+1];
1010 for (int i=0; i<m_nb_cells[idomain]; i++)
1013 if (number_of_connections[inew].find(i)!=number_of_connections[inew].end())
1014 new_index[i+1]=new_index[i]+number_of_connections[inew][i];
1016 new_index[i+1]=new_index[i];
1019 if (new_index[m_nb_cells[idomain]]-1 > 0)
1020 new_value=new int[new_index[m_nb_cells[idomain]]-1];
1026 // INTERP_KERNEL::HashMultiMap<int,int>::iterator iter=cell_corresp[inew].begin();
1028 for (int i=0; i<m_nb_cells[idomain]; i++)
1030 // for (int j=new_index[i];j<new_index[i+1];j++,value_i++,iter++)
1031 // new_value[value_i]=iter->second;
1033 typedef INTERP_KERNEL::HashMultiMap<int,int>::iterator mmiter;
1034 pair<mmiter,mmiter> range=cell_corresp[inew].equal_range(i+1);
1035 for (mmiter it=range.first; it!=range.second; it++)
1037 new_value[value_i]=it->second;
1043 corr[inew] = new MEDMEM::MEDSKYLINEARRAY(m_nb_cells[idomain], new_index[m_nb_cells[idomain]]-1, new_index,new_value,true);
1048 if (new_value !=0) delete[]new_value;
1054 // for (int inew=0; inew<m_nb_domain; inew++)
1055 // if (m_nb_cells[inew]>0)
1056 // delete[] number_of_connections[inew];
1060 //================================================================================
1062 * \brief Return max global face number
1064 //================================================================================
1066 int ParallelTopology::getMaxGlobalFace() const
1069 TGlob2LocsMap::const_iterator g_l_l = m_face_glob_to_loc.begin();
1070 for ( ; g_l_l != m_face_glob_to_loc.end(); ++g_l_l )
1071 if ( g_l_l->first > max )
1076 void ParallelTopology::recreateFaceMapping(const TGeom2FacesByDomian& face_map)
1078 m_face_glob_to_loc.clear();
1079 for (int i=0; i<m_nb_domain;i++)
1080 m_face_loc_to_glob[i].clear();
1082 for (int idomain=0; idomain<m_nb_domain; idomain++)
1085 TGeom2Faces::const_iterator iter = face_map[idomain].begin();
1086 for (; iter != face_map[idomain].end(); iter++)
1088 for (unsigned i=0; i< (iter->second).size(); i++)
1090 MEDSPLITTER_FaceModel* face = (iter->second)[i];
1091 //cout << "global :"<<face->getGlobal()<<" - "<<ilocal<<endl;
1092 m_face_glob_to_loc.insert(make_pair(face->getGlobal(),make_pair(idomain,ilocal)));
1093 m_face_loc_to_glob[idomain].push_back(face->getGlobal());
1097 m_nb_faces[idomain] =ilocal-1;
1101 //================================================================================
1103 * \brief Recreating cell and node mapping after send-reveive and fusion of domain meshes
1105 //================================================================================
1107 void ParallelTopology::recreateMappingAfterFusion(const vector<MEDMEM::MESH*>& meshes)
1109 const char* LOC = "ParallelTopology::recreateMappingAfterFusion(): ";
1111 m_glob_to_loc.clear();
1112 m_node_glob_to_loc.clear();
1113 m_face_glob_to_loc.clear();
1115 for (int idomain=0; idomain<m_nb_domain; idomain++)
1117 m_nb_cells[idomain] = m_nb_nodes[idomain] = m_nb_faces[idomain] = 0;
1118 m_loc_to_glob[idomain].clear();
1119 m_node_loc_to_glob[idomain].clear();
1120 m_face_loc_to_glob[idomain].clear();
1122 if ( !meshes[idomain]->getCoordinateptr() ) continue; // empty domian
1124 //creating cell maps
1126 m_nb_cells[idomain]=meshes[idomain]->getNumberOfElements(MED_EN::MED_CELL,
1127 MED_EN::MED_ALL_ELEMENTS);
1128 if ( m_cell_loc_to_glob_fuse[idomain].size() != m_nb_cells[idomain] )
1129 throw MED_EXCEPTION(MEDMEM::STRING(LOC)<<" invalid nb of fused cells");
1131 m_loc_to_glob[idomain].swap( m_cell_loc_to_glob_fuse[idomain] );
1133 for (int i=0; i<m_nb_cells[idomain]; i++)
1135 int global=m_loc_to_glob[idomain][i];
1136 m_glob_to_loc[global]=make_pair(idomain,i+1);
1139 //creating node maps
1141 m_nb_nodes[idomain]=meshes[idomain]->getNumberOfNodes();
1142 m_node_loc_to_glob[idomain] = ((MEDMEM::MeshFuse*)meshes[idomain])->getNodeNumbers();
1143 if ( m_node_loc_to_glob[idomain].size() != m_nb_nodes[idomain] )
1144 throw MED_EXCEPTION(MEDMEM::STRING(LOC)<<" invalid nb of fused nodes");
1146 // setting mappings for all nodes
1147 for (int inode=0; inode<m_nb_nodes[idomain]; inode++)
1149 int global_number=m_node_loc_to_glob[idomain][inode];
1150 m_node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode+1)));
1154 //creating dimension d-1 component mappings
1156 MED_EN::medEntityMesh constituent_entity =
1157 (m_mesh_dimension == 3 ? MED_EN::MED_FACE : MED_EN::MED_EDGE );
1158 m_nb_faces[idomain] = meshes[idomain]->getNumberOfElements(constituent_entity,
1159 MED_EN::MED_ALL_ELEMENTS);
1160 if ( m_face_loc_to_glob_fuse[idomain].size() != m_nb_faces[idomain] )
1161 throw MED_EXCEPTION(MEDMEM::STRING(LOC)<<" invalid nb of fused faces of domain "<< idomain
1162 << ": expect " << m_nb_faces[idomain]
1163 << " but have " << m_face_loc_to_glob_fuse[idomain].size());
1165 m_face_loc_to_glob[idomain].swap( m_face_loc_to_glob_fuse[idomain] );
1167 for (int iface=0; iface<m_nb_faces[idomain]; iface++)
1169 int global_number=m_face_loc_to_glob[idomain][iface];
1170 m_face_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,iface+1)));