1 // Copyright (C) 2007-2024 CEA, EDF
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, or (at your option) any later version.
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 "MEDPARTITIONER_MeshCollection.hxx"
22 #include "MEDPARTITIONER_ConnectZone.hxx"
23 #include "MEDPARTITIONER_Graph.hxx"
24 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
25 #include "MEDPARTITIONER_MeshCollectionMedAsciiDriver.hxx"
26 #include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx"
27 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
28 #include "MEDPARTITIONER_ParallelTopology.hxx"
29 #include "MEDPARTITIONER_Topology.hxx"
30 #include "MEDPARTITIONER_UserGraph.hxx"
31 #include "MEDPARTITIONER_Utils.hxx"
34 #include "MEDPARTITIONER_JointFinder.hxx"
38 #include "MEDCouplingFieldDouble.hxx"
39 #include "MEDCouplingMemArray.hxx"
40 #include "MEDCouplingNormalizedUnstructuredMesh.hxx"
41 #include "MEDCouplingSkyLineArray.hxx"
42 #include "MEDCouplingUMesh.hxx"
43 #include "MEDLoader.hxx"
44 #include "MEDLoaderBase.hxx"
50 #ifdef MED_ENABLE_PARMETIS
51 #include "MEDPARTITIONER_ParMetisGraph.hxx"
53 #ifdef MED_ENABLE_METIS
54 #include "MEDPARTITIONER_MetisGraph.hxx"
56 #ifdef MED_ENABLE_SCOTCH
57 #include "MEDPARTITIONER_ScotchGraph.hxx"
67 MEDPARTITIONER::MeshCollection::MeshCollection()
69 _owns_topology(false),
71 _domain_selector( 0 ),
72 _i_non_empty_mesh(-1),
73 _driver_type(MEDPARTITIONER::MedXml),
74 _subdomain_boundary_creates( MyGlobals::_Create_Boundary_Faces ),
75 _family_splitting(false),
76 _create_empty_groups(false),
81 /*!constructor creating a new mesh collection (mesh series + topology)
82 *from an old collection and a new topology
84 * On output, the constructor has built the meshes corresponding to the new mesh collection.
85 * The new topology has been updated so that face and node mappings are included.
86 * The families have been cast to their projections in the new topology.
88 * \param initial_collection collection from which the data (coordinates, connectivity) are taken
89 * \param topology topology containing the cell mappings
92 MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection,
94 bool family_splitting,
95 bool create_empty_groups)
96 : _topology(topology),
97 _owns_topology(false),
99 _domain_selector( initialCollection._domain_selector ),
100 _i_non_empty_mesh(-1),
101 _name(initialCollection._name),
102 _driver_type(MEDPARTITIONER::MedXml),
103 _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
104 _family_splitting(family_splitting),
105 _create_empty_groups(create_empty_groups),
108 std::vector<std::vector<std::vector<mcIdType> > > new2oldIds(initialCollection.getTopology()->nbDomain());
109 std::vector<MEDCoupling::DataArrayIdType*> o2nRenumber;
111 castCellMeshes(initialCollection, new2oldIds, o2nRenumber );
113 //defining the name for the collection and the underlying meshes
114 setName(initialCollection.getName());
121 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
122 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
124 if (MyGlobals::_Is0verbose)
125 std::cout<<"treating faces"<<std::endl;
126 NodeMapping nodeMapping;
127 //nodeMapping contains the mapping between old nodes and new nodes
128 // (iolddomain,ioldnode)->(inewdomain,inewnode)
129 createNodeMapping(initialCollection, nodeMapping);
130 std::vector<std::vector<std::vector<mcIdType> > > new2oldFaceIds;
131 castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
137 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
138 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
140 if (MyGlobals::_Is0verbose)
142 if (isParallelMode())
143 std::cout << "ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
145 std::cout << "NOT ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
147 if (MyGlobals::_Is0verbose>10)
148 std::cout<<"treating cell and face families"<<std::endl;
150 castIntField(initialCollection.getMesh(),
152 initialCollection.getCellFamilyIds(),
154 castIntField(initialCollection.getFaceMesh(),
156 initialCollection.getFaceFamilyIds(),
161 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
162 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
164 if (MyGlobals::_Is0verbose)
165 std::cout << "treating groups" << std::endl;
166 _family_info=initialCollection.getFamilyInfo();
167 _group_info=initialCollection.getGroupInfo();
170 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
171 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
173 if (MyGlobals::_Is0verbose)
174 std::cout << "treating fields" << std::endl;
175 castAllFields(initialCollection,"cellFieldDouble");
176 if (_i_non_empty_mesh<0)
178 for (int i=0; i<getNbOfGlobalMeshes(); i++)
182 _i_non_empty_mesh=i; //first existing one local
188 // find faces common with neighbor domains and put them in groups
189 buildBoundaryFaces();
191 //building the connect zones necessary for writing joints
192 buildConnectZones( nodeMapping, o2nRenumber, initialCollection.getTopology()->nbDomain() );
194 // delete o2nRenumber
195 for ( size_t i = 0; i < o2nRenumber.size(); ++i )
196 if ( o2nRenumber[i] )
197 o2nRenumber[i]->decrRef();
201 Creates the meshes using the topology underlying he mesh collection and the mesh data
202 coming from the ancient collection
203 \param initialCollection collection from which the data is extracted to create the new meshes
204 \param [out] o2nRenumber returns for each new domain a permutation array returned by sortCellsInMEDFileFrmt()
207 void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
208 std::vector<std::vector<std::vector<mcIdType> > >& new2oldIds,
209 std::vector<MEDCoupling::DataArrayIdType*> & o2nRenumber)
211 if (MyGlobals::_Verbose>10)
212 std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
214 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes");
216 int nbNewDomain=_topology->nbDomain();
217 int nbOldDomain=initialCollection.getTopology()->nbDomain();
219 _mesh.resize(nbNewDomain);
220 o2nRenumber.resize(nbNewDomain,0);
221 int rank=MyGlobals::_Rank;
222 //splitting the initial domains into smaller bits
223 std::vector<std::vector<MEDCoupling::MEDCouplingUMesh*> > splitMeshes;
224 splitMeshes.resize(nbNewDomain);
225 for (int inew=0; inew<nbNewDomain; inew++)
227 splitMeshes[inew].resize(nbOldDomain, (MEDCoupling::MEDCouplingUMesh*)0);
230 for (int iold=0; iold<nbOldDomain; iold++)
232 if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
234 mcIdType size=(initialCollection._mesh)[iold]->getNumberOfCells();
235 std::vector<mcIdType> globalids(size);
236 initialCollection.getTopology()->getCellList(iold, &globalids[0]);
237 std::vector<mcIdType> ilocalnew(size); //local
238 std::vector<int> ipnew(size); //idomain old
239 _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
241 new2oldIds[iold].resize(nbNewDomain);
242 for (int i=0; i<(int)ilocalnew.size(); i++)
244 new2oldIds[iold][ipnew[i]].push_back(i);
246 for (int inew=0; inew<nbNewDomain; inew++)
248 splitMeshes[inew][iold]=(MEDCoupling::MEDCouplingUMesh*)
249 (initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
250 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
252 if (MyGlobals::_Verbose>400)
253 std::cout<< "proc " << rank << " : a splitMesh iold inew NbCells " << iold << " " << inew << " "
254 << splitMeshes[inew][iold]->getNumberOfCells() << std::endl;
259 if (isParallelMode())
261 //if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
262 for (int iold=0; iold<nbOldDomain; iold++)
263 for(int inew=0; inew<nbNewDomain; inew++)
265 if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
268 if(initialCollection._domain_selector->isMyDomain(iold))
269 _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew));
271 if (_domain_selector->isMyDomain(inew))
272 _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
278 //fusing the split meshes
279 if (MyGlobals::_Verbose>200)
280 std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
281 for (int inew=0; inew<nbNewDomain ;inew++)
283 std::vector<const MEDCoupling::MEDCouplingUMesh*> meshes;
285 for (int i=0; i<(int)splitMeshes[inew].size(); i++)
286 if (splitMeshes[inew][i]!=0)
287 if (splitMeshes[inew][i]->getNumberOfCells()>0)
288 meshes.push_back(splitMeshes[inew][i]);
290 if (!isParallelMode()||_domain_selector->isMyDomain(inew))
292 if (meshes.size()==0)
294 _mesh[inew]=CreateEmptyMEDCouplingUMesh();
295 std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
299 _mesh[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshes(meshes);
300 o2nRenumber[inew]=_mesh[inew]->sortCellsInMEDFileFrmt();
302 mcIdType nbNodesMerged;
305 MEDCoupling::DataArrayIdType* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
306 array->decrRef(); // array is not used in this case
308 _mesh[inew]->zipCoords();
311 for (int i=0;i<(int)splitMeshes[inew].size();i++)
312 if (splitMeshes[inew][i]!=0)
313 splitMeshes[inew][i]->decrRef();
315 if (MyGlobals::_Verbose>300)
316 std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
320 \param initialCollection source mesh collection
321 \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection
323 void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
326 using std::make_pair;
327 // NodeMapping reverseNodeMapping;
328 for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
332 BBTreeOfDim* tree = 0;
334 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
336 // std::map<pair<double,pair<double, double> >, int > nodeClassifier;
337 MEDCoupling::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
338 double* coordsPtr=coords->getPointer();
339 dim = (int)coords->getNumberOfComponents();
340 mcIdType nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
341 bbox=new double[nvertices*2*dim];
343 for (int i=0; i<nvertices*dim;i++)
345 bbox[i*2]=coordsPtr[i]-1e-8;
346 bbox[i*2+1]=coordsPtr[i]+1e-8;
348 tree=new BBTreeOfDim( dim, bbox,0,0,nvertices,1e-9);
351 for (int inew=0; inew<_topology->nbDomain(); inew++)
354 //sending meshes for parallel computation
355 if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))
356 _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
357 else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
359 MEDCoupling::MEDCouplingUMesh* mesh;
360 _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
361 MEDCoupling::DataArrayDouble* coords = mesh->getCoords();
362 for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
364 double* coordsPtr=coords->getPointer()+inode*dim;
365 vector<mcIdType> elems;
366 tree->getElementsAroundPoint(coordsPtr,elems);
367 if (elems.size()==0) continue;
368 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
372 else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
374 if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
377 MEDCoupling::DataArrayDouble* coords = getMesh(inew)->getCoords();
378 for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
380 double* coordsPtr=coords->getPointer()+inode*dim;
381 vector<mcIdType> elems;
382 tree->getElementsAroundPoint(coordsPtr,elems);
383 if (elems.size()==0) continue;
384 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
388 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
397 void getNodeIds(MEDCoupling::MEDCouplingUMesh& meshOne, MEDCoupling::MEDCouplingUMesh& meshTwo, std::vector<mcIdType>& nodeIds)
400 using MEDPARTITIONER::BBTreeOfDim;
401 //if (!&meshOne || !&meshTwo) return; //empty or not existing
403 BBTreeOfDim* tree = 0;
404 mcIdType nv1=meshOne.getNumberOfNodes();
405 MEDCoupling::DataArrayDouble* coords=meshOne.getCoords();
406 int dim = (int)coords->getNumberOfComponents();
408 bbox=new double[nv1*2*dim];
409 double* coordsPtr=coords->getPointer();
410 for (int i=0; i<nv1*dim; i++)
412 bbox[i*2]=coordsPtr[i]-1e-8;
413 bbox[i*2+1]=coordsPtr[i]+1e-8;
415 tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9);
417 mcIdType nv2=meshTwo.getNumberOfNodes();
418 nodeIds.resize(nv2,-1);
419 coords=meshTwo.getCoords();
420 for (int inode=0; inode<nv2; inode++)
422 double* coordsPtr2=coords->getPointer()+inode*dim;
423 vector<mcIdType> elems;
424 tree->getElementsAroundPoint(coordsPtr2,elems);
425 if (elems.size()==0) continue;
426 nodeIds[inode]=elems[0];
433 creates the face meshes on the new domains from the faces on the old domain and the node mapping
434 faces at the interface are duplicated
436 void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
437 const std::multimap<std::pair<int,mcIdType>, std::pair<int,mcIdType> >& nodeMapping,
438 std::vector<std::vector<std::vector<mcIdType> > >& new2oldIds)
440 //splitMeshes structure will contain the partition of
441 //the old faces on the new ones
442 //splitMeshes[4][2] contains the faces from old domain 2
443 //that have to be added to domain 4
449 using std::make_pair;
451 if (MyGlobals::_Verbose>10)
452 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
454 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
456 int nbNewDomain=_topology->nbDomain();
457 int nbOldDomain=initialCollection.getTopology()->nbDomain();
459 vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
460 vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
462 vector< vector<MEDCoupling::MEDCouplingUMesh*> > splitMeshes;
464 splitMeshes.resize(nbNewDomain);
465 for (int inew=0; inew<nbNewDomain; inew++)
467 splitMeshes[inew].resize(nbOldDomain, (MEDCoupling::MEDCouplingUMesh*)0);
469 new2oldIds.resize(nbOldDomain);
470 for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
472 //init null pointer for empty meshes
473 for (int inew=0; inew<nbNewDomain; inew++)
475 for (int iold=0; iold<nbOldDomain; iold++)
477 splitMeshes[inew][iold]=0; //null for empty meshes
478 new2oldIds[iold][inew].clear();
482 //loop over the old domains to analyse the faces and decide
483 //on which new domain they belong
484 for (int iold=0; iold<nbOldDomain; iold++)
486 if (isParallelMode() && !initialCollection._domain_selector->isMyDomain(iold)) continue;
487 if (MyGlobals::_Verbose>400)
488 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes iodDomain "<<iold<<std::endl;
489 //initial face mesh known : in my domain
490 if (meshesCastFrom[iold] != 0)
492 for (mcIdType ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
494 vector<mcIdType> nodes;
495 meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
497 map <int,mcIdType> faces;
499 //analysis of element ielem
500 //counters are set for the element
501 //for each source node, the mapping is interrogated and the domain counters
502 //are incremented for each target node
503 //the face is considered as going to target domains if the counter of the domain
504 //is equal to the number of nodes
505 for (std::size_t inode=0; inode<nodes.size(); inode++)
507 typedef multimap<pair<int,mcIdType>,pair<int,mcIdType> >::const_iterator MI;
508 mcIdType mynode=nodes[inode];
510 pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
511 for (MI iter=myRange.first; iter!=myRange.second; iter++)
513 int inew=iter->second.first;
514 if (faces.find(inew)==faces.end())
521 for (map<int,mcIdType>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
523 if (iter->second==ToIdType(nodes.size()))
524 //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
525 //it is not sure here...
526 //done before writeMedfile on option?... see filterFaceOnCell()
527 new2oldIds[iold][iter->first].push_back(ielem);
531 //creating the splitMeshes from the face ids
532 for (int inew=0; inew<nbNewDomain; inew++)
534 if (meshesCastFrom[iold]->getNumberOfCells() > 0)
536 splitMeshes[inew][iold]=
537 (MEDCoupling::MEDCouplingUMesh*)
538 ( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
539 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
542 splitMeshes[inew][iold]->zipCoords();
546 splitMeshes[inew][iold]=CreateEmptyMEDCouplingUMesh();
552 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
558 if (isParallelMode())
560 MEDCoupling::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
561 for (int iold=0; iold<nbOldDomain; iold++)
562 for (int inew=0; inew<nbNewDomain; inew++)
564 if (initialCollection._domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
566 if (splitMeshes[inew][iold] != 0)
568 _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew));
569 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes send " <<inew<<" "<<iold<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
573 _domain_selector->sendMesh(*(empty), _domain_selector->getProcessorID(inew));
574 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes sen0 " <<inew<<" "<<iold<<std::endl;
577 if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
578 _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
580 //if (splitMeshes[inew][iold])
581 // nb=splitMeshes[inew][iold]->getNumberOfCells();
582 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
588 //fusing the split meshes
589 if (MyGlobals::_Verbose>200)
590 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes fusing" << std::endl;
591 meshesCastTo.resize(nbNewDomain);
592 for (int inew=0; inew<nbNewDomain; inew++)
594 vector<MEDCoupling::MEDCouplingUMesh*> myMeshes;
595 for (int iold=0; iold<nbOldDomain; iold++)
597 MEDCoupling::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
599 if (umesh->getNumberOfCells()>0)
600 myMeshes.push_back(umesh);
604 MEDCoupling::MEDCouplingUMesh *bndMesh = 0;
605 if ( _subdomain_boundary_creates &&
607 _mesh[inew]->getNumberOfCells()>0 )
610 ((MEDCoupling::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
611 if (bndMesh->getNumberOfCells()>0)
612 myMeshes.push_back( bndMesh );
615 if (myMeshes.size()>0)
617 // Need to merge faces that are both in the initial set of faces and in the boundary mesh
618 // which is the last element in myMeshes:
621 for (int iold2 = 0; iold2 < (int)myMeshes.size()-1; iold2++)
622 myMeshes[iold2]->tryToShareSameCoordsPermute(*bndMesh, 1.e-10);
623 vector<const MEDCoupling::MEDCouplingUMesh*> myMeshes_c;
624 for (auto & mesh: myMeshes) myMeshes_c.push_back(mesh);
625 meshesCastTo[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshesOnSameCoords(myMeshes_c);
626 MCAuto<DataArrayIdType> tmp = meshesCastTo[inew]->zipConnectivityTraducer(2);
630 vector<const MEDCoupling::MEDCouplingUMesh*> myMeshes_c;
631 for (auto & mesh: myMeshes) myMeshes_c.push_back(mesh);
632 meshesCastTo[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshes(myMeshes_c);
634 meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
638 MEDCoupling::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
639 meshesCastTo[inew]=empty;
641 for (int iold=0; iold<nbOldDomain; iold++)
642 if (splitMeshes[inew][iold]!=0)
643 splitMeshes[inew][iold]->decrRef();
647 if (MyGlobals::_Verbose>300)
648 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
653 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom,
654 std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo,
655 std::vector<MEDCoupling::DataArrayIdType*>& arrayFrom,
656 std::string nameArrayTo)
660 std::size_t ioldMax=meshesCastFrom.size();
661 std::size_t inewMax=meshesCastTo.size();
664 //preparing bounding box trees for accelerating source-target node identifications
665 if (MyGlobals::_Verbose>99)
666 std::cout<<"making accelerating structures"<<std::endl;
667 std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
668 std::vector<MEDCoupling::DataArrayDouble*>bbox(ioldMax);
669 for (unsigned int iold =0; iold< ioldMax; iold++)
670 if (isParallelMode() && _domain_selector->isMyDomain(iold))
672 MEDCoupling::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->computeCellCenterOfMass();
673 bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
674 acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
675 sourceCoords->decrRef();
677 // send-recv operations
679 for (unsigned int inew=0; inew<inewMax; inew++)
681 for (unsigned int iold=0; iold<ioldMax; iold++)
683 //sending arrays for distant domains
684 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
687 _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
689 mcIdType size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
690 vector<mcIdType> sendIds;
691 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
692 if (size>0) //no empty
694 sendIds.resize(size);
695 std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
700 sendIds.resize(size);
702 SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
704 //receiving arrays from distant domains
705 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
708 vector<mcIdType> recvIds;
709 MEDCoupling::MEDCouplingUMesh* recvMesh;
710 _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
712 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
713 RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
714 remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
715 recvMesh->decrRef(); //cww is it correct?
721 //localc ontributions and aggregation
722 for (unsigned int inew=0; inew<inewMax; inew++)
724 for (unsigned int iold=0; iold<ioldMax; iold++)
725 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
727 remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
730 for (unsigned int iold=0; iold<ioldMax;iold++)
731 if (isParallelMode() && _domain_selector->isMyDomain(iold))
733 bbox[iold]->decrRef();
734 delete acceleratingStructures[iold];
738 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
739 const MEDCoupling::MEDCouplingUMesh& sourceMesh,
740 const MEDCoupling::MEDCouplingUMesh& targetMesh,
741 const mcIdType* fromArray,
742 std::string nameArrayTo,
743 const BBTreeOfDim* myTree)
746 if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
747 MEDCoupling::DataArrayDouble* targetCoords=targetMesh.computeCellCenterOfMass();
748 const double* tc=targetCoords->getConstPointer();
749 mcIdType targetSize=targetMesh.getNumberOfCells();
750 mcIdType sourceSize=sourceMesh.getNumberOfCells();
751 if (MyGlobals::_Verbose>200)
752 std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
753 std::vector<mcIdType> ccI;
755 str=nameArrayTo+"_toArray";
756 cle=Cle1ToStr(str,inew);
759 const BBTreeOfDim* tree;
760 bool cleantree=false;
761 MEDCoupling::DataArrayDouble* sourceBBox=0;
762 int dim = (int)targetCoords->getNumberOfComponents();
765 sourceBBox=sourceMesh.computeCellCenterOfMass()->computeBBoxPerTuple(1e-8);
766 tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
770 //first time iold : create and initiate
771 if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
773 if (MyGlobals::_Is0verbose>100)
774 std::cout << "create " << cle << " size " << targetSize << std::endl;
775 MEDCoupling::DataArrayIdType* p=MEDCoupling::DataArrayIdType::New();
776 p->alloc(targetSize,1);
778 toArray=p->getPointer();
779 _map_dataarray_int[cle]=p;
781 else //other times iold: refind and complete
783 toArray=_map_dataarray_int.find(cle)->second->getPointer();
786 std::map< mcIdType, int > isource2nb; // count coincident elements
787 std::map< mcIdType, int>::iterator i2nb;
789 for (mcIdType itargetnode=0; itargetnode<targetSize; itargetnode++)
791 std::vector<mcIdType> intersectingElems;
792 tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
793 if (intersectingElems.size()!=0)
795 mcIdType isourcenode=intersectingElems[0];
796 if ( intersectingElems.size() > 1 )
798 i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
799 isourcenode = intersectingElems[ i2nb->second++ ];
801 if ( isourcenode < sourceSize ) // protection from invalid elements
803 toArray[itargetnode]=fromArray[isourcenode];
804 ccI.push_back(itargetnode);
805 ccI.push_back(isourcenode);
809 if (MyGlobals::_Verbose>200)
810 std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
811 //memories intersection for future same job on fields (if no existing cle=no intersection)
812 str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
813 if (MyGlobals::_Verbose>700)
814 std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
816 _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
818 targetCoords->decrRef();
819 if (cleantree) delete tree;
820 if (sourceBBox !=0) sourceBBox->decrRef();
823 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
825 if (nameArrayTo!="cellFieldDouble")
826 throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
828 std::string nameTo="typeData=6"; //resume the type of field casted
829 // send-recv operations
830 std::size_t ioldMax=initialCollection.getMesh().size();
831 std::size_t inewMax=this->getMesh().size();
832 std::size_t iFieldMax=initialCollection.getFieldDescriptions().size();
833 if (MyGlobals::_Verbose>10)
834 std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
835 //see collection.prepareFieldDescriptions()
836 for (std::size_t ifield=0; ifield<iFieldMax; ifield++)
838 std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
839 if (descriptionField.find(nameTo)==std::string::npos)
840 continue; //only nameTo accepted in Fields name description
842 for (unsigned int inew=0; inew<inewMax; inew++)
844 for (unsigned int iold=0; iold<ioldMax; iold++)
846 //sending arrays for distant domains
847 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
849 int target=_domain_selector->getProcessorID(inew);
850 MEDCoupling::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
851 if (MyGlobals::_Verbose>10)
852 std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
853 SendDataArrayDouble(field, target);
855 //receiving arrays from distant domains
856 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
858 int source=_domain_selector->getProcessorID(iold);
860 if (MyGlobals::_Verbose>10)
861 std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
862 MEDCoupling::DataArrayDouble* field=RecvDataArrayDouble(source);
863 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
868 //local contributions and aggregation
869 for (unsigned int inew=0; inew<inewMax; inew++)
871 for (unsigned int iold=0; iold<ioldMax; iold++)
872 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
874 MEDCoupling::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
875 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
881 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
882 MEDCoupling::DataArrayDouble* fromArray,
883 std::string nameArrayTo,
884 std::string descriptionField)
885 //here we use 'cellFamily_ccI inew iold' set in remapIntField
887 if (nameArrayTo!="cellFieldDouble")
888 throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
889 std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
891 std::map<std::string,MEDCoupling::DataArrayIdType*>::iterator it1;
892 it1=_map_dataarray_int.find(key);
893 if (it1==_map_dataarray_int.end())
895 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
896 std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
899 //create ccI in remapIntField
900 MEDCoupling::DataArrayIdType *ccI=it1->second;
901 if (MyGlobals::_Verbose>300)
902 std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
904 mcIdType nbcell=this->getMesh()[inew]->getNumberOfCells();
905 std::size_t nbcomp=fromArray->getNumberOfComponents();
906 int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
908 std::string tag="inewFieldDouble="+IntToStr(inew);
909 key=descriptionField+SerializeFromString(tag);
910 mcIdType fromArrayNbOfElem=fromArray->getNbOfElems();
911 mcIdType fromArrayNbOfComp=ToIdType(fromArray->getNumberOfComponents());
912 mcIdType fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
914 if (MyGlobals::_Verbose>1000)
916 std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
917 " fromArray nbOfElems " << fromArrayNbOfElem <<
918 " nbTuples " << fromArray->getNumberOfTuples() <<
919 " nbcells " << fromArrayNbOfCell <<
920 " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
923 MEDCoupling::DataArrayDouble* field=0;
924 std::map<std::string,MEDCoupling::DataArrayDouble*>::iterator it2;
925 it2=_map_dataarray_double.find(key);
926 if (it2==_map_dataarray_double.end())
928 if (MyGlobals::_Verbose>300)
929 std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
930 field=MEDCoupling::DataArrayDouble::New();
931 _map_dataarray_double[key]=field;
932 field->alloc(nbcell*nbPtGauss,nbcomp);
933 field->fillWithZero();
938 if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
940 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
941 " trying remap of field double on cells : \n" << descriptionField << std::endl;
948 field->setPartOfValuesAdv(fromArray,ccI);
952 //replaced by setPartOfValuesAdv if nbPtGauss==1
953 mcIdType iMax=ccI->getNbOfElems();
954 mcIdType* pccI=ccI->getPointer();
955 double* pField=field->getPointer();
956 double* pFrom=fromArray->getPointer();
957 mcIdType itarget, isource, delta=ToIdType(nbPtGauss*nbcomp);
958 for (mcIdType i=0; i<iMax; i=i+2) //cell
962 if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
963 throw INTERP_KERNEL::Exception("Error field override");
964 mcIdType ita=itarget*delta;
965 mcIdType iso=isource*delta;
966 for (mcIdType k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
973 using namespace MEDCoupling;
974 //================================================================================
976 * \brief Sort correspondence ids of one domain and permute ids of the other accordingly
977 * \param [in,out] ids1 - ids of one domain
978 * \param [in,out] ids2 - ids of another domain
979 * \param [in] delta - a delta to change all ids
980 * \param [in] removeEqual - to remove equal ids
981 * \return DataArrayInt* - array of ids joined back
983 //================================================================================
985 DataArrayIdType* sortCorrespondences( DataArrayIdType* ids1,
986 DataArrayIdType* ids2,
988 bool removeEqual = false)
991 MCAuto< DataArrayIdType > renumN2O = ids1->buildPermArrPerLevel();
992 ids1->renumberInPlaceR( renumN2O->begin() );
993 ids2->renumberInPlaceR( renumN2O->begin() );
997 ids1 = ids1->buildUnique();
998 ids2 = ids2->buildUnique();
1002 mcIdType * id = ids1->getPointer();
1003 for ( ; id < ids1->end(); ++id )
1005 id = ids2->getPointer();
1006 for ( ; id < ids2->end(); ++id )
1011 DataArrayIdType* ids12 = DataArrayIdType::Meld( ids1, ids2 ); // two components
1012 ids12->rearrange( 1 ); // make one component
1016 //================================================================================
1018 * \brief Renumber ids according to mesh->sortCellsInMEDFileFrmt()
1019 * \param [in,out] ids - cell ids to renumber
1020 * \param [in] o2nRenumber - renumbering array in "Old to New" mode
1022 //================================================================================
1024 void renumber( DataArrayIdType* ids, const DataArrayIdType* o2nRenumber )
1026 if ( !ids || !o2nRenumber )
1028 mcIdType * id = ids->getPointer();
1029 const mcIdType * o2n = o2nRenumber->getConstPointer();
1030 for ( ; id < ids->end(); ++id )
1037 //================================================================================
1039 * \brief Fill up ConnectZone's stored in _topology with nodal correspondences
1040 * \param [in] nodeMapping - mapping between old nodes and new nodes
1041 * (iolddomain,ioldnode)->(inewdomain,inewnode)
1042 * \param [in] o2nRenumber - renumbering array returned by mesh->sortCellsInMEDFileFrmt()
1044 * \param [in] nbInitialDomains - nb of old domains
1046 //================================================================================
1048 void MEDPARTITIONER::MeshCollection::buildConnectZones( const NodeMapping& nodeMapping,
1049 const std::vector<MEDCoupling::DataArrayIdType*> & o2nRenumber,
1050 int nbInitialDomains)
1052 if ( !MyGlobals::_Create_Joints || _topology->nbDomain() < 2 )
1055 if ( MyGlobals::_World_Size > 1 )
1057 _topology->getCZ().clear();
1058 return; // not implemented for parallel mode
1061 // at construction, _topology creates cell correspondences basing on Graph information,
1063 // 1) add node correspondences,
1064 // 2) split cell correspondences by cell geometry types
1065 // 3) sort ids to be in ascending order
1067 const int nb_domains = _topology->nbDomain();
1069 // ==================================
1070 // 1) add node correspondences
1071 // ==================================
1073 std::vector< std::vector< std::vector< mcIdType > > > nodeCorresp( nb_domains );
1074 for ( int idomain = 0; idomain < nb_domains; ++idomain )
1076 nodeCorresp[ idomain ].resize( nb_domains );
1079 NodeMapping::const_iterator nmIt1, nmIt2 = nodeMapping.begin();
1080 for ( nmIt1 = nmIt2; nmIt1 != nodeMapping.end(); nmIt1 = nmIt2 )
1082 // look for an "old" node mapped into several "new" nodes in different domains
1084 while ( ++nmIt2 != nodeMapping.end() && nmIt2->first == nmIt1->first )
1085 nbSameOld += ( nmIt2->second != nmIt1->second );
1087 if ( nbSameOld > 0 )
1089 NodeMapping::const_iterator nmEnd = nmIt2;
1090 for ( ; true; ++nmIt1 )
1093 if ( ++nmIt2 == nmEnd )
1095 int dom1 = nmIt1->second.first;
1096 mcIdType node1 = nmIt1->second.second;
1097 for ( ; nmIt2 != nmEnd; ++nmIt2 )
1099 int dom2 = nmIt2->second.first;
1100 mcIdType node2 = nmIt2->second.second;
1103 nodeCorresp[ dom1 ][ dom2 ].push_back( node1 );
1104 nodeCorresp[ dom1 ][ dom2 ].push_back( node2 );
1105 nodeCorresp[ dom2 ][ dom1 ].push_back( node2 );
1106 nodeCorresp[ dom2 ][ dom1 ].push_back( node1 );
1113 // add nodeCorresp to czVec
1115 std::vector<MEDPARTITIONER::ConnectZone*>& czVec = _topology->getCZ();
1117 for ( int idomain = 0; idomain < nb_domains; ++idomain )
1119 for ( int idomainNear = 0; idomainNear < nb_domains; ++idomainNear )
1121 std::vector< mcIdType > & corresp = nodeCorresp[ idomain ][ idomainNear ];
1122 if ( corresp.empty() )
1125 MEDPARTITIONER::ConnectZone* cz = 0;
1126 for ( size_t i = 0; i < czVec.size() && !cz; ++i )
1128 czVec[i]->getLocalDomainNumber () == idomain &&
1129 czVec[i]->getDistantDomainNumber() == idomainNear )
1134 cz = new MEDPARTITIONER::ConnectZone();
1135 cz->setName( "Nodal Connect Zone defined by MEDPARTITIONER" );
1136 cz->setLocalDomainNumber ( idomain );
1137 cz->setDistantDomainNumber( idomainNear );
1138 czVec.push_back(cz);
1141 cz->setNodeCorresp( &corresp[0], ToIdType( corresp.size()/2 ));
1145 // ==========================================================
1146 // 2) split cell correspondences by cell geometry types
1147 // ==========================================================
1149 for ( size_t i = 0; i < czVec.size(); ++i )
1151 MEDPARTITIONER::ConnectZone* cz = czVec[i];
1153 cz->getEntityCorrespNumber( 0,0 ) == 0 ||
1154 cz->getLocalDomainNumber () > (int)_mesh.size() ||
1155 cz->getDistantDomainNumber() > (int)_mesh.size() )
1157 MEDCoupling::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber () ];
1158 MEDCoupling::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ];
1160 // separate ids of two domains
1161 const MEDCoupling::MEDCouplingSkyLineArray *corrArray = cz->getEntityCorresp( 0, 0 );
1162 const DataArrayIdType* ids12 = corrArray->getValuesArray();
1163 MCAuto<DataArrayIdType> ids1, ids2, ids12Sorted;
1164 ids1 = ids12->selectByTupleIdSafeSlice( 0, corrArray->getLength(), 2 );
1165 ids2 = ids12->selectByTupleIdSafeSlice( 1, corrArray->getLength(), 2 );
1167 // renumber cells according to mesh->sortCellsInMEDFileFrmt()
1168 renumber( ids1, o2nRenumber[ cz->getLocalDomainNumber() ]);
1169 renumber( ids2, o2nRenumber[ cz->getDistantDomainNumber() ]);
1171 // check nb cell types
1172 std::set<INTERP_KERNEL::NormalizedCellType> types1, types2;
1173 types1 = mesh1->getTypesOfPart( ids1->begin(), ids1->end() );
1174 types2 = mesh2->getTypesOfPart( ids2->begin(), ids2->end() );
1175 if ( types1.size() < 1 || types2.size() < 1 )
1176 continue; // parallel mode?
1178 MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1179 for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1181 czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
1182 czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1185 if ( types1.size() == 1 && types2.size() == 1 ) // split not needed, only sort
1187 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1188 cz->setEntityCorresp( *types1.begin(), *types2.begin(),
1189 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1191 if ( cz21 )// set 2->1 correspondence
1193 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1194 cz21->setEntityCorresp( *types2.begin(), *types1.begin(),
1195 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1198 else // split and sort
1200 typedef std::pair< std::vector< mcIdType >, std::vector< mcIdType > > T2Vecs;
1201 T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ];
1204 const mcIdType nbIds = ids1->getNbOfElems();
1205 const mcIdType * p1 = ids1->begin(), * p2 = ids2->begin();
1206 for ( mcIdType i = 0; i < nbIds; ++i )
1208 t1 = mesh1->getTypeOfCell( p1[ i ]);
1209 t2 = mesh2->getTypeOfCell( p2[ i ]);
1210 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1211 ids.first .push_back( p1[ i ]);
1212 ids.second.push_back( p1[ i ]);
1215 const int maxType = int( INTERP_KERNEL::NORM_MAXTYPE );
1216 for ( t1 = 0; t1 < maxType; ++t1 )
1217 for ( t2 = 0; t2 < maxType; ++t2 )
1219 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1220 if ( ids.first.empty() ) continue;
1221 p1 = & ids.first[0];
1222 p2 = & ids.second[0];
1223 ids1->desallocate();
1224 ids1->pushBackValsSilent( p1, p1+ids.first.size() );
1225 ids2->desallocate();
1226 ids2->pushBackValsSilent( p2, p2+ids.first.size() );
1227 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1229 cz->setEntityCorresp( t1, t2,
1230 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1232 if ( cz21 )// set 2->1 correspondence
1234 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1235 cz21->setEntityCorresp( t2, t1,
1236 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1242 cz->setEntityCorresp( 0, 0, 0, 0 ); // erase ids computed by _topology
1244 cz21->setEntityCorresp( 0, 0, 0, 0 );
1249 // ==========================================
1250 // 3) sort node ids to be in ascending order
1251 // ==========================================
1253 const bool removeEqual = ( nbInitialDomains > 1 );
1255 for ( size_t i = 0; i < czVec.size(); ++i )
1257 MEDPARTITIONER::ConnectZone* cz = czVec[i];
1258 if ( !cz || cz->getNodeNumber() < 1 )
1260 if ( cz->getDistantDomainNumber() < cz->getLocalDomainNumber() )
1261 continue; // treat a pair of domains once
1263 MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1264 for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1266 czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
1267 czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1270 // separate ids of two domains
1271 const MEDCoupling::MEDCouplingSkyLineArray *corrArray = cz->getNodeCorresp();
1272 const DataArrayIdType *ids12 = corrArray->getValuesArray();
1273 MCAuto<DataArrayIdType> ids1, ids2, ids12Sorted;
1274 ids1 = ids12->selectByTupleIdSafeSlice( 0, corrArray->getLength(), 2 );
1275 ids2 = ids12->selectByTupleIdSafeSlice( 1, corrArray->getLength(), 2 );
1277 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/0, removeEqual );
1278 cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1280 if ( cz21 )// set 2->1 correspondence
1282 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0, false );
1283 cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1288 //================================================================================
1290 * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
1291 * group (where "n" and "p" are domain IDs)
1293 //================================================================================
1295 void MEDPARTITIONER::MeshCollection::buildBoundaryFaces()
1297 if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates )
1300 if ( getMeshDimension() < 2 )
1303 using MEDCoupling::MEDCouplingUMesh;
1304 using MEDCoupling::DataArrayDouble;
1305 using MEDCoupling::DataArrayInt;
1307 std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
1308 std::size_t nbMeshes = faceMeshes.size();
1310 //preparing bounding box trees for accelerating search of coincident faces
1311 std::vector<BBTreeOfDim* > bbTrees(nbMeshes);
1312 std::vector<DataArrayDouble*>bbox (nbMeshes);
1313 for (unsigned int inew = 0; inew < nbMeshes-1; inew++)
1314 if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
1316 DataArrayDouble* bcCoords = faceMeshes[inew]->computeCellCenterOfMass();
1317 bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
1318 bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
1319 bbox[inew]->getConstPointer(),0,0,
1320 bbox[inew]->getNumberOfTuples());
1321 bcCoords->decrRef();
1324 // loop on domains to find joint faces between them
1325 for (unsigned int inew1 = 0; inew1 < nbMeshes; inew1++ )
1327 for (unsigned int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
1329 MEDCouplingUMesh* mesh1 = 0;
1330 MEDCouplingUMesh* mesh2 = 0;
1331 //MEDCouplingUMesh* recvMesh = 0;
1332 bool mesh1Here = true, mesh2Here = true;
1333 if (isParallelMode())
1336 mesh1Here = _domain_selector->isMyDomain(inew1);
1337 mesh2Here = _domain_selector->isMyDomain(inew2);
1338 if ( !mesh1Here && mesh2Here )
1340 //send mesh2 to domain of mesh1
1341 _domain_selector->sendMesh(*faceMeshes[inew2],
1342 _domain_selector->getProcessorID(inew1));
1344 else if ( mesh1Here && !mesh2Here )
1346 //receiving mesh2 from a distant domain
1347 _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
1348 if ( faceMeshes[ inew2 ] )
1349 faceMeshes[ inew2 ]->decrRef();
1350 faceMeshes[ inew2 ] = mesh2;
1354 if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
1355 if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
1357 // find coincident faces
1358 std::vector< mcIdType > faces1, faces2;
1359 if ( mesh1 && mesh2 )
1361 const DataArrayDouble* coords2 = mesh2->computeCellCenterOfMass();
1362 const double* c2 = coords2->getConstPointer();
1363 const std::size_t dim = coords2->getNumberOfComponents();
1364 const mcIdType nbFaces2 = mesh2->getNumberOfCells();
1365 const mcIdType nbFaces1 = mesh1->getNumberOfCells();
1367 for (mcIdType i2 = 0; i2 < nbFaces2; i2++)
1369 std::vector<mcIdType> coincFaces;
1370 bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
1371 if (coincFaces.size()!=0)
1373 mcIdType i1 = coincFaces[0];
1374 // if ( coincFaces.size() > 1 )
1376 // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
1377 // i1 = coincFaces[ i2nb->second++ ];
1379 if ( i1 < nbFaces1 ) // protection from invalid elements
1381 faces1.push_back( i1 );
1382 faces2.push_back( i2 );
1389 if ( isParallelMode())
1392 if ( mesh1Here && !mesh2Here )
1394 //send faces2 to domain of recvMesh
1395 SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
1397 else if ( !mesh1Here && mesh2Here )
1399 //receiving ids of faces from a domain of mesh1
1400 RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
1405 // recvMesh->decrRef();
1407 // Create group "JOINT_inew1_inew2_Faces" and corresponding families
1408 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1410 createJointGroup( is2nd ? faces2 : faces1,
1411 inew1 , inew2, is2nd );
1414 } // loop on the 2nd domains (inew2)
1415 } // loop on the 1st domains (inew1)
1418 // delete bounding box trees
1419 for (unsigned int inew = 0; inew < nbMeshes-1; inew++)
1420 if (isParallelMode() && _domain_selector->isMyDomain(inew))
1422 bbox[inew]->decrRef();
1423 delete bbTrees[inew];
1427 //================================================================================
1429 * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
1430 * \param faces - face ids to include into the group
1431 * \param inew1 - index of the 1st domain
1432 * \param inew2 - index of the 2nd domain
1433 * \param is2nd - in which (1st or 2nd) domain to create the group
1435 //================================================================================
1437 void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< mcIdType >& faces,
1442 // get the name of JOINT group
1443 std::string groupName;
1445 std::ostringstream oss;
1447 << (is2nd ? inew2 : inew1 ) << "_"
1448 << (is2nd ? inew1 : inew2 ) << "_"
1449 << ( getMeshDimension()==2 ? "Edge" : "Face" );
1450 groupName = oss.str();
1453 // remove existing "JOINT_*" group
1454 _group_info.erase( groupName );
1456 // get family IDs array
1457 mcIdType* famIDs = 0;
1458 int inew = (is2nd ? inew2 : inew1 );
1459 mcIdType totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
1460 std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
1461 if ( !_map_dataarray_int.count(cle) )
1463 if ( totalNbFaces > 0 )
1465 MEDCoupling::DataArrayIdType* p=MEDCoupling::DataArrayIdType::New();
1466 p->alloc( totalNbFaces, 1 );
1468 famIDs = p->getPointer();
1469 _map_dataarray_int[cle]=p;
1474 famIDs = _map_dataarray_int.find(cle)->second->getPointer();
1476 // find a family ID of an existing JOINT group
1477 mcIdType familyID = 0;
1478 std::map<std::string, mcIdType>::iterator name2id = _family_info.find( groupName );
1479 if ( name2id != _family_info.end() )
1480 familyID = name2id->second;
1482 // remove faces from the familyID-the family
1483 if ( familyID != 0 && famIDs )
1484 for ( mcIdType i = 0; i < totalNbFaces; ++i )
1485 if ( famIDs[i] == familyID )
1488 if ( faces.empty() )
1491 if ( familyID == 0 ) // generate a family ID for JOINT group
1493 std::set< mcIdType > familyIDs;
1494 for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
1495 familyIDs.insert( name2id->second );
1496 // find the next free family ID
1497 int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
1500 if ( !familyIDs.count( ++familyID ))
1503 while ( freeIdCount > 0 );
1506 // push faces to familyID-th group
1507 if ( faces.back() >= totalNbFaces )
1508 throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
1509 for ( size_t i = 0; i < faces.size(); ++i )
1510 famIDs[ faces[i] ] = familyID;
1512 // register JOINT group and family
1513 _family_info[ groupName ] = familyID; // name of the group and family is same
1514 _group_info [ groupName ].push_back( groupName );
1517 /*! constructing the MESH collection from a distributed file
1519 * \param filename name of the master file containing the list of all the MED files
1521 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
1523 _owns_topology(true),
1525 _domain_selector( 0 ),
1526 _i_non_empty_mesh(-1),
1527 _driver_type(MEDPARTITIONER::Undefined),
1528 _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
1529 _family_splitting(false),
1530 _create_empty_groups(false),
1535 _driver=new MeshCollectionMedXmlDriver(this);
1536 _driver->read (filename.c_str());
1537 _driver_type = MedXml;
1540 { // Handle all exceptions
1541 if ( _driver ) delete _driver; _driver=0;
1544 _driver=new MeshCollectionMedAsciiDriver(this);
1545 _driver->read (filename.c_str());
1546 _driver_type=MedAscii;
1552 throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
1555 for ( unsigned int idomain = 0; idomain < _mesh.size(); ++idomain )
1556 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1557 _i_non_empty_mesh = idomain;
1560 /*! Constructing the MESH collection from selected domains of a distributed file
1562 * \param filename - name of the master file containing the list of all the MED files
1563 * \param domainSelector - selector of domains to load
1565 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
1567 _owns_topology(true),
1569 _domain_selector( &domainSelector ),
1570 _i_non_empty_mesh(-1),
1571 _driver_type(MEDPARTITIONER::Undefined),
1572 _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
1573 _family_splitting(false),
1574 _create_empty_groups(false),
1577 std::string myfile=filename;
1578 std::size_t found=myfile.find(".xml");
1579 if (found!=std::string::npos) //file .xml
1583 _driver=new MeshCollectionMedXmlDriver(this);
1584 _driver->read ( filename.c_str(), _domain_selector );
1585 _driver_type = MedXml;
1588 { // Handle all exceptions
1590 throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
1595 found=myfile.find(".med");
1596 if (found!=std::string::npos) //file .med single
1598 //make a temporary file .xml and retry MedXmlDriver
1600 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
1602 <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
1603 <description what=\"\" when=\"\"/>\n \
1605 <mesh name=\"$meshName\"/>\n \
1608 <subdomain number=\"1\"/>\n \
1609 <global_numbering present=\"no\"/>\n \
1612 <subfile id=\"1\">\n \
1613 <name>$fileName</name>\n \
1614 <machine>localhost</machine>\n \
1618 <mesh name=\"$meshName\">\n \
1619 <chunk subdomain=\"1\">\n \
1620 <name>$meshName</name>\n \
1625 std::vector<std::string> meshNames=GetMeshNames(myfile);
1626 xml.replace(xml.find("$fileName"),9,myfile);
1627 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1628 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1629 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1630 std::string nameFileXml(myfile);
1631 nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1632 std::string nameFileXmlDN,nameFileXmlBN;
1633 MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1634 nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1635 if (_domain_selector->rank()==0) //only on to write it
1637 std::ofstream f(nameFileXml.c_str());
1642 if (MyGlobals::_World_Size>1)
1643 MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1647 _driver=new MeshCollectionMedXmlDriver(this);
1648 _driver->read ( nameFileXml.c_str(), _domain_selector );
1649 _driver_type = MedXml;
1652 { // Handle all exceptions
1653 delete _driver; _driver=0;
1654 throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1661 _driver=new MeshCollectionMedAsciiDriver(this);
1662 _driver->read ( filename.c_str(), _domain_selector );
1663 _driver_type=MedAscii;
1669 throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1673 // find non-empty domain mesh
1674 for ( unsigned int idomain = 0; idomain < _mesh.size(); ++idomain )
1675 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1676 _i_non_empty_mesh = idomain;
1680 //check for all proc/file compatibility of _field_descriptions
1682 _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1684 _field_descriptions=MyGlobals::_Field_Descriptions;
1687 catch(INTERP_KERNEL::Exception& e)
1689 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1690 throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1695 //check for all proc/file compatibility of _family_info
1696 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1697 _family_info=DevectorizeToMapOfStringInt(v2);
1699 catch(INTERP_KERNEL::Exception& e)
1701 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1702 throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1707 //check for all proc/file compatibility of _group_info
1708 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1709 _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1711 catch(INTERP_KERNEL::Exception& e)
1713 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1714 throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1719 /*! constructing the MESH collection from a sequential MED-file
1721 * \param filename MED file
1722 * \param meshname name of the mesh that is to be read
1724 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1726 _owns_topology(true),
1728 _domain_selector( 0 ),
1729 _i_non_empty_mesh(-1),
1731 _driver_type(MEDPARTITIONER::MedXml),
1732 _subdomain_boundary_creates(MyGlobals::_Create_Boundary_Faces),
1733 _family_splitting(false),
1734 _create_empty_groups(false),
1737 try // avoid memory leak in case of inexistent filename
1739 retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1745 throw INTERP_KERNEL::Exception("problem reading .med files");
1747 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1748 _i_non_empty_mesh = 0;
1751 MEDPARTITIONER::MeshCollection::~MeshCollection()
1753 for (std::size_t i=0; i<_mesh.size();i++)
1754 if (_mesh[i]!=0) _mesh[i]->decrRef();
1756 for (std::size_t i=0; i<_cell_family_ids.size();i++)
1757 if (_cell_family_ids[i]!=0)
1758 _cell_family_ids[i]->decrRef();
1760 for (std::size_t i=0; i<_face_mesh.size();i++)
1761 if (_face_mesh[i]!=0)
1762 _face_mesh[i]->decrRef();
1764 for (std::size_t i=0; i<_face_family_ids.size();i++)
1765 if (_face_family_ids[i]!=0)
1766 _face_family_ids[i]->decrRef();
1768 for (std::map<std::string, MEDCoupling::DataArrayIdType*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1769 if ((*it).second!=0)
1770 (*it).second->decrRef();
1772 for (std::map<std::string, MEDCoupling::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1773 if ((*it).second!=0)
1774 (*it).second->decrRef();
1777 if (_topology!=0 && _owns_topology)
1780 delete _joint_finder;
1784 /*! constructing the MESH collection from a file
1786 * The method creates as many MED-files as there are domains in the
1787 * collection. It also creates a master file that lists all the MED files.
1788 * The MED files created in this manner contain joints that describe the
1789 * connectivity between subdomains.
1791 * \param filename name of the master file that will contain the list of the MED files
1794 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1796 //suppresses link with driver so that it can be changed for writing
1799 retrieveDriver()->write (filename.c_str(), _domain_selector);
1802 /*! creates or gets the link to the collection driver
1804 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1808 switch(_driver_type)
1811 _driver=new MeshCollectionMedXmlDriver(this);
1814 _driver=new MeshCollectionMedAsciiDriver(this);
1817 throw INTERP_KERNEL::Exception("Unrecognized driver");
1824 /*! gets an existing driver
1827 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1833 * retrieves the mesh dimension
1835 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1837 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1840 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1843 for (size_t i=0; i<_mesh.size(); i++)
1850 mcIdType MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1853 for (size_t i=0; i<_mesh.size(); i++)
1855 if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1860 mcIdType MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1863 for (size_t i=0; i<_face_mesh.size(); i++)
1865 if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1870 std::vector<MEDCoupling::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1875 std::vector<MEDCoupling::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1880 MEDCoupling::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1882 return _mesh[idomain];
1885 MEDCoupling::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1887 return _face_mesh[idomain];
1890 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1893 return _topology->getCZ();
1895 static std::vector<MEDPARTITIONER::ConnectZone*> noCZ;
1899 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1904 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo, bool takeOwneship)
1908 throw INTERP_KERNEL::Exception("topology is already set");
1913 _owns_topology = takeOwneship;
1917 /*! Method creating the cell graph in serial mode
1919 * \param array returns the pointer to the structure that contains the graph
1920 * \param edgeweight returns the pointer to the table that contains the edgeweights
1921 * (only used if indivisible regions are required)
1923 void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1928 using std::make_pair;
1931 if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
1932 const MEDCoupling::MEDCouplingUMesh* mesh=_mesh[0];
1933 if (MyGlobals::_Verbose>50)
1934 std::cout<<"getting nodal connectivity"<<std::endl;
1936 //looking for reverse nodal connectivity i global numbering
1937 if (isParallelMode() && !_domain_selector->isMyDomain(0))
1939 vector<mcIdType> value;
1940 vector<mcIdType> index(1,0);
1942 array = MEDCoupling::MEDCouplingSkyLineArray::New(index,value);
1945 array=mesh->generateGraph();
1947 /*! Method creating the cell graph in multidomain mode
1949 * \param array returns the pointer to the structure that contains the graph
1950 * \param edgeweight returns the pointer to the table that contains the edgeweights
1951 * (only used if indivisible regions are required)
1953 void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1955 using std::multimap;
1957 using std::make_pair;
1960 std::multimap< mcIdType, mcIdType > node2cell;
1961 std::map< pair<mcIdType,mcIdType>, mcIdType > cell2cellcounter;
1962 std::multimap<mcIdType,mcIdType> cell2cell;
1964 std::vector<std::vector<std::multimap<mcIdType,mcIdType> > > commonDistantNodes;
1965 int nbdomain=_topology->nbDomain();
1967 if (isParallelMode())
1969 _joint_finder=new JointFinder(*this);
1970 _joint_finder->findCommonDistantNodes();
1971 commonDistantNodes=_joint_finder->getDistantNodeCell();
1974 if (MyGlobals::_Verbose>500)
1975 _joint_finder->print();
1978 if (MyGlobals::_Verbose>50)
1979 std::cout<<"getting nodal connectivity"<<std::endl;
1980 //looking for reverse nodal connectivity i global numbering
1982 for (int idomain=0; idomain<nbdomain; idomain++)
1984 if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1986 meshDim = _mesh[idomain]->getMeshDimension();
1988 MEDCoupling::DataArrayIdType* index=MEDCoupling::DataArrayIdType::New();
1989 MEDCoupling::DataArrayIdType* revConn=MEDCoupling::DataArrayIdType::New();
1990 mcIdType nbNodes=_mesh[idomain]->getNumberOfNodes();
1991 _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1992 //problem saturation over 1 000 000 nodes for 1 proc
1993 if (MyGlobals::_Verbose>100)
1994 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1995 mcIdType* index_ptr=index->getPointer();
1996 mcIdType* revConnPtr=revConn->getPointer();
1997 for (mcIdType i=0; i<nbNodes; i++)
1999 for (mcIdType icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
2001 mcIdType globalNode=_topology->convertNodeToGlobal(idomain,i);
2002 mcIdType globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
2003 node2cell.insert(make_pair(globalNode, globalCell));
2009 for (int iother=0; iother<nbdomain; iother++)
2011 std::multimap<mcIdType,mcIdType>::iterator it;
2012 int isource=idomain;
2014 for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
2015 it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
2017 mcIdType globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
2018 mcIdType globalCell=(*it).second;
2019 node2cell.insert(make_pair(globalNode, globalCell));
2025 //creating graph arcs (cell to cell relations)
2026 //arcs are stored in terms of (index,value) notation
2029 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
2030 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
2032 //warning here one node have less than or equal effective number of cell with it
2033 //but cell could have more than effective nodes
2034 //because other equals nodes in other domain (with other global inode)
2035 if (MyGlobals::_Verbose>50)
2036 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
2038 for (mcIdType inode=0;inode<_topology->nbNodes();inode++)
2040 typedef multimap<mcIdType,mcIdType>::const_iterator MI;
2041 std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
2042 for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
2043 for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
2045 mcIdType icell1=cell1->second;
2046 mcIdType icell2=cell2->second;
2047 if (icell1>icell2) std::swap(icell1,icell2);
2048 std::map<pair<mcIdType,mcIdType>,mcIdType>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2049 if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2050 else (it->second)++;
2053 // for (mcIdType icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
2055 // typedef multimap<int,mcIdType>::const_iterator MI;
2056 // std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
2057 // for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
2059 // std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
2060 // for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
2062 // mcIdType icell2=cell2->second;
2063 // std::map<pair<int,mcIdType>,mcIdType>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2064 // if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2065 // else (it->second)++;
2071 //converting the counter to a multimap structure
2072 for (std::map<pair<mcIdType,mcIdType>,mcIdType>::const_iterator it=cell2cellcounter.begin();
2073 it!=cell2cellcounter.end();
2075 if (it->second>=meshDim)
2077 cell2cell.insert(std::make_pair(it->first.first,it->first.second));
2078 cell2cell.insert(std::make_pair(it->first.second, it->first.first));
2082 if (MyGlobals::_Verbose>50)
2083 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
2084 //filling up index and value to create skylinearray structure
2085 std::vector <mcIdType> index,value;
2089 for (int idomain=0; idomain<nbdomain; idomain++)
2091 if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
2092 mcIdType nbCells=_mesh[idomain]->getNumberOfCells();
2093 for (mcIdType icell=0; icell<nbCells; icell++)
2096 mcIdType globalCell=_topology->convertCellToGlobal(idomain,icell);
2097 multimap<mcIdType,mcIdType>::iterator it;
2098 pair<multimap<mcIdType,mcIdType>::iterator,multimap<mcIdType,mcIdType>::iterator> ret;
2099 ret=cell2cell.equal_range(globalCell);
2100 for (it=ret.first; it!=ret.second; ++it)
2102 mcIdType ival=(*it).second; //no adding one existing yet
2103 for (mcIdType i=idep ; i<idep+size ; i++)
2112 value.push_back(ival);
2116 idep=index[index.size()-1]+size;
2117 index.push_back(idep);
2121 array=MEDCoupling::MEDCouplingSkyLineArray::New(index,value);
2123 if (MyGlobals::_Verbose>100)
2125 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
2126 index.size()-1 << " " << value.size() << std::endl;
2127 std::size_t max=index.size()>15?15:index.size();
2130 for (std::size_t i=0; i<max; ++i)
2131 std::cout<<index[i]<<" ";
2132 std::cout << "... " << index[index.size()-1] << std::endl;
2133 for (std::size_t i=0; i<max; ++i)
2134 std::cout<< value[i] << " ";
2135 mcIdType ll=index[index.size()-1]-1;
2136 std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
2143 /*! Creates the partition corresponding to the cell graph and the partition number
2145 * \param nbdomain number of subdomains for the newly created graph
2147 * returns a topology based on the new graph
2149 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
2150 Graph::splitter_type split,
2151 const std::string& options_string,
2152 int *user_edge_weights,
2153 int *user_vertices_weights)
2155 if (MyGlobals::_Verbose>10)
2156 std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
2159 throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
2160 MEDCoupling::MEDCouplingSkyLineArray* array=0;
2163 if (_topology->nbDomain()>1 || isParallelMode())
2164 buildParallelCellGraph(array,edgeweights);
2166 buildCellGraph(array,edgeweights);
2168 Graph* cellGraph = 0;
2172 if ( isParallelMode() && MyGlobals::_World_Size > 1 )
2174 #ifdef MED_ENABLE_PARMETIS
2175 if (MyGlobals::_Verbose>10)
2176 std::cout << "ParMETISGraph" << std::endl;
2177 cellGraph=new ParMETISGraph(array,edgeweights);
2182 #ifdef MED_ENABLE_METIS
2183 if (MyGlobals::_Verbose>10)
2184 std::cout << "METISGraph" << std::endl;
2185 cellGraph=new METISGraph(array,edgeweights);
2189 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available (or you're not running in //). Check your products, please.");
2193 #ifdef MED_ENABLE_SCOTCH
2194 if (MyGlobals::_Verbose>10)
2195 std::cout << "SCOTCHGraph" << std::endl;
2196 cellGraph=new SCOTCHGraph(array,edgeweights);
2198 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
2203 //!user-defined weights
2204 if (user_edge_weights!=0)
2205 cellGraph->setEdgesWeights(user_edge_weights);
2206 if (user_vertices_weights!=0)
2207 cellGraph->setVerticesWeights(user_vertices_weights);
2209 if (MyGlobals::_Is0verbose>10)
2210 std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
2211 cellGraph->partGraph(nbdomain, options_string, _domain_selector);
2213 if (MyGlobals::_Is0verbose>10)
2214 std::cout << "building new topology" << std::endl;
2215 //cellGraph is a shared pointer
2216 Topology *topology=0;
2217 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2219 delete [] edgeweights;
2221 if (MyGlobals::_Verbose>11)
2222 std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
2226 /*! Creates a topology for a partition specified by the user
2228 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
2230 * returns a topology based on the new partition
2232 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
2234 MEDCoupling::MEDCouplingSkyLineArray* array=0;
2237 if ( _topology->nbDomain()>1)
2238 buildParallelCellGraph(array,edgeweights);
2240 buildCellGraph(array,edgeweights);
2243 std::set<int> domains;
2244 for (mcIdType i=0; i<_topology->nbCells(); i++)
2246 domains.insert(partition[i]);
2248 cellGraph=new UserGraph(array, partition, _topology->nbCells());
2250 //cellGraph is a shared pointer
2251 Topology *topology=0;
2252 int nbdomain=(int)domains.size();
2253 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2254 // if (array!=0) delete array;
2259 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
2261 for (int i=0; i<_topology->nbDomain(); i++)
2263 std::ostringstream oss;
2265 if (!isParallelMode() || _domain_selector->isMyDomain(i))
2266 _mesh[i]->setName(oss.str());
2270 MEDCoupling::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
2271 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
2272 //something like MCAuto<MEDCouplingFieldDouble> f2=ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
2274 int rank=MyGlobals::_Rank;
2275 std::string tag="ioldFieldDouble="+IntToStr(iold);
2276 std::string descriptionIold=descriptionField+SerializeFromString(tag);
2277 if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
2279 if (MyGlobals::_Verbose>300)
2280 std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
2281 MEDCoupling::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
2284 if (MyGlobals::_Verbose>200)
2285 std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
2286 std::string description, fileName, meshName, fieldName;
2287 int typeField, DT, IT, entity;
2288 fileName=MyGlobals::_File_Names[iold];
2289 if (MyGlobals::_Verbose>10)
2290 std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
2291 FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
2292 meshName=MyGlobals::_Mesh_Names[iold];
2294 MCAuto<MEDCoupling::MEDCouplingField> f2Tmp(ReadField((MEDCoupling::TypeOfField) typeField, fileName, meshName, 0, fieldName, DT, IT));
2295 MCAuto<MEDCoupling::MEDCouplingFieldDouble> f2(MEDCoupling::DynamicCast<MEDCouplingField,MEDCouplingFieldDouble>(f2Tmp));
2297 MEDCoupling::DataArrayDouble* res=f2->getArray();
2298 //to know names of components
2299 std::vector<std::string> browse=BrowseFieldDouble(f2);
2300 std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
2301 if (MyGlobals::_Verbose>10)
2302 std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
2303 MyGlobals::_General_Informations.push_back(localFieldInformation);
2305 _map_dataarray_double[descriptionIold]=res;
2309 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
2310 //to have unique valid fields names/pointers/descriptions for partitionning
2311 //filter _field_descriptions to be in all procs compliant and equal
2313 std::size_t nbfiles=MyGlobals::_File_Names.size(); //nb domains
2316 nbfiles=_topology->nbDomain();
2318 std::vector<std::string> r2;
2319 //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
2320 for (std::size_t i=0; i<_field_descriptions.size(); i++)
2322 std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
2323 for (std::size_t ii=0; ii<r1.size(); ii++)
2324 r2.push_back(r1[ii]);
2326 //here vector(procs*fields) of serialised vector(description) data
2327 _field_descriptions=r2;
2328 std::size_t nbfields=_field_descriptions.size(); //on all domains
2329 if ((nbfields%nbfiles)!=0)
2331 if (MyGlobals::_Rank==0)
2333 std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
2334 << "fileMedNames :" << std::endl
2335 << ReprVectorOfString(MyGlobals::_File_Names)
2336 << "field_descriptions :" << std::endl
2337 << ReprVectorOfString(MyGlobals::_Field_Descriptions);
2339 throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
2341 _field_descriptions.resize(nbfields/nbfiles);
2342 for (std::size_t i=0; i<_field_descriptions.size(); i++)
2344 std::string str=_field_descriptions[i];
2345 str=EraseTagSerialized(str,"idomain=");
2346 str=EraseTagSerialized(str,"fileName=");
2347 _field_descriptions[i]=str;
2351 //returns true if inodes of a face are in inodes of a cell
2352 bool isFaceOncell(std::vector< mcIdType >& inodesFace, std::vector< mcIdType >& inodesCell)
2355 std::size_t nbok=inodesFace.size();
2356 for (std::size_t i=0; i<nbok; i++)
2358 mcIdType ii=inodesFace[i];
2360 std::cout << "isFaceOncell problem inodeface<0" << std::endl;
2361 for (std::size_t j=0; j<inodesCell.size(); j++)
2363 if (ii==inodesCell[j])
2366 break; //inode of face found
2370 break; //inode of face not found do not continue...
2372 return (ires==nbok);
2375 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
2377 for (int inew=0; inew<_topology->nbDomain(); inew++)
2379 if (!isParallelMode() || _domain_selector->isMyDomain(inew))
2381 if (MyGlobals::_Verbose>200)
2382 std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
2383 MEDCoupling::MEDCouplingUMesh* mcel=_mesh[inew];
2384 MEDCoupling::MEDCouplingUMesh* mfac=_face_mesh[inew];
2386 //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
2387 std::vector<mcIdType> nodeIds;
2388 getNodeIds(*mcel, *mfac, nodeIds);
2389 if (nodeIds.size()==0)
2390 continue; //one empty mesh nothing to do
2392 MEDCoupling::DataArrayIdType *revNodalCel=MEDCoupling::DataArrayIdType::New();
2393 MEDCoupling::DataArrayIdType *revNodalIndxCel=MEDCoupling::DataArrayIdType::New();
2394 mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
2395 mcIdType *revC=revNodalCel->getPointer();
2396 mcIdType *revIndxC=revNodalIndxCel->getPointer();
2398 std::vector< mcIdType > faceOnCell;
2399 std::vector< mcIdType > faceNotOnCell;
2400 mcIdType nbface=mfac->getNumberOfCells();
2401 for (mcIdType iface=0; iface<nbface; iface++)
2404 std::vector< mcIdType > inodesFace;
2405 mfac->getNodeIdsOfCell(iface, inodesFace);
2406 mcIdType nbnodFace=ToIdType(inodesFace.size());
2407 if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
2408 continue; // invalid node ids
2409 //set inodesFace in mcel
2411 for (int i=0; i<nbnodFace; i++)
2412 nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
2413 if ( nbok != nbnodFace )
2415 mcIdType inod=inodesFace[0];
2418 std::cout << "filterFaceOnCell problem 1" << std::endl;
2421 mcIdType nbcell=revIndxC[inod+1]-revIndxC[inod];
2422 for (mcIdType j=0; j<nbcell; j++) //look for each cell with inod
2424 mcIdType icel=revC[revIndxC[inod]+j];
2425 std::vector< mcIdType > inodesCell;
2426 mcel->getNodeIdsOfCell(icel, inodesCell);
2427 ok=isFaceOncell(inodesFace, inodesCell);
2432 faceOnCell.push_back(iface);
2436 faceNotOnCell.push_back(iface);
2437 if (MyGlobals::_Is0verbose>300)
2438 std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
2442 revNodalCel->decrRef();
2443 revNodalIndxCel->decrRef();
2445 // std::string keyy;
2446 // keyy=Cle1ToStr("filterFaceOnCell",inew);
2447 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
2448 // keyy=Cle1ToStr("filterNotFaceOnCell",inew);
2449 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
2451 // filter the face mesh
2452 if ( faceOnCell.empty() )
2453 _face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
2455 _face_mesh[inew] = (MEDCoupling::MEDCouplingUMesh *)
2456 mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
2459 // filter the face families
2460 std::string key = Cle1ToStr("faceFamily_toArray",inew);
2461 if ( getMapDataArrayInt().count( key ))
2463 MEDCoupling::DataArrayIdType * & fam = getMapDataArrayInt()[ key ];
2464 MEDCoupling::DataArrayIdType * famFilter = MEDCoupling::DataArrayIdType::New();
2465 famFilter->alloc(faceOnCell.size(),1);
2466 mcIdType* pfamFilter = famFilter->getPointer();
2467 mcIdType* pfam = fam->getPointer();
2468 for ( size_t i=0; i<faceOnCell.size(); i++ )
2469 pfamFilter[i]=pfam[faceOnCell[i]];