X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDPartitioner%2FMEDPARTITIONER_MeshCollection.cxx;h=d0ad9d4cbbce204e9d10d8d6f646b9096100eb74;hb=5f6d7861933f9c7596059d7a847e58624e177216;hp=0016d8595af7a4dcba9dcd90dfa1b9f295043779;hpb=10f37bf6f33a762626d7f1093b2f5450c1688667;p=tools%2Fmedcoupling.git diff --git a/src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx b/src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx index 0016d8595..d0ad9d4cb 100644 --- a/src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx +++ b/src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2012 CEA/DEN, EDF R&D +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -18,40 +18,41 @@ // #include "MEDPARTITIONER_MeshCollection.hxx" + +#include "MEDPARTITIONER_ConnectZone.hxx" +#include "MEDPARTITIONER_Graph.hxx" #include "MEDPARTITIONER_MeshCollectionDriver.hxx" -#include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx" #include "MEDPARTITIONER_MeshCollectionMedAsciiDriver.hxx" +#include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx" #include "MEDPARTITIONER_ParaDomainSelector.hxx" -#include "MEDPARTITIONER_Topology.hxx" #include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_Topology.hxx" +#include "MEDPARTITIONER_UserGraph.hxx" +#include "MEDPARTITIONER_Utils.hxx" -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI #include "MEDPARTITIONER_JointFinder.hxx" #endif -#include "MEDPARTITIONER_Graph.hxx" -#include "MEDPARTITIONER_UserGraph.hxx" -#include "MEDPARTITIONER_Utils.hxx" - -#include "MEDLoaderBase.hxx" -#include "MEDLoader.hxx" +#include "MEDCouplingAutoRefCountObjectPtr.hxx" +#include "MEDCouplingFieldDouble.hxx" #include "MEDCouplingMemArray.hxx" -#include "MEDCouplingUMesh.hxx" #include "MEDCouplingNormalizedUnstructuredMesh.hxx" -#include "MEDCouplingFieldDouble.hxx" -#include "PointLocator3DIntersectorP0P0.hxx" - -#include "MEDCouplingAutoRefCountObjectPtr.hxx" -#include "BBTree.txx" +#include "MEDCouplingSkyLineArray.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDLoader.hxx" +#include "MEDLoaderBase.hxx" -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI #include #endif -#if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS) +#ifdef MED_ENABLE_PARMETIS +#include "MEDPARTITIONER_ParMetisGraph.hxx" +#endif +#ifdef MED_ENABLE_METIS #include "MEDPARTITIONER_MetisGraph.hxx" #endif - #ifdef MED_ENABLE_SCOTCH #include "MEDPARTITIONER_ScotchGraph.hxx" #endif @@ -70,7 +71,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection() _domain_selector( 0 ), _i_non_empty_mesh(-1), _driver_type(MEDPARTITIONER::MedXml), - _subdomain_boundary_creates(false), + _subdomain_boundary_creates( MyGlobals::_Creates_Boundary_Faces ), _family_splitting(false), _create_empty_groups(false), _joint_finder(0) @@ -99,13 +100,15 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection _i_non_empty_mesh(-1), _name(initialCollection._name), _driver_type(MEDPARTITIONER::MedXml), - _subdomain_boundary_creates(false), + _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces), _family_splitting(family_splitting), _create_empty_groups(create_empty_groups), _joint_finder(0) { std::vector > > new2oldIds(initialCollection.getTopology()->nbDomain()); - castCellMeshes(initialCollection, new2oldIds); + std::vector o2nRenumber; + + castCellMeshes(initialCollection, new2oldIds, o2nRenumber ); //defining the name for the collection and the underlying meshes setName(initialCollection.getName()); @@ -114,7 +117,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection //treating faces ///////////////// -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1) MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages #endif @@ -130,7 +133,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection //////////////////// //treating families //////////////////// -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1) MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages #endif @@ -145,16 +148,16 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection std::cout<<"treating cell and face families"<getMesh(), - initialCollection.getCellFamilyIds(), - "cellFamily"); + this->getMesh(), + initialCollection.getCellFamilyIds(), + "cellFamily"); castIntField(initialCollection.getFaceMesh(), - this->getFaceMesh(), - initialCollection.getFaceFamilyIds(), - "faceFamily"); + this->getFaceMesh(), + initialCollection.getFaceFamilyIds(), + "faceFamily"); //treating groups -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1) MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages #endif @@ -162,8 +165,8 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection std::cout << "treating groups" << std::endl; _family_info=initialCollection.getFamilyInfo(); _group_info=initialCollection.getGroupInfo(); - -#ifdef HAVE_MPI2 + +#ifdef HAVE_MPI if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1) MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages #endif @@ -172,7 +175,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection castAllFields(initialCollection,"cellFieldDouble"); if (_i_non_empty_mesh<0) { - for (int i=0; i<_mesh.size(); i++) + for (size_t i=0; i<_mesh.size(); i++) { if (_mesh[i]) { @@ -182,16 +185,28 @@ MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection } } + // find faces common with neighbor domains and put them in groups + buildBoundaryFaces(); + + //building the connect zones necessary for writing joints + buildConnectZones( nodeMapping, o2nRenumber, initialCollection.getTopology()->nbDomain() ); + + // delete o2nRenumber + for ( size_t i = 0; i < o2nRenumber.size(); ++i ) + if ( o2nRenumber[i] ) + o2nRenumber[i]->decrRef(); } /*! - Creates the meshes using the topology underlying he mesh collection and the mesh data + Creates the meshes using the topology underlying he mesh collection and the mesh data coming from the ancient collection \param initialCollection collection from which the data is extracted to create the new meshes + \param [out] o2nRenumber returns for each new domain a permutation array returned by sortCellsInMEDFileFrmt() */ void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection, - std::vector > >& new2oldIds) + std::vector > >& new2oldIds, + std::vector & o2nRenumber) { if (MyGlobals::_Verbose>10) std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl; @@ -202,6 +217,7 @@ void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialColle int nbOldDomain=initialCollection.getTopology()->nbDomain(); _mesh.resize(nbNewDomain); + o2nRenumber.resize(nbNewDomain,0); int rank=MyGlobals::_Rank; //splitting the initial domains into smaller bits std::vector > splitMeshes; @@ -239,7 +255,7 @@ void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialColle } } } -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI if (isParallelMode()) { //if (MyGlobals::_Verbose>300) std::cout<<"proc "< meshes; - + for (int i=0; i<(int)splitMeshes[inew].size(); i++) - if (splitMeshes[inew][i]!=0) + if (splitMeshes[inew][i]!=0) if (splitMeshes[inew][i]->getNumberOfCells()>0) meshes.push_back(splitMeshes[inew][i]); if (!isParallelMode()||_domain_selector->isMyDomain(inew)) { - if (meshes.size()==0) + if (meshes.size()==0) { _mesh[inew]=CreateEmptyMEDCouplingUMesh(); std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl; @@ -281,6 +297,7 @@ void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialColle else { _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes); + o2nRenumber[inew]=_mesh[inew]->sortCellsInMEDFileFrmt(); bool areNodesMerged; int nbNodesMerged; if (meshes.size()>1) @@ -289,19 +306,18 @@ void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialColle array->decrRef(); // array is not used in this case } _mesh[inew]->zipCoords(); - } } for (int i=0;i<(int)splitMeshes[inew].size();i++) if (splitMeshes[inew][i]!=0) splitMeshes[inew][i]->decrRef(); - } + } if (MyGlobals::_Verbose>300) std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl; } /*! - \param initialCollection source mesh collection + \param initialCollection source mesh collection \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection */ void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping) @@ -313,26 +329,28 @@ void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialC { double* bbox; - BBTree<3>* tree; + BBTreeOfDim* tree = 0; + int dim = 3; if (!isParallelMode() || (_domain_selector->isMyDomain(iold))) { // std::map >, int > nodeClassifier; - int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes(); - bbox=new double[nvertices*2*3]; ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords(); double* coordsPtr=coords->getPointer(); + dim = coords->getNumberOfComponents(); + int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes(); + bbox=new double[nvertices*2*dim]; - for (int i=0; i(bbox,0,0,nvertices,1e-9); + tree=new BBTreeOfDim( dim, bbox,0,0,nvertices,1e-9); } for (int inew=0; inew<_topology->nbDomain(); inew++) { -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI //sending meshes for parallel computation if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold)) _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold)); @@ -343,7 +361,7 @@ void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialC ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords(); for (int inode=0; inodegetNumberOfNodes();inode++) { - double* coordsPtr=coords->getPointer()+inode*3; + double* coordsPtr=coords->getPointer()+inode*dim; vector elems; tree->getElementsAroundPoint(coordsPtr,elems); if (elems.size()==0) continue; @@ -353,52 +371,55 @@ void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialC } else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold))) #else - if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold))) + if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold))) #endif - { - ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords(); - for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++) - { - double* coordsPtr=coords->getPointer()+inode*3; - vector elems; - tree->getElementsAroundPoint(coordsPtr,elems); - if (elems.size()==0) continue; - nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode))); - } - } + { + ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords(); + for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++) + { + double* coordsPtr=coords->getPointer()+inode*dim; + vector elems; + tree->getElementsAroundPoint(coordsPtr,elems); + if (elems.size()==0) continue; + nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode))); + } + } } if (!isParallelMode() || (_domain_selector->isMyDomain(iold))) { delete tree; delete[] bbox; } - } + } } void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector& nodeIds) { using std::vector; + using MEDPARTITIONER::BBTreeOfDim; if (!&meshOne || !&meshTwo) return; //empty or not existing double* bbox; - BBTree<3>* tree; + BBTreeOfDim* tree = 0; int nv1=meshOne.getNumberOfNodes(); - bbox=new double[nv1*6]; ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords(); + int dim = coords->getNumberOfComponents(); + + bbox=new double[nv1*2*dim]; double* coordsPtr=coords->getPointer(); - for (int i=0; i(bbox,0,0,nv1,1e-9); - + tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9); + int nv2=meshTwo.getNumberOfNodes(); nodeIds.resize(nv2,-1); coords=meshTwo.getCoords(); for (int inode=0; inodegetPointer()+inode*3; + double* coordsPtr2=coords->getPointer()+inode*dim; vector elems; tree->getElementsAroundPoint(coordsPtr2,elems); if (elems.size()==0) continue; @@ -416,31 +437,30 @@ void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialColle const std::multimap, std::pair >& nodeMapping, std::vector > >& new2oldIds) { - //splitMeshes structure will contain the partition of //the old faces on the new ones //splitMeshes[4][2] contains the faces from old domain 2 //that have to be added to domain 4 - + using std::vector; using std::map; using std::multimap; using std::pair; using std::make_pair; - + if (MyGlobals::_Verbose>10) std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl; if (_topology==0) throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes"); - + int nbNewDomain=_topology->nbDomain(); int nbOldDomain=initialCollection.getTopology()->nbDomain(); - + vector& meshesCastFrom=initialCollection.getFaceMesh(); vector& meshesCastTo=this->getFaceMesh(); - + vector< vector > splitMeshes; - + splitMeshes.resize(nbNewDomain); for (int inew=0; inewisMyDomain(iold) && _domain_selector->isMyDomain(inew)) _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold)); - int nb=0; - if (splitMeshes[inew][iold]) - nb=splitMeshes[inew][iold]->getNumberOfCells(); + //int nb=0; + //if (splitMeshes[inew][iold]) + // nb=splitMeshes[inew][iold]->getNumberOfCells(); //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<getNumberOfCells()<decrRef(); @@ -579,10 +599,22 @@ void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialColle if (umesh->getNumberOfCells()>0) myMeshes.push_back(umesh); } - + + ParaMEDMEM::MEDCouplingUMesh *bndMesh = 0; + if ( _subdomain_boundary_creates && + _mesh[inew] && + _mesh[inew]->getNumberOfCells()>0 ) + { + bndMesh = + ((ParaMEDMEM::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true)); + if (bndMesh->getNumberOfCells()>0) + myMeshes.push_back( bndMesh ); + } + if (myMeshes.size()>0) { meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes); + meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef(); } else { @@ -592,6 +624,8 @@ void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialColle for (int iold=0; iolddecrRef(); + if ( bndMesh ) + bndMesh->decrRef(); } if (MyGlobals::_Verbose>300) std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl; @@ -600,9 +634,9 @@ void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialColle void MEDPARTITIONER::MeshCollection::castIntField(std::vector& meshesCastFrom, - std::vector& meshesCastTo, - std::vector& arrayFrom, - std::string nameArrayTo) + std::vector& meshesCastTo, + std::vector& arrayFrom, + std::string nameArrayTo) { using std::vector; @@ -613,18 +647,18 @@ void MEDPARTITIONER::MeshCollection::castIntField(std::vector99) std::cout<<"making accelerating structures"<* > acceleratingStructures(ioldMax); + std::vector acceleratingStructures(ioldMax); std::vectorbbox(ioldMax); for (int iold =0; iold< ioldMax; iold++) if (isParallelMode() && _domain_selector->isMyDomain(iold)) { ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner(); bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6); - acceleratingStructures[iold]=new BBTree<3,int> (bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples()); + acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples()); + sourceCoords->decrRef(); } - // send-recv operations -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI for (int inew=0; inew* myTree) + const BBTreeOfDim* myTree) { if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist @@ -697,21 +731,22 @@ void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold, const double* tc=targetCoords->getConstPointer(); int targetSize=targetMesh.getNumberOfCells(); int sourceSize=sourceMesh.getNumberOfCells(); -if (MyGlobals::_Verbose>200) -std::cout<<"remap vers target de taille "<200) + std::cout<<"remap vers target de taille "< ccI; std::string str,cle; str=nameArrayTo+"_toArray"; cle=Cle1ToStr(str,inew); int* toArray; - const BBTree<3>* tree; + const BBTreeOfDim* tree; bool cleantree=false; ParaMEDMEM::DataArrayDouble* sourceBBox=0; + int dim = targetCoords->getNumberOfComponents(); if (myTree==0) { sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8); - tree=new BBTree<3> (sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10); + tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10); cleantree=true; } else tree=myTree; @@ -730,17 +765,28 @@ std::cout<<"remap vers target de taille "<second->getPointer(); } - + + std::map< int, int > isource2nb; // count coincident elements + std::map::iterator i2nb; + for (int itargetnode=0; itargetnode intersectingElems; - tree->getElementsAroundPoint(tc+itargetnode*3,intersectingElems); // to be changed in 2D + tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems); if (intersectingElems.size()!=0) { int isourcenode=intersectingElems[0]; - toArray[itargetnode]=fromArray[isourcenode]; - ccI.push_back(itargetnode); - ccI.push_back(isourcenode); + if ( intersectingElems.size() > 1 ) + { + i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first; + isourcenode = intersectingElems[ i2nb->second++ ]; + } + if ( isourcenode < sourceSize ) // protection from invalid elements + { + toArray[itargetnode]=fromArray[isourcenode]; + ccI.push_back(itargetnode); + ccI.push_back(isourcenode); + } } } if (MyGlobals::_Verbose>200) @@ -751,12 +797,11 @@ std::cout<<"remap vers target de taille "<decrRef(); if (cleantree) delete tree; if (sourceBBox !=0) sourceBBox->decrRef(); - - } +} void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo) { @@ -776,7 +821,7 @@ void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollec std::string descriptionField=initialCollection.getFieldDescriptions()[ifield]; if (descriptionField.find(nameTo)==std::string::npos) continue; //only nameTo accepted in Fields name description -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI for (int inew=0; inew renumN2O = ids1->buildPermArrPerLevel(); + ids1->renumberInPlaceR( renumN2O->begin() ); + ids2->renumberInPlaceR( renumN2O->begin() ); + + if ( removeEqual ) + { + ids1 = ids1->buildUnique(); + ids2 = ids2->buildUnique(); + } + if ( delta != 0 ) + { + int * id = ids1->getPointer(); + for ( ; id < ids1->end(); ++id ) + ++(*id); + id = ids2->getPointer(); + for ( ; id < ids2->end(); ++id ) + ++(*id); + } + + // join + DataArrayInt* ids12 = DataArrayInt::Meld( ids1, ids2 ); // two components + ids12->rearrange( 1 ); // make one component + return ids12; + } + + //================================================================================ + /*! + * \brief Renumber ids according to mesh->sortCellsInMEDFileFrmt() + * \param [in,out] ids - cell ids to renumber + * \param [in] o2nRenumber - renumbering array in "Old to New" mode + */ + //================================================================================ + + void renumber( DataArrayInt* ids, const DataArrayInt* o2nRenumber ) + { + if ( !ids || !o2nRenumber ) + return; + int * id = ids->getPointer(); + const int * o2n = o2nRenumber->getConstPointer(); + for ( ; id < ids->end(); ++id ) + { + *id = o2n[ *id ]; + } + } +} + +//================================================================================ +/*! + * \brief Fill up ConnectZone's stored in _topology with nodal correspondences + * \param [in] nodeMapping - mapping between old nodes and new nodes + * (iolddomain,ioldnode)->(inewdomain,inewnode) + * \param [in] o2nRenumber - renumbering array returned by mesh->sortCellsInMEDFileFrmt() + * per a new domain + * \param [in] nbInitialDomains - nb of old domains + */ +//================================================================================ + +void MEDPARTITIONER::MeshCollection::buildConnectZones( const NodeMapping& nodeMapping, + const std::vector & o2nRenumber, + int nbInitialDomains) +{ + if ( !MyGlobals::_Create_Joints || _topology->nbDomain() < 2 ) + return; + + if ( MyGlobals::_World_Size > 1 ) + { + _topology->getCZ().clear(); + return; // not implemented for parallel mode + } + + // at construction, _topology creates cell correspondences basing on Graph information, + // and here we + // 1) add node correspondences, + // 2) split cell correspondences by cell geometry types + // 3) sort ids to be in ascending order + + const int nb_domains = _topology->nbDomain(); + + // ================================== + // 1) add node correspondences + // ================================== + + std::vector< std::vector< std::vector< int > > > nodeCorresp( nb_domains ); + for ( int idomain = 0; idomain < nb_domains; ++idomain ) + { + nodeCorresp[ idomain ].resize( nb_domains ); + } + + NodeMapping::const_iterator nmIt1, nmIt2 = nodeMapping.begin(); + for ( nmIt1 = nmIt2; nmIt1 != nodeMapping.end(); nmIt1 = nmIt2 ) + { + // look for an "old" node mapped into several "new" nodes in different domains + int nbSameOld = 0; + while ( ++nmIt2 != nodeMapping.end() && nmIt2->first == nmIt1->first ) + nbSameOld += ( nmIt2->second != nmIt1->second ); + + if ( nbSameOld > 0 ) + { + NodeMapping::const_iterator nmEnd = nmIt2; + for ( ; true; ++nmIt1 ) + { + nmIt2 = nmIt1; + if ( ++nmIt2 == nmEnd ) + break; + int dom1 = nmIt1->second.first; + int node1 = nmIt1->second.second; + for ( ; nmIt2 != nmEnd; ++nmIt2 ) + { + int dom2 = nmIt2->second.first; + int node2 = nmIt2->second.second; + if ( dom1 != dom2 ) + { + nodeCorresp[ dom1 ][ dom2 ].push_back( node1 ); + nodeCorresp[ dom1 ][ dom2 ].push_back( node2 ); + nodeCorresp[ dom2 ][ dom1 ].push_back( node2 ); + nodeCorresp[ dom2 ][ dom1 ].push_back( node1 ); + } + } + } + } + } + + // add nodeCorresp to czVec + + std::vector& czVec = _topology->getCZ(); + + for ( int idomain = 0; idomain < nb_domains; ++idomain ) + { + for ( int idomainNear = 0; idomainNear < nb_domains; ++idomainNear ) + { + std::vector< int > & corresp = nodeCorresp[ idomain ][ idomainNear ]; + if ( corresp.empty() ) + continue; + + MEDPARTITIONER::ConnectZone* cz = 0; + for ( size_t i = 0; i < czVec.size() && !cz; ++i ) + if ( czVec[i] && + czVec[i]->getLocalDomainNumber () == idomain && + czVec[i]->getDistantDomainNumber() == idomainNear ) + cz = czVec[i]; + + if ( !cz ) + { + cz = new MEDPARTITIONER::ConnectZone(); + cz->setName( "Nodal Connect Zone defined by MEDPARTITIONER" ); + cz->setLocalDomainNumber ( idomain ); + cz->setDistantDomainNumber( idomainNear ); + czVec.push_back(cz); + } + + cz->setNodeCorresp( &corresp[0], corresp.size()/2 ); + } + } + + // ========================================================== + // 2) split cell correspondences by cell geometry types + // ========================================================== + + for ( size_t i = 0; i < czVec.size(); ++i ) + { + MEDPARTITIONER::ConnectZone* cz = czVec[i]; + if ( !cz || + cz->getEntityCorrespNumber( 0,0 ) == 0 || + cz->getLocalDomainNumber () > (int)_mesh.size() || + cz->getDistantDomainNumber() > (int)_mesh.size() ) + continue; + ParaMEDMEM::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber () ]; + ParaMEDMEM::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ]; + + // separate ids of two domains + const ParaMEDMEM::MEDCouplingSkyLineArray *corrArray = cz->getEntityCorresp( 0, 0 ); + const DataArrayInt* ids12 = corrArray->getValueArray(); + MEDCouplingAutoRefCountObjectPtr ids1, ids2, ids12Sorted; + ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 ); + ids2 = ids12->selectByTupleId2( 1, corrArray->getLength(), 2 ); + + // renumber cells according to mesh->sortCellsInMEDFileFrmt() + renumber( ids1, o2nRenumber[ cz->getLocalDomainNumber() ]); + renumber( ids2, o2nRenumber[ cz->getDistantDomainNumber() ]); + + // check nb cell types + std::set types1, types2; + types1 = mesh1->getTypesOfPart( ids1->begin(), ids1->end() ); + types2 = mesh2->getTypesOfPart( ids2->begin(), ids2->end() ); + if ( types1.size() < 1 || types2.size() < 1 ) + continue; // parallel mode? + + MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1 + for ( size_t j = 0; j < czVec.size() && !cz21; ++j ) + if ( czVec[j] && + czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() && + czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() ) + cz21 = czVec[j]; + + if ( types1.size() == 1 && types2.size() == 1 ) // split not needed, only sort + { + ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 ); + cz->setEntityCorresp( *types1.begin(), *types2.begin(), + ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 ); + + if ( cz21 )// set 2->1 correspondence + { + ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 ); + cz21->setEntityCorresp( *types2.begin(), *types1.begin(), + ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 ); + } + } + else // split and sort + { + typedef std::pair< std::vector< int >, std::vector< int > > T2Vecs; + T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ]; + int t1, t2; + + const int nbIds = ids1->getNbOfElems(); + const int * p1 = ids1->begin(), * p2 = ids2->begin(); + for ( int i = 0; i < nbIds; ++i ) + { + t1 = mesh1->getTypeOfCell( p1[ i ]); + t2 = mesh2->getTypeOfCell( p2[ i ]); + T2Vecs & ids = idsByType[ t1 ][ t2 ]; + ids.first .push_back( p1[ i ]); + ids.second.push_back( p1[ i ]); + } + + const int maxType = int( INTERP_KERNEL::NORM_MAXTYPE ); + for ( t1 = 0; t1 < maxType; ++t1 ) + for ( t2 = 0; t2 < maxType; ++t2 ) + { + T2Vecs & ids = idsByType[ t1 ][ t2 ]; + if ( ids.first.empty() ) continue; + p1 = & ids.first[0]; + p2 = & ids.second[0]; + ids1->desallocate(); + ids1->pushBackValsSilent( p1, p1+ids.first.size() ); + ids2->desallocate(); + ids2->pushBackValsSilent( p2, p2+ids.first.size() ); + ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 ); + + cz->setEntityCorresp( t1, t2, + ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 ); + + if ( cz21 )// set 2->1 correspondence + { + ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 ); + cz21->setEntityCorresp( t2, t1, + ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 ); + break; + } + } + }// split and sort + + cz->setEntityCorresp( 0, 0, 0, 0 ); // erase ids computed by _topology + if ( cz21 ) + cz21->setEntityCorresp( 0, 0, 0, 0 ); + + } // loop on czVec + + + // ========================================== + // 3) sort node ids to be in ascending order + // ========================================== + + const bool removeEqual = ( nbInitialDomains > 1 ); + + for ( size_t i = 0; i < czVec.size(); ++i ) + { + MEDPARTITIONER::ConnectZone* cz = czVec[i]; + if ( !cz || cz->getNodeNumber() < 1 ) + continue; + if ( cz->getDistantDomainNumber() < cz->getLocalDomainNumber() ) + continue; // treat a pair of domains once + + MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1 + for ( size_t j = 0; j < czVec.size() && !cz21; ++j ) + if ( czVec[j] && + czVec[j]->getLocalDomainNumber () == cz->getDistantDomainNumber() && + czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() ) + cz21 = czVec[j]; + + // separate ids of two domains + const ParaMEDMEM::MEDCouplingSkyLineArray *corrArray = cz->getNodeCorresp(); + const DataArrayInt *ids12 = corrArray->getValueArray(); + MEDCouplingAutoRefCountObjectPtr ids1, ids2, ids12Sorted; + ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 ); + ids2 = ids12->selectByTupleId2( 1, corrArray->getLength(), 2 ); + + ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/0, removeEqual ); + cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 ); + + if ( cz21 )// set 2->1 correspondence + { + ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0, false ); + cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 ); + } + } +} + +//================================================================================ +/*! + * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces" + * group (where "n" and "p" are domain IDs) + */ +//================================================================================ + +void MEDPARTITIONER::MeshCollection::buildBoundaryFaces() +{ + if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates ) + return; + + if ( getMeshDimension() < 2 ) + return; + + using ParaMEDMEM::MEDCouplingUMesh; + using ParaMEDMEM::DataArrayDouble; + using ParaMEDMEM::DataArrayInt; + + std::vector& faceMeshes = getFaceMesh(); + int nbMeshes = faceMeshes.size(); + + //preparing bounding box trees for accelerating search of coincident faces + std::vector bbTrees(nbMeshes); + std::vectorbbox (nbMeshes); + for (int inew = 0; inew < nbMeshes-1; inew++) + if ( !isParallelMode() || _domain_selector->isMyDomain(inew) ) + { + DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner(); + bbox [inew] = bcCoords->computeBBoxPerTuple(1.e-6); + bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(), + bbox[inew]->getConstPointer(),0,0, + bbox[inew]->getNumberOfTuples()); + bcCoords->decrRef(); + } + + // loop on domains to find joint faces between them + for (int inew1 = 0; inew1 < nbMeshes; inew1++ ) + { + for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ ) + { + MEDCouplingUMesh* mesh1 = 0; + MEDCouplingUMesh* mesh2 = 0; + //MEDCouplingUMesh* recvMesh = 0; + bool mesh1Here = true, mesh2Here = true; + if (isParallelMode()) + { +#ifdef HAVE_MPI + mesh1Here = _domain_selector->isMyDomain(inew1); + mesh2Here = _domain_selector->isMyDomain(inew2); + if ( !mesh1Here && mesh2Here ) + { + //send mesh2 to domain of mesh1 + _domain_selector->sendMesh(*faceMeshes[inew2], + _domain_selector->getProcessorID(inew1)); + } + else if ( mesh1Here && !mesh2Here ) + { + //receiving mesh2 from a distant domain + _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2)); + if ( faceMeshes[ inew2 ] ) + faceMeshes[ inew2 ]->decrRef(); + faceMeshes[ inew2 ] = mesh2; + } +#endif + } + if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ]; + if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ]; + + // find coincident faces + std::vector< int > faces1, faces2; + if ( mesh1 && mesh2 ) + { + const DataArrayDouble* coords2 = mesh2->getBarycenterAndOwner(); + const double* c2 = coords2->getConstPointer(); + const int dim = coords2->getNumberOfComponents(); + const int nbFaces2 = mesh2->getNumberOfCells(); + const int nbFaces1 = mesh1->getNumberOfCells(); + + for (int i2 = 0; i2 < nbFaces2; i2++) + { + std::vector coincFaces; + bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces ); + if (coincFaces.size()!=0) + { + int i1 = coincFaces[0]; + // if ( coincFaces.size() > 1 ) + // { + // i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first; + // i1 = coincFaces[ i2nb->second++ ]; + // } + if ( i1 < nbFaces1 ) // protection from invalid elements + { + faces1.push_back( i1 ); + faces2.push_back( i2 ); + } + } + } + coords2->decrRef(); + } + + if ( isParallelMode()) + { +#ifdef HAVE_MPI + if ( mesh1Here && !mesh2Here ) + { + //send faces2 to domain of recvMesh + SendIntVec(faces2, _domain_selector->getProcessorID(inew2)); + } + else if ( !mesh1Here && mesh2Here ) + { + //receiving ids of faces from a domain of mesh1 + RecvIntVec(faces2, _domain_selector->getProcessorID(inew1)); + } +#endif + } + // if ( recvMesh ) + // recvMesh->decrRef(); + + // Create group "JOINT_inew1_inew2_Faces" and corresponding families + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + createJointGroup( is2nd ? faces2 : faces1, + inew1 , inew2, is2nd ); + } + + } // loop on the 2nd domains (inew2) + } // loop on the 1st domains (inew1) + + + // delete bounding box trees + for (int inew = 0; inew < nbMeshes-1; inew++) + if (isParallelMode() && _domain_selector->isMyDomain(inew)) + { + bbox[inew]->decrRef(); + delete bbTrees[inew]; + } +} + +//================================================================================ +/*! + * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families + * \param faces - face ids to include into the group + * \param inew1 - index of the 1st domain + * \param inew2 - index of the 2nd domain + * \param is2nd - in which (1st or 2nd) domain to create the group + */ +//================================================================================ + +void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces, + const int inew1, + const int inew2, + const bool is2nd ) +{ + // get the name of JOINT group + std::string groupName; + { + std::ostringstream oss; + oss << "JOINT_" + << (is2nd ? inew2 : inew1 ) << "_" + << (is2nd ? inew1 : inew2 ) << "_" + << ( getMeshDimension()==2 ? "Edge" : "Face" ); + groupName = oss.str(); + } + + // remove existing "JOINT_*" group + _group_info.erase( groupName ); + + // get family IDs array + int* famIDs = 0; + int inew = (is2nd ? inew2 : inew1 ); + int totalNbFaces = _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0; + std::string cle = Cle1ToStr( "faceFamily_toArray", inew ); + if ( !_map_dataarray_int.count(cle) ) + { + if ( totalNbFaces > 0 ) + { + ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New(); + p->alloc( totalNbFaces, 1 ); + p->fillWithZero(); + famIDs = p->getPointer(); + _map_dataarray_int[cle]=p; + } + } + else + { + famIDs = _map_dataarray_int.find(cle)->second->getPointer(); + } + // find a family ID of an existing JOINT group + int familyID = 0; + std::map::iterator name2id = _family_info.find( groupName ); + if ( name2id != _family_info.end() ) + familyID = name2id->second; + + // remove faces from the familyID-the family + if ( familyID != 0 && famIDs ) + for ( int i = 0; i < totalNbFaces; ++i ) + if ( famIDs[i] == familyID ) + famIDs[i] = 0; + + if ( faces.empty() ) + return; + + if ( familyID == 0 ) // generate a family ID for JOINT group + { + std::set< int > familyIDs; + for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id ) + familyIDs.insert( name2id->second ); + // find the next free family ID + int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd; + do + { + if ( !familyIDs.count( ++familyID )) + --freeIdCount; + } + while ( freeIdCount > 0 ); + } + + // push faces to familyID-th group + if ( faces.back() >= totalNbFaces ) + throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID"); + for ( size_t i = 0; i < faces.size(); ++i ) + famIDs[ faces[i] ] = familyID; + + // register JOINT group and family + _family_info[ groupName ] = familyID; // name of the group and family is same + _group_info [ groupName ].push_back( groupName ); +} + /*! constructing the MESH collection from a distributed file * * \param filename name of the master file containing the list of all the MED files @@ -917,7 +1508,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename) _domain_selector( 0 ), _i_non_empty_mesh(-1), _driver_type(MEDPARTITIONER::Undefined), - _subdomain_boundary_creates(false), + _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces), _family_splitting(false), _create_empty_groups(false), _joint_finder(0) @@ -961,7 +1552,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, Para _domain_selector( &domainSelector ), _i_non_empty_mesh(-1), _driver_type(MEDPARTITIONER::Undefined), - _subdomain_boundary_creates(false), + _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces), _family_splitting(false), _create_empty_groups(false), _joint_finder(0) @@ -1014,7 +1605,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, Para \n \ \n \ \n"; - std::vector meshNames=MEDLoader::GetMeshNames(myfile.c_str()); + std::vector meshNames=MEDLoader::GetMeshNames(myfile); xml.replace(xml.find("$fileName"),9,myfile); xml.replace(xml.find("$meshName"),9,meshNames[0]); xml.replace(xml.find("$meshName"),9,meshNames[0]); @@ -1030,7 +1621,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, Para f<1) MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml #endif @@ -1070,7 +1661,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, Para try { //check for all proc/file compatibility of _field_descriptions -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions); #else _field_descriptions=MyGlobals::_Field_Descriptions; @@ -1081,7 +1672,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, Para std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl; throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields"); } -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI try { //check for all proc/file compatibility of _family_info @@ -1121,7 +1712,7 @@ MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, cons _i_non_empty_mesh(-1), _name(meshname), _driver_type(MEDPARTITIONER::MedXml), - _subdomain_boundary_creates(false), + _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces), _family_splitting(false), _create_empty_groups(false), _joint_finder(0) @@ -1168,7 +1759,7 @@ MEDPARTITIONER::MeshCollection::~MeshCollection() delete _driver; if (_topology!=0 && _owns_topology) delete _topology; -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI delete _joint_finder; #endif } @@ -1185,9 +1776,6 @@ MEDPARTITIONER::MeshCollection::~MeshCollection() */ void MEDPARTITIONER::MeshCollection::write(const std::string& filename) { - //building the connect zones necessary for writing joints - // if (_topology->nbDomain()>1) - // buildConnectZones(); //suppresses link with driver so that it can be changed for writing delete _driver; _driver=0; @@ -1235,7 +1823,7 @@ int MEDPARTITIONER::MeshCollection::getMeshDimension() const int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const { int nb=0; - for (int i=0; i<_mesh.size(); i++) + for (size_t i=0; i<_mesh.size(); i++) { if (_mesh[i]) nb++; } @@ -1245,7 +1833,7 @@ int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const { int nb=0; - for (int i=0; i<_mesh.size(); i++) + for (size_t i=0; i<_mesh.size(); i++) { if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells(); } @@ -1255,7 +1843,7 @@ int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const { int nb=0; - for (int i=0; i<_face_mesh.size(); i++) + for (size_t i=0; i<_face_mesh.size(); i++) { if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells(); } @@ -1284,7 +1872,11 @@ ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int id std::vector& MEDPARTITIONER::MeshCollection::getCZ() { - return _connect_zones; + if ( _topology ) + return _topology->getCZ(); + + static std::vector noCZ; + return noCZ; } MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const @@ -1292,23 +1884,56 @@ MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const return _topology; } -void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo) +void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo, bool takeOwneship) { if (_topology!=0) { throw INTERP_KERNEL::Exception("topology is already set"); } else - _topology = topo; + { + _topology = topo; + _owns_topology = takeOwneship; + } } -/*! Method creating the cell graph +/*! Method creating the cell graph in serial mode + * + * \param array returns the pointer to the structure that contains the graph + * \param edgeweight returns the pointer to the table that contains the edgeweights + * (only used if indivisible regions are required) + */ +void MEDPARTITIONER::MeshCollection::buildCellGraph(ParaMEDMEM::MEDCouplingSkyLineArray* & array, int *& edgeweights ) +{ + + using std::map; + using std::vector; + using std::make_pair; + using std::pair; + + if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only"); + const ParaMEDMEM::MEDCouplingUMesh* mesh=_mesh[0]; + if (MyGlobals::_Verbose>50) + std::cout<<"getting nodal connectivity"<isMyDomain(0)) + { + vector value; + vector index(1,0); + + array=new ParaMEDMEM::MEDCouplingSkyLineArray(index,value); + return; + } + array=mesh->generateGraph(); +} +/*! Method creating the cell graph in multidomain mode * * \param array returns the pointer to the structure that contains the graph * \param edgeweight returns the pointer to the table that contains the edgeweights * (only used if indivisible regions are required) */ -void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights ) +void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(ParaMEDMEM::MEDCouplingSkyLineArray* & array, int *& edgeweights ) { using std::multimap; using std::vector; @@ -1321,7 +1946,7 @@ void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray std::vector > > commonDistantNodes; int nbdomain=_topology->nbDomain(); -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI if (isParallelMode()) { _joint_finder=new JointFinder(*this); @@ -1336,10 +1961,12 @@ void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray if (MyGlobals::_Verbose>50) std::cout<<"getting nodal connectivity"<isMyDomain(idomain)) continue; + meshDim = _mesh[idomain]->getMeshDimension(); ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New(); ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New(); @@ -1361,7 +1988,7 @@ void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray } revConn->decrRef(); index->decrRef(); -#ifdef HAVE_MPI2 +#ifdef HAVE_MPI for (int iother=0; iother::iterator it; @@ -1428,9 +2055,9 @@ void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray for (std::map,int>::const_iterator it=cell2cellcounter.begin(); it!=cell2cellcounter.end(); it++) - if (it->second>=3) + if (it->second>=meshDim) { - cell2cell.insert(std::make_pair(it->first.first,it->first.second)); //should be adapted for 2D! + cell2cell.insert(std::make_pair(it->first.first,it->first.second)); cell2cell.insert(std::make_pair(it->first.second, it->first.first)); } @@ -1473,8 +2100,8 @@ void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray index.push_back(idep); } } - - array=new MEDPARTITIONER::SkyLineArray(index,value); + + array=new ParaMEDMEM::MEDCouplingSkyLineArray(index,value); if (MyGlobals::_Verbose>100) { @@ -1503,32 +2130,48 @@ void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray * returns a topology based on the new graph */ MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain, - Graph::splitter_type split, + Graph::splitter_type split, const std::string& options_string, int *user_edge_weights, int *user_vertices_weights) { if (MyGlobals::_Verbose>10) std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl; - + if (nbdomain <1) throw INTERP_KERNEL::Exception("Number of subdomains must be > 0"); - MEDPARTITIONER::SkyLineArray* array=0; + ParaMEDMEM::MEDCouplingSkyLineArray* array=0; int* edgeweights=0; - buildCellGraph(array,edgeweights); - - Graph* cellGraph; + + if (_topology->nbDomain()>1 || isParallelMode()) + buildParallelCellGraph(array,edgeweights); + else + buildCellGraph(array,edgeweights); + + Graph* cellGraph = 0; switch (split) { case Graph::METIS: -#if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS) - if (MyGlobals::_Verbose>10) - std::cout << "METISGraph" << std::endl; - cellGraph=new METISGraph(array,edgeweights); -#else - throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please."); + if ( isParallelMode() && MyGlobals::_World_Size > 1 ) + { +#ifdef MED_ENABLE_PARMETIS + if (MyGlobals::_Verbose>10) + std::cout << "ParMETISGraph" << std::endl; + cellGraph=new ParMETISGraph(array,edgeweights); +#endif + } + if ( !cellGraph ) + { +#ifdef MED_ENABLE_METIS + if (MyGlobals::_Verbose>10) + std::cout << "METISGraph" << std::endl; + cellGraph=new METISGraph(array,edgeweights); #endif + } + if ( !cellGraph ) + throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please."); break; + case Graph::SCOTCH: #ifdef MED_ENABLE_SCOTCH if (MyGlobals::_Verbose>10) @@ -1571,10 +2214,14 @@ MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nb */ MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition) { - MEDPARTITIONER::SkyLineArray* array=0; + ParaMEDMEM::MEDCouplingSkyLineArray* array=0; int* edgeweights=0; - buildCellGraph(array,edgeweights); + if ( _topology->nbDomain()>1) + buildParallelCellGraph(array,edgeweights); + else + buildCellGraph(array,edgeweights); + Graph* cellGraph; std::set domains; for (int i=0; i<_topology->nbCells(); i++) @@ -1599,13 +2246,13 @@ void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name) std::ostringstream oss; oss<isMyDomain(i)) - _mesh[i]->setName(oss.str().c_str()); + _mesh[i]->setName(oss.str()); } } ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold) //getField look for and read it if not done, and assume decrRef() in ~MeshCollection; -//something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name.c_str(),f1->getMesh()->getName(),0,f1->getName(),0,1); +//something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1); { int rank=MyGlobals::_Rank; std::string tag="ioldFieldDouble="+IntToStr(iold); @@ -1628,7 +2275,7 @@ ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::strin meshName=MyGlobals::_Mesh_Names[iold]; ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField, - fileName.c_str(), meshName.c_str(), 0, fieldName.c_str(), DT, IT); + fileName, meshName, 0, fieldName, DT, IT); ParaMEDMEM::DataArrayDouble* res=f2->getArray(); //to know names of components @@ -1648,6 +2295,10 @@ void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions() //filter _field_descriptions to be in all procs compliant and equal { int nbfiles=MyGlobals::_File_Names.size(); //nb domains + if (nbfiles==0) + { + nbfiles=_topology->nbDomain(); + } std::vector r2; //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data for (int i=0; i<(int)_field_descriptions.size(); i++) @@ -1709,7 +2360,7 @@ void MEDPARTITIONER::MeshCollection::filterFaceOnCell() { for (int inew=0; inew<_topology->nbDomain(); inew++) { - if (isParallelMode() && _domain_selector->isMyDomain(inew)) + if (!isParallelMode() || _domain_selector->isMyDomain(inew)) { if (MyGlobals::_Verbose>200) std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl; @@ -1721,13 +2372,13 @@ void MEDPARTITIONER::MeshCollection::filterFaceOnCell() getNodeIds(*mcel, *mfac, nodeIds); if (nodeIds.size()==0) continue; //one empty mesh nothing to do - + ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New(); ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New(); mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel); int *revC=revNodalCel->getPointer(); int *revIndxC=revNodalIndxCel->getPointer(); - + std::vector< int > faceOnCell; std::vector< int > faceNotOnCell; int nbface=mfac->getNumberOfCells(); @@ -1737,11 +2388,20 @@ void MEDPARTITIONER::MeshCollection::filterFaceOnCell() std::vector< int > inodesFace; mfac->getNodeIdsOfCell(iface, inodesFace); int nbnodFace=inodesFace.size(); + if ( nbnodFace != mfac->getNumberOfNodesInCell( iface )) + continue; // invalid node ids //set inodesFace in mcel - for (int i=0; i= 0 ); + if ( nbok != nbnodFace ) + continue; int inod=inodesFace[0]; if (inod<0) - std::cout << "filterFaceOnCell problem 1" << std::endl; + { + std::cout << "filterFaceOnCell problem 1" << std::endl; + continue; + } int nbcell=revIndxC[inod+1]-revIndxC[inod]; for (int j=0; jdecrRef(); revNodalIndxCel->decrRef(); - - std::string keyy; - keyy=Cle1ToStr("filterFaceOnCell",inew); - _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell); - keyy=Cle1ToStr("filterNotFaceOnCell",inew); - _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell); + + // std::string keyy; + // keyy=Cle1ToStr("filterFaceOnCell",inew); + // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell); + // keyy=Cle1ToStr("filterNotFaceOnCell",inew); + // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell); + + // filter the face mesh + if ( faceOnCell.empty() ) + _face_mesh[inew] = CreateEmptyMEDCouplingUMesh(); + else + _face_mesh[inew] = (ParaMEDMEM::MEDCouplingUMesh *) + mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true); + mfac->decrRef(); + + // filter the face families + std::string key = Cle1ToStr("faceFamily_toArray",inew); + if ( getMapDataArrayInt().count( key )) + { + ParaMEDMEM::DataArrayInt * & fam = getMapDataArrayInt()[ key ]; + ParaMEDMEM::DataArrayInt * famFilter = ParaMEDMEM::DataArrayInt::New(); + famFilter->alloc(faceOnCell.size(),1); + int* pfamFilter = famFilter->getPointer(); + int* pfam = fam->getPointer(); + for ( size_t i=0; idecrRef(); + fam = famFilter; + } } } }