1 // Copyright (C) 2007-2015 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"
37 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
38 #include "MEDCouplingFieldDouble.hxx"
39 #include "MEDCouplingMemArray.hxx"
40 #include "MEDCouplingNormalizedUnstructuredMesh.hxx"
41 #include "MEDCouplingUMesh.hxx"
42 #include "MEDLoader.hxx"
43 #include "MEDLoaderBase.hxx"
49 #ifdef MED_ENABLE_PARMETIS
50 #include "MEDPARTITIONER_ParMetisGraph.hxx"
52 #ifdef MED_ENABLE_METIS
53 #include "MEDPARTITIONER_MetisGraph.hxx"
55 #ifdef MED_ENABLE_SCOTCH
56 #include "MEDPARTITIONER_ScotchGraph.hxx"
66 MEDPARTITIONER::MeshCollection::MeshCollection()
68 _owns_topology(false),
70 _domain_selector( 0 ),
71 _i_non_empty_mesh(-1),
72 _driver_type(MEDPARTITIONER::MedXml),
73 _subdomain_boundary_creates( MyGlobals::_Creates_Boundary_Faces ),
74 _family_splitting(false),
75 _create_empty_groups(false),
80 /*!constructor creating a new mesh collection (mesh series + topology)
81 *from an old collection and a new topology
83 * On output, the constructor has built the meshes corresponding to the new mesh collection.
84 * The new topology has been updated so that face and node mappings are included.
85 * The families have been cast to their projections in the new topology.
87 * \param initial_collection collection from which the data (coordinates, connectivity) are taken
88 * \param topology topology containing the cell mappings
91 MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection,
93 bool family_splitting,
94 bool create_empty_groups)
95 : _topology(topology),
96 _owns_topology(false),
98 _domain_selector( initialCollection._domain_selector ),
99 _i_non_empty_mesh(-1),
100 _name(initialCollection._name),
101 _driver_type(MEDPARTITIONER::MedXml),
102 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
103 _family_splitting(family_splitting),
104 _create_empty_groups(create_empty_groups),
107 std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
108 std::vector<ParaMEDMEM::DataArrayInt*> o2nRenumber;
110 castCellMeshes(initialCollection, new2oldIds, o2nRenumber );
112 //defining the name for the collection and the underlying meshes
113 setName(initialCollection.getName());
120 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
121 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
123 if (MyGlobals::_Is0verbose)
124 std::cout<<"treating faces"<<std::endl;
125 NodeMapping nodeMapping;
126 //nodeMapping contains the mapping between old nodes and new nodes
127 // (iolddomain,ioldnode)->(inewdomain,inewnode)
128 createNodeMapping(initialCollection, nodeMapping);
129 std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
130 castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
136 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
137 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
139 if (MyGlobals::_Is0verbose)
141 if (isParallelMode())
142 std::cout << "ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
144 std::cout << "NOT ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
146 if (MyGlobals::_Is0verbose>10)
147 std::cout<<"treating cell and face families"<<std::endl;
149 castIntField(initialCollection.getMesh(),
151 initialCollection.getCellFamilyIds(),
153 castIntField(initialCollection.getFaceMesh(),
155 initialCollection.getFaceFamilyIds(),
160 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
161 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
163 if (MyGlobals::_Is0verbose)
164 std::cout << "treating groups" << std::endl;
165 _family_info=initialCollection.getFamilyInfo();
166 _group_info=initialCollection.getGroupInfo();
169 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
170 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
172 if (MyGlobals::_Is0verbose)
173 std::cout << "treating fields" << std::endl;
174 castAllFields(initialCollection,"cellFieldDouble");
175 if (_i_non_empty_mesh<0)
177 for (size_t i=0; i<_mesh.size(); i++)
181 _i_non_empty_mesh=i; //first existing one local
187 // find faces common with neighbor domains and put them in groups
188 buildBoundaryFaces();
190 //building the connect zones necessary for writing joints
191 buildConnectZones( nodeMapping, o2nRenumber, initialCollection.getTopology()->nbDomain() );
193 // delete o2nRenumber
194 for ( size_t i = 0; i < o2nRenumber.size(); ++i )
195 if ( o2nRenumber[i] )
196 o2nRenumber[i]->decrRef();
200 Creates the meshes using the topology underlying he mesh collection and the mesh data
201 coming from the ancient collection
202 \param initialCollection collection from which the data is extracted to create the new meshes
203 \param [out] o2nRenumber returns for each new domain a permutation array returned by sortCellsInMEDFileFrmt()
206 void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
207 std::vector<std::vector<std::vector<int> > >& new2oldIds,
208 std::vector<ParaMEDMEM::DataArrayInt*> & o2nRenumber)
210 if (MyGlobals::_Verbose>10)
211 std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
213 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes");
215 int nbNewDomain=_topology->nbDomain();
216 int nbOldDomain=initialCollection.getTopology()->nbDomain();
218 _mesh.resize(nbNewDomain);
219 o2nRenumber.resize(nbNewDomain,0);
220 int rank=MyGlobals::_Rank;
221 //splitting the initial domains into smaller bits
222 std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
223 splitMeshes.resize(nbNewDomain);
224 for (int inew=0; inew<nbNewDomain; inew++)
226 splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
229 for (int iold=0; iold<nbOldDomain; iold++)
231 if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
233 int size=(initialCollection._mesh)[iold]->getNumberOfCells();
234 std::vector<int> globalids(size);
235 initialCollection.getTopology()->getCellList(iold, &globalids[0]);
236 std::vector<int> ilocalnew(size); //local
237 std::vector<int> ipnew(size); //idomain old
238 _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
240 new2oldIds[iold].resize(nbNewDomain);
241 for (int i=0; i<(int)ilocalnew.size(); i++)
243 new2oldIds[iold][ipnew[i]].push_back(i);
245 for (int inew=0; inew<nbNewDomain; inew++)
247 splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)
248 (initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
249 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
251 if (MyGlobals::_Verbose>400)
252 std::cout<< "proc " << rank << " : a splitMesh iold inew NbCells " << iold << " " << inew << " "
253 << splitMeshes[inew][iold]->getNumberOfCells() << std::endl;
258 if (isParallelMode())
260 //if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
261 for (int iold=0; iold<nbOldDomain; iold++)
262 for(int inew=0; inew<nbNewDomain; inew++)
264 if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
267 if(initialCollection._domain_selector->isMyDomain(iold))
268 _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew));
270 if (_domain_selector->isMyDomain(inew))
271 _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
277 //fusing the split meshes
278 if (MyGlobals::_Verbose>200)
279 std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
280 for (int inew=0; inew<nbNewDomain ;inew++)
282 std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
284 for (int i=0; i<(int)splitMeshes[inew].size(); i++)
285 if (splitMeshes[inew][i]!=0)
286 if (splitMeshes[inew][i]->getNumberOfCells()>0)
287 meshes.push_back(splitMeshes[inew][i]);
289 if (!isParallelMode()||_domain_selector->isMyDomain(inew))
291 if (meshes.size()==0)
293 _mesh[inew]=CreateEmptyMEDCouplingUMesh();
294 std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
298 _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes);
299 o2nRenumber[inew]=_mesh[inew]->sortCellsInMEDFileFrmt();
304 ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
305 array->decrRef(); // array is not used in this case
307 _mesh[inew]->zipCoords();
310 for (int i=0;i<(int)splitMeshes[inew].size();i++)
311 if (splitMeshes[inew][i]!=0)
312 splitMeshes[inew][i]->decrRef();
314 if (MyGlobals::_Verbose>300)
315 std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
319 \param initialCollection source mesh collection
320 \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection
322 void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
325 using std::make_pair;
326 // NodeMapping reverseNodeMapping;
327 for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
331 BBTreeOfDim* tree = 0;
333 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
335 // std::map<pair<double,pair<double, double> >, int > nodeClassifier;
336 ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
337 double* coordsPtr=coords->getPointer();
338 dim = coords->getNumberOfComponents();
339 int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
340 bbox=new double[nvertices*2*dim];
342 for (int i=0; i<nvertices*dim;i++)
344 bbox[i*2]=coordsPtr[i]-1e-8;
345 bbox[i*2+1]=coordsPtr[i]+1e-8;
347 tree=new BBTreeOfDim( dim, bbox,0,0,nvertices,1e-9);
350 for (int inew=0; inew<_topology->nbDomain(); inew++)
353 //sending meshes for parallel computation
354 if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))
355 _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
356 else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
358 ParaMEDMEM::MEDCouplingUMesh* mesh;
359 _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
360 ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
361 for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
363 double* coordsPtr=coords->getPointer()+inode*dim;
365 tree->getElementsAroundPoint(coordsPtr,elems);
366 if (elems.size()==0) continue;
367 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
371 else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
373 if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
376 ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
377 for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
379 double* coordsPtr=coords->getPointer()+inode*dim;
381 tree->getElementsAroundPoint(coordsPtr,elems);
382 if (elems.size()==0) continue;
383 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
387 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
396 void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
399 using MEDPARTITIONER::BBTreeOfDim;
400 if (!&meshOne || !&meshTwo) return; //empty or not existing
402 BBTreeOfDim* tree = 0;
403 int nv1=meshOne.getNumberOfNodes();
404 ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
405 int dim = coords->getNumberOfComponents();
407 bbox=new double[nv1*2*dim];
408 double* coordsPtr=coords->getPointer();
409 for (int i=0; i<nv1*dim; i++)
411 bbox[i*2]=coordsPtr[i]-1e-8;
412 bbox[i*2+1]=coordsPtr[i]+1e-8;
414 tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9);
416 int nv2=meshTwo.getNumberOfNodes();
417 nodeIds.resize(nv2,-1);
418 coords=meshTwo.getCoords();
419 for (int inode=0; inode<nv2; inode++)
421 double* coordsPtr2=coords->getPointer()+inode*dim;
423 tree->getElementsAroundPoint(coordsPtr2,elems);
424 if (elems.size()==0) continue;
425 nodeIds[inode]=elems[0];
432 creates the face meshes on the new domains from the faces on the old domain and the node mapping
433 faces at the interface are duplicated
435 void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
436 const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
437 std::vector<std::vector<std::vector<int> > >& new2oldIds)
439 //splitMeshes structure will contain the partition of
440 //the old faces on the new ones
441 //splitMeshes[4][2] contains the faces from old domain 2
442 //that have to be added to domain 4
448 using std::make_pair;
450 if (MyGlobals::_Verbose>10)
451 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
453 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
455 int nbNewDomain=_topology->nbDomain();
456 int nbOldDomain=initialCollection.getTopology()->nbDomain();
458 vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
459 vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
461 vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
463 splitMeshes.resize(nbNewDomain);
464 for (int inew=0; inew<nbNewDomain; inew++)
466 splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
468 new2oldIds.resize(nbOldDomain);
469 for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
471 //init null pointer for empty meshes
472 for (int inew=0; inew<nbNewDomain; inew++)
474 for (int iold=0; iold<nbOldDomain; iold++)
476 splitMeshes[inew][iold]=0; //null for empty meshes
477 new2oldIds[iold][inew].clear();
481 //loop over the old domains to analyse the faces and decide
482 //on which new domain they belong
483 for (int iold=0; iold<nbOldDomain; iold++)
485 if (isParallelMode() && !initialCollection._domain_selector->isMyDomain(iold)) continue;
486 if (MyGlobals::_Verbose>400)
487 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes iodDomain "<<iold<<std::endl;
488 //initial face mesh known : in my domain
489 if (meshesCastFrom[iold] != 0)
491 for (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
494 meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
498 //analysis of element ielem
499 //counters are set for the element
500 //for each source node, the mapping is interrogated and the domain counters
501 //are incremented for each target node
502 //the face is considered as going to target domains if the counter of the domain
503 //is equal to the number of nodes
504 for (int inode=0; inode<(int)nodes.size(); inode++)
506 typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
507 int mynode=nodes[inode];
509 pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
510 for (MI iter=myRange.first; iter!=myRange.second; iter++)
512 int inew=iter->second.first;
513 if (faces.find(inew)==faces.end())
520 for (map<int,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
522 if (iter->second==(int)nodes.size())
523 //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
524 //it is not sure here...
525 //done before writeMedfile on option?... see filterFaceOnCell()
526 new2oldIds[iold][iter->first].push_back(ielem);
530 //creating the splitMeshes from the face ids
531 for (int inew=0; inew<nbNewDomain; inew++)
533 if (meshesCastFrom[iold]->getNumberOfCells() > 0)
535 splitMeshes[inew][iold]=
536 (ParaMEDMEM::MEDCouplingUMesh*)
537 ( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
538 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
541 splitMeshes[inew][iold]->zipCoords();
545 splitMeshes[inew][iold]=CreateEmptyMEDCouplingUMesh();
551 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
557 if (isParallelMode())
559 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
560 for (int iold=0; iold<nbOldDomain; iold++)
561 for (int inew=0; inew<nbNewDomain; inew++)
563 if (initialCollection._domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
565 if (splitMeshes[inew][iold] != 0)
567 _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew));
568 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes send " <<inew<<" "<<iold<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
572 _domain_selector->sendMesh(*(empty), _domain_selector->getProcessorID(inew));
573 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes sen0 " <<inew<<" "<<iold<<std::endl;
576 if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
577 _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
579 //if (splitMeshes[inew][iold])
580 // nb=splitMeshes[inew][iold]->getNumberOfCells();
581 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
587 //fusing the split meshes
588 if (MyGlobals::_Verbose>200)
589 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes fusing" << std::endl;
590 meshesCastTo.resize(nbNewDomain);
591 for (int inew=0; inew<nbNewDomain; inew++)
593 vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
594 for (int iold=0; iold<nbOldDomain; iold++)
596 ParaMEDMEM::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
598 if (umesh->getNumberOfCells()>0)
599 myMeshes.push_back(umesh);
602 ParaMEDMEM::MEDCouplingUMesh *bndMesh = 0;
603 if ( _subdomain_boundary_creates &&
605 _mesh[inew]->getNumberOfCells()>0 )
608 ((ParaMEDMEM::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
609 if (bndMesh->getNumberOfCells()>0)
610 myMeshes.push_back( bndMesh );
613 if (myMeshes.size()>0)
615 meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
616 meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
620 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
621 meshesCastTo[inew]=empty;
623 for (int iold=0; iold<nbOldDomain; iold++)
624 if (splitMeshes[inew][iold]!=0)
625 splitMeshes[inew][iold]->decrRef();
629 if (MyGlobals::_Verbose>300)
630 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
635 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
636 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
637 std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
638 std::string nameArrayTo)
642 int ioldMax=meshesCastFrom.size();
643 int inewMax=meshesCastTo.size();
646 //preparing bounding box trees for accelerating source-target node identifications
647 if (MyGlobals::_Verbose>99)
648 std::cout<<"making accelerating structures"<<std::endl;
649 std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
650 std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
651 for (int iold =0; iold< ioldMax; iold++)
652 if (isParallelMode() && _domain_selector->isMyDomain(iold))
654 ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
655 bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
656 acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
657 sourceCoords->decrRef();
659 // send-recv operations
661 for (int inew=0; inew<inewMax; inew++)
663 for (int iold=0; iold<ioldMax; iold++)
665 //sending arrays for distant domains
666 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
669 _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
671 int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
673 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
674 if (size>0) //no empty
676 sendIds.resize(size);
677 std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
682 sendIds.resize(size);
684 SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
686 //receiving arrays from distant domains
687 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
691 ParaMEDMEM::MEDCouplingUMesh* recvMesh;
692 _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
694 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
695 RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
696 remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
697 recvMesh->decrRef(); //cww is it correct?
703 //local contributions and aggregation
704 for (int inew=0; inew<inewMax; inew++)
706 for (int iold=0; iold<ioldMax; iold++)
707 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
709 remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
712 for (int iold=0; iold<ioldMax;iold++)
713 if (isParallelMode() && _domain_selector->isMyDomain(iold))
715 bbox[iold]->decrRef();
716 delete acceleratingStructures[iold];
720 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
721 const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
722 const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
723 const int* fromArray,
724 std::string nameArrayTo,
725 const BBTreeOfDim* myTree)
728 if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
729 ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
730 const double* tc=targetCoords->getConstPointer();
731 int targetSize=targetMesh.getNumberOfCells();
732 int sourceSize=sourceMesh.getNumberOfCells();
733 if (MyGlobals::_Verbose>200)
734 std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
735 std::vector<int> ccI;
737 str=nameArrayTo+"_toArray";
738 cle=Cle1ToStr(str,inew);
741 const BBTreeOfDim* tree;
742 bool cleantree=false;
743 ParaMEDMEM::DataArrayDouble* sourceBBox=0;
744 int dim = targetCoords->getNumberOfComponents();
747 sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
748 tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
752 //first time iold : create and initiate
753 if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
755 if (MyGlobals::_Is0verbose>100)
756 std::cout << "create " << cle << " size " << targetSize << std::endl;
757 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
758 p->alloc(targetSize,1);
760 toArray=p->getPointer();
761 _map_dataarray_int[cle]=p;
763 else //other times iold: refind and complete
765 toArray=_map_dataarray_int.find(cle)->second->getPointer();
768 std::map< int, int > isource2nb; // count coincident elements
769 std::map<int,int>::iterator i2nb;
771 for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
773 std::vector<int> intersectingElems;
774 tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
775 if (intersectingElems.size()!=0)
777 int isourcenode=intersectingElems[0];
778 if ( intersectingElems.size() > 1 )
780 i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
781 isourcenode = intersectingElems[ i2nb->second++ ];
783 if ( isourcenode < sourceSize ) // protection from invalid elements
785 toArray[itargetnode]=fromArray[isourcenode];
786 ccI.push_back(itargetnode);
787 ccI.push_back(isourcenode);
791 if (MyGlobals::_Verbose>200)
792 std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
793 //memories intersection for future same job on fields (if no existing cle=no intersection)
794 str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
795 if (MyGlobals::_Verbose>700)
796 std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
798 _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
800 targetCoords->decrRef();
801 if (cleantree) delete tree;
802 if (sourceBBox !=0) sourceBBox->decrRef();
805 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
807 if (nameArrayTo!="cellFieldDouble")
808 throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
810 std::string nameTo="typeData=6"; //resume the type of field casted
811 // send-recv operations
812 int ioldMax=initialCollection.getMesh().size();
813 int inewMax=this->getMesh().size();
814 int iFieldMax=initialCollection.getFieldDescriptions().size();
815 if (MyGlobals::_Verbose>10)
816 std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
817 //see collection.prepareFieldDescriptions()
818 for (int ifield=0; ifield<iFieldMax; ifield++)
820 std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
821 if (descriptionField.find(nameTo)==std::string::npos)
822 continue; //only nameTo accepted in Fields name description
824 for (int inew=0; inew<inewMax; inew++)
826 for (int iold=0; iold<ioldMax; iold++)
828 //sending arrays for distant domains
829 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
831 int target=_domain_selector->getProcessorID(inew);
832 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
833 if (MyGlobals::_Verbose>10)
834 std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
835 SendDataArrayDouble(field, target);
837 //receiving arrays from distant domains
838 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
840 int source=_domain_selector->getProcessorID(iold);
842 if (MyGlobals::_Verbose>10)
843 std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
844 ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
845 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
850 //local contributions and aggregation
851 for (int inew=0; inew<inewMax; inew++)
853 for (int iold=0; iold<ioldMax; iold++)
854 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
856 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
857 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
863 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
864 ParaMEDMEM::DataArrayDouble* fromArray,
865 std::string nameArrayTo,
866 std::string descriptionField)
867 //here we use 'cellFamily_ccI inew iold' set in remapIntField
869 if (nameArrayTo!="cellFieldDouble")
870 throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
871 std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
873 std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
874 it1=_map_dataarray_int.find(key);
875 if (it1==_map_dataarray_int.end())
877 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
878 std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
881 //create ccI in remapIntField
882 ParaMEDMEM::DataArrayInt *ccI=it1->second;
883 if (MyGlobals::_Verbose>300)
884 std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
886 int nbcell=this->getMesh()[inew]->getNumberOfCells();
887 int nbcomp=fromArray->getNumberOfComponents();
888 int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
890 std::string tag="inewFieldDouble="+IntToStr(inew);
891 key=descriptionField+SerializeFromString(tag);
892 int fromArrayNbOfElem=fromArray->getNbOfElems();
893 int fromArrayNbOfComp=fromArray->getNumberOfComponents();
894 int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
896 if (MyGlobals::_Verbose>1000)
898 std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
899 " fromArray nbOfElems " << fromArrayNbOfElem <<
900 " nbTuples " << fromArray->getNumberOfTuples() <<
901 " nbcells " << fromArrayNbOfCell <<
902 " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
905 ParaMEDMEM::DataArrayDouble* field=0;
906 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
907 it2=_map_dataarray_double.find(key);
908 if (it2==_map_dataarray_double.end())
910 if (MyGlobals::_Verbose>300)
911 std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
912 field=ParaMEDMEM::DataArrayDouble::New();
913 _map_dataarray_double[key]=field;
914 field->alloc(nbcell*nbPtGauss,nbcomp);
915 field->fillWithZero();
920 if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
922 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
923 " trying remap of field double on cells : \n" << descriptionField << std::endl;
930 field->setPartOfValuesAdv(fromArray,ccI);
934 //replaced by setPartOfValuesAdv if nbPtGauss==1
935 int iMax=ccI->getNbOfElems();
936 int* pccI=ccI->getPointer();
937 double* pField=field->getPointer();
938 double* pFrom=fromArray->getPointer();
939 int itarget, isource, delta=nbPtGauss*nbcomp;
940 for (int i=0; i<iMax; i=i+2) //cell
944 if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
945 throw INTERP_KERNEL::Exception("Error field override");
946 int ita=itarget*delta;
947 int iso=isource*delta;
948 for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
955 using namespace ParaMEDMEM;
956 //================================================================================
958 * \brief Sort correspondence ids of one domain and permute ids of the other accordingly
959 * \param [in,out] ids1 - ids of one domain
960 * \param [in,out] ids2 - ids of another domain
961 * \param [in] delta - a delta to change all ids
962 * \param [in] removeEqual - to remove equal ids
963 * \return DataArrayInt* - array of ids joined back
965 //================================================================================
967 DataArrayInt* sortCorrespondences( DataArrayInt* ids1,
970 bool removeEqual = false)
973 MEDCouplingAutoRefCountObjectPtr< DataArrayInt > renumN2O = ids1->buildPermArrPerLevel();
974 ids1->renumberInPlaceR( renumN2O->begin() );
975 ids2->renumberInPlaceR( renumN2O->begin() );
979 ids1 = ids1->buildUnique();
980 ids2 = ids2->buildUnique();
984 int * id = ids1->getPointer();
985 for ( ; id < ids1->end(); ++id )
987 id = ids2->getPointer();
988 for ( ; id < ids2->end(); ++id )
993 DataArrayInt* ids12 = DataArrayInt::Meld( ids1, ids2 ); // two components
994 ids12->rearrange( 1 ); // make one component
998 //================================================================================
1000 * \brief Renumber ids according to mesh->sortCellsInMEDFileFrmt()
1001 * \param [in,out] ids - cell ids to renumber
1002 * \param [in] o2nRenumber - renumbering array in "Old to New" mode
1004 //================================================================================
1006 void renumber( DataArrayInt* ids, const DataArrayInt* o2nRenumber )
1008 if ( !ids || !o2nRenumber )
1010 int * id = ids->getPointer();
1011 const int * o2n = o2nRenumber->getConstPointer();
1012 for ( ; id < ids->end(); ++id )
1019 //================================================================================
1021 * \brief Fill up ConnectZone's stored in _topology with nodal correspondences
1022 * \param [in] nodeMapping - mapping between old nodes and new nodes
1023 * (iolddomain,ioldnode)->(inewdomain,inewnode)
1024 * \param [in] o2nRenumber - renumbering array returned by mesh->sortCellsInMEDFileFrmt()
1026 * \param [in] nbInitialDomains - nb of old domains
1028 //================================================================================
1030 void MEDPARTITIONER::MeshCollection::buildConnectZones( const NodeMapping& nodeMapping,
1031 const std::vector<ParaMEDMEM::DataArrayInt*> & o2nRenumber,
1032 int nbInitialDomains)
1034 if ( !MyGlobals::_Create_Joints || _topology->nbDomain() < 2 )
1037 if ( MyGlobals::_World_Size > 1 )
1039 _topology->getCZ().clear();
1040 return; // not implemented for parallel mode
1043 // at construction, _topology creates cell correspondences basing on Graph information,
1045 // 1) add node correspondences,
1046 // 2) split cell correspondences by cell geometry types
1047 // 3) sort ids to be in ascending order
1049 const int nb_domains = _topology->nbDomain();
1051 // ==================================
1052 // 1) add node correspondences
1053 // ==================================
1055 std::vector< std::vector< std::vector< int > > > nodeCorresp( nb_domains );
1056 for ( int idomain = 0; idomain < nb_domains; ++idomain )
1058 nodeCorresp[ idomain ].resize( nb_domains );
1061 NodeMapping::const_iterator nmIt1, nmIt2 = nodeMapping.begin();
1062 for ( nmIt1 = nmIt2; nmIt1 != nodeMapping.end(); nmIt1 = nmIt2 )
1064 // look for an "old" node mapped into several "new" nodes in different domains
1066 while ( ++nmIt2 != nodeMapping.end() && nmIt2->first == nmIt1->first )
1067 nbSameOld += ( nmIt2->second != nmIt1->second );
1069 if ( nbSameOld > 0 )
1071 NodeMapping::const_iterator nmEnd = nmIt2;
1072 for ( ; true; ++nmIt1 )
1075 if ( ++nmIt2 == nmEnd )
1077 int dom1 = nmIt1->second.first;
1078 int node1 = nmIt1->second.second;
1079 for ( ; nmIt2 != nmEnd; ++nmIt2 )
1081 int dom2 = nmIt2->second.first;
1082 int node2 = nmIt2->second.second;
1085 nodeCorresp[ dom1 ][ dom2 ].push_back( node1 );
1086 nodeCorresp[ dom1 ][ dom2 ].push_back( node2 );
1087 nodeCorresp[ dom2 ][ dom1 ].push_back( node2 );
1088 nodeCorresp[ dom2 ][ dom1 ].push_back( node1 );
1095 // add nodeCorresp to czVec
1097 std::vector<MEDPARTITIONER::ConnectZone*>& czVec = _topology->getCZ();
1099 for ( int idomain = 0; idomain < nb_domains; ++idomain )
1101 for ( int idomainNear = 0; idomainNear < nb_domains; ++idomainNear )
1103 std::vector< int > & corresp = nodeCorresp[ idomain ][ idomainNear ];
1104 if ( corresp.empty() )
1107 MEDPARTITIONER::ConnectZone* cz = 0;
1108 for ( size_t i = 0; i < czVec.size() && !cz; ++i )
1110 czVec[i]->getLocalDomainNumber () == idomain &&
1111 czVec[i]->getDistantDomainNumber() == idomainNear )
1116 cz = new MEDPARTITIONER::ConnectZone();
1117 cz->setName( "Nodal Connect Zone defined by MEDPARTITIONER" );
1118 cz->setLocalDomainNumber ( idomain );
1119 cz->setDistantDomainNumber( idomainNear );
1120 czVec.push_back(cz);
1123 cz->setNodeCorresp( &corresp[0], corresp.size()/2 );
1127 // ==========================================================
1128 // 2) split cell correspondences by cell geometry types
1129 // ==========================================================
1131 for ( size_t i = 0; i < czVec.size(); ++i )
1133 MEDPARTITIONER::ConnectZone* cz = czVec[i];
1135 cz->getEntityCorrespNumber( 0,0 ) == 0 ||
1136 cz->getLocalDomainNumber () > (int)_mesh.size() ||
1137 cz->getDistantDomainNumber() > (int)_mesh.size() )
1139 ParaMEDMEM::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber () ];
1140 ParaMEDMEM::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ];
1142 // separate ids of two domains
1143 const ParaMEDMEM::MEDCouplingSkyLineArray *corrArray = cz->getEntityCorresp( 0, 0 );
1144 const DataArrayInt* ids12 = corrArray->getValueArray();
1145 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids1, ids2, ids12Sorted;
1146 ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 );
1147 ids2 = ids12->selectByTupleId2( 1, corrArray->getLength(), 2 );
1149 // renumber cells according to mesh->sortCellsInMEDFileFrmt()
1150 renumber( ids1, o2nRenumber[ cz->getLocalDomainNumber() ]);
1151 renumber( ids2, o2nRenumber[ cz->getDistantDomainNumber() ]);
1153 // check nb cell types
1154 std::set<INTERP_KERNEL::NormalizedCellType> types1, types2;
1155 types1 = mesh1->getTypesOfPart( ids1->begin(), ids1->end() );
1156 types2 = mesh2->getTypesOfPart( ids2->begin(), ids2->end() );
1157 if ( types1.size() < 1 || types2.size() < 1 )
1158 continue; // parallel mode?
1160 MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1161 for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1163 czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
1164 czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1167 if ( types1.size() == 1 && types2.size() == 1 ) // split not needed, only sort
1169 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1170 cz->setEntityCorresp( *types1.begin(), *types2.begin(),
1171 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1173 if ( cz21 )// set 2->1 correspondence
1175 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1176 cz21->setEntityCorresp( *types2.begin(), *types1.begin(),
1177 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1180 else // split and sort
1182 typedef std::pair< std::vector< int >, std::vector< int > > T2Vecs;
1183 T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ];
1186 const int nbIds = ids1->getNbOfElems();
1187 const int * p1 = ids1->begin(), * p2 = ids2->begin();
1188 for ( int i = 0; i < nbIds; ++i )
1190 t1 = mesh1->getTypeOfCell( p1[ i ]);
1191 t2 = mesh2->getTypeOfCell( p2[ i ]);
1192 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1193 ids.first .push_back( p1[ i ]);
1194 ids.second.push_back( p1[ i ]);
1197 const int maxType = int( INTERP_KERNEL::NORM_MAXTYPE );
1198 for ( t1 = 0; t1 < maxType; ++t1 )
1199 for ( t2 = 0; t2 < maxType; ++t2 )
1201 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1202 if ( ids.first.empty() ) continue;
1203 p1 = & ids.first[0];
1204 p2 = & ids.second[0];
1205 ids1->desallocate();
1206 ids1->pushBackValsSilent( p1, p1+ids.first.size() );
1207 ids2->desallocate();
1208 ids2->pushBackValsSilent( p2, p2+ids.first.size() );
1209 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1211 cz->setEntityCorresp( t1, t2,
1212 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1214 if ( cz21 )// set 2->1 correspondence
1216 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1217 cz21->setEntityCorresp( t2, t1,
1218 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1224 cz->setEntityCorresp( 0, 0, 0, 0 ); // erase ids computed by _topology
1226 cz21->setEntityCorresp( 0, 0, 0, 0 );
1231 // ==========================================
1232 // 3) sort node ids to be in ascending order
1233 // ==========================================
1235 const bool removeEqual = ( nbInitialDomains > 1 );
1237 for ( size_t i = 0; i < czVec.size(); ++i )
1239 MEDPARTITIONER::ConnectZone* cz = czVec[i];
1240 if ( !cz || cz->getNodeNumber() < 1 )
1242 if ( cz->getDistantDomainNumber() < cz->getLocalDomainNumber() )
1243 continue; // treat a pair of domains once
1245 MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1246 for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1248 czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() &&
1249 czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1252 // separate ids of two domains
1253 const ParaMEDMEM::MEDCouplingSkyLineArray *corrArray = cz->getNodeCorresp();
1254 const DataArrayInt *ids12 = corrArray->getValueArray();
1255 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids1, ids2, ids12Sorted;
1256 ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 );
1257 ids2 = ids12->selectByTupleId2( 1, corrArray->getLength(), 2 );
1259 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/0, removeEqual );
1260 cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1262 if ( cz21 )// set 2->1 correspondence
1264 ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0, false );
1265 cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1270 //================================================================================
1272 * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
1273 * group (where "n" and "p" are domain IDs)
1275 //================================================================================
1277 void MEDPARTITIONER::MeshCollection::buildBoundaryFaces()
1279 if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates )
1282 if ( getMeshDimension() < 2 )
1285 using ParaMEDMEM::MEDCouplingUMesh;
1286 using ParaMEDMEM::DataArrayDouble;
1287 using ParaMEDMEM::DataArrayInt;
1289 std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
1290 int nbMeshes = faceMeshes.size();
1292 //preparing bounding box trees for accelerating search of coincident faces
1293 std::vector<BBTreeOfDim* > bbTrees(nbMeshes);
1294 std::vector<DataArrayDouble*>bbox (nbMeshes);
1295 for (int inew = 0; inew < nbMeshes-1; inew++)
1296 if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
1298 DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner();
1299 bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
1300 bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
1301 bbox[inew]->getConstPointer(),0,0,
1302 bbox[inew]->getNumberOfTuples());
1303 bcCoords->decrRef();
1306 // loop on domains to find joint faces between them
1307 for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
1309 for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
1311 MEDCouplingUMesh* mesh1 = 0;
1312 MEDCouplingUMesh* mesh2 = 0;
1313 //MEDCouplingUMesh* recvMesh = 0;
1314 bool mesh1Here = true, mesh2Here = true;
1315 if (isParallelMode())
1318 mesh1Here = _domain_selector->isMyDomain(inew1);
1319 mesh2Here = _domain_selector->isMyDomain(inew2);
1320 if ( !mesh1Here && mesh2Here )
1322 //send mesh2 to domain of mesh1
1323 _domain_selector->sendMesh(*faceMeshes[inew2],
1324 _domain_selector->getProcessorID(inew1));
1326 else if ( mesh1Here && !mesh2Here )
1328 //receiving mesh2 from a distant domain
1329 _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
1330 if ( faceMeshes[ inew2 ] )
1331 faceMeshes[ inew2 ]->decrRef();
1332 faceMeshes[ inew2 ] = mesh2;
1336 if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
1337 if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
1339 // find coincident faces
1340 std::vector< int > faces1, faces2;
1341 if ( mesh1 && mesh2 )
1343 const DataArrayDouble* coords2 = mesh2->getBarycenterAndOwner();
1344 const double* c2 = coords2->getConstPointer();
1345 const int dim = coords2->getNumberOfComponents();
1346 const int nbFaces2 = mesh2->getNumberOfCells();
1347 const int nbFaces1 = mesh1->getNumberOfCells();
1349 for (int i2 = 0; i2 < nbFaces2; i2++)
1351 std::vector<int> coincFaces;
1352 bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
1353 if (coincFaces.size()!=0)
1355 int i1 = coincFaces[0];
1356 // if ( coincFaces.size() > 1 )
1358 // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
1359 // i1 = coincFaces[ i2nb->second++ ];
1361 if ( i1 < nbFaces1 ) // protection from invalid elements
1363 faces1.push_back( i1 );
1364 faces2.push_back( i2 );
1371 if ( isParallelMode())
1374 if ( mesh1Here && !mesh2Here )
1376 //send faces2 to domain of recvMesh
1377 SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
1379 else if ( !mesh1Here && mesh2Here )
1381 //receiving ids of faces from a domain of mesh1
1382 RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
1387 // recvMesh->decrRef();
1389 // Create group "JOINT_inew1_inew2_Faces" and corresponding families
1390 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1392 createJointGroup( is2nd ? faces2 : faces1,
1393 inew1 , inew2, is2nd );
1396 } // loop on the 2nd domains (inew2)
1397 } // loop on the 1st domains (inew1)
1400 // delete bounding box trees
1401 for (int inew = 0; inew < nbMeshes-1; inew++)
1402 if (isParallelMode() && _domain_selector->isMyDomain(inew))
1404 bbox[inew]->decrRef();
1405 delete bbTrees[inew];
1409 //================================================================================
1411 * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
1412 * \param faces - face ids to include into the group
1413 * \param inew1 - index of the 1st domain
1414 * \param inew2 - index of the 2nd domain
1415 * \param is2nd - in which (1st or 2nd) domain to create the group
1417 //================================================================================
1419 void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
1424 // get the name of JOINT group
1425 std::string groupName;
1427 std::ostringstream oss;
1429 << (is2nd ? inew2 : inew1 ) << "_"
1430 << (is2nd ? inew1 : inew2 ) << "_"
1431 << ( getMeshDimension()==2 ? "Edge" : "Face" );
1432 groupName = oss.str();
1435 // remove existing "JOINT_*" group
1436 _group_info.erase( groupName );
1438 // get family IDs array
1440 int inew = (is2nd ? inew2 : inew1 );
1441 int totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
1442 std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
1443 if ( !_map_dataarray_int.count(cle) )
1445 if ( totalNbFaces > 0 )
1447 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
1448 p->alloc( totalNbFaces, 1 );
1450 famIDs = p->getPointer();
1451 _map_dataarray_int[cle]=p;
1456 famIDs = _map_dataarray_int.find(cle)->second->getPointer();
1458 // find a family ID of an existing JOINT group
1460 std::map<std::string, int>::iterator name2id = _family_info.find( groupName );
1461 if ( name2id != _family_info.end() )
1462 familyID = name2id->second;
1464 // remove faces from the familyID-the family
1465 if ( familyID != 0 && famIDs )
1466 for ( int i = 0; i < totalNbFaces; ++i )
1467 if ( famIDs[i] == familyID )
1470 if ( faces.empty() )
1473 if ( familyID == 0 ) // generate a family ID for JOINT group
1475 std::set< int > familyIDs;
1476 for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
1477 familyIDs.insert( name2id->second );
1478 // find the next free family ID
1479 int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
1482 if ( !familyIDs.count( ++familyID ))
1485 while ( freeIdCount > 0 );
1488 // push faces to familyID-th group
1489 if ( faces.back() >= totalNbFaces )
1490 throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
1491 for ( size_t i = 0; i < faces.size(); ++i )
1492 famIDs[ faces[i] ] = familyID;
1494 // register JOINT group and family
1495 _family_info[ groupName ] = familyID; // name of the group and family is same
1496 _group_info [ groupName ].push_back( groupName );
1499 /*! constructing the MESH collection from a distributed file
1501 * \param filename name of the master file containing the list of all the MED files
1503 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
1505 _owns_topology(true),
1507 _domain_selector( 0 ),
1508 _i_non_empty_mesh(-1),
1509 _driver_type(MEDPARTITIONER::Undefined),
1510 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1511 _family_splitting(false),
1512 _create_empty_groups(false),
1517 _driver=new MeshCollectionMedXmlDriver(this);
1518 _driver->read (filename.c_str());
1519 _driver_type = MedXml;
1522 { // Handle all exceptions
1523 if ( _driver ) delete _driver; _driver=0;
1526 _driver=new MeshCollectionMedAsciiDriver(this);
1527 _driver->read (filename.c_str());
1528 _driver_type=MedAscii;
1534 throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
1537 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1538 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1539 _i_non_empty_mesh = idomain;
1542 /*! Constructing the MESH collection from selected domains of a distributed file
1544 * \param filename - name of the master file containing the list of all the MED files
1545 * \param domainSelector - selector of domains to load
1547 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
1549 _owns_topology(true),
1551 _domain_selector( &domainSelector ),
1552 _i_non_empty_mesh(-1),
1553 _driver_type(MEDPARTITIONER::Undefined),
1554 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1555 _family_splitting(false),
1556 _create_empty_groups(false),
1559 std::string myfile=filename;
1560 std::size_t found=myfile.find(".xml");
1561 if (found!=std::string::npos) //file .xml
1565 _driver=new MeshCollectionMedXmlDriver(this);
1566 _driver->read ( filename.c_str(), _domain_selector );
1567 _driver_type = MedXml;
1570 { // Handle all exceptions
1572 throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
1577 found=myfile.find(".med");
1578 if (found!=std::string::npos) //file .med single
1580 //make a temporary file .xml and retry MedXmlDriver
1582 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
1584 <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
1585 <description what=\"\" when=\"\"/>\n \
1587 <mesh name=\"$meshName\"/>\n \
1590 <subdomain number=\"1\"/>\n \
1591 <global_numbering present=\"no\"/>\n \
1594 <subfile id=\"1\">\n \
1595 <name>$fileName</name>\n \
1596 <machine>localhost</machine>\n \
1600 <mesh name=\"$meshName\">\n \
1601 <chunk subdomain=\"1\">\n \
1602 <name>$meshName</name>\n \
1607 std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile);
1608 xml.replace(xml.find("$fileName"),9,myfile);
1609 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1610 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1611 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1612 std::string nameFileXml(myfile);
1613 nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1614 std::string nameFileXmlDN,nameFileXmlBN;
1615 MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1616 nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1617 if (_domain_selector->rank()==0) //only on to write it
1619 std::ofstream f(nameFileXml.c_str());
1624 if (MyGlobals::_World_Size>1)
1625 MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1629 _driver=new MeshCollectionMedXmlDriver(this);
1630 _driver->read ( nameFileXml.c_str(), _domain_selector );
1631 _driver_type = MedXml;
1634 { // Handle all exceptions
1635 delete _driver; _driver=0;
1636 throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1643 _driver=new MeshCollectionMedAsciiDriver(this);
1644 _driver->read ( filename.c_str(), _domain_selector );
1645 _driver_type=MedAscii;
1651 throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1655 // find non-empty domain mesh
1656 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1657 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1658 _i_non_empty_mesh = idomain;
1662 //check for all proc/file compatibility of _field_descriptions
1664 _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1666 _field_descriptions=MyGlobals::_Field_Descriptions;
1669 catch(INTERP_KERNEL::Exception& e)
1671 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1672 throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1677 //check for all proc/file compatibility of _family_info
1678 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1679 _family_info=DevectorizeToMapOfStringInt(v2);
1681 catch(INTERP_KERNEL::Exception& e)
1683 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1684 throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1689 //check for all proc/file compatibility of _group_info
1690 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1691 _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1693 catch(INTERP_KERNEL::Exception& e)
1695 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1696 throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1701 /*! constructing the MESH collection from a sequential MED-file
1703 * \param filename MED file
1704 * \param meshname name of the mesh that is to be read
1706 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1708 _owns_topology(true),
1710 _domain_selector( 0 ),
1711 _i_non_empty_mesh(-1),
1713 _driver_type(MEDPARTITIONER::MedXml),
1714 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1715 _family_splitting(false),
1716 _create_empty_groups(false),
1719 try // avoid memory leak in case of inexistent filename
1721 retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1727 throw INTERP_KERNEL::Exception("problem reading .med files");
1729 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1730 _i_non_empty_mesh = 0;
1733 MEDPARTITIONER::MeshCollection::~MeshCollection()
1735 for (int i=0; i<(int)_mesh.size();i++)
1736 if (_mesh[i]!=0) _mesh[i]->decrRef();
1738 for (int i=0; i<(int)_cell_family_ids.size();i++)
1739 if (_cell_family_ids[i]!=0)
1740 _cell_family_ids[i]->decrRef();
1742 for (int i=0; i<(int)_face_mesh.size();i++)
1743 if (_face_mesh[i]!=0)
1744 _face_mesh[i]->decrRef();
1746 for (int i=0; i<(int)_face_family_ids.size();i++)
1747 if (_face_family_ids[i]!=0)
1748 _face_family_ids[i]->decrRef();
1750 for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1751 if ((*it).second!=0)
1752 (*it).second->decrRef();
1754 for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1755 if ((*it).second!=0)
1756 (*it).second->decrRef();
1759 if (_topology!=0 && _owns_topology)
1762 delete _joint_finder;
1766 /*! constructing the MESH collection from a file
1768 * The method creates as many MED-files as there are domains in the
1769 * collection. It also creates a master file that lists all the MED files.
1770 * The MED files created in ths manner contain joints that describe the
1771 * connectivity between subdomains.
1773 * \param filename name of the master file that will contain the list of the MED files
1776 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1778 //suppresses link with driver so that it can be changed for writing
1781 retrieveDriver()->write (filename.c_str(), _domain_selector);
1784 /*! creates or gets the link to the collection driver
1786 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1790 switch(_driver_type)
1793 _driver=new MeshCollectionMedXmlDriver(this);
1796 _driver=new MeshCollectionMedAsciiDriver(this);
1799 throw INTERP_KERNEL::Exception("Unrecognized driver");
1806 /*! gets an existing driver
1809 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1815 * retrieves the mesh dimension
1817 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1819 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1822 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1825 for (size_t i=0; i<_mesh.size(); i++)
1832 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1835 for (size_t i=0; i<_mesh.size(); i++)
1837 if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1842 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1845 for (size_t i=0; i<_face_mesh.size(); i++)
1847 if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1852 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1857 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1862 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1864 return _mesh[idomain];
1867 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1869 return _face_mesh[idomain];
1872 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1875 return _topology->getCZ();
1877 static std::vector<MEDPARTITIONER::ConnectZone*> noCZ;
1881 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1886 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo, bool takeOwneship)
1890 throw INTERP_KERNEL::Exception("topology is already set");
1895 _owns_topology = takeOwneship;
1899 /*! Method creating the cell graph in serial mode
1901 * \param array returns the pointer to the structure that contains the graph
1902 * \param edgeweight returns the pointer to the table that contains the edgeweights
1903 * (only used if indivisible regions are required)
1905 void MEDPARTITIONER::MeshCollection::buildCellGraph(ParaMEDMEM::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1910 using std::make_pair;
1913 if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
1914 const ParaMEDMEM::MEDCouplingUMesh* mesh=_mesh[0];
1915 if (MyGlobals::_Verbose>50)
1916 std::cout<<"getting nodal connectivity"<<std::endl;
1918 //looking for reverse nodal connectivity i global numbering
1919 if (isParallelMode() && !_domain_selector->isMyDomain(0))
1922 vector<int> index(1,0);
1924 array=new ParaMEDMEM::MEDCouplingSkyLineArray(index,value);
1927 array=mesh->generateGraph();
1929 /*! Method creating the cell graph in multidomain mode
1931 * \param array returns the pointer to the structure that contains the graph
1932 * \param edgeweight returns the pointer to the table that contains the edgeweights
1933 * (only used if indivisible regions are required)
1935 void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(ParaMEDMEM::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1937 using std::multimap;
1939 using std::make_pair;
1942 std::multimap< int, int > node2cell;
1943 std::map< pair<int,int>, int > cell2cellcounter;
1944 std::multimap<int,int> cell2cell;
1946 std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
1947 int nbdomain=_topology->nbDomain();
1949 if (isParallelMode())
1951 _joint_finder=new JointFinder(*this);
1952 _joint_finder->findCommonDistantNodes();
1953 commonDistantNodes=_joint_finder->getDistantNodeCell();
1956 if (MyGlobals::_Verbose>500)
1957 _joint_finder->print();
1960 if (MyGlobals::_Verbose>50)
1961 std::cout<<"getting nodal connectivity"<<std::endl;
1962 //looking for reverse nodal connectivity i global numbering
1964 for (int idomain=0; idomain<nbdomain; idomain++)
1966 if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1968 meshDim = _mesh[idomain]->getMeshDimension();
1970 ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
1971 ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1972 int nbNodes=_mesh[idomain]->getNumberOfNodes();
1973 _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1974 //problem saturation over 1 000 000 nodes for 1 proc
1975 if (MyGlobals::_Verbose>100)
1976 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1977 int* index_ptr=index->getPointer();
1978 int* revConnPtr=revConn->getPointer();
1979 for (int i=0; i<nbNodes; i++)
1981 for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
1983 int globalNode=_topology->convertNodeToGlobal(idomain,i);
1984 int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
1985 node2cell.insert(make_pair(globalNode, globalCell));
1991 for (int iother=0; iother<nbdomain; iother++)
1993 std::multimap<int,int>::iterator it;
1994 int isource=idomain;
1996 for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
1997 it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
1999 int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
2000 int globalCell=(*it).second;
2001 node2cell.insert(make_pair(globalNode, globalCell));
2007 //creating graph arcs (cell to cell relations)
2008 //arcs are stored in terms of (index,value) notation
2011 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
2012 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
2014 //warning here one node have less than or equal effective number of cell with it
2015 //but cell could have more than effective nodes
2016 //because other equals nodes in other domain (with other global inode)
2017 if (MyGlobals::_Verbose>50)
2018 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
2020 for (int inode=0;inode<_topology->nbNodes();inode++)
2022 typedef multimap<int,int>::const_iterator MI;
2023 std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
2024 for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
2025 for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
2027 int icell1=cell1->second;
2028 int icell2=cell2->second;
2029 if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
2030 std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2031 if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2032 else (it->second)++;
2035 // for (int icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
2037 // typedef multimap<int,int>::const_iterator MI;
2038 // std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
2039 // for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
2041 // std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
2042 // for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
2044 // int icell2=cell2->second;
2045 // std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2046 // if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2047 // else (it->second)++;
2053 //converting the counter to a multimap structure
2054 for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
2055 it!=cell2cellcounter.end();
2057 if (it->second>=meshDim)
2059 cell2cell.insert(std::make_pair(it->first.first,it->first.second));
2060 cell2cell.insert(std::make_pair(it->first.second, it->first.first));
2064 if (MyGlobals::_Verbose>50)
2065 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
2066 //filling up index and value to create skylinearray structure
2067 std::vector <int> index,value;
2071 for (int idomain=0; idomain<nbdomain; idomain++)
2073 if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
2074 int nbCells=_mesh[idomain]->getNumberOfCells();
2075 for (int icell=0; icell<nbCells; icell++)
2078 int globalCell=_topology->convertCellToGlobal(idomain,icell);
2079 multimap<int,int>::iterator it;
2080 pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
2081 ret=cell2cell.equal_range(globalCell);
2082 for (it=ret.first; it!=ret.second; ++it)
2084 int ival=(*it).second; //no adding one existing yet
2085 for (int i=idep ; i<idep+size ; i++)
2094 value.push_back(ival);
2098 idep=index[index.size()-1]+size;
2099 index.push_back(idep);
2103 array=new ParaMEDMEM::MEDCouplingSkyLineArray(index,value);
2105 if (MyGlobals::_Verbose>100)
2107 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
2108 index.size()-1 << " " << value.size() << std::endl;
2109 int max=index.size()>15?15:index.size();
2112 for (int i=0; i<max; ++i)
2113 std::cout<<index[i]<<" ";
2114 std::cout << "... " << index[index.size()-1] << std::endl;
2115 for (int i=0; i<max; ++i)
2116 std::cout<< value[i] << " ";
2117 int ll=index[index.size()-1]-1;
2118 std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
2125 /*! Creates the partition corresponding to the cell graph and the partition number
2127 * \param nbdomain number of subdomains for the newly created graph
2129 * returns a topology based on the new graph
2131 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
2132 Graph::splitter_type split,
2133 const std::string& options_string,
2134 int *user_edge_weights,
2135 int *user_vertices_weights)
2137 if (MyGlobals::_Verbose>10)
2138 std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
2141 throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
2142 ParaMEDMEM::MEDCouplingSkyLineArray* array=0;
2145 if (_topology->nbDomain()>1 || isParallelMode())
2146 buildParallelCellGraph(array,edgeweights);
2148 buildCellGraph(array,edgeweights);
2150 Graph* cellGraph = 0;
2154 if ( isParallelMode() && MyGlobals::_World_Size > 1 )
2156 #ifdef MED_ENABLE_PARMETIS
2157 if (MyGlobals::_Verbose>10)
2158 std::cout << "ParMETISGraph" << std::endl;
2159 cellGraph=new ParMETISGraph(array,edgeweights);
2164 #ifdef MED_ENABLE_METIS
2165 if (MyGlobals::_Verbose>10)
2166 std::cout << "METISGraph" << std::endl;
2167 cellGraph=new METISGraph(array,edgeweights);
2171 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
2175 #ifdef MED_ENABLE_SCOTCH
2176 if (MyGlobals::_Verbose>10)
2177 std::cout << "SCOTCHGraph" << std::endl;
2178 cellGraph=new SCOTCHGraph(array,edgeweights);
2180 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
2185 //!user-defined weights
2186 if (user_edge_weights!=0)
2187 cellGraph->setEdgesWeights(user_edge_weights);
2188 if (user_vertices_weights!=0)
2189 cellGraph->setVerticesWeights(user_vertices_weights);
2191 if (MyGlobals::_Is0verbose>10)
2192 std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
2193 cellGraph->partGraph(nbdomain, options_string, _domain_selector);
2195 if (MyGlobals::_Is0verbose>10)
2196 std::cout << "building new topology" << std::endl;
2197 //cellGraph is a shared pointer
2198 Topology *topology=0;
2199 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2201 delete [] edgeweights;
2203 if (MyGlobals::_Verbose>11)
2204 std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
2208 /*! Creates a topology for a partition specified by the user
2210 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
2212 * returns a topology based on the new partition
2214 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
2216 ParaMEDMEM::MEDCouplingSkyLineArray* array=0;
2219 if ( _topology->nbDomain()>1)
2220 buildParallelCellGraph(array,edgeweights);
2222 buildCellGraph(array,edgeweights);
2225 std::set<int> domains;
2226 for (int i=0; i<_topology->nbCells(); i++)
2228 domains.insert(partition[i]);
2230 cellGraph=new UserGraph(array, partition, _topology->nbCells());
2232 //cellGraph is a shared pointer
2233 Topology *topology=0;
2234 int nbdomain=domains.size();
2235 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2236 // if (array!=0) delete array;
2241 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
2243 for (int i=0; i<_topology->nbDomain(); i++)
2245 std::ostringstream oss;
2247 if (!isParallelMode() || _domain_selector->isMyDomain(i))
2248 _mesh[i]->setName(oss.str());
2252 ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
2253 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
2254 //something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
2256 int rank=MyGlobals::_Rank;
2257 std::string tag="ioldFieldDouble="+IntToStr(iold);
2258 std::string descriptionIold=descriptionField+SerializeFromString(tag);
2259 if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
2261 if (MyGlobals::_Verbose>300)
2262 std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
2263 ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
2266 if (MyGlobals::_Verbose>200)
2267 std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
2268 std::string description, fileName, meshName, fieldName;
2269 int typeField, DT, IT, entity;
2270 fileName=MyGlobals::_File_Names[iold];
2271 if (MyGlobals::_Verbose>10)
2272 std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
2273 FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
2274 meshName=MyGlobals::_Mesh_Names[iold];
2276 ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
2277 fileName, meshName, 0, fieldName, DT, IT);
2279 ParaMEDMEM::DataArrayDouble* res=f2->getArray();
2280 //to know names of components
2281 std::vector<std::string> browse=BrowseFieldDouble(f2);
2282 std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
2283 if (MyGlobals::_Verbose>10)
2284 std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
2285 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 ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
2367 ParaMEDMEM::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 ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
2376 ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::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] = (ParaMEDMEM::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 ParaMEDMEM::DataArrayInt * & fam = getMapDataArrayInt()[ key ];
2447 ParaMEDMEM::DataArrayInt * famFilter = ParaMEDMEM::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]];