1 // Copyright (C) 2007-2012 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"
51 #if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS)
52 #include "MEDPARTITIONER_MetisGraph.hxx"
55 #ifdef MED_ENABLE_SCOTCH
56 #include "MEDPARTITIONER_ScotchGraph.hxx"
66 MEDPARTITIONER::MeshCollection::MeshCollection()
68 _owns_topology(false),
70 _domain_selector( 0 ),
71 _i_non_empty_mesh(-1),
72 _driver_type(MEDPARTITIONER::MedXml),
73 _subdomain_boundary_creates(false),
74 _family_splitting(false),
75 _create_empty_groups(false),
80 /*!constructor creating a new mesh collection (mesh series + topology)
81 *from an old collection and a new topology
83 * On output, the constructor has built the meshes corresponding to the new mesh collection.
84 * The new topology has been updated so that face and node mappings are included.
85 * The families have been cast to their projections in the new topology.
87 * \param initial_collection collection from which the data (coordinates, connectivity) are taken
88 * \param topology topology containing the cell mappings
91 MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection,
93 bool family_splitting,
94 bool create_empty_groups)
95 : _topology(topology),
96 _owns_topology(false),
98 _domain_selector( initialCollection._domain_selector ),
99 _i_non_empty_mesh(-1),
100 _name(initialCollection._name),
101 _driver_type(MEDPARTITIONER::MedXml),
102 _subdomain_boundary_creates(false),
103 _family_splitting(family_splitting),
104 _create_empty_groups(create_empty_groups),
107 std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
108 castCellMeshes(initialCollection, new2oldIds);
110 //defining the name for the collection and the underlying meshes
111 setName(initialCollection.getName());
118 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
119 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
121 if (MyGlobals::_Is0verbose)
122 std::cout<<"treating faces"<<std::endl;
123 NodeMapping nodeMapping;
124 //nodeMapping contains the mapping between old nodes and new nodes
125 // (iolddomain,ioldnode)->(inewdomain,inewnode)
126 createNodeMapping(initialCollection, nodeMapping);
127 std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
128 castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
134 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
135 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
137 if (MyGlobals::_Is0verbose)
139 if (isParallelMode())
140 std::cout << "ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
142 std::cout << "NOT ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
144 if (MyGlobals::_Is0verbose>10)
145 std::cout<<"treating cell and face families"<<std::endl;
147 castIntField(initialCollection.getMesh(),
149 initialCollection.getCellFamilyIds(),
151 castIntField(initialCollection.getFaceMesh(),
153 initialCollection.getFaceFamilyIds(),
158 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
159 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
161 if (MyGlobals::_Is0verbose)
162 std::cout << "treating groups" << std::endl;
163 _family_info=initialCollection.getFamilyInfo();
164 _group_info=initialCollection.getGroupInfo();
167 if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
168 MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
170 if (MyGlobals::_Is0verbose)
171 std::cout << "treating fields" << std::endl;
172 castAllFields(initialCollection,"cellFieldDouble");
173 if (_i_non_empty_mesh<0)
175 for (int i=0; i<_mesh.size(); i++)
179 _i_non_empty_mesh=i; //first existing one local
188 Creates the meshes using the topology underlying he mesh collection and the mesh data
189 coming from the ancient collection
190 \param initialCollection collection from which the data is extracted to create the new meshes
193 void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
194 std::vector<std::vector<std::vector<int> > >& new2oldIds)
196 if (MyGlobals::_Verbose>10)
197 std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
199 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes");
201 int nbNewDomain=_topology->nbDomain();
202 int nbOldDomain=initialCollection.getTopology()->nbDomain();
204 _mesh.resize(nbNewDomain);
205 int rank=MyGlobals::_Rank;
206 //splitting the initial domains into smaller bits
207 std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
208 splitMeshes.resize(nbNewDomain);
209 for (int inew=0; inew<nbNewDomain; inew++)
211 splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
214 for (int iold=0; iold<nbOldDomain; iold++)
216 if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
218 int size=(initialCollection._mesh)[iold]->getNumberOfCells();
219 std::vector<int> globalids(size);
220 initialCollection.getTopology()->getCellList(iold, &globalids[0]);
221 std::vector<int> ilocalnew(size); //local
222 std::vector<int> ipnew(size); //idomain old
223 _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
225 new2oldIds[iold].resize(nbNewDomain);
226 for (int i=0; i<(int)ilocalnew.size(); i++)
228 new2oldIds[iold][ipnew[i]].push_back(i);
230 for (int inew=0; inew<nbNewDomain; inew++)
232 splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)
233 (initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
234 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
236 if (MyGlobals::_Verbose>400)
237 std::cout<< "proc " << rank << " : a splitMesh iold inew NbCells " << iold << " " << inew << " "
238 << splitMeshes[inew][iold]->getNumberOfCells() << std::endl;
243 if (isParallelMode())
245 //if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
246 for (int iold=0; iold<nbOldDomain; iold++)
247 for(int inew=0; inew<nbNewDomain; inew++)
249 if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
252 if(initialCollection._domain_selector->isMyDomain(iold))
253 _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew));
255 if (_domain_selector->isMyDomain(inew))
256 _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
262 //fusing the split meshes
263 if (MyGlobals::_Verbose>200)
264 std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
265 for (int inew=0; inew<nbNewDomain ;inew++)
267 std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
269 for (int i=0; i<(int)splitMeshes[inew].size(); i++)
270 if (splitMeshes[inew][i]!=0)
271 if (splitMeshes[inew][i]->getNumberOfCells()>0)
272 meshes.push_back(splitMeshes[inew][i]);
274 if (!isParallelMode()||_domain_selector->isMyDomain(inew))
276 if (meshes.size()==0)
278 _mesh[inew]=CreateEmptyMEDCouplingUMesh();
279 std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
283 _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes);
288 ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
289 array->decrRef(); // array is not used in this case
291 _mesh[inew]->zipCoords();
295 for (int i=0;i<(int)splitMeshes[inew].size();i++)
296 if (splitMeshes[inew][i]!=0)
297 splitMeshes[inew][i]->decrRef();
299 if (MyGlobals::_Verbose>300)
300 std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
304 \param initialCollection source mesh collection
305 \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection
307 void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
310 using std::make_pair;
311 // NodeMapping reverseNodeMapping;
312 for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
317 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
319 // std::map<pair<double,pair<double, double> >, int > nodeClassifier;
320 int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
321 bbox=new double[nvertices*2*3];
322 ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
323 double* coordsPtr=coords->getPointer();
325 for (int i=0; i<nvertices*3;i++)
327 bbox[i*2]=coordsPtr[i]-1e-8;
328 bbox[i*2+1]=coordsPtr[i]+1e-8;
330 tree=new BBTree<3>(bbox,0,0,nvertices,1e-9);
333 for (int inew=0; inew<_topology->nbDomain(); inew++)
336 //sending meshes for parallel computation
337 if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))
338 _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
339 else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
341 ParaMEDMEM::MEDCouplingUMesh* mesh;
342 _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
343 ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
344 for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
346 double* coordsPtr=coords->getPointer()+inode*3;
348 tree->getElementsAroundPoint(coordsPtr,elems);
349 if (elems.size()==0) continue;
350 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
354 else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
356 if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
359 ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
360 for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
362 double* coordsPtr=coords->getPointer()+inode*3;
364 tree->getElementsAroundPoint(coordsPtr,elems);
365 if (elems.size()==0) continue;
366 nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
370 if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
379 void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
382 if (!&meshOne || !&meshTwo) return; //empty or not existing
385 int nv1=meshOne.getNumberOfNodes();
386 bbox=new double[nv1*6];
387 ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
388 double* coordsPtr=coords->getPointer();
389 for (int i=0; i<nv1*3; i++)
391 bbox[i*2]=coordsPtr[i]-1e-8;
392 bbox[i*2+1]=coordsPtr[i]+1e-8;
394 tree=new BBTree<3>(bbox,0,0,nv1,1e-9);
396 int nv2=meshTwo.getNumberOfNodes();
397 nodeIds.resize(nv2,-1);
398 coords=meshTwo.getCoords();
399 for (int inode=0; inode<nv2; inode++)
401 double* coordsPtr2=coords->getPointer()+inode*3;
403 tree->getElementsAroundPoint(coordsPtr2,elems);
404 if (elems.size()==0) continue;
405 nodeIds[inode]=elems[0];
412 creates the face meshes on the new domains from the faces on the old domain and the node mapping
413 faces at the interface are duplicated
415 void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
416 const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
417 std::vector<std::vector<std::vector<int> > >& new2oldIds)
420 //splitMeshes structure will contain the partition of
421 //the old faces on the new ones
422 //splitMeshes[4][2] contains the faces from old domain 2
423 //that have to be added to domain 4
429 using std::make_pair;
431 if (MyGlobals::_Verbose>10)
432 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
434 throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
436 int nbNewDomain=_topology->nbDomain();
437 int nbOldDomain=initialCollection.getTopology()->nbDomain();
439 vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
440 vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
442 vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
444 splitMeshes.resize(nbNewDomain);
445 for (int inew=0; inew<nbNewDomain; inew++)
447 splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
449 new2oldIds.resize(nbOldDomain);
450 for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
452 //init null pointer for empty meshes
453 for (int inew=0; inew<nbNewDomain; inew++)
455 for (int iold=0; iold<nbOldDomain; iold++)
457 splitMeshes[inew][iold]=0; //null for empty meshes
458 new2oldIds[iold][inew].clear();
462 //loop over the old domains to analyse the faces and decide
463 //on which new domain they belong
464 for (int iold=0; iold<nbOldDomain; iold++)
466 if (isParallelMode() && !initialCollection._domain_selector->isMyDomain(iold)) continue;
467 if (MyGlobals::_Verbose>400)
468 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes iodDomain "<<iold<<std::endl;
469 //initial face mesh known : in my domain
470 if (meshesCastFrom[iold] != 0)
472 for (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
475 meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
479 //analysis of element ielem
480 //counters are set for the element
481 //for each source node, the mapping is interrogated and the domain counters
482 //are incremented for each target node
483 //the face is considered as going to target domains if the counter of the domain
484 //is equal to the number of nodes
485 for (int inode=0; inode<(int)nodes.size(); inode++)
487 typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
488 int mynode=nodes[inode];
490 pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
491 for (MI iter=myRange.first; iter!=myRange.second; iter++)
493 int inew=iter->second.first;
494 if (faces.find(inew)==faces.end())
501 for (map<int,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
503 if (iter->second==(int)nodes.size())
504 //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
505 //it is not sure here...
506 //done before writeMedfile on option?... see filterFaceOnCell()
507 new2oldIds[iold][iter->first].push_back(ielem);
511 //creating the splitMeshes from the face ids
512 for (int inew=0; inew<nbNewDomain; inew++)
514 if (meshesCastFrom[iold]->getNumberOfCells() > 0)
516 splitMeshes[inew][iold]=
517 (ParaMEDMEM::MEDCouplingUMesh*)
518 ( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
519 &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
522 splitMeshes[inew][iold]->zipCoords();
526 splitMeshes[inew][iold]=CreateEmptyMEDCouplingUMesh();
532 std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
538 if (isParallelMode())
540 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
541 for (int iold=0; iold<nbOldDomain; iold++)
542 for (int inew=0; inew<nbNewDomain; inew++)
544 if (initialCollection._domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
546 if (splitMeshes[inew][iold] != 0)
548 _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew));
549 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes send " <<inew<<" "<<iold<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
553 _domain_selector->sendMesh(*(empty), _domain_selector->getProcessorID(inew));
554 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes sen0 " <<inew<<" "<<iold<<std::endl;
557 if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
558 _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
560 if (splitMeshes[inew][iold])
561 nb=splitMeshes[inew][iold]->getNumberOfCells();
562 //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
568 //fusing the split meshes
569 if (MyGlobals::_Verbose>200)
570 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes fusing" << std::endl;
571 meshesCastTo.resize(nbNewDomain);
572 for (int inew=0; inew<nbNewDomain; inew++)
574 vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
575 for (int iold=0; iold<nbOldDomain; iold++)
577 ParaMEDMEM::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
579 if (umesh->getNumberOfCells()>0)
580 myMeshes.push_back(umesh);
583 if (myMeshes.size()>0)
585 meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
589 ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
590 meshesCastTo[inew]=empty;
592 for (int iold=0; iold<nbOldDomain; iold++)
593 if (splitMeshes[inew][iold]!=0)
594 splitMeshes[inew][iold]->decrRef();
596 if (MyGlobals::_Verbose>300)
597 std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
602 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
603 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
604 std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
605 std::string nameArrayTo)
609 int ioldMax=meshesCastFrom.size();
610 int inewMax=meshesCastTo.size();
613 //preparing bounding box trees for accelerating source-target node identifications
614 if (MyGlobals::_Verbose>99)
615 std::cout<<"making accelerating structures"<<std::endl;
616 std::vector<BBTree<3,int>* > acceleratingStructures(ioldMax);
617 std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
618 for (int iold =0; iold< ioldMax; iold++)
619 if (isParallelMode() && _domain_selector->isMyDomain(iold))
621 ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
622 bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
623 acceleratingStructures[iold]=new BBTree<3,int> (bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
626 // send-recv operations
628 for (int inew=0; inew<inewMax; inew++)
630 for (int iold=0; iold<ioldMax; iold++)
632 //sending arrays for distant domains
633 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
636 _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
638 int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
640 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
641 if (size>0) //no empty
643 sendIds.resize(size);
644 std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
649 sendIds.resize(size);
651 SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
653 //receiving arrays from distant domains
654 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
658 ParaMEDMEM::MEDCouplingUMesh* recvMesh;
659 _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
661 if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
662 RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
663 remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
664 recvMesh->decrRef(); //cww is it correct?
670 //local contributions and aggregation
671 for (int inew=0; inew<inewMax; inew++)
673 for (int iold=0; iold<ioldMax; iold++)
674 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
676 remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
679 for (int iold=0; iold<ioldMax;iold++)
680 if (isParallelMode() && _domain_selector->isMyDomain(iold))
682 bbox[iold]->decrRef();
683 delete acceleratingStructures[iold];
687 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
688 const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
689 const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
690 const int* fromArray,
691 std::string nameArrayTo,
692 const BBTree<3,int>* myTree)
695 if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
696 ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
697 const double* tc=targetCoords->getConstPointer();
698 int targetSize=targetMesh.getNumberOfCells();
699 int sourceSize=sourceMesh.getNumberOfCells();
700 if (MyGlobals::_Verbose>200)
701 std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
702 std::vector<int> ccI;
704 str=nameArrayTo+"_toArray";
705 cle=Cle1ToStr(str,inew);
708 const BBTree<3>* tree;
709 bool cleantree=false;
710 ParaMEDMEM::DataArrayDouble* sourceBBox=0;
713 sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
714 tree=new BBTree<3> (sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
718 //first time iold : create and initiate
719 if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
721 if (MyGlobals::_Is0verbose>100)
722 std::cout << "create " << cle << " size " << targetSize << std::endl;
723 ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
724 p->alloc(targetSize,1);
726 toArray=p->getPointer();
727 _map_dataarray_int[cle]=p;
729 else //other times iold: refind and complete
731 toArray=_map_dataarray_int.find(cle)->second->getPointer();
734 for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
736 std::vector<int> intersectingElems;
737 tree->getElementsAroundPoint(tc+itargetnode*3,intersectingElems); // to be changed in 2D
738 if (intersectingElems.size()!=0)
740 int isourcenode=intersectingElems[0];
741 toArray[itargetnode]=fromArray[isourcenode];
742 ccI.push_back(itargetnode);
743 ccI.push_back(isourcenode);
746 if (MyGlobals::_Verbose>200)
747 std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
748 //memories intersection for future same job on fields (if no existing cle=no intersection)
749 str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
750 if (MyGlobals::_Verbose>700)
751 std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
753 _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
755 targetCoords->decrRef();
756 if (cleantree) delete tree;
757 if (sourceBBox !=0) sourceBBox->decrRef();
761 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
763 if (nameArrayTo!="cellFieldDouble")
764 throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
766 std::string nameTo="typeData=6"; //resume the type of field casted
767 // send-recv operations
768 int ioldMax=initialCollection.getMesh().size();
769 int inewMax=this->getMesh().size();
770 int iFieldMax=initialCollection.getFieldDescriptions().size();
771 if (MyGlobals::_Verbose>10)
772 std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
773 //see collection.prepareFieldDescriptions()
774 for (int ifield=0; ifield<iFieldMax; ifield++)
776 std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
777 if (descriptionField.find(nameTo)==std::string::npos)
778 continue; //only nameTo accepted in Fields name description
780 for (int inew=0; inew<inewMax; inew++)
782 for (int iold=0; iold<ioldMax; iold++)
784 //sending arrays for distant domains
785 if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
787 int target=_domain_selector->getProcessorID(inew);
788 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
789 if (MyGlobals::_Verbose>10)
790 std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
791 SendDataArrayDouble(field, target);
793 //receiving arrays from distant domains
794 if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
796 int source=_domain_selector->getProcessorID(iold);
798 if (MyGlobals::_Verbose>10)
799 std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
800 ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
801 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
806 //local contributions and aggregation
807 for (int inew=0; inew<inewMax; inew++)
809 for (int iold=0; iold<ioldMax; iold++)
810 if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
812 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
813 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
819 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
820 ParaMEDMEM::DataArrayDouble* fromArray,
821 std::string nameArrayTo,
822 std::string descriptionField)
823 //here we use 'cellFamily_ccI inew iold' set in remapIntField
825 if (nameArrayTo!="cellFieldDouble")
826 throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
827 std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
829 std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
830 it1=_map_dataarray_int.find(key);
831 if (it1==_map_dataarray_int.end())
833 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
834 std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
837 //create ccI in remapIntField
838 ParaMEDMEM::DataArrayInt *ccI=it1->second;
839 if (MyGlobals::_Verbose>300)
840 std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
842 int nbcell=this->getMesh()[inew]->getNumberOfCells();
843 int nbcomp=fromArray->getNumberOfComponents();
844 int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
846 std::string tag="inewFieldDouble="+IntToStr(inew);
847 key=descriptionField+SerializeFromString(tag);
848 int fromArrayNbOfElem=fromArray->getNbOfElems();
849 int fromArrayNbOfComp=fromArray->getNumberOfComponents();
850 int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
852 if (MyGlobals::_Verbose>1000)
854 std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
855 " fromArray nbOfElems " << fromArrayNbOfElem <<
856 " nbTuples " << fromArray->getNumberOfTuples() <<
857 " nbcells " << fromArrayNbOfCell <<
858 " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
861 ParaMEDMEM::DataArrayDouble* field=0;
862 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
863 it2=_map_dataarray_double.find(key);
864 if (it2==_map_dataarray_double.end())
866 if (MyGlobals::_Verbose>300)
867 std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
868 field=ParaMEDMEM::DataArrayDouble::New();
869 _map_dataarray_double[key]=field;
870 field->alloc(nbcell*nbPtGauss,nbcomp);
871 field->fillWithZero();
876 if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
878 std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
879 " trying remap of field double on cells : \n" << descriptionField << std::endl;
886 field->setPartOfValuesAdv(fromArray,ccI);
890 //replaced by setPartOfValuesAdv if nbPtGauss==1
891 int iMax=ccI->getNbOfElems();
892 int* pccI=ccI->getPointer();
893 double* pField=field->getPointer();
894 double* pFrom=fromArray->getPointer();
895 int itarget, isource, delta=nbPtGauss*nbcomp;
896 for (int i=0; i<iMax; i=i+2) //cell
900 if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
901 throw INTERP_KERNEL::Exception("Error field override");
902 int ita=itarget*delta;
903 int iso=isource*delta;
904 for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
909 /*! constructing the MESH collection from a distributed file
911 * \param filename name of the master file containing the list of all the MED files
913 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
915 _owns_topology(true),
917 _domain_selector( 0 ),
918 _i_non_empty_mesh(-1),
919 _driver_type(MEDPARTITIONER::Undefined),
920 _subdomain_boundary_creates(false),
921 _family_splitting(false),
922 _create_empty_groups(false),
927 _driver=new MeshCollectionMedXmlDriver(this);
928 _driver->read (filename.c_str());
929 _driver_type = MedXml;
932 { // Handle all exceptions
933 if ( _driver ) delete _driver; _driver=0;
936 _driver=new MeshCollectionMedAsciiDriver(this);
937 _driver->read (filename.c_str());
938 _driver_type=MedAscii;
944 throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
947 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
948 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
949 _i_non_empty_mesh = idomain;
952 /*! Constructing the MESH collection from selected domains of a distributed file
954 * \param filename - name of the master file containing the list of all the MED files
955 * \param domainSelector - selector of domains to load
957 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
959 _owns_topology(true),
961 _domain_selector( &domainSelector ),
962 _i_non_empty_mesh(-1),
963 _driver_type(MEDPARTITIONER::Undefined),
964 _subdomain_boundary_creates(false),
965 _family_splitting(false),
966 _create_empty_groups(false),
969 std::string myfile=filename;
970 std::size_t found=myfile.find(".xml");
971 if (found!=std::string::npos) //file .xml
975 _driver=new MeshCollectionMedXmlDriver(this);
976 _driver->read ( filename.c_str(), _domain_selector );
977 _driver_type = MedXml;
980 { // Handle all exceptions
982 throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
987 found=myfile.find(".med");
988 if (found!=std::string::npos) //file .med single
990 //make a temporary file .xml and retry MedXmlDriver
992 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
994 <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
995 <description what=\"\" when=\"\"/>\n \
997 <mesh name=\"$meshName\"/>\n \
1000 <subdomain number=\"1\"/>\n \
1001 <global_numbering present=\"no\"/>\n \
1004 <subfile id=\"1\">\n \
1005 <name>$fileName</name>\n \
1006 <machine>localhost</machine>\n \
1010 <mesh name=\"$meshName\">\n \
1011 <chunk subdomain=\"1\">\n \
1012 <name>$meshName</name>\n \
1017 std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile.c_str());
1018 xml.replace(xml.find("$fileName"),9,myfile);
1019 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1020 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1021 xml.replace(xml.find("$meshName"),9,meshNames[0]);
1022 std::string nameFileXml(myfile);
1023 nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1024 std::string nameFileXmlDN,nameFileXmlBN;
1025 MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1026 nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1027 if (_domain_selector->rank()==0) //only on to write it
1029 std::ofstream f(nameFileXml.c_str());
1034 if (MyGlobals::_World_Size>1)
1035 MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1039 _driver=new MeshCollectionMedXmlDriver(this);
1040 _driver->read ( nameFileXml.c_str(), _domain_selector );
1041 _driver_type = MedXml;
1044 { // Handle all exceptions
1045 delete _driver; _driver=0;
1046 throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1053 _driver=new MeshCollectionMedAsciiDriver(this);
1054 _driver->read ( filename.c_str(), _domain_selector );
1055 _driver_type=MedAscii;
1061 throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1065 // find non-empty domain mesh
1066 for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1067 if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1068 _i_non_empty_mesh = idomain;
1072 //check for all proc/file compatibility of _field_descriptions
1074 _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1076 _field_descriptions=MyGlobals::_Field_Descriptions;
1079 catch(INTERP_KERNEL::Exception& e)
1081 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1082 throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1087 //check for all proc/file compatibility of _family_info
1088 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1089 _family_info=DevectorizeToMapOfStringInt(v2);
1091 catch(INTERP_KERNEL::Exception& e)
1093 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1094 throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1099 //check for all proc/file compatibility of _group_info
1100 std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1101 _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1103 catch(INTERP_KERNEL::Exception& e)
1105 std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1106 throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1111 /*! constructing the MESH collection from a sequential MED-file
1113 * \param filename MED file
1114 * \param meshname name of the mesh that is to be read
1116 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1118 _owns_topology(true),
1120 _domain_selector( 0 ),
1121 _i_non_empty_mesh(-1),
1123 _driver_type(MEDPARTITIONER::MedXml),
1124 _subdomain_boundary_creates(false),
1125 _family_splitting(false),
1126 _create_empty_groups(false),
1129 try // avoid memory leak in case of inexistent filename
1131 retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1137 throw INTERP_KERNEL::Exception("problem reading .med files");
1139 if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1140 _i_non_empty_mesh = 0;
1143 MEDPARTITIONER::MeshCollection::~MeshCollection()
1145 for (int i=0; i<(int)_mesh.size();i++)
1146 if (_mesh[i]!=0) _mesh[i]->decrRef();
1148 for (int i=0; i<(int)_cell_family_ids.size();i++)
1149 if (_cell_family_ids[i]!=0)
1150 _cell_family_ids[i]->decrRef();
1152 for (int i=0; i<(int)_face_mesh.size();i++)
1153 if (_face_mesh[i]!=0)
1154 _face_mesh[i]->decrRef();
1156 for (int i=0; i<(int)_face_family_ids.size();i++)
1157 if (_face_family_ids[i]!=0)
1158 _face_family_ids[i]->decrRef();
1160 for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1161 if ((*it).second!=0)
1162 (*it).second->decrRef();
1164 for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1165 if ((*it).second!=0)
1166 (*it).second->decrRef();
1169 if (_topology!=0 && _owns_topology)
1172 delete _joint_finder;
1176 /*! constructing the MESH collection from a file
1178 * The method creates as many MED-files as there are domains in the
1179 * collection. It also creates a master file that lists all the MED files.
1180 * The MED files created in ths manner contain joints that describe the
1181 * connectivity between subdomains.
1183 * \param filename name of the master file that will contain the list of the MED files
1186 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1188 //building the connect zones necessary for writing joints
1189 // if (_topology->nbDomain()>1)
1190 // buildConnectZones();
1191 //suppresses link with driver so that it can be changed for writing
1194 retrieveDriver()->write (filename.c_str(), _domain_selector);
1197 /*! creates or gets the link to the collection driver
1199 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1203 switch(_driver_type)
1206 _driver=new MeshCollectionMedXmlDriver(this);
1209 _driver=new MeshCollectionMedAsciiDriver(this);
1212 throw INTERP_KERNEL::Exception("Unrecognized driver");
1219 /*! gets an existing driver
1222 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1228 * retrieves the mesh dimension
1230 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1232 return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1235 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1238 for (int i=0; i<_mesh.size(); i++)
1245 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1248 for (int i=0; i<_mesh.size(); i++)
1250 if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1255 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1258 for (int i=0; i<_face_mesh.size(); i++)
1260 if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1265 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1270 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1275 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1277 return _mesh[idomain];
1280 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1282 return _face_mesh[idomain];
1285 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1287 return _connect_zones;
1290 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1295 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo)
1299 throw INTERP_KERNEL::Exception("topology is already set");
1305 /*! Method creating the cell graph
1307 * \param array returns the pointer to the structure that contains the graph
1308 * \param edgeweight returns the pointer to the table that contains the edgeweights
1309 * (only used if indivisible regions are required)
1311 void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
1313 using std::multimap;
1315 using std::make_pair;
1318 std::multimap< int, int > node2cell;
1319 std::map< pair<int,int>, int > cell2cellcounter;
1320 std::multimap<int,int> cell2cell;
1322 std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
1323 int nbdomain=_topology->nbDomain();
1325 if (isParallelMode())
1327 _joint_finder=new JointFinder(*this);
1328 _joint_finder->findCommonDistantNodes();
1329 commonDistantNodes=_joint_finder->getDistantNodeCell();
1332 if (MyGlobals::_Verbose>500)
1333 _joint_finder->print();
1336 if (MyGlobals::_Verbose>50)
1337 std::cout<<"getting nodal connectivity"<<std::endl;
1338 //looking for reverse nodal connectivity i global numbering
1339 for (int idomain=0; idomain<nbdomain; idomain++)
1341 if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1344 ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
1345 ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1346 int nbNodes=_mesh[idomain]->getNumberOfNodes();
1347 _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1348 //problem saturation over 1 000 000 nodes for 1 proc
1349 if (MyGlobals::_Verbose>100)
1350 std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1351 int* index_ptr=index->getPointer();
1352 int* revConnPtr=revConn->getPointer();
1353 for (int i=0; i<nbNodes; i++)
1355 for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
1357 int globalNode=_topology->convertNodeToGlobal(idomain,i);
1358 int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
1359 node2cell.insert(make_pair(globalNode, globalCell));
1365 for (int iother=0; iother<nbdomain; iother++)
1367 std::multimap<int,int>::iterator it;
1368 int isource=idomain;
1370 for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
1371 it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
1373 int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
1374 int globalCell=(*it).second;
1375 node2cell.insert(make_pair(globalNode, globalCell));
1381 //creating graph arcs (cell to cell relations)
1382 //arcs are stored in terms of (index,value) notation
1385 // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
1386 // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
1388 //warning here one node have less than or equal effective number of cell with it
1389 //but cell could have more than effective nodes
1390 //because other equals nodes in other domain (with other global inode)
1391 if (MyGlobals::_Verbose>50)
1392 std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
1394 for (int inode=0;inode<_topology->nbNodes();inode++)
1396 typedef multimap<int,int>::const_iterator MI;
1397 std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
1398 for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
1399 for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
1401 int icell1=cell1->second;
1402 int icell2=cell2->second;
1403 if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
1404 std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
1405 if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
1406 else (it->second)++;
1409 // for (int icell1=0; icell1<_topology->nbCells(); icell1++) //on all nodes
1411 // typedef multimap<int,int>::const_iterator MI;
1412 // std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
1413 // for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++) //on nodes with icell
1415 // std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
1416 // for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++) //on one of these cell
1418 // int icell2=cell2->second;
1419 // std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
1420 // if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
1421 // else (it->second)++;
1427 //converting the counter to a multimap structure
1428 for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
1429 it!=cell2cellcounter.end();
1433 cell2cell.insert(std::make_pair(it->first.first,it->first.second)); //should be adapted for 2D!
1434 cell2cell.insert(std::make_pair(it->first.second, it->first.first));
1438 if (MyGlobals::_Verbose>50)
1439 std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
1440 //filling up index and value to create skylinearray structure
1441 std::vector <int> index,value;
1445 for (int idomain=0; idomain<nbdomain; idomain++)
1447 if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
1448 int nbCells=_mesh[idomain]->getNumberOfCells();
1449 for (int icell=0; icell<nbCells; icell++)
1452 int globalCell=_topology->convertCellToGlobal(idomain,icell);
1453 multimap<int,int>::iterator it;
1454 pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
1455 ret=cell2cell.equal_range(globalCell);
1456 for (it=ret.first; it!=ret.second; ++it)
1458 int ival=(*it).second; //no adding one existing yet
1459 for (int i=idep ; i<idep+size ; i++)
1468 value.push_back(ival);
1472 idep=index[index.size()-1]+size;
1473 index.push_back(idep);
1477 array=new MEDPARTITIONER::SkyLineArray(index,value);
1479 if (MyGlobals::_Verbose>100)
1481 std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
1482 index.size()-1 << " " << value.size() << std::endl;
1483 int max=index.size()>15?15:index.size();
1486 for (int i=0; i<max; ++i)
1487 std::cout<<index[i]<<" ";
1488 std::cout << "... " << index[index.size()-1] << std::endl;
1489 for (int i=0; i<max; ++i)
1490 std::cout<< value[i] << " ";
1491 int ll=index[index.size()-1]-1;
1492 std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
1499 /*! Creates the partition corresponding to the cell graph and the partition number
1501 * \param nbdomain number of subdomains for the newly created graph
1503 * returns a topology based on the new graph
1505 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
1506 Graph::splitter_type split,
1507 const std::string& options_string,
1508 int *user_edge_weights,
1509 int *user_vertices_weights)
1511 if (MyGlobals::_Verbose>10)
1512 std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
1515 throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
1516 MEDPARTITIONER::SkyLineArray* array=0;
1518 buildCellGraph(array,edgeweights);
1524 #if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS)
1525 if (MyGlobals::_Verbose>10)
1526 std::cout << "METISGraph" << std::endl;
1527 cellGraph=new METISGraph(array,edgeweights);
1529 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
1533 #ifdef MED_ENABLE_SCOTCH
1534 if (MyGlobals::_Verbose>10)
1535 std::cout << "SCOTCHGraph" << std::endl;
1536 cellGraph=new SCOTCHGraph(array,edgeweights);
1538 throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
1543 //!user-defined weights
1544 if (user_edge_weights!=0)
1545 cellGraph->setEdgesWeights(user_edge_weights);
1546 if (user_vertices_weights!=0)
1547 cellGraph->setVerticesWeights(user_vertices_weights);
1549 if (MyGlobals::_Is0verbose>10)
1550 std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
1551 cellGraph->partGraph(nbdomain, options_string, _domain_selector);
1553 if (MyGlobals::_Is0verbose>10)
1554 std::cout << "building new topology" << std::endl;
1555 //cellGraph is a shared pointer
1556 Topology *topology=0;
1557 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
1559 delete [] edgeweights;
1561 if (MyGlobals::_Verbose>11)
1562 std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
1566 /*! Creates a topology for a partition specified by the user
1568 * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
1570 * returns a topology based on the new partition
1572 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
1574 MEDPARTITIONER::SkyLineArray* array=0;
1577 buildCellGraph(array,edgeweights);
1579 std::set<int> domains;
1580 for (int i=0; i<_topology->nbCells(); i++)
1582 domains.insert(partition[i]);
1584 cellGraph=new UserGraph(array, partition, _topology->nbCells());
1586 //cellGraph is a shared pointer
1587 Topology *topology=0;
1588 int nbdomain=domains.size();
1589 topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
1590 // if (array!=0) delete array;
1595 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
1597 for (int i=0; i<_topology->nbDomain(); i++)
1599 std::ostringstream oss;
1601 if (!isParallelMode() || _domain_selector->isMyDomain(i))
1602 _mesh[i]->setName(oss.str().c_str());
1606 ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
1607 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
1608 //something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name.c_str(),f1->getMesh()->getName(),0,f1->getName(),0,1);
1610 int rank=MyGlobals::_Rank;
1611 std::string tag="ioldFieldDouble="+IntToStr(iold);
1612 std::string descriptionIold=descriptionField+SerializeFromString(tag);
1613 if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
1615 if (MyGlobals::_Verbose>300)
1616 std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
1617 ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
1620 if (MyGlobals::_Verbose>200)
1621 std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
1622 std::string description, fileName, meshName, fieldName;
1623 int typeField, DT, IT, entity;
1624 fileName=MyGlobals::_File_Names[iold];
1625 if (MyGlobals::_Verbose>10)
1626 std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
1627 FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
1628 meshName=MyGlobals::_Mesh_Names[iold];
1630 ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
1631 fileName.c_str(), meshName.c_str(), 0, fieldName.c_str(), DT, IT);
1633 ParaMEDMEM::DataArrayDouble* res=f2->getArray();
1634 //to know names of components
1635 std::vector<std::string> browse=BrowseFieldDouble(f2);
1636 std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
1637 if (MyGlobals::_Verbose>10)
1638 std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
1639 MyGlobals::_General_Informations.push_back(localFieldInformation);
1642 _map_dataarray_double[descriptionIold]=res;
1646 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
1647 //to have unique valid fields names/pointers/descriptions for partitionning
1648 //filter _field_descriptions to be in all procs compliant and equal
1650 int nbfiles=MyGlobals::_File_Names.size(); //nb domains
1651 std::vector<std::string> r2;
1652 //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
1653 for (int i=0; i<(int)_field_descriptions.size(); i++)
1655 std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
1656 for (int ii=0; ii<(int)r1.size(); ii++)
1657 r2.push_back(r1[ii]);
1659 //here vector(procs*fields) of serialised vector(description) data
1660 _field_descriptions=r2;
1661 int nbfields=_field_descriptions.size(); //on all domains
1662 if ((nbfields%nbfiles)!=0)
1664 if (MyGlobals::_Rank==0)
1666 std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
1667 << "fileMedNames :" << std::endl
1668 << ReprVectorOfString(MyGlobals::_File_Names)
1669 << "field_descriptions :" << std::endl
1670 << ReprVectorOfString(MyGlobals::_Field_Descriptions);
1672 throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
1674 _field_descriptions.resize(nbfields/nbfiles);
1675 for (int i=0; i<(int)_field_descriptions.size(); i++)
1677 std::string str=_field_descriptions[i];
1678 str=EraseTagSerialized(str,"idomain=");
1679 str=EraseTagSerialized(str,"fileName=");
1680 _field_descriptions[i]=str;
1684 //returns true if inodes of a face are in inodes of a cell
1685 bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >& inodesCell)
1688 int nbok=inodesFace.size();
1689 for (int i=0; i<nbok; i++)
1691 int ii=inodesFace[i];
1693 std::cout << "isFaceOncell problem inodeface<0" << std::endl;
1694 for (int j=0; j<(int)inodesCell.size(); j++)
1696 if (ii==inodesCell[j])
1699 break; //inode of face found
1703 break; //inode of face not found do not continue...
1705 return (ires==nbok);
1708 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
1710 for (int inew=0; inew<_topology->nbDomain(); inew++)
1712 if (isParallelMode() && _domain_selector->isMyDomain(inew))
1714 if (MyGlobals::_Verbose>200)
1715 std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
1716 ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
1717 ParaMEDMEM::MEDCouplingUMesh* mfac=_face_mesh[inew];
1719 //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
1720 std::vector<int> nodeIds;
1721 getNodeIds(*mcel, *mfac, nodeIds);
1722 if (nodeIds.size()==0)
1723 continue; //one empty mesh nothing to do
1725 ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
1726 ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New();
1727 mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
1728 int *revC=revNodalCel->getPointer();
1729 int *revIndxC=revNodalIndxCel->getPointer();
1731 std::vector< int > faceOnCell;
1732 std::vector< int > faceNotOnCell;
1733 int nbface=mfac->getNumberOfCells();
1734 for (int iface=0; iface<nbface; iface++)
1737 std::vector< int > inodesFace;
1738 mfac->getNodeIdsOfCell(iface, inodesFace);
1739 int nbnodFace=inodesFace.size();
1740 //set inodesFace in mcel
1741 for (int i=0; i<nbnodFace; i++) inodesFace[i]=nodeIds[inodesFace[i]];
1742 int inod=inodesFace[0];
1744 std::cout << "filterFaceOnCell problem 1" << std::endl;
1745 int nbcell=revIndxC[inod+1]-revIndxC[inod];
1746 for (int j=0; j<nbcell; j++) //look for each cell with inod
1748 int icel=revC[revIndxC[inod]+j];
1749 std::vector< int > inodesCell;
1750 mcel->getNodeIdsOfCell(icel, inodesCell);
1751 ok=isFaceOncell(inodesFace, inodesCell);
1756 faceOnCell.push_back(iface);
1760 faceNotOnCell.push_back(iface);
1761 if (MyGlobals::_Is0verbose>300)
1762 std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
1766 revNodalCel->decrRef();
1767 revNodalIndxCel->decrRef();
1770 keyy=Cle1ToStr("filterFaceOnCell",inew);
1771 _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
1772 keyy=Cle1ToStr("filterNotFaceOnCell",inew);
1773 _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);