1 // Copyright (C) 2007-2014 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"
21 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
22 #include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx"
23 #include "MEDPARTITIONER_MeshCollectionMedAsciiDriver.hxx"
24 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
25 #include "MEDPARTITIONER_Topology.hxx"
26 #include "MEDPARTITIONER_ParallelTopology.hxx"
29 #include "MEDPARTITIONER_JointFinder.hxx"
32 #include "MEDPARTITIONER_Graph.hxx"
33 #include "MEDPARTITIONER_UserGraph.hxx"
34 #include "MEDPARTITIONER_Utils.hxx"
36 #include "MEDLoaderBase.hxx"
37 #include "MEDLoader.hxx"
38 #include "MEDCouplingMemArray.hxx"
39 #include "MEDCouplingUMesh.hxx"
40 #include "MEDCouplingNormalizedUnstructuredMesh.hxx"
41 #include "MEDCouplingFieldDouble.hxx"
42 #include "PointLocator3DIntersectorP0P0.hxx"
44 #include "MEDCouplingAutoRefCountObjectPtr.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 castCellMeshes(initialCollection, new2oldIds);
111 //defining the name for the collection and the underlying meshes
112 setName(initialCollection.getName());
119 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
120 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
122 if (MyGlobals::_Is0verbose)
123 std::cout<<"treating faces"<<std::endl;
124 NodeMapping nodeMapping;
125 //nodeMapping contains the mapping between old nodes and new nodes
126 // (iolddomain,ioldnode)->(inewdomain,inewnode)
127 createNodeMapping(initialCollection, nodeMapping);
128 std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
129 castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
135 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
136 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
138 if (MyGlobals::_Is0verbose)
140 if (isParallelMode())
141 std::cout << "ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
143 std::cout << "NOT ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
145 if (MyGlobals::_Is0verbose>10)
146 std::cout<<"treating cell and face families"<<std::endl;
148 castIntField(initialCollection.getMesh(),
150 initialCollection.getCellFamilyIds(),
152 castIntField(initialCollection.getFaceMesh(),
154 initialCollection.getFaceFamilyIds(),
159 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
160 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
162 if (MyGlobals::_Is0verbose)
163 std::cout << "treating groups" << std::endl;
164 _family_info=initialCollection.getFamilyInfo();
165 _group_info=initialCollection.getGroupInfo();
168 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
169 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
171 if (MyGlobals::_Is0verbose)
172 std::cout << "treating fields" << std::endl;
173 castAllFields(initialCollection,"cellFieldDouble");
174 if (_i_non_empty_mesh<0)
176 for (int i=0; i<_mesh.size(); i++)
180 _i_non_empty_mesh=i; //first existing one local
189 Creates the meshes using the topology underlying he mesh collection and the mesh data
190 coming from the ancient collection
191 \param initialCollection collection from which the data is extracted to create the new meshes
194 void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
195 std::vector<std::vector<std::vector<int> > >& new2oldIds)
197 if (MyGlobals::_Verbose>10)
198 std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
200 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes");
202 int nbNewDomain=_topology->nbDomain();
203 int nbOldDomain=initialCollection.getTopology()->nbDomain();
205 _mesh.resize(nbNewDomain);
206 int rank=MyGlobals::_Rank;
207 //splitting the initial domains into smaller bits
208 std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
209 splitMeshes.resize(nbNewDomain);
210 for (int inew=0; inew<nbNewDomain; inew++)
212 splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
215 for (int iold=0; iold<nbOldDomain; iold++)
217 if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
219 int size=(initialCollection._mesh)[iold]->getNumberOfCells();
220 std::vector<int> globalids(size);
221 initialCollection.getTopology()->getCellList(iold, &globalids[0]);
222 std::vector<int> ilocalnew(size); //local
223 std::vector<int> ipnew(size); //idomain old
224 _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
226 new2oldIds[iold].resize(nbNewDomain);
227 for (int i=0; i<(int)ilocalnew.size(); i++)
229 new2oldIds[iold][ipnew[i]].push_back(i);
231 for (int inew=0; inew<nbNewDomain; inew++)
233 splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)
234 (initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
235 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
237 if (MyGlobals::_Verbose>400)
238 std::cout<< "proc " << rank << " : a splitMesh iold inew NbCells " << iold << " " << inew << " "
239 << splitMeshes[inew][iold]->getNumberOfCells() << std::endl;
244 if (isParallelMode())
246 //if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
247 for (int iold=0; iold<nbOldDomain; iold++)
248 for(int inew=0; inew<nbNewDomain; inew++)
250 if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
253 if(initialCollection._domain_selector->isMyDomain(iold))
254 _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew));
256 if (_domain_selector->isMyDomain(inew))
257 _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
263 //fusing the split meshes
264 if (MyGlobals::_Verbose>200)
265 std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
266 for (int inew=0; inew<nbNewDomain ;inew++)
268 std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
270 for (int i=0; i<(int)splitMeshes[inew].size(); i++)
271 if (splitMeshes[inew][i]!=0)
272 if (splitMeshes[inew][i]->getNumberOfCells()>0)
273 meshes.push_back(splitMeshes[inew][i]);
275 if (!isParallelMode()||_domain_selector->isMyDomain(inew))
277 if (meshes.size()==0)
279 _mesh[inew]=CreateEmptyMEDCouplingUMesh();
280 std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
284 _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes);
289 ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
290 array->decrRef(); // array is not used in this case
292 _mesh[inew]->zipCoords();
296 for (int i=0;i<(int)splitMeshes[inew].size();i++)
297 if (splitMeshes[inew][i]!=0)
298 splitMeshes[inew][i]->decrRef();
300 if (MyGlobals::_Verbose>300)
301 std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
305 \param initialCollection source mesh collection
306 \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection
308 void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
311 using std::make_pair;
312 // NodeMapping reverseNodeMapping;
313 for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
317 BBTreeOfDim* tree = 0;
319 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
321 // std::map<pair<double,pair<double, double> >, int > nodeClassifier;
322 ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
323 double* coordsPtr=coords->getPointer();
324 dim = coords->getNumberOfComponents();
325 int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
326 bbox=new double[nvertices*2*dim];
328 for (int i=0; i<nvertices*dim;i++)
330 bbox[i*2]=coordsPtr[i]-1e-8;
331 bbox[i*2+1]=coordsPtr[i]+1e-8;
333 tree=new BBTreeOfDim( dim, bbox,0,0,nvertices,1e-9);
336 for (int inew=0; inew<_topology->nbDomain(); inew++)
339 //sending meshes for parallel computation
340 if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))
341 _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
342 else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
344 ParaMEDMEM::MEDCouplingUMesh* mesh;
345 _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
346 ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
347 for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
349 double* coordsPtr=coords->getPointer()+inode*dim;
351 tree->getElementsAroundPoint(coordsPtr,elems);
352 if (elems.size()==0) continue;
353 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
357 else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
359 if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
362 ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
363 for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
365 double* coordsPtr=coords->getPointer()+inode*dim;
367 tree->getElementsAroundPoint(coordsPtr,elems);
368 if (elems.size()==0) continue;
369 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
373 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
382 void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
385 using MEDPARTITIONER::BBTreeOfDim;
386 if (!&meshOne || !&meshTwo) return; //empty or not existing
388 BBTreeOfDim* tree = 0;
389 int nv1=meshOne.getNumberOfNodes();
390 ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
391 int dim = coords->getNumberOfComponents();
393 bbox=new double[nv1*2*dim];
394 double* coordsPtr=coords->getPointer();
395 for (int i=0; i<nv1*dim; i++)
397 bbox[i*2]=coordsPtr[i]-1e-8;
398 bbox[i*2+1]=coordsPtr[i]+1e-8;
400 tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9);
402 int nv2=meshTwo.getNumberOfNodes();
403 nodeIds.resize(nv2,-1);
404 coords=meshTwo.getCoords();
405 for (int inode=0; inode<nv2; inode++)
407 double* coordsPtr2=coords->getPointer()+inode*dim;
409 tree->getElementsAroundPoint(coordsPtr2,elems);
410 if (elems.size()==0) continue;
411 nodeIds[inode]=elems[0];
418 creates the face meshes on the new domains from the faces on the old domain and the node mapping
419 faces at the interface are duplicated
421 void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
422 const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
423 std::vector<std::vector<std::vector<int> > >& new2oldIds)
425 //splitMeshes structure will contain the partition of
426 //the old faces on the new ones
427 //splitMeshes[4][2] contains the faces from old domain 2
428 //that have to be added to domain 4
434 using std::make_pair;
436 if (MyGlobals::_Verbose>10)
437 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
439 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
441 int nbNewDomain=_topology->nbDomain();
442 int nbOldDomain=initialCollection.getTopology()->nbDomain();
444 vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
445 vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
447 vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
449 splitMeshes.resize(nbNewDomain);
450 for (int inew=0; inew<nbNewDomain; inew++)
452 splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
454 new2oldIds.resize(nbOldDomain);
455 for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
457 //init null pointer for empty meshes
458 for (int inew=0; inew<nbNewDomain; inew++)
460 for (int iold=0; iold<nbOldDomain; iold++)
462 splitMeshes[inew][iold]=0; //null for empty meshes
463 new2oldIds[iold][inew].clear();
467 //loop over the old domains to analyse the faces and decide
468 //on which new domain they belong
469 for (int iold=0; iold<nbOldDomain; iold++)
471 if (isParallelMode() && !initialCollection._domain_selector->isMyDomain(iold)) continue;
472 if (MyGlobals::_Verbose>400)
473 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes iodDomain "<<iold<<std::endl;
474 //initial face mesh known : in my domain
475 if (meshesCastFrom[iold] != 0)
477 for (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
480 meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
484 //analysis of element ielem
485 //counters are set for the element
486 //for each source node, the mapping is interrogated and the domain counters
487 //are incremented for each target node
488 //the face is considered as going to target domains if the counter of the domain
489 //is equal to the number of nodes
490 for (int inode=0; inode<(int)nodes.size(); inode++)
492 typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
493 int mynode=nodes[inode];
495 pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
496 for (MI iter=myRange.first; iter!=myRange.second; iter++)
498 int inew=iter->second.first;
499 if (faces.find(inew)==faces.end())
506 for (map<int,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
508 if (iter->second==(int)nodes.size())
509 //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
510 //it is not sure here...
511 //done before writeMedfile on option?... see filterFaceOnCell()
512 new2oldIds[iold][iter->first].push_back(ielem);
516 //creating the splitMeshes from the face ids
517 for (int inew=0; inew<nbNewDomain; inew++)
519 if (meshesCastFrom[iold]->getNumberOfCells() > 0)
521 splitMeshes[inew][iold]=
522 (ParaMEDMEM::MEDCouplingUMesh*)
523 ( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
524 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
527 splitMeshes[inew][iold]->zipCoords();
531 splitMeshes[inew][iold]=CreateEmptyMEDCouplingUMesh();
537 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
543 if (isParallelMode())
545 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
546 for (int iold=0; iold<nbOldDomain; iold++)
547 for (int inew=0; inew<nbNewDomain; inew++)
549 if (initialCollection._domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
551 if (splitMeshes[inew][iold] != 0)
553 _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew));
554 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes send " <<inew<<" "<<iold<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
558 _domain_selector->sendMesh(*(empty), _domain_selector->getProcessorID(inew));
559 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes sen0 " <<inew<<" "<<iold<<std::endl;
562 if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
563 _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
565 //if (splitMeshes[inew][iold])
566 // nb=splitMeshes[inew][iold]->getNumberOfCells();
567 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
573 //fusing the split meshes
574 if (MyGlobals::_Verbose>200)
575 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes fusing" << std::endl;
576 meshesCastTo.resize(nbNewDomain);
577 for (int inew=0; inew<nbNewDomain; inew++)
579 vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
580 for (int iold=0; iold<nbOldDomain; iold++)
582 ParaMEDMEM::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
584 if (umesh->getNumberOfCells()>0)
585 myMeshes.push_back(umesh);
588 ParaMEDMEM::MEDCouplingUMesh *bndMesh = 0;
589 if ( _subdomain_boundary_creates &&
591 _mesh[inew]->getNumberOfCells()>0 )
594 ((ParaMEDMEM::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
595 if (bndMesh->getNumberOfCells()>0)
596 myMeshes.push_back( bndMesh );
599 if (myMeshes.size()>0)
601 meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
602 meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
606 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
607 meshesCastTo[inew]=empty;
609 for (int iold=0; iold<nbOldDomain; iold++)
610 if (splitMeshes[inew][iold]!=0)
611 splitMeshes[inew][iold]->decrRef();
615 if (MyGlobals::_Verbose>300)
616 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
621 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
622 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
623 std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
624 std::string nameArrayTo)
628 int ioldMax=meshesCastFrom.size();
629 int inewMax=meshesCastTo.size();
632 //preparing bounding box trees for accelerating source-target node identifications
633 if (MyGlobals::_Verbose>99)
634 std::cout<<"making accelerating structures"<<std::endl;
635 std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
636 std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
637 for (int iold =0; iold< ioldMax; iold++)
638 if (isParallelMode() && _domain_selector->isMyDomain(iold))
640 ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
641 bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
642 acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
643 sourceCoords->decrRef();
645 // send-recv operations
647 for (int inew=0; inew<inewMax; inew++)
649 for (int iold=0; iold<ioldMax; iold++)
651 //sending arrays for distant domains
652 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
655 _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
657 int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
659 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
660 if (size>0) //no empty
662 sendIds.resize(size);
663 std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
668 sendIds.resize(size);
670 SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
672 //receiving arrays from distant domains
673 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
677 ParaMEDMEM::MEDCouplingUMesh* recvMesh;
678 _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
680 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
681 RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
682 remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
683 recvMesh->decrRef(); //cww is it correct?
689 //local contributions and aggregation
690 for (int inew=0; inew<inewMax; inew++)
692 for (int iold=0; iold<ioldMax; iold++)
693 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
695 remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
698 for (int iold=0; iold<ioldMax;iold++)
699 if (isParallelMode() && _domain_selector->isMyDomain(iold))
701 bbox[iold]->decrRef();
702 delete acceleratingStructures[iold];
706 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
707 const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
708 const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
709 const int* fromArray,
710 std::string nameArrayTo,
711 const BBTreeOfDim* myTree)
714 if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
715 ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
716 const double* tc=targetCoords->getConstPointer();
717 int targetSize=targetMesh.getNumberOfCells();
718 int sourceSize=sourceMesh.getNumberOfCells();
719 if (MyGlobals::_Verbose>200)
720 std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
721 std::vector<int> ccI;
723 str=nameArrayTo+"_toArray";
724 cle=Cle1ToStr(str,inew);
727 const BBTreeOfDim* tree;
728 bool cleantree=false;
729 ParaMEDMEM::DataArrayDouble* sourceBBox=0;
730 int dim = targetCoords->getNumberOfComponents();
733 sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
734 tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
738 //first time iold : create and initiate
739 if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
741 if (MyGlobals::_Is0verbose>100)
742 std::cout << "create " << cle << " size " << targetSize << std::endl;
743 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
744 p->alloc(targetSize,1);
746 toArray=p->getPointer();
747 _map_dataarray_int[cle]=p;
749 else //other times iold: refind and complete
751 toArray=_map_dataarray_int.find(cle)->second->getPointer();
754 std::map< int, int > isource2nb; // count coincident elements
755 std::map<int,int>::iterator i2nb;
757 for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
759 std::vector<int> intersectingElems;
760 tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
761 if (intersectingElems.size()!=0)
763 int isourcenode=intersectingElems[0];
764 if ( intersectingElems.size() > 1 )
766 i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
767 isourcenode = intersectingElems[ i2nb->second++ ];
769 if ( isourcenode < sourceSize ) // protection from invalid elements
771 toArray[itargetnode]=fromArray[isourcenode];
772 ccI.push_back(itargetnode);
773 ccI.push_back(isourcenode);
777 if (MyGlobals::_Verbose>200)
778 std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
779 //memories intersection for future same job on fields (if no existing cle=no intersection)
780 str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
781 if (MyGlobals::_Verbose>700)
782 std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
784 _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
786 targetCoords->decrRef();
787 if (cleantree) delete tree;
788 if (sourceBBox !=0) sourceBBox->decrRef();
791 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
793 if (nameArrayTo!="cellFieldDouble")
794 throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
796 std::string nameTo="typeData=6"; //resume the type of field casted
797 // send-recv operations
798 int ioldMax=initialCollection.getMesh().size();
799 int inewMax=this->getMesh().size();
800 int iFieldMax=initialCollection.getFieldDescriptions().size();
801 if (MyGlobals::_Verbose>10)
802 std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
803 //see collection.prepareFieldDescriptions()
804 for (int ifield=0; ifield<iFieldMax; ifield++)
806 std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
807 if (descriptionField.find(nameTo)==std::string::npos)
808 continue; //only nameTo accepted in Fields name description
810 for (int inew=0; inew<inewMax; inew++)
812 for (int iold=0; iold<ioldMax; iold++)
814 //sending arrays for distant domains
815 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
817 int target=_domain_selector->getProcessorID(inew);
818 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
819 if (MyGlobals::_Verbose>10)
820 std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
821 SendDataArrayDouble(field, target);
823 //receiving arrays from distant domains
824 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
826 int source=_domain_selector->getProcessorID(iold);
828 if (MyGlobals::_Verbose>10)
829 std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
830 ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
831 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
836 //local contributions and aggregation
837 for (int inew=0; inew<inewMax; inew++)
839 for (int iold=0; iold<ioldMax; iold++)
840 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
842 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
843 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
849 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
850 ParaMEDMEM::DataArrayDouble* fromArray,
851 std::string nameArrayTo,
852 std::string descriptionField)
853 //here we use 'cellFamily_ccI inew iold' set in remapIntField
855 if (nameArrayTo!="cellFieldDouble")
856 throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
857 std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
859 std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
860 it1=_map_dataarray_int.find(key);
861 if (it1==_map_dataarray_int.end())
863 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
864 std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
867 //create ccI in remapIntField
868 ParaMEDMEM::DataArrayInt *ccI=it1->second;
869 if (MyGlobals::_Verbose>300)
870 std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
872 int nbcell=this->getMesh()[inew]->getNumberOfCells();
873 int nbcomp=fromArray->getNumberOfComponents();
874 int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
876 std::string tag="inewFieldDouble="+IntToStr(inew);
877 key=descriptionField+SerializeFromString(tag);
878 int fromArrayNbOfElem=fromArray->getNbOfElems();
879 int fromArrayNbOfComp=fromArray->getNumberOfComponents();
880 int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
882 if (MyGlobals::_Verbose>1000)
884 std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
885 " fromArray nbOfElems " << fromArrayNbOfElem <<
886 " nbTuples " << fromArray->getNumberOfTuples() <<
887 " nbcells " << fromArrayNbOfCell <<
888 " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
891 ParaMEDMEM::DataArrayDouble* field=0;
892 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
893 it2=_map_dataarray_double.find(key);
894 if (it2==_map_dataarray_double.end())
896 if (MyGlobals::_Verbose>300)
897 std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
898 field=ParaMEDMEM::DataArrayDouble::New();
899 _map_dataarray_double[key]=field;
900 field->alloc(nbcell*nbPtGauss,nbcomp);
901 field->fillWithZero();
906 if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
908 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
909 " trying remap of field double on cells : \n" << descriptionField << std::endl;
916 field->setPartOfValuesAdv(fromArray,ccI);
920 //replaced by setPartOfValuesAdv if nbPtGauss==1
921 int iMax=ccI->getNbOfElems();
922 int* pccI=ccI->getPointer();
923 double* pField=field->getPointer();
924 double* pFrom=fromArray->getPointer();
925 int itarget, isource, delta=nbPtGauss*nbcomp;
926 for (int i=0; i<iMax; i=i+2) //cell
930 if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
931 throw INTERP_KERNEL::Exception("Error field override");
932 int ita=itarget*delta;
933 int iso=isource*delta;
934 for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
939 //================================================================================
941 * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
942 * group (where "n" and "p" are domain IDs)
944 //================================================================================
946 void MEDPARTITIONER::MeshCollection::buildConnectZones()
948 if ( getMeshDimension() < 2 )
951 using ParaMEDMEM::MEDCouplingUMesh;
952 using ParaMEDMEM::DataArrayDouble;
953 using ParaMEDMEM::DataArrayInt;
955 std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
956 int nbMeshes = faceMeshes.size();
958 //preparing bounding box trees for accelerating search of coincident faces
959 std::vector<BBTreeOfDim* > bbTrees(nbMeshes);
960 std::vector<DataArrayDouble*>bbox (nbMeshes);
961 for (int inew = 0; inew < nbMeshes-1; inew++)
962 if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
964 DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner();
965 bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
966 bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
967 bbox[inew]->getConstPointer(),0,0,
968 bbox[inew]->getNumberOfTuples());
972 // loop on domains to find joint faces between them
973 for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
975 for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
977 MEDCouplingUMesh* mesh1 = 0;
978 MEDCouplingUMesh* mesh2 = 0;
979 //MEDCouplingUMesh* recvMesh = 0;
980 bool mesh1Here = true, mesh2Here = true;
981 if (isParallelMode())
984 mesh1Here = _domain_selector->isMyDomain(inew1);
985 mesh2Here = _domain_selector->isMyDomain(inew2);
986 if ( !mesh1Here && mesh2Here )
988 //send mesh2 to domain of mesh1
989 _domain_selector->sendMesh(*faceMeshes[inew2],
990 _domain_selector->getProcessorID(inew1));
992 else if ( mesh1Here && !mesh2Here )
994 //receiving mesh2 from a distant domain
995 _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
996 if ( faceMeshes[ inew2 ] )
997 faceMeshes[ inew2 ]->decrRef();
998 faceMeshes[ inew2 ] = mesh2;
1002 if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
1003 if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
1005 // find coincident faces
1006 std::vector< int > faces1, faces2;
1007 if ( mesh1 && mesh2 )
1009 const DataArrayDouble* coords2 = mesh2->getBarycenterAndOwner();
1010 const double* c2 = coords2->getConstPointer();
1011 const int dim = coords2->getNumberOfComponents();
1012 const int nbFaces2 = mesh2->getNumberOfCells();
1013 const int nbFaces1 = mesh1->getNumberOfCells();
1015 for (int i2 = 0; i2 < nbFaces2; i2++)
1017 std::vector<int> coincFaces;
1018 bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
1019 if (coincFaces.size()!=0)
1021 int i1 = coincFaces[0];
1022 // if ( coincFaces.size() > 1 )
1024 // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
1025 // i1 = coincFaces[ i2nb->second++ ];
1027 if ( i1 < nbFaces1 ) // protection from invalid elements
1029 faces1.push_back( i1 );
1030 faces2.push_back( i2 );
1037 if ( isParallelMode())
1040 if ( mesh1Here && !mesh2Here )
1042 //send faces2 to domain of recvMesh
1043 SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
1045 else if ( !mesh1Here && mesh2Here )
1047 //receiving ids of faces from a domain of mesh1
1048 RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
1053 // recvMesh->decrRef();
1055 // Create group "JOINT_inew1_inew2_Faces" and corresponding families
1056 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1058 createJointGroup( is2nd ? faces2 : faces1,
1059 inew1 , inew2, is2nd );
1062 } // loop on the 2nd domains (inew2)
1063 } // loop on the 1st domains (inew1)
1066 // delete bounding box trees
1067 for (int inew = 0; inew < nbMeshes-1; inew++)
1068 if (isParallelMode() && _domain_selector->isMyDomain(inew))
1070 bbox[inew]->decrRef();
1071 delete bbTrees[inew];
1075 //================================================================================
1077 * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
1078 * \param faces - face ids to include into the group
1079 * \param inew1 - index of the 1st domain
1080 * \param inew2 - index of the 2nd domain
1081 * \param is2nd - in which (1st or 2nd) domain to create the group
1083 //================================================================================
1085 void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
1090 // get the name of JOINT group
1091 std::string groupName;
1093 std::ostringstream oss;
1095 << (is2nd ? inew2 : inew1 ) << "_"
1096 << (is2nd ? inew1 : inew2 ) << "_"
1097 << ( getMeshDimension()==2 ? "Edge" : "Face" );
1098 groupName = oss.str();
1101 // remove existing "JOINT_*" group
1102 _group_info.erase( groupName );
1104 // get family IDs array
1106 int inew = (is2nd ? inew2 : inew1 );
1107 int totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
1108 std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
1109 if ( !_map_dataarray_int.count(cle) )
1111 if ( totalNbFaces > 0 )
1113 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
1114 p->alloc( totalNbFaces, 1 );
1116 famIDs = p->getPointer();
1117 _map_dataarray_int[cle]=p;
1122 famIDs = _map_dataarray_int.find(cle)->second->getPointer();
1124 // find a family ID of an existing JOINT group
1126 std::map<std::string, int>::iterator name2id = _family_info.find( groupName );
1127 if ( name2id != _family_info.end() )
1128 familyID = name2id->second;
1130 // remove faces from the familyID-the family
1131 if ( familyID != 0 && famIDs )
1132 for ( size_t i = 0; i < totalNbFaces; ++i )
1133 if ( famIDs[i] == familyID )
1136 if ( faces.empty() )
1139 if ( familyID == 0 ) // generate a family ID for JOINT group
1141 std::set< int > familyIDs;
1142 for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
1143 familyIDs.insert( name2id->second );
1144 // find the next free family ID
1145 int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
1148 if ( !familyIDs.count( ++familyID ))
1151 while ( freeIdCount > 0 );
1154 // push faces to familyID-th group
1155 if ( faces.back() >= totalNbFaces )
1156 throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
1157 for ( size_t i = 0; i < faces.size(); ++i )
1158 famIDs[ faces[i] ] = familyID;
1160 // register JOINT group and family
1161 _family_info[ groupName ] = familyID; // name of the group and family is same
1162 _group_info [ groupName ].push_back( groupName );
1165 /*! constructing the MESH collection from a distributed file
1167 * \param filename name of the master file containing the list of all the MED files
1169 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
1171 _owns_topology(true),
1173 _domain_selector( 0 ),
1174 _i_non_empty_mesh(-1),
1175 _driver_type(MEDPARTITIONER::Undefined),
1176 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1177 _family_splitting(false),
1178 _create_empty_groups(false),
1183 _driver=new MeshCollectionMedXmlDriver(this);
1184 _driver->read (filename.c_str());
1185 _driver_type = MedXml;
1188 { // Handle all exceptions
1189 if ( _driver ) delete _driver; _driver=0;
1192 _driver=new MeshCollectionMedAsciiDriver(this);
1193 _driver->read (filename.c_str());
1194 _driver_type=MedAscii;
1200 throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
1203 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1204 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1205 _i_non_empty_mesh = idomain;
1208 /*! Constructing the MESH collection from selected domains of a distributed file
1210 * \param filename - name of the master file containing the list of all the MED files
1211 * \param domainSelector - selector of domains to load
1213 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
1215 _owns_topology(true),
1217 _domain_selector( &domainSelector ),
1218 _i_non_empty_mesh(-1),
1219 _driver_type(MEDPARTITIONER::Undefined),
1220 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1221 _family_splitting(false),
1222 _create_empty_groups(false),
1225 std::string myfile=filename;
1226 std::size_t found=myfile.find(".xml");
1227 if (found!=std::string::npos) //file .xml
1231 _driver=new MeshCollectionMedXmlDriver(this);
1232 _driver->read ( filename.c_str(), _domain_selector );
1233 _driver_type = MedXml;
1236 { // Handle all exceptions
1238 throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
1243 found=myfile.find(".med");
1244 if (found!=std::string::npos) //file .med single
1246 //make a temporary file .xml and retry MedXmlDriver
1248 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
1250 <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
1251 <description what=\"\" when=\"\"/>\n \
1253 <mesh name=\"$meshName\"/>\n \
1256 <subdomain number=\"1\"/>\n \
1257 <global_numbering present=\"no\"/>\n \
1260 <subfile id=\"1\">\n \
1261 <name>$fileName</name>\n \
1262 <machine>localhost</machine>\n \
1266 <mesh name=\"$meshName\">\n \
1267 <chunk subdomain=\"1\">\n \
1268 <name>$meshName</name>\n \
1273 std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile);
1274 xml.replace(xml.find("$fileName"),9,myfile);
1275 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1276 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1277 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1278 std::string nameFileXml(myfile);
1279 nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1280 std::string nameFileXmlDN,nameFileXmlBN;
1281 MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1282 nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1283 if (_domain_selector->rank()==0) //only on to write it
1285 std::ofstream f(nameFileXml.c_str());
1290 if (MyGlobals::_World_Size>1)
1291 MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1295 _driver=new MeshCollectionMedXmlDriver(this);
1296 _driver->read ( nameFileXml.c_str(), _domain_selector );
1297 _driver_type = MedXml;
1300 { // Handle all exceptions
1301 delete _driver; _driver=0;
1302 throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1309 _driver=new MeshCollectionMedAsciiDriver(this);
1310 _driver->read ( filename.c_str(), _domain_selector );
1311 _driver_type=MedAscii;
1317 throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1321 // find non-empty domain mesh
1322 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1323 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1324 _i_non_empty_mesh = idomain;
1328 //check for all proc/file compatibility of _field_descriptions
1330 _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1332 _field_descriptions=MyGlobals::_Field_Descriptions;
1335 catch(INTERP_KERNEL::Exception& e)
1337 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1338 throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1343 //check for all proc/file compatibility of _family_info
1344 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1345 _family_info=DevectorizeToMapOfStringInt(v2);
1347 catch(INTERP_KERNEL::Exception& e)
1349 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1350 throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1355 //check for all proc/file compatibility of _group_info
1356 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1357 _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1359 catch(INTERP_KERNEL::Exception& e)
1361 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1362 throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1367 /*! constructing the MESH collection from a sequential MED-file
1369 * \param filename MED file
1370 * \param meshname name of the mesh that is to be read
1372 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1374 _owns_topology(true),
1376 _domain_selector( 0 ),
1377 _i_non_empty_mesh(-1),
1379 _driver_type(MEDPARTITIONER::MedXml),
1380 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1381 _family_splitting(false),
1382 _create_empty_groups(false),
1385 try // avoid memory leak in case of inexistent filename
1387 retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1393 throw INTERP_KERNEL::Exception("problem reading .med files");
1395 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1396 _i_non_empty_mesh = 0;
1399 MEDPARTITIONER::MeshCollection::~MeshCollection()
1401 for (int i=0; i<(int)_mesh.size();i++)
1402 if (_mesh[i]!=0) _mesh[i]->decrRef();
1404 for (int i=0; i<(int)_cell_family_ids.size();i++)
1405 if (_cell_family_ids[i]!=0)
1406 _cell_family_ids[i]->decrRef();
1408 for (int i=0; i<(int)_face_mesh.size();i++)
1409 if (_face_mesh[i]!=0)
1410 _face_mesh[i]->decrRef();
1412 for (int i=0; i<(int)_face_family_ids.size();i++)
1413 if (_face_family_ids[i]!=0)
1414 _face_family_ids[i]->decrRef();
1416 for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1417 if ((*it).second!=0)
1418 (*it).second->decrRef();
1420 for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1421 if ((*it).second!=0)
1422 (*it).second->decrRef();
1425 if (_topology!=0 && _owns_topology)
1428 delete _joint_finder;
1432 /*! constructing the MESH collection from a file
1434 * The method creates as many MED-files as there are domains in the
1435 * collection. It also creates a master file that lists all the MED files.
1436 * The MED files created in ths manner contain joints that describe the
1437 * connectivity between subdomains.
1439 * \param filename name of the master file that will contain the list of the MED files
1442 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1444 //building the connect zones necessary for writing joints
1445 if (_topology->nbDomain()>1 && _subdomain_boundary_creates )
1446 buildConnectZones();
1447 //suppresses link with driver so that it can be changed for writing
1450 retrieveDriver()->write (filename.c_str(), _domain_selector);
1453 /*! creates or gets the link to the collection driver
1455 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1459 switch(_driver_type)
1462 _driver=new MeshCollectionMedXmlDriver(this);
1465 _driver=new MeshCollectionMedAsciiDriver(this);
1468 throw INTERP_KERNEL::Exception("Unrecognized driver");
1475 /*! gets an existing driver
1478 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1484 * retrieves the mesh dimension
1486 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1488 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1491 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1494 for (int i=0; i<_mesh.size(); i++)
1501 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1504 for (int i=0; i<_mesh.size(); i++)
1506 if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1511 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1514 for (int i=0; i<_face_mesh.size(); i++)
1516 if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1521 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1526 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1531 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1533 return _mesh[idomain];
1536 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1538 return _face_mesh[idomain];
1541 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1543 return _connect_zones;
1546 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1551 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo)
1555 throw INTERP_KERNEL::Exception("topology is already set");
1561 /*! Method creating the cell graph in serial mode
1563 * \param array returns the pointer to the structure that contains the graph
1564 * \param edgeweight returns the pointer to the table that contains the edgeweights
1565 * (only used if indivisible regions are required)
1567 void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
1572 using std::make_pair;
1575 if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
1576 const ParaMEDMEM::MEDCouplingUMesh* mesh=_mesh[0];
1577 if (MyGlobals::_Verbose>50)
1578 std::cout<<"getting nodal connectivity"<<std::endl;
1580 //looking for reverse nodal connectivity i global numbering
1581 if (isParallelMode() && !_domain_selector->isMyDomain(0))
1584 vector<int> index(1,0);
1586 array=new MEDPARTITIONER::SkyLineArray(index,value);
1590 int meshDim = mesh->getMeshDimension();
1592 ParaMEDMEM::DataArrayInt* indexr=ParaMEDMEM::DataArrayInt::New();
1593 ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1594 int nbNodes=mesh->getNumberOfNodes();
1595 mesh->getReverseNodalConnectivity(revConn,indexr);
1596 //problem saturation over 1 000 000 nodes for 1 proc
1597 if (MyGlobals::_Verbose>100)
1598 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1599 const int* indexr_ptr=indexr->getConstPointer();
1600 const int* revConn_ptr=revConn->getConstPointer();
1602 const ParaMEDMEM::DataArrayInt* index;
1603 const ParaMEDMEM::DataArrayInt* conn;
1604 conn=mesh->getNodalConnectivity();
1605 index=mesh->getNodalConnectivityIndex();
1606 int nbCells=mesh->getNumberOfCells();
1608 if (MyGlobals::_Verbose>100)
1609 std::cout << "proc " << MyGlobals::_Rank << " : getNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1610 const int* index_ptr=index->getConstPointer();
1611 const int* conn_ptr=conn->getConstPointer();
1613 //creating graph arcs (cell to cell relations)
1614 //arcs are stored in terms of (index,value) notation
1617 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
1618 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
1620 //warning here one node have less than or equal effective number of cell with it
1621 //but cell could have more than effective nodes
1622 //because other equals nodes in other domain (with other global inode)
1623 if (MyGlobals::_Verbose>50)
1624 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
1626 vector <int> cell2cell_index(nbCells+1,0);
1627 vector <int> cell2cell;
1628 cell2cell.reserve(3*nbCells);
1630 for (int icell=0; icell<nbCells;icell++)
1632 map<int,int > counter;
1633 for (int iconn=index_ptr[icell]; iconn<index_ptr[icell+1];iconn++)
1635 int inode=conn_ptr[iconn];
1636 for (int iconnr=indexr_ptr[inode]; iconnr<indexr_ptr[inode+1];iconnr++)
1638 int icell2=revConn_ptr[iconnr];
1639 map<int,int>::iterator iter=counter.find(icell2);
1640 if (iter!=counter.end()) (iter->second)++;
1641 else counter.insert(make_pair(icell2,1));
1644 for (map<int,int>::const_iterator iter=counter.begin();
1645 iter!=counter.end();
1647 if (iter->second >= meshDim)
1649 cell2cell_index[icell+1]++;
1650 cell2cell.push_back(iter->first);
1658 cell2cell_index[0]=0;
1659 for (int icell=0; icell<nbCells;icell++)
1660 cell2cell_index[icell+1]=cell2cell_index[icell]+cell2cell_index[icell+1];
1663 if (MyGlobals::_Verbose>50)
1664 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
1666 //filling up index and value to create skylinearray structure
1667 array=new MEDPARTITIONER::SkyLineArray(cell2cell_index,cell2cell);
1669 if (MyGlobals::_Verbose>100)
1671 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
1672 cell2cell_index.size()-1 << " " << cell2cell.size() << std::endl;
1673 int max=cell2cell_index.size()>15?15:cell2cell_index.size();
1674 if (cell2cell_index.size()>1)
1676 for (int i=0; i<max; ++i)
1677 std::cout<<cell2cell_index[i]<<" ";
1678 std::cout << "... " << cell2cell_index[cell2cell_index.size()-1] << std::endl;
1679 for (int i=0; i<max; ++i)
1680 std::cout<< cell2cell[i] << " ";
1681 int ll=cell2cell_index[cell2cell_index.size()-1]-1;
1682 std::cout << "... (" << ll << ") " << cell2cell[ll-1] << " " << cell2cell[ll] << std::endl;
1687 /*! Method creating the cell graph in multidomain mode
1689 * \param array returns the pointer to the structure that contains the graph
1690 * \param edgeweight returns the pointer to the table that contains the edgeweights
1691 * (only used if indivisible regions are required)
1693 void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
1695 using std::multimap;
1697 using std::make_pair;
1700 std::multimap< int, int > node2cell;
1701 std::map< pair<int,int>, int > cell2cellcounter;
1702 std::multimap<int,int> cell2cell;
1704 std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
1705 int nbdomain=_topology->nbDomain();
1707 if (isParallelMode())
1709 _joint_finder=new JointFinder(*this);
1710 _joint_finder->findCommonDistantNodes();
1711 commonDistantNodes=_joint_finder->getDistantNodeCell();
1714 if (MyGlobals::_Verbose>500)
1715 _joint_finder->print();
1718 if (MyGlobals::_Verbose>50)
1719 std::cout<<"getting nodal connectivity"<<std::endl;
1720 //looking for reverse nodal connectivity i global numbering
1722 for (int idomain=0; idomain<nbdomain; idomain++)
1724 if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1726 meshDim = _mesh[idomain]->getMeshDimension();
1728 ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
1729 ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1730 int nbNodes=_mesh[idomain]->getNumberOfNodes();
1731 _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1732 //problem saturation over 1 000 000 nodes for 1 proc
1733 if (MyGlobals::_Verbose>100)
1734 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1735 int* index_ptr=index->getPointer();
1736 int* revConnPtr=revConn->getPointer();
1737 for (int i=0; i<nbNodes; i++)
1739 for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
1741 int globalNode=_topology->convertNodeToGlobal(idomain,i);
1742 int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
1743 node2cell.insert(make_pair(globalNode, globalCell));
1749 for (int iother=0; iother<nbdomain; iother++)
1751 std::multimap<int,int>::iterator it;
1752 int isource=idomain;
1754 for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
1755 it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
1757 int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
1758 int globalCell=(*it).second;
1759 node2cell.insert(make_pair(globalNode, globalCell));
1765 //creating graph arcs (cell to cell relations)
1766 //arcs are stored in terms of (index,value) notation
1769 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
1770 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
1772 //warning here one node have less than or equal effective number of cell with it
1773 //but cell could have more than effective nodes
1774 //because other equals nodes in other domain (with other global inode)
1775 if (MyGlobals::_Verbose>50)
1776 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
1778 for (int inode=0;inode<_topology->nbNodes();inode++)
1780 typedef multimap<int,int>::const_iterator MI;
1781 std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
1782 for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
1783 for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
1785 int icell1=cell1->second;
1786 int icell2=cell2->second;
1787 if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
1788 std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
1789 if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
1790 else (it->second)++;
1793 // for (int icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
1795 // typedef multimap<int,int>::const_iterator MI;
1796 // std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
1797 // for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
1799 // std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
1800 // for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
1802 // int icell2=cell2->second;
1803 // std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
1804 // if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
1805 // else (it->second)++;
1811 //converting the counter to a multimap structure
1812 for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
1813 it!=cell2cellcounter.end();
1815 if (it->second>=meshDim)
1817 cell2cell.insert(std::make_pair(it->first.first,it->first.second));
1818 cell2cell.insert(std::make_pair(it->first.second, it->first.first));
1822 if (MyGlobals::_Verbose>50)
1823 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
1824 //filling up index and value to create skylinearray structure
1825 std::vector <int> index,value;
1829 for (int idomain=0; idomain<nbdomain; idomain++)
1831 if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
1832 int nbCells=_mesh[idomain]->getNumberOfCells();
1833 for (int icell=0; icell<nbCells; icell++)
1836 int globalCell=_topology->convertCellToGlobal(idomain,icell);
1837 multimap<int,int>::iterator it;
1838 pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
1839 ret=cell2cell.equal_range(globalCell);
1840 for (it=ret.first; it!=ret.second; ++it)
1842 int ival=(*it).second; //no adding one existing yet
1843 for (int i=idep ; i<idep+size ; i++)
1852 value.push_back(ival);
1856 idep=index[index.size()-1]+size;
1857 index.push_back(idep);
1861 array=new MEDPARTITIONER::SkyLineArray(index,value);
1863 if (MyGlobals::_Verbose>100)
1865 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
1866 index.size()-1 << " " << value.size() << std::endl;
1867 int max=index.size()>15?15:index.size();
1870 for (int i=0; i<max; ++i)
1871 std::cout<<index[i]<<" ";
1872 std::cout << "... " << index[index.size()-1] << std::endl;
1873 for (int i=0; i<max; ++i)
1874 std::cout<< value[i] << " ";
1875 int ll=index[index.size()-1]-1;
1876 std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
1883 /*! Creates the partition corresponding to the cell graph and the partition number
1885 * \param nbdomain number of subdomains for the newly created graph
1887 * returns a topology based on the new graph
1889 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
1890 Graph::splitter_type split,
1891 const std::string& options_string,
1892 int *user_edge_weights,
1893 int *user_vertices_weights)
1895 if (MyGlobals::_Verbose>10)
1896 std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
1899 throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
1900 MEDPARTITIONER::SkyLineArray* array=0;
1903 if (_topology->nbDomain()>1 || isParallelMode())
1904 buildParallelCellGraph(array,edgeweights);
1906 buildCellGraph(array,edgeweights);
1908 Graph* cellGraph = 0;
1912 if ( isParallelMode() && MyGlobals::_World_Size > 1 )
1914 #ifdef MED_ENABLE_PARMETIS
1915 if (MyGlobals::_Verbose>10)
1916 std::cout << "ParMETISGraph" << std::endl;
1917 cellGraph=new ParMETISGraph(array,edgeweights);
1922 #ifdef MED_ENABLE_METIS
1923 if (MyGlobals::_Verbose>10)
1924 std::cout << "METISGraph" << std::endl;
1925 cellGraph=new METISGraph(array,edgeweights);
1929 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
1933 #ifdef MED_ENABLE_SCOTCH
1934 if (MyGlobals::_Verbose>10)
1935 std::cout << "SCOTCHGraph" << std::endl;
1936 cellGraph=new SCOTCHGraph(array,edgeweights);
1938 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
1943 //!user-defined weights
1944 if (user_edge_weights!=0)
1945 cellGraph->setEdgesWeights(user_edge_weights);
1946 if (user_vertices_weights!=0)
1947 cellGraph->setVerticesWeights(user_vertices_weights);
1949 if (MyGlobals::_Is0verbose>10)
1950 std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
1951 cellGraph->partGraph(nbdomain, options_string, _domain_selector);
1953 if (MyGlobals::_Is0verbose>10)
1954 std::cout << "building new topology" << std::endl;
1955 //cellGraph is a shared pointer
1956 Topology *topology=0;
1957 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
1959 delete [] edgeweights;
1961 if (MyGlobals::_Verbose>11)
1962 std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
1966 /*! Creates a topology for a partition specified by the user
1968 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
1970 * returns a topology based on the new partition
1972 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
1974 MEDPARTITIONER::SkyLineArray* array=0;
1977 if ( _topology->nbDomain()>1)
1978 buildParallelCellGraph(array,edgeweights);
1980 buildCellGraph(array,edgeweights);
1983 std::set<int> domains;
1984 for (int i=0; i<_topology->nbCells(); i++)
1986 domains.insert(partition[i]);
1988 cellGraph=new UserGraph(array, partition, _topology->nbCells());
1990 //cellGraph is a shared pointer
1991 Topology *topology=0;
1992 int nbdomain=domains.size();
1993 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
1994 // if (array!=0) delete array;
1999 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
2001 for (int i=0; i<_topology->nbDomain(); i++)
2003 std::ostringstream oss;
2005 if (!isParallelMode() || _domain_selector->isMyDomain(i))
2006 _mesh[i]->setName(oss.str());
2010 ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
2011 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
2012 //something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
2014 int rank=MyGlobals::_Rank;
2015 std::string tag="ioldFieldDouble="+IntToStr(iold);
2016 std::string descriptionIold=descriptionField+SerializeFromString(tag);
2017 if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
2019 if (MyGlobals::_Verbose>300)
2020 std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
2021 ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
2024 if (MyGlobals::_Verbose>200)
2025 std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
2026 std::string description, fileName, meshName, fieldName;
2027 int typeField, DT, IT, entity;
2028 fileName=MyGlobals::_File_Names[iold];
2029 if (MyGlobals::_Verbose>10)
2030 std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
2031 FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
2032 meshName=MyGlobals::_Mesh_Names[iold];
2034 ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
2035 fileName, meshName, 0, fieldName, DT, IT);
2037 ParaMEDMEM::DataArrayDouble* res=f2->getArray();
2038 //to know names of components
2039 std::vector<std::string> browse=BrowseFieldDouble(f2);
2040 std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
2041 if (MyGlobals::_Verbose>10)
2042 std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
2043 MyGlobals::_General_Informations.push_back(localFieldInformation);
2046 _map_dataarray_double[descriptionIold]=res;
2050 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
2051 //to have unique valid fields names/pointers/descriptions for partitionning
2052 //filter _field_descriptions to be in all procs compliant and equal
2054 int nbfiles=MyGlobals::_File_Names.size(); //nb domains
2055 std::vector<std::string> r2;
2056 //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
2057 for (int i=0; i<(int)_field_descriptions.size(); i++)
2059 std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
2060 for (int ii=0; ii<(int)r1.size(); ii++)
2061 r2.push_back(r1[ii]);
2063 //here vector(procs*fields) of serialised vector(description) data
2064 _field_descriptions=r2;
2065 int nbfields=_field_descriptions.size(); //on all domains
2066 if ((nbfields%nbfiles)!=0)
2068 if (MyGlobals::_Rank==0)
2070 std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
2071 << "fileMedNames :" << std::endl
2072 << ReprVectorOfString(MyGlobals::_File_Names)
2073 << "field_descriptions :" << std::endl
2074 << ReprVectorOfString(MyGlobals::_Field_Descriptions);
2076 throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
2078 _field_descriptions.resize(nbfields/nbfiles);
2079 for (int i=0; i<(int)_field_descriptions.size(); i++)
2081 std::string str=_field_descriptions[i];
2082 str=EraseTagSerialized(str,"idomain=");
2083 str=EraseTagSerialized(str,"fileName=");
2084 _field_descriptions[i]=str;
2088 //returns true if inodes of a face are in inodes of a cell
2089 bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >& inodesCell)
2092 int nbok=inodesFace.size();
2093 for (int i=0; i<nbok; i++)
2095 int ii=inodesFace[i];
2097 std::cout << "isFaceOncell problem inodeface<0" << std::endl;
2098 for (int j=0; j<(int)inodesCell.size(); j++)
2100 if (ii==inodesCell[j])
2103 break; //inode of face found
2107 break; //inode of face not found do not continue...
2109 return (ires==nbok);
2112 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
2114 for (int inew=0; inew<_topology->nbDomain(); inew++)
2116 if (!isParallelMode() || _domain_selector->isMyDomain(inew))
2118 if (MyGlobals::_Verbose>200)
2119 std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
2120 ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
2121 ParaMEDMEM::MEDCouplingUMesh* mfac=_face_mesh[inew];
2123 //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
2124 std::vector<int> nodeIds;
2125 getNodeIds(*mcel, *mfac, nodeIds);
2126 if (nodeIds.size()==0)
2127 continue; //one empty mesh nothing to do
2129 ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
2130 ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New();
2131 mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
2132 int *revC=revNodalCel->getPointer();
2133 int *revIndxC=revNodalIndxCel->getPointer();
2135 std::vector< int > faceOnCell;
2136 std::vector< int > faceNotOnCell;
2137 int nbface=mfac->getNumberOfCells();
2138 for (int iface=0; iface<nbface; iface++)
2141 std::vector< int > inodesFace;
2142 mfac->getNodeIdsOfCell(iface, inodesFace);
2143 int nbnodFace=inodesFace.size();
2144 if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
2145 continue; // invalid node ids
2146 //set inodesFace in mcel
2148 for (int i=0; i<nbnodFace; i++)
2149 nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
2150 if ( nbok != nbnodFace )
2152 int inod=inodesFace[0];
2155 std::cout << "filterFaceOnCell problem 1" << std::endl;
2158 int nbcell=revIndxC[inod+1]-revIndxC[inod];
2159 for (int j=0; j<nbcell; j++) //look for each cell with inod
2161 int icel=revC[revIndxC[inod]+j];
2162 std::vector< int > inodesCell;
2163 mcel->getNodeIdsOfCell(icel, inodesCell);
2164 ok=isFaceOncell(inodesFace, inodesCell);
2169 faceOnCell.push_back(iface);
2173 faceNotOnCell.push_back(iface);
2174 if (MyGlobals::_Is0verbose>300)
2175 std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
2179 revNodalCel->decrRef();
2180 revNodalIndxCel->decrRef();
2182 // std::string keyy;
2183 // keyy=Cle1ToStr("filterFaceOnCell",inew);
2184 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
2185 // keyy=Cle1ToStr("filterNotFaceOnCell",inew);
2186 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
2188 // filter the face mesh
2189 if ( faceOnCell.empty() )
2190 _face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
2192 _face_mesh[inew] = (ParaMEDMEM::MEDCouplingUMesh *)
2193 mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
2196 // filter the face families
2197 std::string key = Cle1ToStr("faceFamily_toArray",inew);
2198 if ( getMapDataArrayInt().count( key ))
2200 ParaMEDMEM::DataArrayInt * & fam = getMapDataArrayInt()[ key ];
2201 ParaMEDMEM::DataArrayInt * famFilter = ParaMEDMEM::DataArrayInt::New();
2202 famFilter->alloc(faceOnCell.size(),1);
2203 int* pfamFilter = famFilter->getPointer();
2204 int* pfam = fam->getPointer();
2205 for ( size_t i=0; i<faceOnCell.size(); i++ )
2206 pfamFilter[i]=pfam[faceOnCell[i]];