1 // Copyright (C) 2007-2016 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, 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::_Creates_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::_Creates_Boundary_Faces),
104 _family_splitting(family_splitting),
105 _create_empty_groups(create_empty_groups),
108 std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
109 std::vector<MEDCoupling::DataArrayInt*> 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<int> > > 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 (size_t i=0; i<_mesh.size(); 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<int> > >& new2oldIds,
209 std::vector<MEDCoupling::DataArrayInt*> & 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 int size=(initialCollection._mesh)[iold]->getNumberOfCells();
235 std::vector<int> globalids(size);
236 initialCollection.getTopology()->getCellList(iold, &globalids[0]);
237 std::vector<int> 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();
305 MEDCoupling::DataArrayInt* 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 = coords->getNumberOfComponents();
340 int 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;
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;
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<int>& nodeIds)
400 using MEDPARTITIONER::BBTreeOfDim;
401 //if (!&meshOne || !&meshTwo) return; //empty or not existing
403 BBTreeOfDim* tree = 0;
404 int nv1=meshOne.getNumberOfNodes();
405 MEDCoupling::DataArrayDouble* coords=meshOne.getCoords();
406 int dim = 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 int 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;
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,int>, std::pair<int,int> >& nodeMapping,
438 std::vector<std::vector<std::vector<int> > >& 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 (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
495 meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
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 (int inode=0; inode<(int)nodes.size(); inode++)
507 typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
508 int 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,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
523 if (iter->second==(int)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<const 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);
603 MEDCoupling::MEDCouplingUMesh *bndMesh = 0;
604 if ( _subdomain_boundary_creates &&
606 _mesh[inew]->getNumberOfCells()>0 )
609 ((MEDCoupling::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
610 if (bndMesh->getNumberOfCells()>0)
611 myMeshes.push_back( bndMesh );
614 if (myMeshes.size()>0)
616 meshesCastTo[inew]=MEDCoupling::MEDCouplingUMesh::MergeUMeshes(myMeshes);
617 meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
621 MEDCoupling::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
622 meshesCastTo[inew]=empty;
624 for (int iold=0; iold<nbOldDomain; iold++)
625 if (splitMeshes[inew][iold]!=0)
626 splitMeshes[inew][iold]->decrRef();
630 if (MyGlobals::_Verbose>300)
631 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
636 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastFrom,
637 std::vector<MEDCoupling::MEDCouplingUMesh*>& meshesCastTo,
638 std::vector<MEDCoupling::DataArrayInt*>& arrayFrom,
639 std::string nameArrayTo)
643 int ioldMax=meshesCastFrom.size();
644 int inewMax=meshesCastTo.size();
647 //preparing bounding box trees for accelerating source-target node identifications
648 if (MyGlobals::_Verbose>99)
649 std::cout<<"making accelerating structures"<<std::endl;
650 std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
651 std::vector<MEDCoupling::DataArrayDouble*>bbox(ioldMax);
652 for (int iold =0; iold< ioldMax; iold++)
653 if (isParallelMode() && _domain_selector->isMyDomain(iold))
655 MEDCoupling::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->computeCellCenterOfMass();
656 bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
657 acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
658 sourceCoords->decrRef();
660 // send-recv operations
662 for (int inew=0; inew<inewMax; inew++)
664 for (int iold=0; iold<ioldMax; iold++)
666 //sending arrays for distant domains
667 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
670 _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
672 int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
674 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
675 if (size>0) //no empty
677 sendIds.resize(size);
678 std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
683 sendIds.resize(size);
685 SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
687 //receiving arrays from distant domains
688 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
692 MEDCoupling::MEDCouplingUMesh* recvMesh;
693 _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
695 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
696 RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
697 remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
698 recvMesh->decrRef(); //cww is it correct?
704 //local contributions and aggregation
705 for (int inew=0; inew<inewMax; inew++)
707 for (int iold=0; iold<ioldMax; iold++)
708 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
710 remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
713 for (int iold=0; iold<ioldMax;iold++)
714 if (isParallelMode() && _domain_selector->isMyDomain(iold))
716 bbox[iold]->decrRef();
717 delete acceleratingStructures[iold];
721 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
722 const MEDCoupling::MEDCouplingUMesh& sourceMesh,
723 const MEDCoupling::MEDCouplingUMesh& targetMesh,
724 const int* fromArray,
725 std::string nameArrayTo,
726 const BBTreeOfDim* myTree)
729 if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
730 MEDCoupling::DataArrayDouble* targetCoords=targetMesh.computeCellCenterOfMass();
731 const double* tc=targetCoords->getConstPointer();
732 int targetSize=targetMesh.getNumberOfCells();
733 int sourceSize=sourceMesh.getNumberOfCells();
734 if (MyGlobals::_Verbose>200)
735 std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
736 std::vector<int> ccI;
738 str=nameArrayTo+"_toArray";
739 cle=Cle1ToStr(str,inew);
742 const BBTreeOfDim* tree;
743 bool cleantree=false;
744 MEDCoupling::DataArrayDouble* sourceBBox=0;
745 int dim = targetCoords->getNumberOfComponents();
748 sourceBBox=sourceMesh.computeCellCenterOfMass()->computeBBoxPerTuple(1e-8);
749 tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
753 //first time iold : create and initiate
754 if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
756 if (MyGlobals::_Is0verbose>100)
757 std::cout << "create " << cle << " size " << targetSize << std::endl;
758 MEDCoupling::DataArrayInt* p=MEDCoupling::DataArrayInt::New();
759 p->alloc(targetSize,1);
761 toArray=p->getPointer();
762 _map_dataarray_int[cle]=p;
764 else //other times iold: refind and complete
766 toArray=_map_dataarray_int.find(cle)->second->getPointer();
769 std::map< int, int > isource2nb; // count coincident elements
770 std::map<int,int>::iterator i2nb;
772 for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
774 std::vector<int> intersectingElems;
775 tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
776 if (intersectingElems.size()!=0)
778 int isourcenode=intersectingElems[0];
779 if ( intersectingElems.size() > 1 )
781 i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
782 isourcenode = intersectingElems[ i2nb->second++ ];
784 if ( isourcenode < sourceSize ) // protection from invalid elements
786 toArray[itargetnode]=fromArray[isourcenode];
787 ccI.push_back(itargetnode);
788 ccI.push_back(isourcenode);
792 if (MyGlobals::_Verbose>200)
793 std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
794 //memories intersection for future same job on fields (if no existing cle=no intersection)
795 str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
796 if (MyGlobals::_Verbose>700)
797 std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
799 _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
801 targetCoords->decrRef();
802 if (cleantree) delete tree;
803 if (sourceBBox !=0) sourceBBox->decrRef();
806 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
808 if (nameArrayTo!="cellFieldDouble")
809 throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
811 std::string nameTo="typeData=6"; //resume the type of field casted
812 // send-recv operations
813 int ioldMax=initialCollection.getMesh().size();
814 int inewMax=this->getMesh().size();
815 int iFieldMax=initialCollection.getFieldDescriptions().size();
816 if (MyGlobals::_Verbose>10)
817 std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
818 //see collection.prepareFieldDescriptions()
819 for (int ifield=0; ifield<iFieldMax; ifield++)
821 std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
822 if (descriptionField.find(nameTo)==std::string::npos)
823 continue; //only nameTo accepted in Fields name description
825 for (int inew=0; inew<inewMax; inew++)
827 for (int iold=0; iold<ioldMax; iold++)
829 //sending arrays for distant domains
830 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
832 int target=_domain_selector->getProcessorID(inew);
833 MEDCoupling::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
834 if (MyGlobals::_Verbose>10)
835 std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
836 SendDataArrayDouble(field, target);
838 //receiving arrays from distant domains
839 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
841 int source=_domain_selector->getProcessorID(iold);
843 if (MyGlobals::_Verbose>10)
844 std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
845 MEDCoupling::DataArrayDouble* field=RecvDataArrayDouble(source);
846 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
851 //local contributions and aggregation
852 for (int inew=0; inew<inewMax; inew++)
854 for (int iold=0; iold<ioldMax; iold++)
855 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
857 MEDCoupling::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
858 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
864 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
865 MEDCoupling::DataArrayDouble* fromArray,
866 std::string nameArrayTo,
867 std::string descriptionField)
868 //here we use 'cellFamily_ccI inew iold' set in remapIntField
870 if (nameArrayTo!="cellFieldDouble")
871 throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
872 std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
874 std::map<std::string,MEDCoupling::DataArrayInt*>::iterator it1;
875 it1=_map_dataarray_int.find(key);
876 if (it1==_map_dataarray_int.end())
878 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
879 std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
882 //create ccI in remapIntField
883 MEDCoupling::DataArrayInt *ccI=it1->second;
884 if (MyGlobals::_Verbose>300)
885 std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
887 int nbcell=this->getMesh()[inew]->getNumberOfCells();
888 int nbcomp=fromArray->getNumberOfComponents();
889 int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
891 std::string tag="inewFieldDouble="+IntToStr(inew);
892 key=descriptionField+SerializeFromString(tag);
893 int fromArrayNbOfElem=fromArray->getNbOfElems();
894 int fromArrayNbOfComp=fromArray->getNumberOfComponents();
895 int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
897 if (MyGlobals::_Verbose>1000)
899 std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
900 " fromArray nbOfElems " << fromArrayNbOfElem <<
901 " nbTuples " << fromArray->getNumberOfTuples() <<
902 " nbcells " << fromArrayNbOfCell <<
903 " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
906 MEDCoupling::DataArrayDouble* field=0;
907 std::map<std::string,MEDCoupling::DataArrayDouble*>::iterator it2;
908 it2=_map_dataarray_double.find(key);
909 if (it2==_map_dataarray_double.end())
911 if (MyGlobals::_Verbose>300)
912 std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
913 field=MEDCoupling::DataArrayDouble::New();
914 _map_dataarray_double[key]=field;
915 field->alloc(nbcell*nbPtGauss,nbcomp);
916 field->fillWithZero();
921 if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
923 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
924 " trying remap of field double on cells : \n" << descriptionField << std::endl;
931 field->setPartOfValuesAdv(fromArray,ccI);
935 //replaced by setPartOfValuesAdv if nbPtGauss==1
936 int iMax=ccI->getNbOfElems();
937 int* pccI=ccI->getPointer();
938 double* pField=field->getPointer();
939 double* pFrom=fromArray->getPointer();
940 int itarget, isource, delta=nbPtGauss*nbcomp;
941 for (int i=0; i<iMax; i=i+2) //cell
945 if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
946 throw INTERP_KERNEL::Exception("Error field override");
947 int ita=itarget*delta;
948 int iso=isource*delta;
949 for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
956 using namespace MEDCoupling;
957 //================================================================================
959 * \brief Sort correspondence ids of one domain and permute ids of the other accordingly
960 * \param [in,out] ids1 - ids of one domain
961 * \param [in,out] ids2 - ids of another domain
962 * \param [in] delta - a delta to change all ids
963 * \param [in] removeEqual - to remove equal ids
964 * \return DataArrayInt* - array of ids joined back
966 //================================================================================
968 DataArrayInt* sortCorrespondences( DataArrayInt* ids1,
971 bool removeEqual = false)
974 MCAuto< DataArrayInt > renumN2O = ids1->buildPermArrPerLevel();
975 ids1->renumberInPlaceR( renumN2O->begin() );
976 ids2->renumberInPlaceR( renumN2O->begin() );
980 ids1 = ids1->buildUnique();
981 ids2 = ids2->buildUnique();
985 int * id = ids1->getPointer();
986 for ( ; id < ids1->end(); ++id )
988 id = ids2->getPointer();
989 for ( ; id < ids2->end(); ++id )
994 DataArrayInt* ids12 = DataArrayInt::Meld( ids1, ids2 ); // two components
995 ids12->rearrange( 1 ); // make one component
999 //================================================================================
1001 * \brief Renumber ids according to mesh->sortCellsInMEDFileFrmt()
1002 * \param [in,out] ids - cell ids to renumber
1003 * \param [in] o2nRenumber - renumbering array in "Old to New" mode
1005 //================================================================================
1007 void renumber( DataArrayInt* ids, const DataArrayInt* o2nRenumber )
1009 if ( !ids || !o2nRenumber )
1011 int * id = ids->getPointer();
1012 const int * o2n = o2nRenumber->getConstPointer();
1013 for ( ; id < ids->end(); ++id )
1020 //================================================================================
1022 * \brief Fill up ConnectZone's stored in _topology with nodal correspondences
1023 * \param [in] nodeMapping - mapping between old nodes and new nodes
1024 * (iolddomain,ioldnode)->(inewdomain,inewnode)
1025 * \param [in] o2nRenumber - renumbering array returned by mesh->sortCellsInMEDFileFrmt()
1027 * \param [in] nbInitialDomains - nb of old domains
1029 //================================================================================
1031 void MEDPARTITIONER::MeshCollection::buildConnectZones( const NodeMapping& nodeMapping,
1032 const std::vector<MEDCoupling::DataArrayInt*> & o2nRenumber,
1033 int nbInitialDomains)
1035 if ( !MyGlobals::_Create_Joints || _topology->nbDomain() < 2 )
1038 if ( MyGlobals::_World_Size > 1 )
1040 _topology->getCZ().clear();
1041 return; // not implemented for parallel mode
1044 // at construction, _topology creates cell correspondences basing on Graph information,
1046 // 1) add node correspondences,
1047 // 2) split cell correspondences by cell geometry types
1048 // 3) sort ids to be in ascending order
1050 const int nb_domains = _topology->nbDomain();
1052 // ==================================
1053 // 1) add node correspondences
1054 // ==================================
1056 std::vector< std::vector< std::vector< int > > > nodeCorresp( nb_domains );
1057 for ( int idomain = 0; idomain < nb_domains; ++idomain )
1059 nodeCorresp[ idomain ].resize( nb_domains );
1062 NodeMapping::const_iterator nmIt1, nmIt2 = nodeMapping.begin();
1063 for ( nmIt1 = nmIt2; nmIt1 != nodeMapping.end(); nmIt1 = nmIt2 )
1065 // look for an "old" node mapped into several "new" nodes in different domains
1067 while ( ++nmIt2 != nodeMapping.end() && nmIt2->first == nmIt1->first )
1068 nbSameOld += ( nmIt2->second != nmIt1->second );
1070 if ( nbSameOld > 0 )
1072 NodeMapping::const_iterator nmEnd = nmIt2;
1073 for ( ; true; ++nmIt1 )
1076 if ( ++nmIt2 == nmEnd )
1078 int dom1 = nmIt1->second.first;
1079 int node1 = nmIt1->second.second;
1080 for ( ; nmIt2 != nmEnd; ++nmIt2 )
1082 int dom2 = nmIt2->second.first;
1083 int node2 = nmIt2->second.second;
1086 nodeCorresp[ dom1 ][ dom2 ].push_back( node1 );
1087 nodeCorresp[ dom1 ][ dom2 ].push_back( node2 );
1088 nodeCorresp[ dom2 ][ dom1 ].push_back( node2 );
1089 nodeCorresp[ dom2 ][ dom1 ].push_back( node1 );
1096 // add nodeCorresp to czVec
1098 std::vector<MEDPARTITIONER::ConnectZone*>& czVec = _topology->getCZ();
1100 for ( int idomain = 0; idomain < nb_domains; ++idomain )
1102 for ( int idomainNear = 0; idomainNear < nb_domains; ++idomainNear )
1104 std::vector< int > & corresp = nodeCorresp[ idomain ][ idomainNear ];
1105 if ( corresp.empty() )
1108 MEDPARTITIONER::ConnectZone* cz = 0;
1109 for ( size_t i = 0; i < czVec.size() && !cz; ++i )
1111 czVec[i]->getLocalDomainNumber () == idomain &&
1112 czVec[i]->getDistantDomainNumber() == idomainNear )
1117 cz = new MEDPARTITIONER::ConnectZone();
1118 cz->setName( "Nodal Connect Zone defined by MEDPARTITIONER" );
1119 cz->setLocalDomainNumber ( idomain );
1120 cz->setDistantDomainNumber( idomainNear );
1121 czVec.push_back(cz);
1124 cz->setNodeCorresp( &corresp[0], corresp.size()/2 );
1128 // ==========================================================
1129 // 2) split cell correspondences by cell geometry types
1130 // ==========================================================
1132 for ( size_t i = 0; i < czVec.size(); ++i )
1134 MEDPARTITIONER::ConnectZone* cz = czVec[i];
1136 cz->getEntityCorrespNumber( 0,0 ) == 0 ||
1137 cz->getLocalDomainNumber () > (int)_mesh.size() ||
1138 cz->getDistantDomainNumber() > (int)_mesh.size() )
1140 MEDCoupling::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber () ];
1141 MEDCoupling::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ];
1143 // separate ids of two domains
1144 const MEDCoupling::MEDCouplingSkyLineArray *corrArray = cz->getEntityCorresp( 0, 0 );
1145 const DataArrayInt* ids12 = corrArray->getValuesArray();
1146 MCAuto<DataArrayInt> ids1, ids2, ids12Sorted;
1147 ids1 = ids12->selectByTupleIdSafeSlice( 0, corrArray->getLength(), 2 );
1148 ids2 = ids12->selectByTupleIdSafeSlice( 1, corrArray->getLength(), 2 );
1150 // renumber cells according to mesh->sortCellsInMEDFileFrmt()
1151 renumber( ids1, o2nRenumber[ cz->getLocalDomainNumber() ]);
1152 renumber( ids2, o2nRenumber[ cz->getDistantDomainNumber() ]);
1154 // check nb cell types
1155 std::set<INTERP_KERNEL::NormalizedCellType> types1, types2;
1156 types1 = mesh1->getTypesOfPart( ids1->begin(), ids1->end() );
1157 types2 = mesh2->getTypesOfPart( ids2->begin(), ids2->end() );
1158 if ( types1.size() < 1 || types2.size() < 1 )
1159 continue; // parallel mode?
1161 MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1162 for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1164 czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
1165 czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1168 if ( types1.size() == 1 && types2.size() == 1 ) // split not needed, only sort
1170 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1171 cz->setEntityCorresp( *types1.begin(), *types2.begin(),
1172 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1174 if ( cz21 )// set 2->1 correspondence
1176 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1177 cz21->setEntityCorresp( *types2.begin(), *types1.begin(),
1178 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1181 else // split and sort
1183 typedef std::pair< std::vector< int >, std::vector< int > > T2Vecs;
1184 T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ];
1187 const int nbIds = ids1->getNbOfElems();
1188 const int * p1 = ids1->begin(), * p2 = ids2->begin();
1189 for ( int i = 0; i < nbIds; ++i )
1191 t1 = mesh1->getTypeOfCell( p1[ i ]);
1192 t2 = mesh2->getTypeOfCell( p2[ i ]);
1193 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1194 ids.first .push_back( p1[ i ]);
1195 ids.second.push_back( p1[ i ]);
1198 const int maxType = int( INTERP_KERNEL::NORM_MAXTYPE );
1199 for ( t1 = 0; t1 < maxType; ++t1 )
1200 for ( t2 = 0; t2 < maxType; ++t2 )
1202 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1203 if ( ids.first.empty() ) continue;
1204 p1 = & ids.first[0];
1205 p2 = & ids.second[0];
1206 ids1->desallocate();
1207 ids1->pushBackValsSilent( p1, p1+ids.first.size() );
1208 ids2->desallocate();
1209 ids2->pushBackValsSilent( p2, p2+ids.first.size() );
1210 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1212 cz->setEntityCorresp( t1, t2,
1213 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1215 if ( cz21 )// set 2->1 correspondence
1217 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1218 cz21->setEntityCorresp( t2, t1,
1219 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1225 cz->setEntityCorresp( 0, 0, 0, 0 ); // erase ids computed by _topology
1227 cz21->setEntityCorresp( 0, 0, 0, 0 );
1232 // ==========================================
1233 // 3) sort node ids to be in ascending order
1234 // ==========================================
1236 const bool removeEqual = ( nbInitialDomains > 1 );
1238 for ( size_t i = 0; i < czVec.size(); ++i )
1240 MEDPARTITIONER::ConnectZone* cz = czVec[i];
1241 if ( !cz || cz->getNodeNumber() < 1 )
1243 if ( cz->getDistantDomainNumber() < cz->getLocalDomainNumber() )
1244 continue; // treat a pair of domains once
1246 MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1247 for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1249 czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
1250 czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1253 // separate ids of two domains
1254 const MEDCoupling::MEDCouplingSkyLineArray *corrArray = cz->getNodeCorresp();
1255 const DataArrayInt *ids12 = corrArray->getValuesArray();
1256 MCAuto<DataArrayInt> ids1, ids2, ids12Sorted;
1257 ids1 = ids12->selectByTupleIdSafeSlice( 0, corrArray->getLength(), 2 );
1258 ids2 = ids12->selectByTupleIdSafeSlice( 1, corrArray->getLength(), 2 );
1260 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/0, removeEqual );
1261 cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1263 if ( cz21 )// set 2->1 correspondence
1265 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0, false );
1266 cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1271 //================================================================================
1273 * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
1274 * group (where "n" and "p" are domain IDs)
1276 //================================================================================
1278 void MEDPARTITIONER::MeshCollection::buildBoundaryFaces()
1280 if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates )
1283 if ( getMeshDimension() < 2 )
1286 using MEDCoupling::MEDCouplingUMesh;
1287 using MEDCoupling::DataArrayDouble;
1288 using MEDCoupling::DataArrayInt;
1290 std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
1291 int nbMeshes = faceMeshes.size();
1293 //preparing bounding box trees for accelerating search of coincident faces
1294 std::vector<BBTreeOfDim* > bbTrees(nbMeshes);
1295 std::vector<DataArrayDouble*>bbox (nbMeshes);
1296 for (int inew = 0; inew < nbMeshes-1; inew++)
1297 if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
1299 DataArrayDouble* bcCoords = faceMeshes[inew]->computeCellCenterOfMass();
1300 bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
1301 bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
1302 bbox[inew]->getConstPointer(),0,0,
1303 bbox[inew]->getNumberOfTuples());
1304 bcCoords->decrRef();
1307 // loop on domains to find joint faces between them
1308 for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
1310 for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
1312 MEDCouplingUMesh* mesh1 = 0;
1313 MEDCouplingUMesh* mesh2 = 0;
1314 //MEDCouplingUMesh* recvMesh = 0;
1315 bool mesh1Here = true, mesh2Here = true;
1316 if (isParallelMode())
1319 mesh1Here = _domain_selector->isMyDomain(inew1);
1320 mesh2Here = _domain_selector->isMyDomain(inew2);
1321 if ( !mesh1Here && mesh2Here )
1323 //send mesh2 to domain of mesh1
1324 _domain_selector->sendMesh(*faceMeshes[inew2],
1325 _domain_selector->getProcessorID(inew1));
1327 else if ( mesh1Here && !mesh2Here )
1329 //receiving mesh2 from a distant domain
1330 _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
1331 if ( faceMeshes[ inew2 ] )
1332 faceMeshes[ inew2 ]->decrRef();
1333 faceMeshes[ inew2 ] = mesh2;
1337 if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
1338 if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
1340 // find coincident faces
1341 std::vector< int > faces1, faces2;
1342 if ( mesh1 && mesh2 )
1344 const DataArrayDouble* coords2 = mesh2->computeCellCenterOfMass();
1345 const double* c2 = coords2->getConstPointer();
1346 const int dim = coords2->getNumberOfComponents();
1347 const int nbFaces2 = mesh2->getNumberOfCells();
1348 const int nbFaces1 = mesh1->getNumberOfCells();
1350 for (int i2 = 0; i2 < nbFaces2; i2++)
1352 std::vector<int> coincFaces;
1353 bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
1354 if (coincFaces.size()!=0)
1356 int i1 = coincFaces[0];
1357 // if ( coincFaces.size() > 1 )
1359 // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
1360 // i1 = coincFaces[ i2nb->second++ ];
1362 if ( i1 < nbFaces1 ) // protection from invalid elements
1364 faces1.push_back( i1 );
1365 faces2.push_back( i2 );
1372 if ( isParallelMode())
1375 if ( mesh1Here && !mesh2Here )
1377 //send faces2 to domain of recvMesh
1378 SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
1380 else if ( !mesh1Here && mesh2Here )
1382 //receiving ids of faces from a domain of mesh1
1383 RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
1388 // recvMesh->decrRef();
1390 // Create group "JOINT_inew1_inew2_Faces" and corresponding families
1391 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1393 createJointGroup( is2nd ? faces2 : faces1,
1394 inew1 , inew2, is2nd );
1397 } // loop on the 2nd domains (inew2)
1398 } // loop on the 1st domains (inew1)
1401 // delete bounding box trees
1402 for (int inew = 0; inew < nbMeshes-1; inew++)
1403 if (isParallelMode() && _domain_selector->isMyDomain(inew))
1405 bbox[inew]->decrRef();
1406 delete bbTrees[inew];
1410 //================================================================================
1412 * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
1413 * \param faces - face ids to include into the group
1414 * \param inew1 - index of the 1st domain
1415 * \param inew2 - index of the 2nd domain
1416 * \param is2nd - in which (1st or 2nd) domain to create the group
1418 //================================================================================
1420 void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
1425 // get the name of JOINT group
1426 std::string groupName;
1428 std::ostringstream oss;
1430 << (is2nd ? inew2 : inew1 ) << "_"
1431 << (is2nd ? inew1 : inew2 ) << "_"
1432 << ( getMeshDimension()==2 ? "Edge" : "Face" );
1433 groupName = oss.str();
1436 // remove existing "JOINT_*" group
1437 _group_info.erase( groupName );
1439 // get family IDs array
1441 int inew = (is2nd ? inew2 : inew1 );
1442 int totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
1443 std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
1444 if ( !_map_dataarray_int.count(cle) )
1446 if ( totalNbFaces > 0 )
1448 MEDCoupling::DataArrayInt* p=MEDCoupling::DataArrayInt::New();
1449 p->alloc( totalNbFaces, 1 );
1451 famIDs = p->getPointer();
1452 _map_dataarray_int[cle]=p;
1457 famIDs = _map_dataarray_int.find(cle)->second->getPointer();
1459 // find a family ID of an existing JOINT group
1461 std::map<std::string, int>::iterator name2id = _family_info.find( groupName );
1462 if ( name2id != _family_info.end() )
1463 familyID = name2id->second;
1465 // remove faces from the familyID-the family
1466 if ( familyID != 0 && famIDs )
1467 for ( int i = 0; i < totalNbFaces; ++i )
1468 if ( famIDs[i] == familyID )
1471 if ( faces.empty() )
1474 if ( familyID == 0 ) // generate a family ID for JOINT group
1476 std::set< int > familyIDs;
1477 for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
1478 familyIDs.insert( name2id->second );
1479 // find the next free family ID
1480 int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
1483 if ( !familyIDs.count( ++familyID ))
1486 while ( freeIdCount > 0 );
1489 // push faces to familyID-th group
1490 if ( faces.back() >= totalNbFaces )
1491 throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
1492 for ( size_t i = 0; i < faces.size(); ++i )
1493 famIDs[ faces[i] ] = familyID;
1495 // register JOINT group and family
1496 _family_info[ groupName ] = familyID; // name of the group and family is same
1497 _group_info [ groupName ].push_back( groupName );
1500 /*! constructing the MESH collection from a distributed file
1502 * \param filename name of the master file containing the list of all the MED files
1504 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
1506 _owns_topology(true),
1508 _domain_selector( 0 ),
1509 _i_non_empty_mesh(-1),
1510 _driver_type(MEDPARTITIONER::Undefined),
1511 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1512 _family_splitting(false),
1513 _create_empty_groups(false),
1518 _driver=new MeshCollectionMedXmlDriver(this);
1519 _driver->read (filename.c_str());
1520 _driver_type = MedXml;
1523 { // Handle all exceptions
1524 if ( _driver ) delete _driver; _driver=0;
1527 _driver=new MeshCollectionMedAsciiDriver(this);
1528 _driver->read (filename.c_str());
1529 _driver_type=MedAscii;
1535 throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
1538 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1539 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1540 _i_non_empty_mesh = idomain;
1543 /*! Constructing the MESH collection from selected domains of a distributed file
1545 * \param filename - name of the master file containing the list of all the MED files
1546 * \param domainSelector - selector of domains to load
1548 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
1550 _owns_topology(true),
1552 _domain_selector( &domainSelector ),
1553 _i_non_empty_mesh(-1),
1554 _driver_type(MEDPARTITIONER::Undefined),
1555 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1556 _family_splitting(false),
1557 _create_empty_groups(false),
1560 std::string myfile=filename;
1561 std::size_t found=myfile.find(".xml");
1562 if (found!=std::string::npos) //file .xml
1566 _driver=new MeshCollectionMedXmlDriver(this);
1567 _driver->read ( filename.c_str(), _domain_selector );
1568 _driver_type = MedXml;
1571 { // Handle all exceptions
1573 throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
1578 found=myfile.find(".med");
1579 if (found!=std::string::npos) //file .med single
1581 //make a temporary file .xml and retry MedXmlDriver
1583 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
1585 <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
1586 <description what=\"\" when=\"\"/>\n \
1588 <mesh name=\"$meshName\"/>\n \
1591 <subdomain number=\"1\"/>\n \
1592 <global_numbering present=\"no\"/>\n \
1595 <subfile id=\"1\">\n \
1596 <name>$fileName</name>\n \
1597 <machine>localhost</machine>\n \
1601 <mesh name=\"$meshName\">\n \
1602 <chunk subdomain=\"1\">\n \
1603 <name>$meshName</name>\n \
1608 std::vector<std::string> meshNames=GetMeshNames(myfile);
1609 xml.replace(xml.find("$fileName"),9,myfile);
1610 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1611 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1612 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1613 std::string nameFileXml(myfile);
1614 nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1615 std::string nameFileXmlDN,nameFileXmlBN;
1616 MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1617 nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1618 if (_domain_selector->rank()==0) //only on to write it
1620 std::ofstream f(nameFileXml.c_str());
1625 if (MyGlobals::_World_Size>1)
1626 MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1630 _driver=new MeshCollectionMedXmlDriver(this);
1631 _driver->read ( nameFileXml.c_str(), _domain_selector );
1632 _driver_type = MedXml;
1635 { // Handle all exceptions
1636 delete _driver; _driver=0;
1637 throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1644 _driver=new MeshCollectionMedAsciiDriver(this);
1645 _driver->read ( filename.c_str(), _domain_selector );
1646 _driver_type=MedAscii;
1652 throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1656 // find non-empty domain mesh
1657 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1658 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1659 _i_non_empty_mesh = idomain;
1663 //check for all proc/file compatibility of _field_descriptions
1665 _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1667 _field_descriptions=MyGlobals::_Field_Descriptions;
1670 catch(INTERP_KERNEL::Exception& e)
1672 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1673 throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1678 //check for all proc/file compatibility of _family_info
1679 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1680 _family_info=DevectorizeToMapOfStringInt(v2);
1682 catch(INTERP_KERNEL::Exception& e)
1684 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1685 throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1690 //check for all proc/file compatibility of _group_info
1691 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1692 _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1694 catch(INTERP_KERNEL::Exception& e)
1696 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1697 throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1702 /*! constructing the MESH collection from a sequential MED-file
1704 * \param filename MED file
1705 * \param meshname name of the mesh that is to be read
1707 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1709 _owns_topology(true),
1711 _domain_selector( 0 ),
1712 _i_non_empty_mesh(-1),
1714 _driver_type(MEDPARTITIONER::MedXml),
1715 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1716 _family_splitting(false),
1717 _create_empty_groups(false),
1720 try // avoid memory leak in case of inexistent filename
1722 retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1728 throw INTERP_KERNEL::Exception("problem reading .med files");
1730 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1731 _i_non_empty_mesh = 0;
1734 MEDPARTITIONER::MeshCollection::~MeshCollection()
1736 for (int i=0; i<(int)_mesh.size();i++)
1737 if (_mesh[i]!=0) _mesh[i]->decrRef();
1739 for (int i=0; i<(int)_cell_family_ids.size();i++)
1740 if (_cell_family_ids[i]!=0)
1741 _cell_family_ids[i]->decrRef();
1743 for (int i=0; i<(int)_face_mesh.size();i++)
1744 if (_face_mesh[i]!=0)
1745 _face_mesh[i]->decrRef();
1747 for (int i=0; i<(int)_face_family_ids.size();i++)
1748 if (_face_family_ids[i]!=0)
1749 _face_family_ids[i]->decrRef();
1751 for (std::map<std::string, MEDCoupling::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1752 if ((*it).second!=0)
1753 (*it).second->decrRef();
1755 for (std::map<std::string, MEDCoupling::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1756 if ((*it).second!=0)
1757 (*it).second->decrRef();
1760 if (_topology!=0 && _owns_topology)
1763 delete _joint_finder;
1767 /*! constructing the MESH collection from a file
1769 * The method creates as many MED-files as there are domains in the
1770 * collection. It also creates a master file that lists all the MED files.
1771 * The MED files created in ths manner contain joints that describe the
1772 * connectivity between subdomains.
1774 * \param filename name of the master file that will contain the list of the MED files
1777 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1779 //suppresses link with driver so that it can be changed for writing
1782 retrieveDriver()->write (filename.c_str(), _domain_selector);
1785 /*! creates or gets the link to the collection driver
1787 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1791 switch(_driver_type)
1794 _driver=new MeshCollectionMedXmlDriver(this);
1797 _driver=new MeshCollectionMedAsciiDriver(this);
1800 throw INTERP_KERNEL::Exception("Unrecognized driver");
1807 /*! gets an existing driver
1810 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1816 * retrieves the mesh dimension
1818 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1820 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1823 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1826 for (size_t i=0; i<_mesh.size(); i++)
1833 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1836 for (size_t i=0; i<_mesh.size(); i++)
1838 if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1843 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1846 for (size_t i=0; i<_face_mesh.size(); i++)
1848 if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1853 std::vector<MEDCoupling::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1858 std::vector<MEDCoupling::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1863 MEDCoupling::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1865 return _mesh[idomain];
1868 MEDCoupling::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1870 return _face_mesh[idomain];
1873 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1876 return _topology->getCZ();
1878 static std::vector<MEDPARTITIONER::ConnectZone*> noCZ;
1882 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1887 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo, bool takeOwneship)
1891 throw INTERP_KERNEL::Exception("topology is already set");
1896 _owns_topology = takeOwneship;
1900 /*! Method creating the cell graph in serial mode
1902 * \param array returns the pointer to the structure that contains the graph
1903 * \param edgeweight returns the pointer to the table that contains the edgeweights
1904 * (only used if indivisible regions are required)
1906 void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1911 using std::make_pair;
1914 if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
1915 const MEDCoupling::MEDCouplingUMesh* mesh=_mesh[0];
1916 if (MyGlobals::_Verbose>50)
1917 std::cout<<"getting nodal connectivity"<<std::endl;
1919 //looking for reverse nodal connectivity i global numbering
1920 if (isParallelMode() && !_domain_selector->isMyDomain(0))
1923 vector<int> index(1,0);
1925 array = MEDCoupling::MEDCouplingSkyLineArray::New(index,value);
1928 array=mesh->generateGraph();
1930 /*! Method creating the cell graph in multidomain mode
1932 * \param array returns the pointer to the structure that contains the graph
1933 * \param edgeweight returns the pointer to the table that contains the edgeweights
1934 * (only used if indivisible regions are required)
1936 void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(MEDCoupling::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1938 using std::multimap;
1940 using std::make_pair;
1943 std::multimap< int, int > node2cell;
1944 std::map< pair<int,int>, int > cell2cellcounter;
1945 std::multimap<int,int> cell2cell;
1947 std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
1948 int nbdomain=_topology->nbDomain();
1950 if (isParallelMode())
1952 _joint_finder=new JointFinder(*this);
1953 _joint_finder->findCommonDistantNodes();
1954 commonDistantNodes=_joint_finder->getDistantNodeCell();
1957 if (MyGlobals::_Verbose>500)
1958 _joint_finder->print();
1961 if (MyGlobals::_Verbose>50)
1962 std::cout<<"getting nodal connectivity"<<std::endl;
1963 //looking for reverse nodal connectivity i global numbering
1965 for (int idomain=0; idomain<nbdomain; idomain++)
1967 if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1969 meshDim = _mesh[idomain]->getMeshDimension();
1971 MEDCoupling::DataArrayInt* index=MEDCoupling::DataArrayInt::New();
1972 MEDCoupling::DataArrayInt* revConn=MEDCoupling::DataArrayInt::New();
1973 int nbNodes=_mesh[idomain]->getNumberOfNodes();
1974 _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1975 //problem saturation over 1 000 000 nodes for 1 proc
1976 if (MyGlobals::_Verbose>100)
1977 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1978 int* index_ptr=index->getPointer();
1979 int* revConnPtr=revConn->getPointer();
1980 for (int i=0; i<nbNodes; i++)
1982 for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
1984 int globalNode=_topology->convertNodeToGlobal(idomain,i);
1985 int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
1986 node2cell.insert(make_pair(globalNode, globalCell));
1992 for (int iother=0; iother<nbdomain; iother++)
1994 std::multimap<int,int>::iterator it;
1995 int isource=idomain;
1997 for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
1998 it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
2000 int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
2001 int globalCell=(*it).second;
2002 node2cell.insert(make_pair(globalNode, globalCell));
2008 //creating graph arcs (cell to cell relations)
2009 //arcs are stored in terms of (index,value) notation
2012 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
2013 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
2015 //warning here one node have less than or equal effective number of cell with it
2016 //but cell could have more than effective nodes
2017 //because other equals nodes in other domain (with other global inode)
2018 if (MyGlobals::_Verbose>50)
2019 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
2021 for (int inode=0;inode<_topology->nbNodes();inode++)
2023 typedef multimap<int,int>::const_iterator MI;
2024 std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
2025 for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
2026 for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
2028 int icell1=cell1->second;
2029 int icell2=cell2->second;
2030 if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
2031 std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2032 if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2033 else (it->second)++;
2036 // for (int icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
2038 // typedef multimap<int,int>::const_iterator MI;
2039 // std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
2040 // for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
2042 // std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
2043 // for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
2045 // int icell2=cell2->second;
2046 // std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2047 // if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2048 // else (it->second)++;
2054 //converting the counter to a multimap structure
2055 for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
2056 it!=cell2cellcounter.end();
2058 if (it->second>=meshDim)
2060 cell2cell.insert(std::make_pair(it->first.first,it->first.second));
2061 cell2cell.insert(std::make_pair(it->first.second, it->first.first));
2065 if (MyGlobals::_Verbose>50)
2066 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
2067 //filling up index and value to create skylinearray structure
2068 std::vector <int> index,value;
2072 for (int idomain=0; idomain<nbdomain; idomain++)
2074 if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
2075 int nbCells=_mesh[idomain]->getNumberOfCells();
2076 for (int icell=0; icell<nbCells; icell++)
2079 int globalCell=_topology->convertCellToGlobal(idomain,icell);
2080 multimap<int,int>::iterator it;
2081 pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
2082 ret=cell2cell.equal_range(globalCell);
2083 for (it=ret.first; it!=ret.second; ++it)
2085 int ival=(*it).second; //no adding one existing yet
2086 for (int i=idep ; i<idep+size ; i++)
2095 value.push_back(ival);
2099 idep=index[index.size()-1]+size;
2100 index.push_back(idep);
2104 array=MEDCoupling::MEDCouplingSkyLineArray::New(index,value);
2106 if (MyGlobals::_Verbose>100)
2108 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
2109 index.size()-1 << " " << value.size() << std::endl;
2110 int max=index.size()>15?15:index.size();
2113 for (int i=0; i<max; ++i)
2114 std::cout<<index[i]<<" ";
2115 std::cout << "... " << index[index.size()-1] << std::endl;
2116 for (int i=0; i<max; ++i)
2117 std::cout<< value[i] << " ";
2118 int ll=index[index.size()-1]-1;
2119 std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
2126 /*! Creates the partition corresponding to the cell graph and the partition number
2128 * \param nbdomain number of subdomains for the newly created graph
2130 * returns a topology based on the new graph
2132 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
2133 Graph::splitter_type split,
2134 const std::string& options_string,
2135 int *user_edge_weights,
2136 int *user_vertices_weights)
2138 if (MyGlobals::_Verbose>10)
2139 std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
2142 throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
2143 MEDCoupling::MEDCouplingSkyLineArray* array=0;
2146 if (_topology->nbDomain()>1 || isParallelMode())
2147 buildParallelCellGraph(array,edgeweights);
2149 buildCellGraph(array,edgeweights);
2151 Graph* cellGraph = 0;
2155 if ( isParallelMode() && MyGlobals::_World_Size > 1 )
2157 #ifdef MED_ENABLE_PARMETIS
2158 if (MyGlobals::_Verbose>10)
2159 std::cout << "ParMETISGraph" << std::endl;
2160 cellGraph=new ParMETISGraph(array,edgeweights);
2165 #ifdef MED_ENABLE_METIS
2166 if (MyGlobals::_Verbose>10)
2167 std::cout << "METISGraph" << std::endl;
2168 cellGraph=new METISGraph(array,edgeweights);
2172 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
2176 #ifdef MED_ENABLE_SCOTCH
2177 if (MyGlobals::_Verbose>10)
2178 std::cout << "SCOTCHGraph" << std::endl;
2179 cellGraph=new SCOTCHGraph(array,edgeweights);
2181 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
2186 //!user-defined weights
2187 if (user_edge_weights!=0)
2188 cellGraph->setEdgesWeights(user_edge_weights);
2189 if (user_vertices_weights!=0)
2190 cellGraph->setVerticesWeights(user_vertices_weights);
2192 if (MyGlobals::_Is0verbose>10)
2193 std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
2194 cellGraph->partGraph(nbdomain, options_string, _domain_selector);
2196 if (MyGlobals::_Is0verbose>10)
2197 std::cout << "building new topology" << std::endl;
2198 //cellGraph is a shared pointer
2199 Topology *topology=0;
2200 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2202 delete [] edgeweights;
2204 if (MyGlobals::_Verbose>11)
2205 std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
2209 /*! Creates a topology for a partition specified by the user
2211 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
2213 * returns a topology based on the new partition
2215 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
2217 MEDCoupling::MEDCouplingSkyLineArray* array=0;
2220 if ( _topology->nbDomain()>1)
2221 buildParallelCellGraph(array,edgeweights);
2223 buildCellGraph(array,edgeweights);
2226 std::set<int> domains;
2227 for (int i=0; i<_topology->nbCells(); i++)
2229 domains.insert(partition[i]);
2231 cellGraph=new UserGraph(array, partition, _topology->nbCells());
2233 //cellGraph is a shared pointer
2234 Topology *topology=0;
2235 int nbdomain=domains.size();
2236 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2237 // if (array!=0) delete array;
2242 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
2244 for (int i=0; i<_topology->nbDomain(); i++)
2246 std::ostringstream oss;
2248 if (!isParallelMode() || _domain_selector->isMyDomain(i))
2249 _mesh[i]->setName(oss.str());
2253 MEDCoupling::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
2254 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
2255 //something like MCAuto<MEDCouplingFieldDouble> f2=ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
2257 int rank=MyGlobals::_Rank;
2258 std::string tag="ioldFieldDouble="+IntToStr(iold);
2259 std::string descriptionIold=descriptionField+SerializeFromString(tag);
2260 if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
2262 if (MyGlobals::_Verbose>300)
2263 std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
2264 MEDCoupling::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
2267 if (MyGlobals::_Verbose>200)
2268 std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
2269 std::string description, fileName, meshName, fieldName;
2270 int typeField, DT, IT, entity;
2271 fileName=MyGlobals::_File_Names[iold];
2272 if (MyGlobals::_Verbose>10)
2273 std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
2274 FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
2275 meshName=MyGlobals::_Mesh_Names[iold];
2277 MEDCoupling::MCAuto<MEDCoupling::MEDCouplingFieldDouble> f2=ReadField((MEDCoupling::TypeOfField) typeField,
2278 fileName, meshName, 0, fieldName, DT, IT);
2280 MEDCoupling::DataArrayDouble* res=f2->getArray();
2281 //to know names of components
2282 std::vector<std::string> browse=BrowseFieldDouble(f2);
2283 std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
2284 if (MyGlobals::_Verbose>10)
2285 std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
2286 MyGlobals::_General_Informations.push_back(localFieldInformation);
2288 _map_dataarray_double[descriptionIold]=res;
2292 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
2293 //to have unique valid fields names/pointers/descriptions for partitionning
2294 //filter _field_descriptions to be in all procs compliant and equal
2296 int nbfiles=MyGlobals::_File_Names.size(); //nb domains
2299 nbfiles=_topology->nbDomain();
2301 std::vector<std::string> r2;
2302 //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
2303 for (int i=0; i<(int)_field_descriptions.size(); i++)
2305 std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
2306 for (int ii=0; ii<(int)r1.size(); ii++)
2307 r2.push_back(r1[ii]);
2309 //here vector(procs*fields) of serialised vector(description) data
2310 _field_descriptions=r2;
2311 int nbfields=_field_descriptions.size(); //on all domains
2312 if ((nbfields%nbfiles)!=0)
2314 if (MyGlobals::_Rank==0)
2316 std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
2317 << "fileMedNames :" << std::endl
2318 << ReprVectorOfString(MyGlobals::_File_Names)
2319 << "field_descriptions :" << std::endl
2320 << ReprVectorOfString(MyGlobals::_Field_Descriptions);
2322 throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
2324 _field_descriptions.resize(nbfields/nbfiles);
2325 for (int i=0; i<(int)_field_descriptions.size(); i++)
2327 std::string str=_field_descriptions[i];
2328 str=EraseTagSerialized(str,"idomain=");
2329 str=EraseTagSerialized(str,"fileName=");
2330 _field_descriptions[i]=str;
2334 //returns true if inodes of a face are in inodes of a cell
2335 bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >& inodesCell)
2338 int nbok=inodesFace.size();
2339 for (int i=0; i<nbok; i++)
2341 int ii=inodesFace[i];
2343 std::cout << "isFaceOncell problem inodeface<0" << std::endl;
2344 for (int j=0; j<(int)inodesCell.size(); j++)
2346 if (ii==inodesCell[j])
2349 break; //inode of face found
2353 break; //inode of face not found do not continue...
2355 return (ires==nbok);
2358 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
2360 for (int inew=0; inew<_topology->nbDomain(); inew++)
2362 if (!isParallelMode() || _domain_selector->isMyDomain(inew))
2364 if (MyGlobals::_Verbose>200)
2365 std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
2366 MEDCoupling::MEDCouplingUMesh* mcel=_mesh[inew];
2367 MEDCoupling::MEDCouplingUMesh* mfac=_face_mesh[inew];
2369 //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
2370 std::vector<int> nodeIds;
2371 getNodeIds(*mcel, *mfac, nodeIds);
2372 if (nodeIds.size()==0)
2373 continue; //one empty mesh nothing to do
2375 MEDCoupling::DataArrayInt *revNodalCel=MEDCoupling::DataArrayInt::New();
2376 MEDCoupling::DataArrayInt *revNodalIndxCel=MEDCoupling::DataArrayInt::New();
2377 mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
2378 int *revC=revNodalCel->getPointer();
2379 int *revIndxC=revNodalIndxCel->getPointer();
2381 std::vector< int > faceOnCell;
2382 std::vector< int > faceNotOnCell;
2383 int nbface=mfac->getNumberOfCells();
2384 for (int iface=0; iface<nbface; iface++)
2387 std::vector< int > inodesFace;
2388 mfac->getNodeIdsOfCell(iface, inodesFace);
2389 int nbnodFace=inodesFace.size();
2390 if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
2391 continue; // invalid node ids
2392 //set inodesFace in mcel
2394 for (int i=0; i<nbnodFace; i++)
2395 nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
2396 if ( nbok != nbnodFace )
2398 int inod=inodesFace[0];
2401 std::cout << "filterFaceOnCell problem 1" << std::endl;
2404 int nbcell=revIndxC[inod+1]-revIndxC[inod];
2405 for (int j=0; j<nbcell; j++) //look for each cell with inod
2407 int icel=revC[revIndxC[inod]+j];
2408 std::vector< int > inodesCell;
2409 mcel->getNodeIdsOfCell(icel, inodesCell);
2410 ok=isFaceOncell(inodesFace, inodesCell);
2415 faceOnCell.push_back(iface);
2419 faceNotOnCell.push_back(iface);
2420 if (MyGlobals::_Is0verbose>300)
2421 std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
2425 revNodalCel->decrRef();
2426 revNodalIndxCel->decrRef();
2428 // std::string keyy;
2429 // keyy=Cle1ToStr("filterFaceOnCell",inew);
2430 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
2431 // keyy=Cle1ToStr("filterNotFaceOnCell",inew);
2432 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
2434 // filter the face mesh
2435 if ( faceOnCell.empty() )
2436 _face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
2438 _face_mesh[inew] = (MEDCoupling::MEDCouplingUMesh *)
2439 mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
2442 // filter the face families
2443 std::string key = Cle1ToStr("faceFamily_toArray",inew);
2444 if ( getMapDataArrayInt().count( key ))
2446 MEDCoupling::DataArrayInt * & fam = getMapDataArrayInt()[ key ];
2447 MEDCoupling::DataArrayInt * famFilter = MEDCoupling::DataArrayInt::New();
2448 famFilter->alloc(faceOnCell.size(),1);
2449 int* pfamFilter = famFilter->getPointer();
2450 int* pfam = fam->getPointer();
2451 for ( size_t i=0; i<faceOnCell.size(); i++ )
2452 pfamFilter[i]=pfam[faceOnCell[i]];