1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
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);
605 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
606 meshesCastTo[inew]=empty;
608 for (int iold=0; iold<nbOldDomain; iold++)
609 if (splitMeshes[inew][iold]!=0)
610 splitMeshes[inew][iold]->decrRef();
614 if (MyGlobals::_Verbose>300)
615 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
620 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
621 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
622 std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
623 std::string nameArrayTo)
627 int ioldMax=meshesCastFrom.size();
628 int inewMax=meshesCastTo.size();
631 //preparing bounding box trees for accelerating source-target node identifications
632 if (MyGlobals::_Verbose>99)
633 std::cout<<"making accelerating structures"<<std::endl;
634 std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
635 std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
636 for (int iold =0; iold< ioldMax; iold++)
637 if (isParallelMode() && _domain_selector->isMyDomain(iold))
639 ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
640 bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
641 acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
642 sourceCoords->decrRef();
644 // send-recv operations
646 for (int inew=0; inew<inewMax; inew++)
648 for (int iold=0; iold<ioldMax; iold++)
650 //sending arrays for distant domains
651 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
654 _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
656 int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
658 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
659 if (size>0) //no empty
661 sendIds.resize(size);
662 std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
667 sendIds.resize(size);
669 SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
671 //receiving arrays from distant domains
672 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
676 ParaMEDMEM::MEDCouplingUMesh* recvMesh;
677 _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
679 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
680 RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
681 remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
682 recvMesh->decrRef(); //cww is it correct?
688 //local contributions and aggregation
689 for (int inew=0; inew<inewMax; inew++)
691 for (int iold=0; iold<ioldMax; iold++)
692 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
694 remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
697 for (int iold=0; iold<ioldMax;iold++)
698 if (isParallelMode() && _domain_selector->isMyDomain(iold))
700 bbox[iold]->decrRef();
701 delete acceleratingStructures[iold];
705 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
706 const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
707 const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
708 const int* fromArray,
709 std::string nameArrayTo,
710 const BBTreeOfDim* myTree)
713 if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
714 ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
715 const double* tc=targetCoords->getConstPointer();
716 int targetSize=targetMesh.getNumberOfCells();
717 int sourceSize=sourceMesh.getNumberOfCells();
718 if (MyGlobals::_Verbose>200)
719 std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
720 std::vector<int> ccI;
722 str=nameArrayTo+"_toArray";
723 cle=Cle1ToStr(str,inew);
726 const BBTreeOfDim* tree;
727 bool cleantree=false;
728 ParaMEDMEM::DataArrayDouble* sourceBBox=0;
729 int dim = targetCoords->getNumberOfComponents();
732 sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
733 tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
737 //first time iold : create and initiate
738 if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
740 if (MyGlobals::_Is0verbose>100)
741 std::cout << "create " << cle << " size " << targetSize << std::endl;
742 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
743 p->alloc(targetSize,1);
745 toArray=p->getPointer();
746 _map_dataarray_int[cle]=p;
748 else //other times iold: refind and complete
750 toArray=_map_dataarray_int.find(cle)->second->getPointer();
753 std::map< int, int > isource2nb; // count coincident elements
754 std::map<int,int>::iterator i2nb;
756 for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
758 std::vector<int> intersectingElems;
759 tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
760 if (intersectingElems.size()!=0)
762 int isourcenode=intersectingElems[0];
763 if ( intersectingElems.size() > 1 )
765 i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
766 isourcenode = intersectingElems[ i2nb->second++ ];
768 if ( isourcenode < sourceSize ) // protection from invalid elements
770 toArray[itargetnode]=fromArray[isourcenode];
771 ccI.push_back(itargetnode);
772 ccI.push_back(isourcenode);
776 if (MyGlobals::_Verbose>200)
777 std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
778 //memories intersection for future same job on fields (if no existing cle=no intersection)
779 str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
780 if (MyGlobals::_Verbose>700)
781 std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
783 _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
785 targetCoords->decrRef();
786 if (cleantree) delete tree;
787 if (sourceBBox !=0) sourceBBox->decrRef();
790 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
792 if (nameArrayTo!="cellFieldDouble")
793 throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
795 std::string nameTo="typeData=6"; //resume the type of field casted
796 // send-recv operations
797 int ioldMax=initialCollection.getMesh().size();
798 int inewMax=this->getMesh().size();
799 int iFieldMax=initialCollection.getFieldDescriptions().size();
800 if (MyGlobals::_Verbose>10)
801 std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
802 //see collection.prepareFieldDescriptions()
803 for (int ifield=0; ifield<iFieldMax; ifield++)
805 std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
806 if (descriptionField.find(nameTo)==std::string::npos)
807 continue; //only nameTo accepted in Fields name description
809 for (int inew=0; inew<inewMax; inew++)
811 for (int iold=0; iold<ioldMax; iold++)
813 //sending arrays for distant domains
814 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
816 int target=_domain_selector->getProcessorID(inew);
817 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
818 if (MyGlobals::_Verbose>10)
819 std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
820 SendDataArrayDouble(field, target);
822 //receiving arrays from distant domains
823 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
825 int source=_domain_selector->getProcessorID(iold);
827 if (MyGlobals::_Verbose>10)
828 std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
829 ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
830 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
835 //local contributions and aggregation
836 for (int inew=0; inew<inewMax; inew++)
838 for (int iold=0; iold<ioldMax; iold++)
839 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
841 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
842 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
848 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
849 ParaMEDMEM::DataArrayDouble* fromArray,
850 std::string nameArrayTo,
851 std::string descriptionField)
852 //here we use 'cellFamily_ccI inew iold' set in remapIntField
854 if (nameArrayTo!="cellFieldDouble")
855 throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
856 std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
858 std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
859 it1=_map_dataarray_int.find(key);
860 if (it1==_map_dataarray_int.end())
862 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
863 std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
866 //create ccI in remapIntField
867 ParaMEDMEM::DataArrayInt *ccI=it1->second;
868 if (MyGlobals::_Verbose>300)
869 std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
871 int nbcell=this->getMesh()[inew]->getNumberOfCells();
872 int nbcomp=fromArray->getNumberOfComponents();
873 int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
875 std::string tag="inewFieldDouble="+IntToStr(inew);
876 key=descriptionField+SerializeFromString(tag);
877 int fromArrayNbOfElem=fromArray->getNbOfElems();
878 int fromArrayNbOfComp=fromArray->getNumberOfComponents();
879 int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
881 if (MyGlobals::_Verbose>1000)
883 std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
884 " fromArray nbOfElems " << fromArrayNbOfElem <<
885 " nbTuples " << fromArray->getNumberOfTuples() <<
886 " nbcells " << fromArrayNbOfCell <<
887 " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
890 ParaMEDMEM::DataArrayDouble* field=0;
891 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
892 it2=_map_dataarray_double.find(key);
893 if (it2==_map_dataarray_double.end())
895 if (MyGlobals::_Verbose>300)
896 std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
897 field=ParaMEDMEM::DataArrayDouble::New();
898 _map_dataarray_double[key]=field;
899 field->alloc(nbcell*nbPtGauss,nbcomp);
900 field->fillWithZero();
905 if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
907 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
908 " trying remap of field double on cells : \n" << descriptionField << std::endl;
915 field->setPartOfValuesAdv(fromArray,ccI);
919 //replaced by setPartOfValuesAdv if nbPtGauss==1
920 int iMax=ccI->getNbOfElems();
921 int* pccI=ccI->getPointer();
922 double* pField=field->getPointer();
923 double* pFrom=fromArray->getPointer();
924 int itarget, isource, delta=nbPtGauss*nbcomp;
925 for (int i=0; i<iMax; i=i+2) //cell
929 if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
930 throw INTERP_KERNEL::Exception("Error field override");
931 int ita=itarget*delta;
932 int iso=isource*delta;
933 for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
938 //================================================================================
940 * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
941 * group (where "n" and "p" are domain IDs)
943 //================================================================================
945 void MEDPARTITIONER::MeshCollection::buildConnectZones()
947 if ( getMeshDimension() < 2 )
950 using ParaMEDMEM::MEDCouplingUMesh;
951 using ParaMEDMEM::DataArrayDouble;
952 using ParaMEDMEM::DataArrayInt;
954 std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
955 int nbMeshes = faceMeshes.size();
957 //preparing bounding box trees for accelerating search of coincident faces
958 std::vector<BBTreeOfDim* > bbTrees(nbMeshes);
959 std::vector<DataArrayDouble*>bbox (nbMeshes);
960 for (int inew = 0; inew < nbMeshes-1; inew++)
961 if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
963 DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner();
964 bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
965 bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
966 bbox[inew]->getConstPointer(),0,0,
967 bbox[inew]->getNumberOfTuples());
971 // loop on domains to find joint faces between them
972 for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
974 for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
976 MEDCouplingUMesh* mesh1 = 0;
977 MEDCouplingUMesh* mesh2 = 0;
978 //MEDCouplingUMesh* recvMesh = 0;
979 bool mesh1Here = true, mesh2Here = true;
980 if (isParallelMode())
983 mesh1Here = _domain_selector->isMyDomain(inew1);
984 mesh2Here = _domain_selector->isMyDomain(inew2);
985 if ( !mesh1Here && mesh2Here )
987 //send mesh2 to domain of mesh1
988 _domain_selector->sendMesh(*faceMeshes[inew2],
989 _domain_selector->getProcessorID(inew1));
991 else if ( mesh1Here && !mesh2Here )
993 //receiving mesh2 from a distant domain
994 _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
995 if ( faceMeshes[ inew2 ] )
996 faceMeshes[ inew2 ]->decrRef();
997 faceMeshes[ inew2 ] = mesh2;
1001 if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
1002 if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
1004 // find coincident faces
1005 std::vector< int > faces1, faces2;
1006 if ( mesh1 && mesh2 )
1008 const DataArrayDouble* coords2 = mesh2->getBarycenterAndOwner();
1009 const double* c2 = coords2->getConstPointer();
1010 const int dim = coords2->getNumberOfComponents();
1011 const int nbFaces2 = mesh2->getNumberOfCells();
1012 const int nbFaces1 = mesh1->getNumberOfCells();
1014 for (int i2 = 0; i2 < nbFaces2; i2++)
1016 std::vector<int> coincFaces;
1017 bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
1018 if (coincFaces.size()!=0)
1020 int i1 = coincFaces[0];
1021 // if ( coincFaces.size() > 1 )
1023 // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
1024 // i1 = coincFaces[ i2nb->second++ ];
1026 if ( i1 < nbFaces1 ) // protection from invalid elements
1028 faces1.push_back( i1 );
1029 faces2.push_back( i2 );
1036 if ( isParallelMode())
1039 if ( mesh1Here && !mesh2Here )
1041 //send faces2 to domain of recvMesh
1042 SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
1044 else if ( !mesh1Here && mesh2Here )
1046 //receiving ids of faces from a domain of mesh1
1047 RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
1052 // recvMesh->decrRef();
1054 // Create group "JOINT_inew1_inew2_Faces" and corresponding families
1055 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1057 createJointGroup( is2nd ? faces2 : faces1,
1058 inew1 , inew2, is2nd );
1061 } // loop on the 2nd domains (inew2)
1062 } // loop on the 1st domains (inew1)
1065 // delete bounding box trees
1066 for (int inew = 0; inew < nbMeshes-1; inew++)
1067 if (isParallelMode() && _domain_selector->isMyDomain(inew))
1069 bbox[inew]->decrRef();
1070 delete bbTrees[inew];
1074 //================================================================================
1076 * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
1077 * \param faces - face ids to include into the group
1078 * \param inew1 - index of the 1st domain
1079 * \param inew2 - index of the 2nd domain
1080 * \param is2nd - in which (1st or 2nd) domain to create the group
1082 //================================================================================
1084 void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
1089 // get the name of JOINT group
1090 std::string groupName;
1092 std::ostringstream oss;
1094 << (is2nd ? inew2 : inew1 ) << "_"
1095 << (is2nd ? inew1 : inew2 ) << "_"
1096 << ( getMeshDimension()==2 ? "Edge" : "Face" );
1097 groupName = oss.str();
1100 // remove existing "JOINT_*" group
1101 _group_info.erase( groupName );
1103 // get family IDs array
1105 int inew = (is2nd ? inew2 : inew1 );
1106 int totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
1107 std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
1108 if ( !_map_dataarray_int.count(cle) )
1110 if ( totalNbFaces > 0 )
1112 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
1113 p->alloc( totalNbFaces, 1 );
1115 famIDs = p->getPointer();
1116 _map_dataarray_int[cle]=p;
1121 famIDs = _map_dataarray_int.find(cle)->second->getPointer();
1123 // find a family ID of an existing JOINT group
1125 std::map<std::string, int>::iterator name2id = _family_info.find( groupName );
1126 if ( name2id != _family_info.end() )
1127 familyID = name2id->second;
1129 // remove faces from the familyID-the family
1130 if ( familyID != 0 && famIDs )
1131 for ( size_t i = 0; i < totalNbFaces; ++i )
1132 if ( famIDs[i] == familyID )
1135 if ( faces.empty() )
1138 if ( familyID == 0 ) // generate a family ID for JOINT group
1140 std::set< int > familyIDs;
1141 for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
1142 familyIDs.insert( name2id->second );
1143 // find the next free family ID
1144 int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
1147 if ( !familyIDs.count( ++familyID ))
1150 while ( freeIdCount > 0 );
1153 // push faces to familyID-th group
1154 if ( faces.back() >= totalNbFaces )
1155 throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
1156 for ( size_t i = 0; i < faces.size(); ++i )
1157 famIDs[ faces[i] ] = familyID;
1159 // register JOINT group and family
1160 _family_info[ groupName ] = familyID; // name of the group and family is same
1161 _group_info [ groupName ].push_back( groupName );
1164 /*! constructing the MESH collection from a distributed file
1166 * \param filename name of the master file containing the list of all the MED files
1168 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
1170 _owns_topology(true),
1172 _domain_selector( 0 ),
1173 _i_non_empty_mesh(-1),
1174 _driver_type(MEDPARTITIONER::Undefined),
1175 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1176 _family_splitting(false),
1177 _create_empty_groups(false),
1182 _driver=new MeshCollectionMedXmlDriver(this);
1183 _driver->read (filename.c_str());
1184 _driver_type = MedXml;
1187 { // Handle all exceptions
1188 if ( _driver ) delete _driver; _driver=0;
1191 _driver=new MeshCollectionMedAsciiDriver(this);
1192 _driver->read (filename.c_str());
1193 _driver_type=MedAscii;
1199 throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
1202 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1203 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1204 _i_non_empty_mesh = idomain;
1207 /*! Constructing the MESH collection from selected domains of a distributed file
1209 * \param filename - name of the master file containing the list of all the MED files
1210 * \param domainSelector - selector of domains to load
1212 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
1214 _owns_topology(true),
1216 _domain_selector( &domainSelector ),
1217 _i_non_empty_mesh(-1),
1218 _driver_type(MEDPARTITIONER::Undefined),
1219 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1220 _family_splitting(false),
1221 _create_empty_groups(false),
1224 std::string myfile=filename;
1225 std::size_t found=myfile.find(".xml");
1226 if (found!=std::string::npos) //file .xml
1230 _driver=new MeshCollectionMedXmlDriver(this);
1231 _driver->read ( filename.c_str(), _domain_selector );
1232 _driver_type = MedXml;
1235 { // Handle all exceptions
1237 throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
1242 found=myfile.find(".med");
1243 if (found!=std::string::npos) //file .med single
1245 //make a temporary file .xml and retry MedXmlDriver
1247 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
1249 <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
1250 <description what=\"\" when=\"\"/>\n \
1252 <mesh name=\"$meshName\"/>\n \
1255 <subdomain number=\"1\"/>\n \
1256 <global_numbering present=\"no\"/>\n \
1259 <subfile id=\"1\">\n \
1260 <name>$fileName</name>\n \
1261 <machine>localhost</machine>\n \
1265 <mesh name=\"$meshName\">\n \
1266 <chunk subdomain=\"1\">\n \
1267 <name>$meshName</name>\n \
1272 std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile.c_str());
1273 xml.replace(xml.find("$fileName"),9,myfile);
1274 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1275 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1276 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1277 std::string nameFileXml(myfile);
1278 nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1279 std::string nameFileXmlDN,nameFileXmlBN;
1280 MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1281 nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1282 if (_domain_selector->rank()==0) //only on to write it
1284 std::ofstream f(nameFileXml.c_str());
1289 if (MyGlobals::_World_Size>1)
1290 MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1294 _driver=new MeshCollectionMedXmlDriver(this);
1295 _driver->read ( nameFileXml.c_str(), _domain_selector );
1296 _driver_type = MedXml;
1299 { // Handle all exceptions
1300 delete _driver; _driver=0;
1301 throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1308 _driver=new MeshCollectionMedAsciiDriver(this);
1309 _driver->read ( filename.c_str(), _domain_selector );
1310 _driver_type=MedAscii;
1316 throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1320 // find non-empty domain mesh
1321 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1322 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1323 _i_non_empty_mesh = idomain;
1327 //check for all proc/file compatibility of _field_descriptions
1329 _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1331 _field_descriptions=MyGlobals::_Field_Descriptions;
1334 catch(INTERP_KERNEL::Exception& e)
1336 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1337 throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1342 //check for all proc/file compatibility of _family_info
1343 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1344 _family_info=DevectorizeToMapOfStringInt(v2);
1346 catch(INTERP_KERNEL::Exception& e)
1348 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1349 throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1354 //check for all proc/file compatibility of _group_info
1355 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1356 _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1358 catch(INTERP_KERNEL::Exception& e)
1360 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1361 throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1366 /*! constructing the MESH collection from a sequential MED-file
1368 * \param filename MED file
1369 * \param meshname name of the mesh that is to be read
1371 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1373 _owns_topology(true),
1375 _domain_selector( 0 ),
1376 _i_non_empty_mesh(-1),
1378 _driver_type(MEDPARTITIONER::MedXml),
1379 _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1380 _family_splitting(false),
1381 _create_empty_groups(false),
1384 try // avoid memory leak in case of inexistent filename
1386 retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1392 throw INTERP_KERNEL::Exception("problem reading .med files");
1394 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1395 _i_non_empty_mesh = 0;
1398 MEDPARTITIONER::MeshCollection::~MeshCollection()
1400 for (int i=0; i<(int)_mesh.size();i++)
1401 if (_mesh[i]!=0) _mesh[i]->decrRef();
1403 for (int i=0; i<(int)_cell_family_ids.size();i++)
1404 if (_cell_family_ids[i]!=0)
1405 _cell_family_ids[i]->decrRef();
1407 for (int i=0; i<(int)_face_mesh.size();i++)
1408 if (_face_mesh[i]!=0)
1409 _face_mesh[i]->decrRef();
1411 for (int i=0; i<(int)_face_family_ids.size();i++)
1412 if (_face_family_ids[i]!=0)
1413 _face_family_ids[i]->decrRef();
1415 for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1416 if ((*it).second!=0)
1417 (*it).second->decrRef();
1419 for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1420 if ((*it).second!=0)
1421 (*it).second->decrRef();
1424 if (_topology!=0 && _owns_topology)
1427 delete _joint_finder;
1431 /*! constructing the MESH collection from a file
1433 * The method creates as many MED-files as there are domains in the
1434 * collection. It also creates a master file that lists all the MED files.
1435 * The MED files created in ths manner contain joints that describe the
1436 * connectivity between subdomains.
1438 * \param filename name of the master file that will contain the list of the MED files
1441 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1443 //building the connect zones necessary for writing joints
1444 if (_topology->nbDomain()>1 && _subdomain_boundary_creates )
1445 buildConnectZones();
1446 //suppresses link with driver so that it can be changed for writing
1449 retrieveDriver()->write (filename.c_str(), _domain_selector);
1452 /*! creates or gets the link to the collection driver
1454 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1458 switch(_driver_type)
1461 _driver=new MeshCollectionMedXmlDriver(this);
1464 _driver=new MeshCollectionMedAsciiDriver(this);
1467 throw INTERP_KERNEL::Exception("Unrecognized driver");
1474 /*! gets an existing driver
1477 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1483 * retrieves the mesh dimension
1485 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1487 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1490 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1493 for (int i=0; i<_mesh.size(); i++)
1500 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1503 for (int i=0; i<_mesh.size(); i++)
1505 if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1510 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1513 for (int i=0; i<_face_mesh.size(); i++)
1515 if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1520 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1525 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1530 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1532 return _mesh[idomain];
1535 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1537 return _face_mesh[idomain];
1540 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1542 return _connect_zones;
1545 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1550 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo)
1554 throw INTERP_KERNEL::Exception("topology is already set");
1560 /*! Method creating the cell graph in serial mode
1562 * \param array returns the pointer to the structure that contains the graph
1563 * \param edgeweight returns the pointer to the table that contains the edgeweights
1564 * (only used if indivisible regions are required)
1566 void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
1571 using std::make_pair;
1574 if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
1575 const ParaMEDMEM::MEDCouplingUMesh* mesh=_mesh[0];
1576 if (MyGlobals::_Verbose>50)
1577 std::cout<<"getting nodal connectivity"<<std::endl;
1579 //looking for reverse nodal connectivity i global numbering
1580 if (isParallelMode() && !_domain_selector->isMyDomain(0))
1583 vector<int> index(1,0);
1585 array=new MEDPARTITIONER::SkyLineArray(index,value);
1589 int meshDim = mesh->getMeshDimension();
1591 ParaMEDMEM::DataArrayInt* indexr=ParaMEDMEM::DataArrayInt::New();
1592 ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1593 int nbNodes=mesh->getNumberOfNodes();
1594 mesh->getReverseNodalConnectivity(revConn,indexr);
1595 //problem saturation over 1 000 000 nodes for 1 proc
1596 if (MyGlobals::_Verbose>100)
1597 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1598 const int* indexr_ptr=indexr->getConstPointer();
1599 const int* revConn_ptr=revConn->getConstPointer();
1601 const ParaMEDMEM::DataArrayInt* index;
1602 const ParaMEDMEM::DataArrayInt* conn;
1603 conn=mesh->getNodalConnectivity();
1604 index=mesh->getNodalConnectivityIndex();
1605 int nbCells=mesh->getNumberOfCells();
1607 if (MyGlobals::_Verbose>100)
1608 std::cout << "proc " << MyGlobals::_Rank << " : getNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1609 const int* index_ptr=index->getConstPointer();
1610 const int* conn_ptr=conn->getConstPointer();
1612 //creating graph arcs (cell to cell relations)
1613 //arcs are stored in terms of (index,value) notation
1616 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
1617 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
1619 //warning here one node have less than or equal effective number of cell with it
1620 //but cell could have more than effective nodes
1621 //because other equals nodes in other domain (with other global inode)
1622 if (MyGlobals::_Verbose>50)
1623 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
1625 vector <int> cell2cell_index(nbCells+1,0);
1626 vector <int> cell2cell;
1627 cell2cell.reserve(3*nbCells);
1629 for (int icell=0; icell<nbCells;icell++)
1631 map<int,int > counter;
1632 for (int iconn=index_ptr[icell]; iconn<index_ptr[icell+1];iconn++)
1634 int inode=conn_ptr[iconn];
1635 for (int iconnr=indexr_ptr[inode]; iconnr<indexr_ptr[inode+1];iconnr++)
1637 int icell2=revConn_ptr[iconnr];
1638 map<int,int>::iterator iter=counter.find(icell2);
1639 if (iter!=counter.end()) (iter->second)++;
1640 else counter.insert(make_pair(icell2,1));
1643 for (map<int,int>::const_iterator iter=counter.begin();
1644 iter!=counter.end();
1646 if (iter->second >= meshDim)
1648 cell2cell_index[icell+1]++;
1649 cell2cell.push_back(iter->first);
1657 cell2cell_index[0]=0;
1658 for (int icell=0; icell<nbCells;icell++)
1659 cell2cell_index[icell+1]=cell2cell_index[icell]+cell2cell_index[icell+1];
1662 if (MyGlobals::_Verbose>50)
1663 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
1665 //filling up index and value to create skylinearray structure
1666 array=new MEDPARTITIONER::SkyLineArray(cell2cell_index,cell2cell);
1668 if (MyGlobals::_Verbose>100)
1670 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
1671 cell2cell_index.size()-1 << " " << cell2cell.size() << std::endl;
1672 int max=cell2cell_index.size()>15?15:cell2cell_index.size();
1673 if (cell2cell_index.size()>1)
1675 for (int i=0; i<max; ++i)
1676 std::cout<<cell2cell_index[i]<<" ";
1677 std::cout << "... " << cell2cell_index[cell2cell_index.size()-1] << std::endl;
1678 for (int i=0; i<max; ++i)
1679 std::cout<< cell2cell[i] << " ";
1680 int ll=cell2cell_index[cell2cell_index.size()-1]-1;
1681 std::cout << "... (" << ll << ") " << cell2cell[ll-1] << " " << cell2cell[ll] << std::endl;
1686 /*! Method creating the cell graph in multidomain mode
1688 * \param array returns the pointer to the structure that contains the graph
1689 * \param edgeweight returns the pointer to the table that contains the edgeweights
1690 * (only used if indivisible regions are required)
1692 void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
1694 using std::multimap;
1696 using std::make_pair;
1699 std::multimap< int, int > node2cell;
1700 std::map< pair<int,int>, int > cell2cellcounter;
1701 std::multimap<int,int> cell2cell;
1703 std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
1704 int nbdomain=_topology->nbDomain();
1706 if (isParallelMode())
1708 _joint_finder=new JointFinder(*this);
1709 _joint_finder->findCommonDistantNodes();
1710 commonDistantNodes=_joint_finder->getDistantNodeCell();
1713 if (MyGlobals::_Verbose>500)
1714 _joint_finder->print();
1717 if (MyGlobals::_Verbose>50)
1718 std::cout<<"getting nodal connectivity"<<std::endl;
1719 //looking for reverse nodal connectivity i global numbering
1721 for (int idomain=0; idomain<nbdomain; idomain++)
1723 if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1725 meshDim = _mesh[idomain]->getMeshDimension();
1727 ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
1728 ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1729 int nbNodes=_mesh[idomain]->getNumberOfNodes();
1730 _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1731 //problem saturation over 1 000 000 nodes for 1 proc
1732 if (MyGlobals::_Verbose>100)
1733 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1734 int* index_ptr=index->getPointer();
1735 int* revConnPtr=revConn->getPointer();
1736 for (int i=0; i<nbNodes; i++)
1738 for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
1740 int globalNode=_topology->convertNodeToGlobal(idomain,i);
1741 int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
1742 node2cell.insert(make_pair(globalNode, globalCell));
1748 for (int iother=0; iother<nbdomain; iother++)
1750 std::multimap<int,int>::iterator it;
1751 int isource=idomain;
1753 for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
1754 it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
1756 int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
1757 int globalCell=(*it).second;
1758 node2cell.insert(make_pair(globalNode, globalCell));
1764 //creating graph arcs (cell to cell relations)
1765 //arcs are stored in terms of (index,value) notation
1768 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
1769 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
1771 //warning here one node have less than or equal effective number of cell with it
1772 //but cell could have more than effective nodes
1773 //because other equals nodes in other domain (with other global inode)
1774 if (MyGlobals::_Verbose>50)
1775 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
1777 for (int inode=0;inode<_topology->nbNodes();inode++)
1779 typedef multimap<int,int>::const_iterator MI;
1780 std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
1781 for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
1782 for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
1784 int icell1=cell1->second;
1785 int icell2=cell2->second;
1786 if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
1787 std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
1788 if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
1789 else (it->second)++;
1792 // for (int icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
1794 // typedef multimap<int,int>::const_iterator MI;
1795 // std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
1796 // for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
1798 // std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
1799 // for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
1801 // int icell2=cell2->second;
1802 // std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
1803 // if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
1804 // else (it->second)++;
1810 //converting the counter to a multimap structure
1811 for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
1812 it!=cell2cellcounter.end();
1814 if (it->second>=meshDim)
1816 cell2cell.insert(std::make_pair(it->first.first,it->first.second));
1817 cell2cell.insert(std::make_pair(it->first.second, it->first.first));
1821 if (MyGlobals::_Verbose>50)
1822 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
1823 //filling up index and value to create skylinearray structure
1824 std::vector <int> index,value;
1828 for (int idomain=0; idomain<nbdomain; idomain++)
1830 if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
1831 int nbCells=_mesh[idomain]->getNumberOfCells();
1832 for (int icell=0; icell<nbCells; icell++)
1835 int globalCell=_topology->convertCellToGlobal(idomain,icell);
1836 multimap<int,int>::iterator it;
1837 pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
1838 ret=cell2cell.equal_range(globalCell);
1839 for (it=ret.first; it!=ret.second; ++it)
1841 int ival=(*it).second; //no adding one existing yet
1842 for (int i=idep ; i<idep+size ; i++)
1851 value.push_back(ival);
1855 idep=index[index.size()-1]+size;
1856 index.push_back(idep);
1860 array=new MEDPARTITIONER::SkyLineArray(index,value);
1862 if (MyGlobals::_Verbose>100)
1864 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
1865 index.size()-1 << " " << value.size() << std::endl;
1866 int max=index.size()>15?15:index.size();
1869 for (int i=0; i<max; ++i)
1870 std::cout<<index[i]<<" ";
1871 std::cout << "... " << index[index.size()-1] << std::endl;
1872 for (int i=0; i<max; ++i)
1873 std::cout<< value[i] << " ";
1874 int ll=index[index.size()-1]-1;
1875 std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
1882 /*! Creates the partition corresponding to the cell graph and the partition number
1884 * \param nbdomain number of subdomains for the newly created graph
1886 * returns a topology based on the new graph
1888 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
1889 Graph::splitter_type split,
1890 const std::string& options_string,
1891 int *user_edge_weights,
1892 int *user_vertices_weights)
1894 if (MyGlobals::_Verbose>10)
1895 std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
1898 throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
1899 MEDPARTITIONER::SkyLineArray* array=0;
1902 if (_topology->nbDomain()>1 || isParallelMode())
1903 buildParallelCellGraph(array,edgeweights);
1905 buildCellGraph(array,edgeweights);
1907 Graph* cellGraph = 0;
1911 if ( isParallelMode() && MyGlobals::_World_Size > 1 )
1913 #ifdef MED_ENABLE_PARMETIS
1914 if (MyGlobals::_Verbose>10)
1915 std::cout << "ParMETISGraph" << std::endl;
1916 cellGraph=new ParMETISGraph(array,edgeweights);
1921 #ifdef MED_ENABLE_METIS
1922 if (MyGlobals::_Verbose>10)
1923 std::cout << "METISGraph" << std::endl;
1924 cellGraph=new METISGraph(array,edgeweights);
1928 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
1932 #ifdef MED_ENABLE_SCOTCH
1933 if (MyGlobals::_Verbose>10)
1934 std::cout << "SCOTCHGraph" << std::endl;
1935 cellGraph=new SCOTCHGraph(array,edgeweights);
1937 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
1942 //!user-defined weights
1943 if (user_edge_weights!=0)
1944 cellGraph->setEdgesWeights(user_edge_weights);
1945 if (user_vertices_weights!=0)
1946 cellGraph->setVerticesWeights(user_vertices_weights);
1948 if (MyGlobals::_Is0verbose>10)
1949 std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
1950 cellGraph->partGraph(nbdomain, options_string, _domain_selector);
1952 if (MyGlobals::_Is0verbose>10)
1953 std::cout << "building new topology" << std::endl;
1954 //cellGraph is a shared pointer
1955 Topology *topology=0;
1956 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
1958 delete [] edgeweights;
1960 if (MyGlobals::_Verbose>11)
1961 std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
1965 /*! Creates a topology for a partition specified by the user
1967 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
1969 * returns a topology based on the new partition
1971 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
1973 MEDPARTITIONER::SkyLineArray* array=0;
1976 if ( _topology->nbDomain()>1)
1977 buildParallelCellGraph(array,edgeweights);
1979 buildCellGraph(array,edgeweights);
1982 std::set<int> domains;
1983 for (int i=0; i<_topology->nbCells(); i++)
1985 domains.insert(partition[i]);
1987 cellGraph=new UserGraph(array, partition, _topology->nbCells());
1989 //cellGraph is a shared pointer
1990 Topology *topology=0;
1991 int nbdomain=domains.size();
1992 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
1993 // if (array!=0) delete array;
1998 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
2000 for (int i=0; i<_topology->nbDomain(); i++)
2002 std::ostringstream oss;
2004 if (!isParallelMode() || _domain_selector->isMyDomain(i))
2005 _mesh[i]->setName(oss.str().c_str());
2009 ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
2010 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
2011 //something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name.c_str(),f1->getMesh()->getName(),0,f1->getName(),0,1);
2013 int rank=MyGlobals::_Rank;
2014 std::string tag="ioldFieldDouble="+IntToStr(iold);
2015 std::string descriptionIold=descriptionField+SerializeFromString(tag);
2016 if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
2018 if (MyGlobals::_Verbose>300)
2019 std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
2020 ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
2023 if (MyGlobals::_Verbose>200)
2024 std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
2025 std::string description, fileName, meshName, fieldName;
2026 int typeField, DT, IT, entity;
2027 fileName=MyGlobals::_File_Names[iold];
2028 if (MyGlobals::_Verbose>10)
2029 std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
2030 FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
2031 meshName=MyGlobals::_Mesh_Names[iold];
2033 ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
2034 fileName.c_str(), meshName.c_str(), 0, fieldName.c_str(), DT, IT);
2036 ParaMEDMEM::DataArrayDouble* res=f2->getArray();
2037 //to know names of components
2038 std::vector<std::string> browse=BrowseFieldDouble(f2);
2039 std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
2040 if (MyGlobals::_Verbose>10)
2041 std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
2042 MyGlobals::_General_Informations.push_back(localFieldInformation);
2045 _map_dataarray_double[descriptionIold]=res;
2049 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
2050 //to have unique valid fields names/pointers/descriptions for partitionning
2051 //filter _field_descriptions to be in all procs compliant and equal
2053 int nbfiles=MyGlobals::_File_Names.size(); //nb domains
2054 std::vector<std::string> r2;
2055 //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
2056 for (int i=0; i<(int)_field_descriptions.size(); i++)
2058 std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
2059 for (int ii=0; ii<(int)r1.size(); ii++)
2060 r2.push_back(r1[ii]);
2062 //here vector(procs*fields) of serialised vector(description) data
2063 _field_descriptions=r2;
2064 int nbfields=_field_descriptions.size(); //on all domains
2065 if ((nbfields%nbfiles)!=0)
2067 if (MyGlobals::_Rank==0)
2069 std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
2070 << "fileMedNames :" << std::endl
2071 << ReprVectorOfString(MyGlobals::_File_Names)
2072 << "field_descriptions :" << std::endl
2073 << ReprVectorOfString(MyGlobals::_Field_Descriptions);
2075 throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
2077 _field_descriptions.resize(nbfields/nbfiles);
2078 for (int i=0; i<(int)_field_descriptions.size(); i++)
2080 std::string str=_field_descriptions[i];
2081 str=EraseTagSerialized(str,"idomain=");
2082 str=EraseTagSerialized(str,"fileName=");
2083 _field_descriptions[i]=str;
2087 //returns true if inodes of a face are in inodes of a cell
2088 bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >& inodesCell)
2091 int nbok=inodesFace.size();
2092 for (int i=0; i<nbok; i++)
2094 int ii=inodesFace[i];
2096 std::cout << "isFaceOncell problem inodeface<0" << std::endl;
2097 for (int j=0; j<(int)inodesCell.size(); j++)
2099 if (ii==inodesCell[j])
2102 break; //inode of face found
2106 break; //inode of face not found do not continue...
2108 return (ires==nbok);
2111 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
2113 for (int inew=0; inew<_topology->nbDomain(); inew++)
2115 if (!isParallelMode() || _domain_selector->isMyDomain(inew))
2117 if (MyGlobals::_Verbose>200)
2118 std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
2119 ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
2120 ParaMEDMEM::MEDCouplingUMesh* mfac=_face_mesh[inew];
2122 //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
2123 std::vector<int> nodeIds;
2124 getNodeIds(*mcel, *mfac, nodeIds);
2125 if (nodeIds.size()==0)
2126 continue; //one empty mesh nothing to do
2128 ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
2129 ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New();
2130 mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
2131 int *revC=revNodalCel->getPointer();
2132 int *revIndxC=revNodalIndxCel->getPointer();
2134 std::vector< int > faceOnCell;
2135 std::vector< int > faceNotOnCell;
2136 int nbface=mfac->getNumberOfCells();
2137 for (int iface=0; iface<nbface; iface++)
2140 std::vector< int > inodesFace;
2141 mfac->getNodeIdsOfCell(iface, inodesFace);
2142 int nbnodFace=inodesFace.size();
2143 if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
2144 continue; // invalid node ids
2145 //set inodesFace in mcel
2147 for (int i=0; i<nbnodFace; i++)
2148 nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
2149 if ( nbok != nbnodFace )
2151 int inod=inodesFace[0];
2154 std::cout << "filterFaceOnCell problem 1" << std::endl;
2157 int nbcell=revIndxC[inod+1]-revIndxC[inod];
2158 for (int j=0; j<nbcell; j++) //look for each cell with inod
2160 int icel=revC[revIndxC[inod]+j];
2161 std::vector< int > inodesCell;
2162 mcel->getNodeIdsOfCell(icel, inodesCell);
2163 ok=isFaceOncell(inodesFace, inodesCell);
2168 faceOnCell.push_back(iface);
2172 faceNotOnCell.push_back(iface);
2173 if (MyGlobals::_Is0verbose>300)
2174 std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
2178 revNodalCel->decrRef();
2179 revNodalIndxCel->decrRef();
2181 // std::string keyy;
2182 // keyy=Cle1ToStr("filterFaceOnCell",inew);
2183 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
2184 // keyy=Cle1ToStr("filterNotFaceOnCell",inew);
2185 // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
2187 // filter the face mesh
2188 if ( faceOnCell.empty() )
2189 _face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
2191 _face_mesh[inew] = (ParaMEDMEM::MEDCouplingUMesh *)
2192 mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
2195 // filter the face families
2196 std::string key = Cle1ToStr("faceFamily_toArray",inew);
2197 if ( getMapDataArrayInt().count( key ))
2199 ParaMEDMEM::DataArrayInt * & fam = getMapDataArrayInt()[ key ];
2200 ParaMEDMEM::DataArrayInt * famFilter = ParaMEDMEM::DataArrayInt::New();
2201 famFilter->alloc(faceOnCell.size(),1);
2202 int* pfamFilter = famFilter->getPointer();
2203 int* pfam = fam->getPointer();
2204 for ( size_t i=0; i<faceOnCell.size(); i++ )
2205 pfamFilter[i]=pfam[faceOnCell[i]];