From: ageay Date: Wed, 29 Feb 2012 11:22:30 +0000 (+0000) Subject: Merge from CVW_ReadyToMergeToV6main from branch BR_MEDPartitioner_dev. X-Git-Tag: V6_main_FINAL~813 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=24eb4e649f41ff8dbf32d1a3d5780306cee2fbcc;p=tools%2Fmedcoupling.git Merge from CVW_ReadyToMergeToV6main from branch BR_MEDPartitioner_dev. --- diff --git a/src/MEDPartitioner/MEDPARTITIONER.hxx b/src/MEDPartitioner/MEDPARTITIONER.hxx new file mode 100755 index 000000000..e18a0da3a --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER.hxx @@ -0,0 +1,38 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// File : MEDPARTITIONER.hxx +// Author : Alexander A. BORODIN +// Module : MED +// Exporting/Importing defines for Windows Platform +// +#ifndef MEDPARTITIONER_HXX_ +#define MEDPARTITIONER_HXX_ + +#ifdef WNT +# if defined MEDPARTITIONER_EXPORTS || defined medsplitter_EXPORTS +# define MEDPARTITIONER_EXPORT __declspec( dllexport ) +# else +# define MEDPARTITIONER_EXPORT __declspec( dllimport ) +# endif +#else +# define MEDPARTITIONER_EXPORT +#endif + +#endif //MEDPARTITIONER_HXX_ diff --git a/src/MEDPartitioner/MEDPARTITIONER_ConnectZone.cxx b/src/MEDPartitioner/MEDPARTITIONER_ConnectZone.cxx new file mode 100644 index 000000000..040410714 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_ConnectZone.cxx @@ -0,0 +1,286 @@ +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +// few Med Memory include files +#include "MEDCouplingUMesh.hxx" +#include "MEDPARTITIONER_SkyLineArray.hxx" +#include "MEDPARTITIONER_ConnectZone.hxx" + +// few STL include files +#include + +using namespace MEDPARTITIONER; + +CONNECTZONE::CONNECTZONE(): + _name("") + ,_description("") + ,_distantDomainNumber(0) + ,_localDomainNumber(0) + ,_nodeCorresp(0) + ,_faceCorresp(0) +{ + _entityCorresp.clear(); +} + +CONNECTZONE::~CONNECTZONE(){ + if (_nodeCorresp !=0) delete _nodeCorresp; + if (_faceCorresp !=0) delete _faceCorresp; + for (std::map < std::pair ,MEDPARTITIONER::MEDSKYLINEARRAY * >::iterator + iter = _entityCorresp.begin(); iter != _entityCorresp.end(); iter++) + { + delete iter->second; + } +} + +CONNECTZONE::CONNECTZONE(const CONNECTZONE & myConnectZone): + _name(myConnectZone._name) + ,_description(myConnectZone._description) + ,_distantDomainNumber(myConnectZone._distantDomainNumber) + ,_localDomainNumber(myConnectZone._localDomainNumber) + ,_nodeCorresp(myConnectZone._nodeCorresp) + ,_faceCorresp(myConnectZone._faceCorresp) + ,_entityCorresp(myConnectZone._entityCorresp) +{ +} +std::string CONNECTZONE::getName() const +{ + return _name; +} +std::string CONNECTZONE::getDescription() const +{ + return _description; +} +int CONNECTZONE::getDistantDomainNumber() const +{ + return _distantDomainNumber; +} +int CONNECTZONE::getLocalDomainNumber() const +{ + return _localDomainNumber; +} + +ParaMEDMEM::MEDCouplingUMesh* CONNECTZONE::getLocalMesh() const +{ + return _localMesh; +} + +ParaMEDMEM::MEDCouplingUMesh * CONNECTZONE::getDistantMesh() const +{ + return _distantMesh; +} + +bool CONNECTZONE::isEntityCorrespPresent(int localEntity, + int distantEntity) const +{ + typedef std::map, MEDPARTITIONER::MEDSKYLINEARRAY*>::const_iterator map_iter; + + for (map_iter iter=_entityCorresp.begin(); iter != _entityCorresp.end(); iter++) + { + if ((iter->first).first==localEntity && (iter->first).second==distantEntity) + return true; + } + return false; +} + +const int * CONNECTZONE::getNodeCorrespIndex() const +{ + return _nodeCorresp->getIndex(); +} + +const int * CONNECTZONE::getNodeCorrespValue() const +{ + return _nodeCorresp->getValue(); +} +int CONNECTZONE::getNodeNumber() const +{ + return _nodeCorresp->getNumberOf(); +} +const int * CONNECTZONE::getFaceCorrespIndex() const +{ + return _faceCorresp->getIndex(); +} + +const int * CONNECTZONE::getFaceCorrespValue() const +{ + return _faceCorresp->getValue(); +} +int CONNECTZONE::getFaceNumber() const +{ + return _faceCorresp->getNumberOf(); +} +const int * CONNECTZONE::getEntityCorrespIndex(int localEntity, + int distantEntity) const +{ + typedef std::map, MEDPARTITIONER::MEDSKYLINEARRAY*>::const_iterator map_iter; + + for (map_iter iter=_entityCorresp.begin(); iter != _entityCorresp.end(); iter++) + { + if ((iter->first).first==localEntity && (iter->first).second==distantEntity) + return iter->second->getIndex(); + } + return 0; +} + +const int * CONNECTZONE::getEntityCorrespValue(int localEntity, + int distantEntity) const +{ + typedef std::map, MEDPARTITIONER::MEDSKYLINEARRAY*>::const_iterator map_iter; + + for (map_iter iter=_entityCorresp.begin(); iter != _entityCorresp.end(); iter++) + { + if ((iter->first).first==localEntity && (iter->first).second==distantEntity) + return iter->second->getValue(); + } + return 0; +} + +int CONNECTZONE::getEntityCorrespNumber(int localEntity, + int distantEntity) const +{ + typedef std::map, MEDPARTITIONER::MEDSKYLINEARRAY*>::const_iterator map_iter; + + for (map_iter iter=_entityCorresp.begin(); iter != _entityCorresp.end(); iter++) + { + if ((iter->first).first==localEntity && (iter->first).second==distantEntity) + return iter->second->getNumberOf(); + } + return 0; +} + + +int CONNECTZONE::getEntityCorrespLength(int localEntity, + int distantEntity) const +{ + typedef std::map, MEDPARTITIONER::MEDSKYLINEARRAY*>::const_iterator map_iter; + + for (map_iter iter=_entityCorresp.begin(); iter != _entityCorresp.end(); iter++) + { + if ((iter->first).first==localEntity && (iter->first).second==distantEntity) + return iter->second->getLength(); + } + return 0; +} + +void CONNECTZONE::setName(std::string name) +{ + _name=name; +} +void CONNECTZONE::setDescription(std::string description) +{ + _description=description; +} +void CONNECTZONE::setDistantDomainNumber(int distantDomainNumber) +{ + _distantDomainNumber=distantDomainNumber; +} +void CONNECTZONE::setLocalDomainNumber(int localDomainNumber) +{ + _localDomainNumber=localDomainNumber; +} +void CONNECTZONE::setLocalMesh(ParaMEDMEM::MEDCouplingUMesh * localMesh) +{ + _localMesh=localMesh; +} + +void CONNECTZONE::setDistantMesh(ParaMEDMEM::MEDCouplingUMesh * distantMesh) +{ + _distantMesh=distantMesh; +} + +/*! transforms an int array containing + * the node-node connections + * to a MEDSKYLINEARRAY + */ +void CONNECTZONE::setNodeCorresp(int * nodeCorresp, int nbnode) +{ + std::vector index(nbnode+1),value(2*nbnode); + for (int i=0; i index(nbface+1),value(2*nbface); + for (int i=0; i index(nbentity+1),value(2*nbentity); + for (int i=0; i +#include + +namespace MEDPARTITIONER { +class MEDPARTITIONER_EXPORT CONNECTZONE +{ +private : + std::string _name; + std::string _description; + int _localDomainNumber; + int _distantDomainNumber; + + ParaMEDMEM::MEDCouplingUMesh * _localMesh; + ParaMEDMEM::MEDCouplingUMesh * _distantMesh; + + MEDPARTITIONER::MEDSKYLINEARRAY * _nodeCorresp; + MEDPARTITIONER::MEDSKYLINEARRAY * _faceCorresp; + + std::map < std::pair , MEDPARTITIONER::MEDSKYLINEARRAY * > _entityCorresp; + +public : + CONNECTZONE(); + ~CONNECTZONE(); + CONNECTZONE(const CONNECTZONE & myConnectZone); + + std::string getName() const ; + std::string getDescription() const ; + int getDistantDomainNumber() const ; + int getLocalDomainNumber() const ; + ParaMEDMEM::MEDCouplingUMesh * getLocalMesh() const ; + ParaMEDMEM::MEDCouplingUMesh * getDistantMesh() const ; + + bool isEntityCorrespPresent(int localEntity,int distantEntity) const; + const int * getNodeCorrespIndex() const; + const int * getNodeCorrespValue() const; + int getNodeNumber() const; + const int * getFaceCorrespIndex() const; + const int * getFaceCorrespValue() const; + int getFaceNumber() const; + const int * getEntityCorrespIndex(int localEntity, + int distantEntity) const; + const int * getEntityCorrespValue(int localEntity, + int distantEntity) const; + int getEntityCorrespNumber(int localEntity, + int distantEntity) const; + int getEntityCorrespLength(int localEntity, + int distantEntity) const; + void setName(std::string name) ; + void setDescription(std::string description) ; + void setDistantDomainNumber(int distantDomainNumber) ; + void setLocalDomainNumber(int distantDomainNumber) ; + void setLocalMesh(ParaMEDMEM::MEDCouplingUMesh * localMesh) ; + void setDistantMesh(ParaMEDMEM::MEDCouplingUMesh * distantMesh) ; + + void setNodeCorresp(int * nodeCorresp, int nbnode); + void setNodeCorresp(MEDPARTITIONER::MEDSKYLINEARRAY* array); + void setFaceCorresp(int * faceCorresp, int nbface); + void setFaceCorresp(MEDPARTITIONER::MEDSKYLINEARRAY* array); + void setEntityCorresp(int localEntity, + int distantEntity, + int * entityCorresp, int nbentity); + void setEntityCorresp(int localEntity, + int distantEntity, + MEDPARTITIONER::MEDSKYLINEARRAY* array); +}; +} +# endif diff --git a/src/MEDPartitioner/MEDPARTITIONER_Graph.cxx b/src/MEDPartitioner/MEDPARTITIONER_Graph.cxx new file mode 100644 index 000000000..66fd5497d --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_Graph.cxx @@ -0,0 +1,47 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#include "MEDPARTITIONER_Graph.hxx" + +using namespace MEDPARTITIONER; + +Graph::Graph(MEDPARTITIONER::MEDSKYLINEARRAY* array, int* edgeweight):_graph(array),_partition(0),_edgeweight(edgeweight),_cellweight(0) +{ +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +Graph::~Graph() +{ + if (_partition) + { + delete _partition; + _partition=0; + } + if (_graph) + { + delete _graph; + _graph=0; + } +} + diff --git a/src/MEDPartitioner/MEDPARTITIONER_Graph.hxx b/src/MEDPartitioner/MEDPARTITIONER_Graph.hxx new file mode 100644 index 000000000..a1967a176 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_Graph.hxx @@ -0,0 +1,67 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef MEDPARTITIONER_GRAPH_HXX_ +#define MEDPARTITIONER_GRAPH_HXX_ + +#include "MEDPARTITIONER.hxx" +#include "MEDPARTITIONER_SkyLineArray.hxx" +#include +namespace MEDPARTITIONER { + +class ParaDomainSelector; +class MEDPARTITIONER_EXPORT Graph +{ + public: + typedef enum {METIS,SCOTCH} splitter_type; + + Graph(){}; + + //creates a graph from a SKYLINEARRAY + Graph(MEDPARTITIONER::MEDSKYLINEARRAY* graph, int* edgeweight=0); + + virtual ~Graph(); + + void setEdgesWeights(int* edgeweight){_edgeweight=edgeweight;} + void setVerticesWeights(int* cellweight){_cellweight=cellweight;} + + //computes partitioning of the graph + virtual void partGraph(int ndomain, const std::string&, ParaDomainSelector* sel=0)=0; + + //! returns the partitioning + const int* getPart() const {return _partition->getValue();} + + //! returns the number of graph vertices (which can correspond to the cells in the mesh!) + int nbVertices() const {return _graph->getNumberOf();} + + const MEDPARTITIONER::MEDSKYLINEARRAY* getGraph() const {return _graph;} + + protected: + + MEDPARTITIONER::MEDSKYLINEARRAY* _graph; + + MEDPARTITIONER::MEDSKYLINEARRAY* _partition; + + int* _edgeweight; + + int* _cellweight; + }; + +} +#endif /*GRAPH_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_JointFinder.cxx b/src/MEDPartitioner/MEDPARTITIONER_JointFinder.cxx new file mode 100644 index 000000000..348764b84 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_JointFinder.cxx @@ -0,0 +1,194 @@ +#include "MEDPARTITIONER_JointFinder.hxx" +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_Topology.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDPARTITIONER_utils.hxx" +#include "BBTree.txx" + + +/*! Method contributing to the distant cell graph + */ +using namespace MEDPARTITIONER; + +JointFinder::JointFinder(const MESHCollection& mc):_mesh_collection(mc), _topology(mc.getTopology()),_domain_selector(mc.getParaDomainSelector()) +{ +} + +JointFinder::~JointFinder() +{ + //if (MyGlobals::_is0verbose>100) cout<<"TODO ~JointFinder"<nbDomain(); + _distant_node_cell.resize(nbdomain); + _node_node.resize(nbdomain); + for (int i=0; inbProcs(); + std::vector* > bbtree(nbdomain,(BBTree<3>*) 0); + std::vector rev(nbdomain,(DataArrayInt*) 0); + std::vector revIndx(nbdomain,(DataArrayInt*) 0); + int meshDim; + int spaceDim; + + //init rev and revIndx and bbtree for my domain (of me:proc n) + for (int mydomain=0; mydomainisMyDomain(mydomain)) continue; + const ParaMEDMEM::MEDCouplingUMesh* myMesh=_mesh_collection.getMesh(mydomain); + meshDim=myMesh->getMeshDimension(); + spaceDim= myMesh->getSpaceDimension(); + rev[mydomain] = ParaMEDMEM::DataArrayInt::New(); + revIndx[mydomain] = ParaMEDMEM::DataArrayInt::New(); + myMesh->getReverseNodalConnectivity(rev[mydomain],revIndx[mydomain]); + double* bbx=new double[2*spaceDim*myMesh->getNumberOfNodes()]; + for (int i=0; igetNumberOfNodes()*spaceDim; i++) + { + const double* coords=myMesh->getCoords()->getConstPointer(); + bbx[2*i]=(coords[i])-1e-12; + bbx[2*i+1]=bbx[2*i]+2e-12; + } + bbtree[mydomain]=new BBTree<3> (bbx,0,0,myMesh->getNumberOfNodes(),-1e-12); + delete[] bbx; + } + + //send my domains to other proc an receive other domains from other proc + for (int isource=0; isourceisMyDomain(isource)&&_domain_selector->isMyDomain(itarget)) continue; + if (_domain_selector->isMyDomain(isource)) + { + //preparing data for treatment on target proc + int targetProc = _domain_selector->getProcessorID(itarget); + + std::vector vec(spaceDim*sourceMesh->getNumberOfNodes()); + //cvw cout<<"\nproc "<<_domain_selector->rank()<<" : numberOfNodes "<getNumberOfNodes()<getCoords()->getConstPointer(),sourceMesh->getCoords()->getConstPointer()+sourceMesh->getNumberOfNodes()*spaceDim,&vec[0]); + sendDoubleVec(vec,targetProc); + + //retrieving target data for storage in commonDistantNodes array + std::vector localCorrespondency; + recvIntVec(localCorrespondency, targetProc); + //cvw cout<<"\nproc "<<_domain_selector->rank()<<" : nodeCellCorrespondency "; + for (int i=0; iisMyDomain(itarget)) + { + //receiving data from source proc + int sourceProc = isource%nbproc; + std::vector recvVec; + recvDoubleVec(recvVec,sourceProc); + std::map commonNodes; // (local nodes, distant nodes) list + //cvw cout<<"\nproc "<<_domain_selector->rank()<<" : commonNodes "; + for (int inode=0; inode<(recvVec.size()/meshDim); inode++) + { + double* bbox=new double[2*spaceDim]; + for (int i=0; i inodes; + bbtree[itarget]->getIntersectingElems(bbox,inodes); + delete[] bbox; + + if (inodes.size()>0) + { + commonNodes.insert(std::make_pair(inodes[0],inode)); + //cvw cout<<" "< nodeCellCorrespondency; + for (std::map::iterator iter=commonNodes.begin(); iter!=commonNodes.end(); iter++) + { + _node_node[itarget][isource].push_back(std::make_pair(iter->first, iter->second));//storing node pairs in a vector + const int* revIndxPtr=revIndx[itarget]->getConstPointer(); + const int* revPtr=rev[itarget]->getConstPointer(); + for (int icell=revIndxPtr[iter->first]; icellfirst+1]; icell++) + { + nodeCellCorrespondency.push_back(iter->second); // + int globalCell=_topology->convertCellToGlobal(itarget,revPtr[icell]); + nodeCellCorrespondency.push_back(globalCell); + //nodeCellCorrespondency.push_back(revPtr[icell]); //need to set at global numerotation + //cout<<"processor "<second<<" cellLoc "<rank()<<" : JointFinder sendIntVec "<<_domain_selector->rank()<decrRef(); + if (revIndx[i]!=0) revIndx[i]->decrRef(); + if (bbtree[i]!=0) delete bbtree[i]; + } + + if (MyGlobals::_verbose>100) + std::cout<<"proc "<<_domain_selector->rank()<<" : end JointFinder::findCommonDistantNodes"< > > & JointFinder::getDistantNodeCell() +{ + return _distant_node_cell; +} + +std::vector > > >& JointFinder::getNodeNode() +{ + return _node_node; +} + +void JointFinder::print() +//it is for debug on small arrays under mpi 2,3 cpus +{ + int nbdomain=_topology->nbDomain(); + //MPI_Barrier(MPI_COMM_WORLD); + if (MyGlobals::_is0verbose>0) + cout<<"\nJointFinder print node-node (nn)iproc|itarget|isource|i|inodefirst-inodesecond\n\n"<< + "JointFinder print distantNode=cell (nc)iproc|itarget|isource|inode=icell\n\n"; + for (int isource=0; isourcerank()<rank()<<" : JointFinder _distant_node_cell itarget/isource/inode=icell"<::iterator it; + for (it=_distant_node_cell[isource][itarget].begin() ; it!=_distant_node_cell[isource][itarget].end(); it++) + { + cout<<" nc"<<_domain_selector->rank()<<"|"< +#include +namespace MEDPARTITIONER { + class MESHCollection; + class ParaDomainSelector; + class Topology; +class JointFinder +{ +public: + JointFinder(const MESHCollection& mc); + ~JointFinder(); + void findCommonDistantNodes(); + void print(); + std::vector > >& getDistantNodeCell(); + std::vector > > >& getNodeNode(); + std::vector > > _distant_node_cell; +private: + const MESHCollection& _mesh_collection; + const ParaDomainSelector* _domain_selector; + const Topology* _topology; + //std::vector > > _distant_node_cell; + std::vector > > > _node_node; +}; +}; +#endif diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx new file mode 100644 index 000000000..eaab9f521 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.cxx @@ -0,0 +1,1971 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingNormalizedUnstructuredMesh.hxx" +#include "MEDCouplingMemArray.hxx" +#include "PointLocator3DIntersectorP0P0.hxx" +#include "BBTree.txx" +#include "MEDPARTITIONER_utils.hxx" + +#include "MEDPARTITIONER_Graph.hxx" + +#include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" + +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_MESHCollectionDriver.hxx" +#include "MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx" +#include "MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx" +#include "MEDPARTITIONER_JointFinder.hxx" +#include "MEDPARTITIONER_UserGraph.hxx" +#include "MEDLoader.hxx" +#include "MEDCouplingFieldDouble.hxx" + +#ifdef HAVE_MPI2 +#include +#endif + +#ifdef ENABLE_METIS +#include "MEDPARTITIONER_METISGraph.hxx" +#endif +#ifdef ENABLE_SCOTCH +#include "MEDPARTITIONER_SCOTCHGraph.hxx" +#endif + +#include +#include +#include + +#include + +#include +#include + +using namespace MEDPARTITIONER; + +#ifndef WNT +using namespace __gnu_cxx; +#else +using namespace std; +#endif + + +MESHCollection::MESHCollection() + : _topology(0), + _owns_topology(false), + _driver(0), + _domain_selector( 0 ), + _i_non_empty_mesh(-1), + _driver_type(MEDPARTITIONER::MedXML), + _subdomain_boundary_creates(false), + _family_splitting(false), + _create_empty_groups(false), + _joint_finder(0) +{ +} + +/*!constructor creating a new mesh collection (mesh series + topology) + *from an old collection and a new topology + * + * On output, the constructor has built the meshes corresponding to the new mesh collection. + * The new topology has been updated so that face and node mappings are included. + * The families have been cast to their projections in the new topology. + * + * \param initial_collection collection from which the data (coordinates, connectivity) are taken + * \param topology topology containing the cell mappings + */ + +MESHCollection::MESHCollection( + MESHCollection& initialCollection, + Topology* topology, + bool family_splitting, + bool create_empty_groups) //cvwat04 + : _name(initialCollection._name), + _topology(topology), + _owns_topology(false), + _driver(0), + _domain_selector( initialCollection._domain_selector ), + _i_non_empty_mesh(-1), + _driver_type(MEDPARTITIONER::MedXML), + _subdomain_boundary_creates(false), + _family_splitting(family_splitting), + _create_empty_groups(create_empty_groups), + _joint_finder(0) +{ + std::vector > > new2oldIds(initialCollection.getTopology()->nbDomain()); + if (MyGlobals::_verbose>10) std::cout<<"proc "<(inewdomain,inewnode) + createNodeMapping(initialCollection, nodeMapping); + //cvw std::cout<<"castMeshes"< > > new2oldFaceIds; + castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds); + + //////////////////// + //treating families + //////////////////// + + if (MyGlobals::_is0verbose) + if (isParallelMode()) std::cout<<"ParallelMode on "<nbDomain()<<" Domains"<nbDomain()<<" Domains"<10) std::cout<<"treating cell and face families"<getMesh(), + initialCollection.getCellFamilyIds(), + "cellFamily"); + castIntField2(initialCollection.getFaceMesh(), + this->getFaceMesh(), + initialCollection.getFaceFamilyIds(), + "faceFamily"); + + ////////////////// + //treating groups + ////////////////// + if (MyGlobals::_is0verbose) std::cout<<"treating groups"< > >& new2oldIds) +{ + if (_topology==0) throw INTERP_KERNEL::Exception(LOCALIZED("Topology has not been defined on call to castCellMeshes")); + + int nbNewDomain=_topology->nbDomain(); + int nbOldDomain=initialCollection.getTopology()->nbDomain(); + + _mesh.resize(nbNewDomain); + int rank=MyGlobals::_rank; + //if (MyGlobals::_verbose>10) std::cout<<"proc "< > splitMeshes; + splitMeshes.resize(nbNewDomain); + for (int inew=0; inewisMyDomain(iold)) + { + int size=(initialCollection._mesh)[iold]->getNumberOfCells(); + std::vector globalids(size); + initialCollection.getTopology()->getCellList(iold, &globalids[0]); + std::vector ilocalnew(size); //local + std::vector ipnew(size); //idomain old + //cvw work locally + _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]); + + new2oldIds[iold].resize(nbNewDomain); + for (int i=0; ibuildPartOfMySelf( + &new2oldIds[iold][inew][0], + &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(), + true); + if (MyGlobals::_verbose>400) + std::cout<<"proc "<getNumberOfCells()<300) std::cout<<"proc "<isMyDomain(iold) && _domain_selector->isMyDomain(inew)) continue; + + if(initialCollection._domain_selector->isMyDomain(iold)) + _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew)); + + if (_domain_selector->isMyDomain(inew)) + _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold)); + + } + } + + //fusing the split meshes + if (MyGlobals::_verbose>200) std::cout<<"proc "< meshes; + + for (int i=0; i< splitMeshes[inew].size();i++) + 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) + { + _mesh[inew]=createEmptyMEDCouplingUMesh(); + //throw INTERP_KERNEL::Exception(LOCALIZED("castCellMeshes fusing : no meshes")); + cout<<"WARNING : castCellMeshes fusing : no meshes try another number of processors"<mergeNodes(1e-12,areNodesMerged,nbNodesMerged); + array->decrRef(); // array is not used in this case + _mesh[inew]->zipCoords(); + } + } + for (int i=0; i< splitMeshes[inew].size(); i++) + if (splitMeshes[inew][i]!=0) splitMeshes[inew][i]->decrRef(); + } + if (MyGlobals::_verbose>300) std::cout<<"proc "<nbDomain();iold++) + { + + double* bbox; + BBTree<3>* tree; + if (!isParallelMode() || (_domain_selector->isMyDomain(iold))) + { + // std::map >, int > nodeClassifier; + int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes(); + bbox=new double[nvertices*6]; + ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords(); + double* coordsPtr=coords->getPointer(); + + for (int i=0; i(bbox,0,0,nvertices,1e-9); + } + + for (int inew=0; inew<_topology->nbDomain(); inew++) //cvwat12 + { + //sending meshes for parallel computation + if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold)) + _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold)); + else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold)) + { + ParaMEDMEM::MEDCouplingUMesh* mesh; + _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew)); + ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords(); + for (int inode=0; inodegetNumberOfNodes();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))); + } + mesh->decrRef(); + } + else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold))) + { + 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))); + } + } + } + if (!isParallelMode() || (_domain_selector->isMyDomain(iold))) + { + delete tree; + delete[] bbox; + } + } + +} + +//getNodeIds(meshCell, meshFace, nodeIds) +//inodeCell=nodeIds[inodeFace] +//(put the biggest mesh in One) +//if no corresponding node then inodeCell==-1 +void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, vector& nodeIds) +{ + using std::vector; + if (!&meshOne || !&meshTwo) return; //empty or not existing + double* bbox; + BBTree<3>* tree; + int nv1=meshOne.getNumberOfNodes(); + bbox=new double[nv1*6]; + ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords(); + double* coordsPtr=coords->getPointer(); + for (int i=0; i(bbox,0,0,nv1,1e-9); + + int nv2=meshTwo.getNumberOfNodes(); + nodeIds.resize(nv2,-1); + coords=meshTwo.getCoords(); + for (int inode=0; inodegetPointer()+inode*3; + vector elems; + tree->getElementsAroundPoint(coordsPtr,elems); + if (elems.size()==0) continue; + nodeIds[inode]=elems[0]; + } + delete tree; + delete[] bbox; +} + +/*! + creates the face meshes on the new domains from the faces on the old domain and the node mapping + faces at the interface are duplicated +*/ +void MESHCollection::castFaceMeshes(MESHCollection& initialCollection, + 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; + + vector& meshesCastFrom=initialCollection.getFaceMesh(); + vector& meshesCastTo=this->getFaceMesh(); + + vector< vector > splitMeshes; + int newSize=_topology->nbDomain(); + int fromSize=meshesCastFrom.size(); + + splitMeshes.resize(newSize); + for (int inew=0;inewisMyDomain(iold)) continue; + if (meshesCastFrom[iold] != 0) + { + for (int ielem=0; ielemgetNumberOfCells(); ielem++) + { + vector nodes; + meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes); + + map faces; + + //analysis of element ielem + //counters are set for the element + //for each source node, the mapping is interrogated and the domain counters + //are incremented for each target node + //the face is considered as going to target domains if the counter of the domain + //is equal to the number of nodes + for (int inode=0; inode,pair >::const_iterator MI; + int mynode=nodes[inode]; + + pair myRange = nodeMapping.equal_range(make_pair(iold,mynode)); + for (MI iter=myRange.first; iter!=myRange.second; iter++) + { + int inew=iter->second.first; + if (faces.find(inew)==faces.end()) + faces[inew]=1; + else + faces[inew]++; + } + } + + for (map::iterator iter=faces.begin(); iter!=faces.end(); iter++) + { + if (iter->second==nodes.size()) + //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]? + //it is not sure here... + //done before writeMedfile on option?... see filterFaceOnCell() + new2oldIds[iold][iter->first].push_back(ielem); + } + } + + //creating the splitMeshes from the face ids + for (int inew=0;inew<_topology->nbDomain();inew++) + { + if (meshesCastFrom[iold]->getNumberOfCells() > 0) //cvw + { + splitMeshes[inew][iold]= + (ParaMEDMEM::MEDCouplingUMesh*) + ( meshesCastFrom[iold]->buildPartOfMySelf( + &new2oldIds[iold][inew][0], + &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(), + true) + ); + splitMeshes[inew][iold]->zipCoords(); + } + else + { + //std::cout<<"one empty mesh from "<isMyDomain(iold)<<" "<< + _domain_selector->isMyDomain(inew)<isMyDomain(iold) && !_domain_selector->isMyDomain(inew)) + if (splitMeshes[inew][iold] != 0) { + //cvw std::cout<<"send NOT empty mesh "<getName()<<" "<sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew)); + } + else { + //std::cout<<"send empty mesh "<sendMesh(*(empty), _domain_selector->getProcessorID(inew)); + } + if (!_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)) + _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold)); + } + empty->decrRef(); + } + + //recollecting the bits of splitMeshes to fuse them into one + if (MyGlobals::_verbose>300) std::cout<<"proc "< myMeshes; + for (int iold=0; ioldgetNumberOfCells()>0) { + myMeshes.push_back(umesh); + } + //else { + // std::cout<<"one empty mesh "<0) + { + meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes); + } + else + { + //std::cout<<"one empty meshes to merge "<setName("emptyMesh"); + //empty->setMeshDimension(3); + //empty->allocateCells(0); + ParaMEDMEM::MEDCouplingUMesh *empty=createEmptyMEDCouplingUMesh(); + meshesCastTo[inew]=empty; + } + // meshesCastTo[inew]->zipCoords(); + for (int iold=0; iolddecrRef(); + } + //if (MyGlobals::_verbose>1) std::cout<<"proc "<setCoords(sourceCoords); + vector c; + vector cI; + tmpMesh->getNodeIdsNearPoints(targetCoords->getConstPointer(),targetMesh.getNumberOfCells(),1e-10,c,cI); + if (cI.size()!= targetMesh.getNumberOfCells()+1) + throw INTERP_KERNEL::Exception(LOCALIZED("Error in source/target projection")); + for (int itargetnode=0; itargetnodedecrRef(); + targetCoords->decrRef(); + tmpMesh->decrRef(); +} + +void MESHCollection::castIntField2( //cvwat08 + std::vector& meshesCastFrom, + std::vector& meshesCastTo, + std::vector& arrayFrom, + std::string nameArrayTo) +{ + using std::vector; + + int ioldMax=meshesCastFrom.size(); + int inewMax=meshesCastTo.size(); + // send-recv operations + for (int inew=0; inewisMyDomain(iold) && !_domain_selector->isMyDomain(inew)) + { + //send mesh + _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew)); + //send vector + int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1! + vectorsendIds; + if (MyGlobals::_verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField sendIntVec size "<0) //no empty + { + sendIds.resize(size); + std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]); + } + else //empty + { + size=0; + sendIds.resize(size); + } + //std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField sendIntVec "<getProcessorID(inew)); + } + //receiving arrays from distant domains + if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)) + { + //receive mesh + vector recvIds; + ParaMEDMEM::MEDCouplingUMesh* recvMesh; + _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold)); + //receive vector + if (MyGlobals::_verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<getProcessorID(iold)); + remapIntField2(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo); + recvMesh->decrRef(); //cww is it correct? + } + } + } + + //local contributions and aggregation + for (int inew=0; inewisMyDomain(iold) && _domain_selector->isMyDomain(inew))) + { + remapIntField2(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo); + } + } +} + +void MESHCollection::remapIntField2(int inew,int iold, + const ParaMEDMEM::MEDCouplingUMesh& sourceMesh, + const ParaMEDMEM::MEDCouplingUMesh& targetMesh, + const int* fromArray, + string nameArrayTo) +//here we store ccI for next use in future call of castAllFields and remapDoubleField2 +{ + using std::vector; + //cout<<"remapIntField2 "<setCoords(sourceCoords); + vector c; + vector cI; + vector ccI; //memorize intersection target<-source(inew,iold) + string str,cle; + str=nameArrayTo+"_toArray"; + cle=cle1ToStr(str,inew); + int* toArray; + int targetSize=targetMesh.getNumberOfCells(); + //first time iold : create and initiate + if (_mapDataArrayInt.find(cle)==_mapDataArrayInt.end()) + { + if (MyGlobals::_is0verbose>100) cout<<"create "<alloc(targetSize,1); + p->fillWithZero(); + toArray=p->getPointer(); + _mapDataArrayInt[cle]=p; + } + else //other times iold: refind and complete + { + toArray=_mapDataArrayInt.find(cle)->second->getPointer(); + } + tmpMesh->getNodeIdsNearPoints(targetCoords->getConstPointer(),targetSize,1e-10,c,cI); + if (cI.size()!=targetSize+1) + throw INTERP_KERNEL::Exception(LOCALIZED("Error in source/target projection")); + for (int itargetnode=0; itargetnode700) cout<<"proc "<decrRef(); + targetCoords->decrRef(); + tmpMesh->decrRef(); +} + +void MESHCollection::castAllFields(MESHCollection& initialCollection, string nameArrayTo) //cvwat08 +{ + using std::vector; + if (nameArrayTo!="cellFieldDouble") + throw INTERP_KERNEL::Exception(LOCALIZED("Error castAllField only on cellFieldDouble")); + + string nameTo="typeData=6"; //resume the type of field casted + // send-recv operations + int ioldMax=initialCollection.getMesh().size(); + int inewMax=this->getMesh().size(); + int iFieldMax=initialCollection.getFieldDescriptions().size(); + if (MyGlobals::_verbose>10) cout<<"castAllFields with:\n"< initialCollection.getFieldDescriptions() + //getFieldDescriptions() is a string description of field coherent and tested and set BEFORE. + //see collection.prepareFieldDescriptions() + for (int ifield=0; ifieldisMyDomain(iold) && !_domain_selector->isMyDomain(inew)) + { + int target=_domain_selector->getProcessorID(inew); + ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold); //cvwat14 + //getField look for and read it if not done, and assume decrRef() in ~MESHCollection; + if (MyGlobals::_verbose>10) + std::cout<<"proc "<<_domain_selector->rank()<<" : castAllFields sendDouble"<isMyDomain(iold) && _domain_selector->isMyDomain(inew)) + { + int source=_domain_selector->getProcessorID(iold); + //receive vector + if (MyGlobals::_verbose>10) + std::cout<<"proc "<<_domain_selector->rank()<<" : castAllFields recvDouble"<isMyDomain(iold) && _domain_selector->isMyDomain(inew))) + { + ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold); //cvwat14 + remapDoubleField3(inew,iold,field,nameArrayTo,descriptionField); + } + } + } +} + +void MESHCollection::remapDoubleField3(int inew, int iold, + ParaMEDMEM::DataArrayDouble* fromArray, + string nameArrayTo, + string descriptionField) +//here we use 'cellFamily_ccI inew iold' set in remapIntField2 +{ + using std::vector; + if (nameArrayTo!="cellFieldDouble") + throw INTERP_KERNEL::Exception(LOCALIZED("Error remapDoubleField3 only on cellFieldDouble")); + string cle=cle2ToStr("cellFamily_ccI",inew,iold); + + map::iterator it1; + it1=_mapDataArrayInt.find(cle); + if (it1==_mapDataArrayInt.end()) + { + cerr<<"proc "<second; + if (MyGlobals::_verbose>300) cout<<"proc "<getNbOfElems()<getMesh()[inew]->getNumberOfCells(); //number of cell of mesh + int nbcomp=fromArray->getNumberOfComponents(); + int nbPtGauss=strToInt(extractFromDescription(descriptionField, "nbPtGauss=")); + + //int nbk=fromArray->getNumberOfTuples(); + + //cle=reprGenericDescription(descriptionField)+" "+intToStr(inew); + string tag="inewFieldDouble="+intToStr(inew); + cle=descriptionField+serializeFromString(tag); + //cout<<"descriptionField in remapDoubleField3 : "<getNbOfElems(); + int fromArrayNbOfComp=fromArray->getNumberOfComponents(); + int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss; + + if (MyGlobals::_verbose>1000) + { + cout<<"proc "<getNumberOfTuples()<< + " nbcells "<getNumberOfComponents()<::iterator it2; + it2=_mapDataArrayDouble.find(cle); + if (it2==_mapDataArrayDouble.end()) + { + if (MyGlobals::_verbose>300) cout<<"proc "<alloc(nbcell*nbPtGauss,nbcomp); + field->fillWithZero(); + } + else + { + field=it2->second; + if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp) + { + cerr<<"proc "<getNbOfElems()<setPartOfValuesAdv(fromArray,ccI); + } + else + { + //replaced by setPartOfValuesAdv if nbPtGauss==1 + int iMax=ccI->getNbOfElems(); + int* pccI=ccI->getPointer(); + double* pField=field->getPointer(); + double* pFrom=fromArray->getPointer(); + int itarget, isource, delta=nbPtGauss*nbcomp; + for (int i=0; i=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell)) + throw INTERP_KERNEL::Exception(LOCALIZED("Error field override")); + int ita=itarget*delta; + int iso=isource*delta; + for (int k=0; ksetCoords(sourceCoords); + vector c; + vector cI; + tmpMesh->getNodeIdsNearPoints(targetCoords->getConstPointer(),targetMesh.getNumberOfCells(),1e-10,c,cI); + if (cI.size()!= targetMesh.getNumberOfCells()+1) + throw INTERP_KERNEL::Exception(LOCALIZED("Error in source/target projection")); + + for (int itargetnode=0; itargetnodedecrRef(); + targetCoords->decrRef(); + tmpMesh->decrRef(); +} +*/ + +/*! constructing the MESH collection from a distributed file + * + * \param filename name of the master file containing the list of all the MED files + */ +MESHCollection::MESHCollection(const std::string& filename) + : _topology(0), + _owns_topology(true), + _driver(0), + _domain_selector( 0 ), + _i_non_empty_mesh(-1), + _driver_type(MEDPARTITIONER::Undefined), + _subdomain_boundary_creates(false), + _family_splitting(false), + _create_empty_groups(false), + _joint_finder(0) +{ + // char filenamechar[256]; + // strcpy(filenamechar,filename.c_str()); + try + { + _driver=new MESHCollectionMedXMLDriver(this); + _driver->read (filename.c_str()); + _driver_type = MedXML; + } + catch(...) + { // Handle all exceptions + if ( _driver ) delete _driver; _driver=0; + try + { + _driver=new MESHCollectionMedAsciiDriver(this); + _driver->read (filename.c_str()); + _driver_type=MedAscii; + } + catch(...) + { + if ( _driver ) delete _driver; _driver=0; + throw INTERP_KERNEL::Exception(LOCALIZED("file does not comply with any recognized format")); + } + } + for ( int idomain = 0; idomain < _mesh.size(); ++idomain ) + if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 ) + _i_non_empty_mesh = idomain; +} + +/*! Constructing the MESH collection from selected domains of a distributed file + * + * \param filename - name of the master file containing the list of all the MED files + * \param domainSelector - selector of domains to load + */ +MESHCollection::MESHCollection(const std::string& filename, ParaDomainSelector& domainSelector) //cvwat01 + : _topology(0), + _owns_topology(true), + _driver(0), + //cvw _domain_selector( domainSelector.nbProcs() > 1 ? & domainSelector : 0 ), + _domain_selector( &domainSelector ), + _i_non_empty_mesh(-1), + _driver_type(MEDPARTITIONER::Undefined), + _subdomain_boundary_creates(false), + _family_splitting(false), + _create_empty_groups(false), + _joint_finder(0) +{ + std::string myfile=filename; + std::size_t found=myfile.find(".xml"); + if (found!=std::string::npos) //file .xml + { + try + { + _driver=new MESHCollectionMedXMLDriver(this); //cvwat02 + _driver->read ( (char*)filename.c_str(), _domain_selector ); + _driver_type = MedXML; + } + catch(...) + { // Handle all exceptions + if ( _driver ) delete _driver; _driver=0; + throw INTERP_KERNEL::Exception(LOCALIZED("file .xml does not comply with any recognized format")); + } + } + else + { + found=myfile.find(".med"); + if (found!=std::string::npos) //file .med single + { + //make a temporary file .xml and retry MedXMLDriver + std::string xml="\ +\n \ +\n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + $fileName\n \ + localhost\n \ + \n \ + \n \ + \n \ + \n \ + \n \ + $meshName\n \ + \n \ + \n \ + \n \ +\n"; + std::vector meshNames=MEDLoader::GetMeshNames(myfile.c_str()); + xml.replace(xml.find("$fileName"),9,myfile); + xml.replace(xml.find("$meshName"),9,meshNames[0]); + xml.replace(xml.find("$meshName"),9,meshNames[0]); + xml.replace(xml.find("$meshName"),9,meshNames[0]); + //std::cout<rank()==0) //only on to write it + { + std::ofstream f(nameFileXml.c_str()); + f<read ( (char*)nameFileXml.c_str(), _domain_selector ); + _driver_type = MedXML; + } + catch(...) + { // Handle all exceptions + if ( _driver ) delete _driver; _driver=0; + throw INTERP_KERNEL::Exception(LOCALIZED("file medpartitioner_xxx.xml does not comply with any recognized format")); + } + } + else //no extension + { + try + { + _driver=new MESHCollectionMedAsciiDriver(this); + _driver->read ( (char*)filename.c_str(), _domain_selector ); + _driver_type=MedAscii; + } + catch(...) + { + if ( _driver ) delete _driver; _driver=0; + throw INTERP_KERNEL::Exception(LOCALIZED("file name with no extension does not comply with any recognized format")); + } + } + } + + /*done in MESHCollectionMedXMLDriver read + if ( isParallelMode() ) + // to know nb of cells on each proc to compute global cell ids from locally global + _domain_selector->gatherNbOf( getMesh() );*/ + + // find non-empty domain mesh + for ( int idomain = 0; idomain < _mesh.size(); ++idomain ) + if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 ) + _i_non_empty_mesh = idomain; + + try + { + //check for all proc/file compatibility of _fieldDescriptions + //*MyGlobals::_fileNames=allgathervVectorOfString(*MyGlobals::_fileNames); + _fieldDescriptions=allgathervVectorOfString(MyGlobals::_fieldDescriptions); //cvwat07 + } + catch(INTERP_KERNEL::Exception& e) + { + cerr<<"proc "< v2=allgathervVectorOfString(vectorizeFromMapOfStringInt(_familyInfo)); + _familyInfo=devectorizeToMapOfStringInt(v2); + } + catch(INTERP_KERNEL::Exception& e) + { + cerr<<"proc "< v2=allgathervVectorOfString( + vectorizeFromMapOfStringVectorOfString(_groupInfo)); + _groupInfo=deleteDuplicatesInMapOfStringVectorOfString( + devectorizeToMapOfStringVectorOfString(v2)); + } + catch(INTERP_KERNEL::Exception& e) + { + cerr<<"proc "< _meshes=MEDLoader::GetMeshNames(filename); + //std::vector< std::string > _fields=MEDLoader::GetAllFieldNamesOnMesh(filename,meshname[0]); + //cout<<"number of fields "<<_fields.size()<readSeq (filename.c_str(),meshname.c_str()); + } + catch (...) + { + if ( _driver ) delete _driver; _driver=0; + throw INTERP_KERNEL::Exception(LOCALIZED("problem reading .med files")); + } + if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 ) + _i_non_empty_mesh = 0; +} + +MESHCollection::~MESHCollection() +{ + for (int i=0; i<_mesh.size();i++) + if (_mesh[i]!=0) _mesh[i]->decrRef(); + + for (int i=0; i<_cellFamilyIds.size();i++) + if (_cellFamilyIds[i]!=0) _cellFamilyIds[i]->decrRef(); + + for (int i=0; i<_faceMesh.size();i++) + if (_faceMesh[i]!=0) _faceMesh[i]->decrRef(); + + for (int i=0; i<_faceFamilyIds.size();i++) + if (_faceFamilyIds[i]!=0) _faceFamilyIds[i]->decrRef(); + + for (map::iterator it=_mapDataArrayInt.begin() ; it!=_mapDataArrayInt.end(); it++ ) + if ((*it).second!=0) (*it).second->decrRef(); + + for (map::iterator it=_mapDataArrayDouble.begin() ; it!=_mapDataArrayDouble.end(); it++ ) + if ((*it).second!=0) (*it).second->decrRef(); + + if (_driver !=0) {delete _driver; _driver=0;} + if (_topology!=0 && _owns_topology) {delete _topology; _topology=0;} + + if (_joint_finder!=0) {delete _joint_finder; _joint_finder=0;} +} + +/*! constructing the MESH collection from a file + * + * The method creates as many MED-files as there are domains in the + * collection. It also creates a master file that lists all the MED files. + * The MED files created in ths manner contain joints that describe the + * connectivity between subdomains. + * + * \param filename name of the master file that will contain the list of the MED files + * + */ +void 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 + if (_driver!=0) delete _driver; + _driver=0; + + //char filenamechar[256]; + // strcpy(filenamechar,filename.c_str()); + retrieveDriver()->write (filename.c_str(), _domain_selector); +} + +/*! creates or gets the link to the collection driver + */ +MESHCollectionDriver* MESHCollection::retrieveDriver() +{ + if (_driver==0) + { + switch(_driver_type) + { + case MedXML: + _driver=new MESHCollectionMedXMLDriver(this); + break; + case MedAscii: + _driver=new MESHCollectionMedAsciiDriver(this); + break; + default: + throw INTERP_KERNEL::Exception(LOCALIZED("Unrecognized driver")); + } + } + return _driver; +} + + +/*! gets an existing driver + * + */ +MESHCollectionDriver* MESHCollection::getDriver() const { + return _driver; +} + +// /*! retrieves the mesh dimension*/ +int MESHCollection::getMeshDimension() const { + return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension(); +} + +std::vector& MESHCollection::getMesh() { + return _mesh; +} + +std::vector& MESHCollection::getFaceMesh() { + return _faceMesh; +} + +ParaMEDMEM::MEDCouplingUMesh* MESHCollection::getMesh(int idomain) const { + return _mesh[idomain]; +} + +ParaMEDMEM::MEDCouplingUMesh* MESHCollection::getFaceMesh(int idomain) { + return _faceMesh[idomain]; +} + +std::vector& MESHCollection::getCZ() { + return _connect_zones; +} + +Topology* MESHCollection::getTopology() const { + return _topology; +} + +void MESHCollection::setTopology(Topology* topo) { + if (_topology!=0) + { + throw INTERP_KERNEL::Exception(LOCALIZED("topology is already set")); + } + else + _topology = topo; +} + +/*! Method creating the cell graph + * + * \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 MESHCollection::buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array, int *& edgeweights ) //cvwat09 +{ + using std::multimap; + using std::vector; + using std::make_pair; + using std::pair; + + multimap< int, int > node2cell; + multimap< int, int > cell2cell; + multimap< int, int > cell2node; + + vector > > commonDistantNodes; + int nbdomain=_topology->nbDomain(); + if (isParallelMode()) + { + _joint_finder=new JointFinder(*this); + _joint_finder->findCommonDistantNodes(); + commonDistantNodes=_joint_finder->getDistantNodeCell(); + } + + if (MyGlobals::_verbose>500) _joint_finder->print(); + + //looking for reverse nodal connectivity i global numbering + for (int idomain=0; idomainisMyDomain(idomain)) continue; + + /*obsolete + int offsetCell=0, offsetNode=0; + if (isParallelMode()) + { + offsetCell=_domain_selector->getDomainCellShift(idomain); + offsetNode=_domain_selector->getDomainNodeShift(idomain); + }*/ + + ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New(); + ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New(); + int nbNodes=_mesh[idomain]->getNumberOfNodes(); + //cout<<"proc "<getReverseNodalConnectivity(revConn,index); + //problem saturation over 1 000 000 nodes for 1 proc + if (MyGlobals::_verbose>100) cout<<"proc "<getPointer(); + int* revConnPtr=revConn->getPointer(); + //if (MyGlobals::_verbose>100) cout<<"proc "<convertNodeToGlobal(idomain,i); + int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]); + node2cell.insert(make_pair(globalNode, globalCell)); + //cvw cout<<" "<decrRef(); + index->decrRef(); + //vector > > dNC=getDistantNodeCell() + for (int iother=0; iother::iterator it; + int isource=idomain; + int itarget=iother; + for (it=_joint_finder->_distant_node_cell[isource][itarget].begin(); + it!=_joint_finder->_distant_node_cell[isource][itarget].end(); it++) + { + int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first); + int globalCell=(*it).second; + node2cell.insert(make_pair(globalNode, globalCell)); + //cout<<"processor "<100) cout<<"proc "<nbNodes()<nbNodes(); inode++) //on all nodes + { + typedef multimap::const_iterator MI; + pair myRange=node2cell.equal_range(inode); + for (MI cell1=myRange.first; cell1!=myRange.second; cell1++) //on cells with inode + { + pair myNodes1=cell2node.equal_range(cell1->second); //nodes of one cell + for (MI cell2=myRange.first; cell2!=myRange.second; cell2++) //on one of these cell + { + if (cell1->second!=cell2->second) //in cells of my domain + { + pair myNodes2=cell2node.equal_range(cell2->second); //on nodes of this cell + int nbcn=0; //number of common nodes between cells: at least 3 for cells + for (MI it1=myNodes1.first; it1!=myNodes1.second; ++it1) + { + for (MI it2=myNodes2.first; it2!=myNodes2.second; ++it2) + { + if ((*it1).second==(*it2).second) + { + nbcn=nbcn+1 ; break; + } + } + } + if (nbcn>=3) //cvw TODO if 2d cells set at 2 + cell2cell.insert(make_pair(cell1->second,cell2->second)); + //note here there is some global numerotation of cell which come from other domain (not mydomain) + //cout<<" "<second<<"!"<second; //cvw + } + } + } + } + + if (MyGlobals::_verbose>100) cout<<"proc "< index,value; + index.push_back(0); + int idep=0; + + for (int idomain=0; idomainisMyDomain(idomain)) continue; + int nbCells=_mesh[idomain]->getNumberOfCells(); + for (int icell=0; icellconvertCellToGlobal(idomain,icell); + multimap::iterator it; + pair::iterator,multimap::iterator> ret; + ret=cell2cell.equal_range(globalCell); + //cout<<" "<100) + { + std::cout<<"\nproc "<<_domain_selector->rank()<<" : end MESHCollection::buildCellGraph "<< + index.size()-1<<" "<1) + { + for (int i=0; i<10; ++i) cout< node2cell; + multimap< int, int > cell2cell; + + vector > > commonDistantNodes; + int nbdomain=_topology->nbDomain(); + if (isParallelMode()) + { + _joint_finder=new JointFinder(*this); + _joint_finder->findCommonDistantNodes(); + commonDistantNodes=_joint_finder->getDistantNodeCell(); + } + + //looking for reverse nodal connectivity i global numbering + for (int idomain=0; idomainisMyDomain(idomain)) continue; + int offset=0, procOffset=0; + if (isParallelMode()) + { + offset=_domain_selector->getDomainCellShift(idomain); + procOffset=_domain_selector->getProcShift(); + } + ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New(); + ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New(); + _mesh[idomain]->getReverseNodalConnectivity(revConn,index); + int* index_ptr=index->getPointer(); + ParaMEDMEM::DataArrayInt* globalRevConn=revConn->deepCpy(); + _topology->convertCellToGlobal(idomain, + revConn->getPointer(), + index_ptr[_mesh[idomain]->getNumberOfNodes()], + globalRevConn->getPointer()); + + int min=std::numeric_limits::max(),max=-1; + for (int i=0; i< index_ptr[_mesh[idomain]->getNumberOfNodes()];i++) + { + if (revConn->getPointer()[i]getPointer()[i]; + if (revConn->getPointer()[i]>max) max=revConn->getPointer()[i]; + } + int* globalNodeIds=new int[_mesh[idomain]->getNumberOfNodes()]; + _topology->getNodeList(idomain,globalNodeIds); + + int* globalRevConnPtr=globalRevConn->getPointer(); + for (int i=0; i<_mesh[idomain]->getNumberOfNodes();i++) + { + for (int icell=index_ptr[i]; icellnbDomain();mydomain++) + { + if (_domain_selector->isMyDomain(mydomain)) continue; + multimap::iterator iter; + for (iter=commonDistantNodes[idomain][mydomain].begin();iter!=commonDistantNodes[idomain][mydomain].end();iter++) + { + int ilocnode=iter->first; + int icell=iter->second; + node2cell.insert(make_pair(globalNodeIds[ilocnode],icell+offset)); + } + } + } + + globalRevConn->decrRef(); + revConn->decrRef(); + index->decrRef(); + delete[] globalNodeIds; + } + + //creating graph arcs (cell to cell relations) + //arcs are stored in terms of (index,value) notation + // 0 3 5 6 6 + // 1 2 3 2 3 3 + // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3) + // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell + + //warning here one node have effective number of cell with it + //but cell could have more than effective nodes + //is it normal??? because other equals nodes in other domain (with other inode)??? + + int mincell,maxcell; + if (isParallelMode()) + { + mincell=_domain_selector->getProcShift(); + maxcell=mincell; + for (int i=0; iisMyDomain(i)) maxcell+=_mesh[i]->getNumberOfCells(); + } + else + { + mincell=0; + maxcell=_topology->nbCells(); + } + for (int inode=0; inode<_topology->nbNodes(); inode++) //on all nodes + { + typedef multimap::const_iterator MI; + pair myRange = node2cell.equal_range(inode); + + + for (MI cell1=myRange.first; cell1!=myRange.second; cell1++) //on cells with inode + { + for (MI cell2 = myRange.first; cell2!=myRange.second; cell2++) //on one of these cell + { + if (cell1->second!=cell2->second && cell1->second>=mincell && cell1->secondsecond,cell2->second)); + } + } + } + } + //cout<<"proc "<::const_iterator iter; + + iter=cell2cell.begin(); + vector index,value; + index.push_back(0); + int idep=0; + while (iter != cell2cell.end()) + { + multimap::const_iterator next_iter = cell2cell.upper_bound(iter->first); + int size=0; + for (multimap::const_iterator current_iter=iter; current_iter!=next_iter; current_iter++) + { + int ival=current_iter->second; //no adding one existing yet + for (int i=idep ; i10) cout<<"proc "< 0")); + MEDPARTITIONER::MEDSKYLINEARRAY* array=0; + int* edgeweights=0; + + //cout<<"Building cell graph... "; + // if ( _domain_selector ) + // buildCellGraphParallel(array,edgeweights); + // else + buildCellGraph(array,edgeweights); //cvwat09 + //MPI_Barrier(MPI_COMM_WORLD); + //cout<<"proc "<10) cout<<"METISGraph"<10) cout<<"SCOTCHGraph"<setEdgesWeights(user_edge_weights); + if (user_vertices_weights!=0) + cellGraph->setVerticesWeights(user_vertices_weights); + + if (MyGlobals::_is0verbose>10) cout<<"partitioning graph on "<partGraph(nbdomain, options_string, _domain_selector); //cvwat10 + + if (MyGlobals::_is0verbose>10) cout<<"building new topology"<11) cout<<"proc "< domains; + for (int i=0; i<_topology->nbCells(); i++) + { + domains.insert(partition[i]); + } + int nbdomain=domains.size(); + + cellGraph=(Graph*)(new UserGraph(array, partition, _topology->nbCells())); + + //cellGraph is a shared pointer + Topology* topology = new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension()); + + // if (array!=0) delete array; + delete cellGraph; + return topology; +} + +void MESHCollection::setDomainNames(const std::string& name) +{ + for (int i=0; i<_topology->nbDomain(); i++) + { + std::ostringstream oss; + oss<isMyDomain(i)) + _mesh[i]->setName(oss.str().c_str()); + } +} + +ParaMEDMEM::DataArrayDouble* 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); + { + int rank=MyGlobals::_rank; + string tag="ioldFieldDouble="+intToStr(iold); + string descriptionIold=descriptionField+serializeFromString(tag); + if (_mapDataArrayDouble.find(descriptionIold)!=_mapDataArrayDouble.end()) + { + if (MyGlobals::_verbose>300) cout<<"proc "<reprZip()<200) cout<<"proc "<10) + cout<<"proc "<getArray(); + //to know names of components + vector browse=browseFieldDouble(f2); + //done yet + //double time=f2->getTime(IT,DT); + //browse.push_back("time="+doubleToStr(time)); + string localFieldInformation=descriptionIold+serializeFromVectorOfString(browse); + if (MyGlobals::_verbose>10) cout<<"proc "<incrRef(); //free field, keep res + f2->decrRef(); + _mapDataArrayDouble[descriptionIold]=res; + + //duplicate it! because f2->decRef!! + //DataArrayDouble* res=f2->getArray()->deepCpy(); + //f2->decrRef(); + //cout<reprZip()< r2; + //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data + for (int i=0; i<_fieldDescriptions.size(); i++) + { + vector r1=deserializeToVectorOfString(_fieldDescriptions[i]); + for (int i=0; i& inodesFace,vector< int >& inodesCell) +{ + int ires=0; + int nbok=inodesFace.size(); + for (int i=0; inbDomain(); inew++) + { + if (isParallelMode() && _domain_selector->isMyDomain(inew)) + { + if (MyGlobals::_verbose>200) + std::cout<<"proc "<getNumberOfCells()< nodeIds; + //cout<<"proc "<getReverseNodalConnectivity(revNodalCel,revNodalIndxCel); + int *revC=revNodalCel->getPointer(); + int *revIndxC=revNodalIndxCel->getPointer(); + + vector< int > faceOnCell; + vector< int > faceNotOnCell; + int nbface=mfac->getNumberOfCells(); + for (int iface=0; iface inodesFace; + mfac->getNodeIdsOfCell(iface, inodesFace); + int nbnodFace=inodesFace.size(); + //set inodesFace in mcel + for (int i=0; i inodesCell; + mcel->getNodeIdsOfCell(icel, inodesCell); + ok=isFaceOncell(inodesFace, inodesCell); + if (ok) break; + } + if (ok) + { + faceOnCell.push_back(iface); + //if (MyGlobals::_is0verbose) cout<<"face on cell "<decrRef(); + revNodalIndxCel->decrRef(); + + string cle; + cle=cle1ToStr("filterFaceOnCell",inew); + _mapDataArrayInt[cle]=createDataArrayIntFromVector(faceOnCell); + cle=cle1ToStr("filterNotFaceOnCell",inew); + _mapDataArrayInt[cle]=createDataArrayIntFromVector(faceNotOnCell); + + /*ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New(); + ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New(); + _mesh[idomain]->getReverseNodalConnectivity(revConn,index); + int* index_ptr=index->getPointer();*/ + + /*if (MyGlobals::_is0verbose) + { + cout<<"proc "<nbDomain(); inew++) + { + if (isParallelMode() && _domain_selector->isMyDomain(inew)) + { + if (MyGlobals::_verbose>1) std::cout<<"proc "<getNumberOfCells()<buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx); + revDesc->decrRef(); + desc->decrRef(); + descIndx->decrRef(); + int nbOfCells=meshDM1->getNumberOfCells(); + const int *revDescIndxC=revDescIndx->getConstPointer(); + std::vector boundaryCells; + for(int i=0; idecrRef(); + bool keepCoords=true; + MEDCouplingUMesh *ret=(MEDCouplingUMesh *)meshDM1->buildPartOfMySelf(&boundaryCells[0],&boundaryCells[0]+boundaryCells.size(),keepCoords); + meshDM1->decrRef(); + //don't know what to do with result yet.. + //_faceMesh[inew]->decrRef(); + //_faceMesh[inew]=ret; + } + } +} +*/ + diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx new file mode 100644 index 000000000..e30a9b639 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollection.hxx @@ -0,0 +1,249 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef __MEDPARTITIONER_MESHCOLLECTION_HXX__ +#define __MEDPARTITIONER_MESHCOLLECTION_HXX__ + +#include "MEDPARTITIONER.hxx" + +#include "MEDPARTITIONER_Graph.hxx" +#include "MEDCouplingUMesh.hxx" + +//#include "MEDPARTITIONER_FaceModel.hxx" +//#include "boost/shared_ptr.hpp" +#include +#include +#include +namespace ParaMEDMEM +{ + class MEDCouplingUMesh; + class DataArrayInt; +} + + +namespace MEDPARTITIONER +{ + + class Topology; + class MESHCollectionDriver; + class ParaDomainSelector; + class MEDSKYLINEARRAY; + class CONNECTZONE; + class JointFinder; + typedef enum{MedAscii, MedXML, Undefined} DriverType; + + typedef std::multimap, std::pair > NodeMapping ; + typedef std::vector > NodeList; + + class MEDPARTITIONER_EXPORT MESHCollection + { + + public: + + //Default constructor + MESHCollection(); + + //Constructing from an existing mesh and a new topology + MESHCollection(MESHCollection&, Topology*, bool family_splitting=false, bool create_empty_groups=false); + + //Constructing the mesh collection from a file + MESHCollection(const std::string& filename); + + //Constructing the mesh collection from a file + MESHCollection(const std::string& filename, ParaDomainSelector& domainSelector); + + //Constructing the mesh collection from a file + MESHCollection(const std::string& filename, const std::string& meshname); + + ~MESHCollection(); + + bool isParallelMode() const { return _domain_selector; } + + //writing to a distributed file + void write(const std::string& filename); + + //getting the driver + MESHCollectionDriver* retrieveDriver(); + MESHCollectionDriver* getDriver() const; + void setDriverType(MEDPARTITIONER::DriverType type) {_driver_type=type;} + + //creation of the cell graph + void buildCellGraph(MEDPARTITIONER::MEDSKYLINEARRAY* & array,int *& edgeweights ); + + //creation and partition of the associated graph + Topology* createPartition(int nbdomain, Graph::splitter_type type = Graph::METIS, + const std::string& ="", int* edgeweights=0, int* verticesweights=0); + + //creation of a user specified partition + Topology* createPartition(const int* partition); + + //getting mesh dimension + int getMeshDimension() const ; + + //getting a reference to mesh vector + std::vector& getMesh(); + std::vector& getFaceMesh(); + std::vector >& getGroupMeshes(); + + ParaMEDMEM::MEDCouplingUMesh* getMesh(int idomain) const; + ParaMEDMEM::MEDCouplingUMesh* getFaceMesh(int idomain); + std::vector& getGroupMeshes(int idomain); + + std::vector& getCellFamilyIds() {return _cellFamilyIds;} + std::vector& getFaceFamilyIds() {return _faceFamilyIds;} + + std::map& getMapDataArrayInt() {return _mapDataArrayInt;} + std::map& getMapDataArrayDouble() {return _mapDataArrayDouble;} + + std::map& getFamilyInfo() {return _familyInfo;} + std::map >& getGroupInfo() {return _groupInfo;} + + ParaMEDMEM::DataArrayDouble* getField(std::string descriptionField, int iold); + std::vector& getFieldDescriptions() {return _fieldDescriptions;} + void prepareFieldDescriptions(); + void filterFaceOnCell(); + + //getting a reference to connect zones vector + std::vector& getCZ(); + + //getting a pointer to topology + Topology* getTopology() const ; + ParaDomainSelector* getParaDomainSelector() const{return _domain_selector;} + //settig a new topology + void setTopology(Topology* topology); + + //getting/setting the name of the global mesh (as opposed + //to the name of a subdomain \a nn, which is name_nn) + std::string getName() const {return _name;} + void setName(const std::string& name){_name=name;} + void setDomainNames(const std::string& name); + + //getting/setting the description of the global mesh + std::string getDescription() const {return _description;} + void setDescription(const std::string& name) { _description=name;} + + //creates the node mapping between an old collection and the present one + void createNodeMapping(MESHCollection& initialCollection, + std::multimap,std::pair >& nodeMapping); + + void castCellMeshes(MESHCollection& initialCollection, + std::vector > >& new2oldIds); + + //creates faces on the new collection + void castFaceMeshes(MESHCollection& initialCollection, + const std::multimap, std::pair >& nodeMapping, + std::vector > >& new2oldIds); + + private: + + void castIntField2(std::vector& meshesCastFrom, + std::vector& meshesCastTo, + std::vector& arrayFrom, + std::string nameArrayTo); + + void castAllFields(MESHCollection& initialCollection, + std::string nameArrayTo); + + void findCommonDistantNodes(std::vector > >& commonDistantNodes); + + void remapIntField(const ParaMEDMEM::MEDCouplingUMesh& sourceMesh, + const ParaMEDMEM::MEDCouplingUMesh& targetMesh, + const int* fromArray, + int* toArray); + + void remapIntField2(int inew, int iold, + const ParaMEDMEM::MEDCouplingUMesh& sourceMesh, + const ParaMEDMEM::MEDCouplingUMesh& targetMesh, + const int* fromArray, + std::string nameArrayTo); + + void remapDoubleField3(int inew, int iold, + ParaMEDMEM::DataArrayDouble* fromArray, + std::string nameArrayTo, + std::string descriptionField); + + //!link to mesh_collection topology + Topology* _topology; + + //!control over topology + bool _owns_topology; + + //!link to graph + //Graph* _cell_graph; + + //! Driver for read/write operations + MESHCollectionDriver* _driver; + + //! Parallelizer - mark of parallel execution mode + ParaDomainSelector* _domain_selector; + + //!links to meshes + std::vector _mesh; + std::vector _faceMesh; + + //!index of a non empty mesh within _mesh (in parallel mode all of meshes can be empty) + int _i_non_empty_mesh; + + //!links to connectzones + std::vector _connect_zones; + + //!family ids storages + std::vector _cellFamilyIds; + std::vector _faceFamilyIds; + + //!DataArrayInt* storages + std::map _mapDataArrayInt; + //!DataArrayDouble* storages + std::map _mapDataArrayDouble; + + std::vector _fieldDescriptions; //fields to be partitioned + + //!group family conversion + std::map _familyInfo; + std::map > _groupInfo; + + //!list of groups that are not to be splitted + std::vector _indivisible_regions; + + //!name of global mesh + std::string _name; + + //!description of global mesh + std::string _description; + + //! specifies the driver associated to the collection + DriverType _driver_type; + + /*! flag specifying that the splitter should create boundary constituent entity + so that they are written in joints*/ + bool _subdomain_boundary_creates; + + /*! flag specifying that families must be preserved by the + splitting*/ + bool _family_splitting; + + /*! flag specifying that groups must be created on all domains, + even if they are empty*/ + bool _create_empty_groups; + + JointFinder* _joint_finder; + }; + +}//of namespace + +#endif /*MESHCOLLECTION_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionDriver.cxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionDriver.cxx new file mode 100644 index 000000000..78b1e5edf --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionDriver.cxx @@ -0,0 +1,557 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" +#include "MEDLoader.hxx" +#include "MEDFileMesh.hxx" + +extern "C" { +#include "med.h" +} +//MEDPARTITIONER includes +#include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_MESHCollectionDriver.hxx" +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" +#include "MEDPARTITIONER_utils.hxx" + +using namespace MEDPARTITIONER; +using namespace std; + +//template inclusion +//#include "MEDPARTITIONER_MESHCollectionDriver.H" + +// med_geometrie_element typmai[MED_NBR_GEOMETRIE_MAILLE+2] = { MED_POINT1, +// MED_SEG2, +// MED_SEG3, +// MED_TRIA3, +// MED_TRIA6, +// MED_QUAD4, +// MED_QUAD8, +// MED_TETRA4, +// MED_TETRA10, +// MED_HEXA8, +// MED_HEXA20, +// MED_PENTA6, +// MED_PENTA15, +// MED_PYRA5, +// MED_PYRA13, +// MED_POLYGONE, +// MED_POLYEDRE }; + +MESHCollectionDriver::MESHCollectionDriver(MESHCollection* collection):_collection(collection) +{ +} + + +/*!reads a unique MED File v>=2.1 + * and mounts the corresponding mesh in memory + *\param filename binary file + *\param meshname mesh name in the MED file + * */ +int MESHCollectionDriver::readSeq(const char* filename, const char* meshname) +{ + cout<<"readSeq"<getMesh()).push_back(mfm->getLevel0Mesh(false)); + (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false)); + + //reading family ids + + ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy()); + ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy()); + (_collection->getCellFamilyIds()).push_back(cellIds); + (_collection->getFaceFamilyIds()).push_back(faceIds); + + //reading groups + (_collection->getFamilyInfo())=mfm->getFamilyInfo(); + (_collection->getGroupInfo())=mfm->getGroupInfo(); + +// (_collection->getGroupMeshes()).resize(groupNames.size()); + +// for (int i=0; i< groupNames.size();i++) +// { +// vector myGroup; +// myGroup.push_back(groupNames[i]); +// (_collection->getGroupMeshes())[i].push_back(MEDLoader::ReadUMeshFromGroups(filename,meshname,-1,myGroup)); +// } + + + (_collection->getCZ()).clear(); + /*cvw + vector cellglobal,nodeglobal,faceglobal; + cellglobal.resize(1); + nodeglobal.resize(1); + faceglobal.resize(1); + cellglobal[0]=0; + nodeglobal[0]=0; + faceglobal[0]=0; + //creation of topology from mesh + //connectzone argument is 0 + ParallelTopology* aPT = new ParallelTopology + ((_collection->getMesh()), (_collection->getCZ()), cellglobal, nodeglobal, faceglobal); + */ + + ParallelTopology* aPT = new ParallelTopology((_collection->getMesh())); + _collection->setTopology(aPT); + _collection->setName(meshname); + _collection->setDomainNames(meshname); + return 0; +} + + +//================================================================================ +/*! + * \brief Return mesh dimension from distributed med file had being read + */ +//================================================================================ + +void MESHCollectionDriver::readSubdomain(vector& cellglobal, //cvwat03 + vector& faceglobal, + vector& nodeglobal, int idomain) +{ + string meshname=MyGlobals::_meshNames[idomain]; + string file=MyGlobals::_fileNames[idomain]; + + //cout << "Reading "< nonEmpty=mfm->getNonEmptyLevels(); + + try + { + (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false); + //reading families groups + ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy()); + (_collection->getCellFamilyIds())[idomain]=cellIds; + } + catch(...) + { + (_collection->getMesh())[idomain]=createEmptyMEDCouplingUMesh(); // or 0 if you want tests; + ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New(); + empty->alloc(0,1); + (_collection->getCellFamilyIds())[idomain]=empty; + cout<<"\nNO Level0Mesh (Cells)\n"; + } + try + { + if (nonEmpty.size()>1 && nonEmpty[1]==-1) + { + (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false); + //reading families groups + ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy()); + (_collection->getFaceFamilyIds())[idomain]=faceIds; + } + else + { + throw "no faces"; + } + } + catch(...) + { + //ParaMEDMEM::MEDCouplingUMesh *umesh=ParaMEDMEM::MEDCouplingUMesh::New(); //empty one + //umesh->setMeshDimension(3); + //umesh->allocateCells(0); + //int nb=umesh->getNumberOfCells(); //no use if no allocateCells(0)! because thrown exception + //cout<<"\nempty mesh"<getFaceMesh())[idomain]=createEmptyMEDCouplingUMesh(); // or 0 if you want test; + ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New(); + (_collection->getFaceFamilyIds())[idomain]=empty; + if (MyGlobals::_verbose>10) cout<<"proc "<getFamilyInfo()=mfm->getFamilyInfo(); + _collection->getGroupInfo()=mfm->getGroupInfo(); + + mfm->decrRef(); + + vector localInformation; + string str; + localInformation.push_back(str+"ioldDomain="+intToStr(idomain)); + localInformation.push_back(str+"meshName="+meshname); + MyGlobals::_generalInformations.push_back(serializeFromVectorOfString(localInformation)); + vector localFields=browseAllFieldsOnMesh(file, meshname, idomain); //cvwat07 + if (localFields.size()>0) + MyGlobals::_fieldDescriptions.push_back(serializeFromVectorOfString(localFields)); + //cout<< "End Reading "< nonEmpty=mfm->getNonEmptyLevels(); + + try + { + (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false); + //reading families groups + ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy()); + (_collection->getCellFamilyIds())[idomain]=cellIds; + } + catch(...) + { + (_collection->getMesh())[idomain]=createEmptyMEDCouplingUMesh(); // or 0 if you want tests; + ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New(); + empty->alloc(0,1); + (_collection->getCellFamilyIds())[idomain]=empty; + cout<<"\nNO Level0Mesh (Cells)\n"; + } + try + { + if (nonEmpty.size()>1 && nonEmpty[1]==-1) + { + (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false); + //reading families groups + ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy()); + (_collection->getFaceFamilyIds())[idomain]=faceIds; + } + else + { + throw "no faces"; + } + } + catch(...) + { + (_collection->getFaceMesh())[idomain]=createEmptyMEDCouplingUMesh(); // or 0 if you want test; + ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New(); + (_collection->getFaceFamilyIds())[idomain]=empty; + if (MyGlobals::_verbose>10) cout<<"proc "<getFamilyInfo()=mfm->getFamilyInfo(); + _collection->getGroupInfo()=mfm->getGroupInfo(); + + mfm->decrRef(); + + vector localInformation; + string str; + localInformation.push_back(str+"ioldDomain="+intToStr(idomain)); + localInformation.push_back(str+"meshName="+meshname); + MyGlobals::_generalInformations.push_back(serializeFromVectorOfString(localInformation)); + vector localFields=browseAllFieldsOnMesh(file, meshname, idomain); //cvwat07 + if (localFields.size()>0) + MyGlobals::_fieldDescriptions.push_back(serializeFromVectorOfString(localFields)); + //cout<< "End Reading "< meshes; + ParaMEDMEM::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain); + ParaMEDMEM::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain); + ParaMEDMEM::MEDCouplingUMesh* faceMeshFilter=0; + + string finalMeshName=extractFromDescription(MyGlobals::_generalInformations[0], "finalMeshName="); + string cleFilter=cle1ToStr("filterFaceOnCell",idomain); + DataArrayInt* filter=0; + if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end()) + { + filter=_collection->getMapDataArrayInt().find(cleFilter)->second; + int* index=filter->getPointer(); + faceMeshFilter=(MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true); + faceMesh=faceMeshFilter; + } + cellMesh->setName(finalMeshName.c_str()); + meshes.push_back(cellMesh); + + //cellMesh->zipCoords(); + //faceMesh->zipCoords(); + + faceMesh->checkCoherency(); + if (faceMesh->getNumberOfCells()>0) + { + faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10); + meshes.push_back(faceMesh); + } + + /*do not work + ParaMEDMEM::MEDFileUMesh* mfm2=ParaMEDMEM::MEDFileUMesh::New(); + MEDFileUMesh* mfm2 = static_cast(cellMesh->getMeshes()->getMeshAtPos(0)); + MEDFileUMesh* mfm2 = ParaMEDMEM::MEDFileUMesh::New(cellMesh); + string fname="FUM_"+distfilename; + mfm2->setMeshAtLevel(0, cellMesh ); + mfm2->setMeshAtLevel(-1, faceMesh ); + mfm2->write(fname.c_str(),0); + mfm2->decrRef(); + */ + + ParaMEDMEM::MEDCouplingUMesh* boundaryMesh=0; + //ParaMEDMEM::MEDCouplingUMesh* boundaryMesh1=0; + //ParaMEDMEM::MEDCouplingUMesh* finalboundaryMesh=0; + if (MyGlobals::_creates_boundary_faces>0) + { + //try to write Boundary meshes + bool keepCoords=false; //TODO or true + boundaryMesh=(MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords); + boundaryMesh->setName("boundaryMesh"); + //cout<<"boundaryMesh "<getNumberOfCells()<checkCoherency(); + //boundaryMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10); + //meshes.push_back(boundaryMesh); + //string boundary="boundary_"+distfilename; + + /*try to find joint do no work + int rang=MyGlobals::_rank; + if (rang==1) (_collection->getParaDomainSelector())->sendMesh(*(boundaryMesh),0); + if (rang==0) + { + (_collection->getParaDomainSelector())->recvMesh(boundaryMesh1,1); + //vector meshes; + //vector corr; + //meshes.push_back(boundaryMesh); + //meshes.push_back(boundaryMesh1); + //need share the same coords + //boundaryMesh1->tryToShareSameCoordsPermute(*boundaryMesh, 1e-10); + //finalboundaryMesh=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,2, corr); + //boundaryMesh=finalboundaryMesh; + + boundaryMesh->zipCoords(); + boundaryMesh1->zipCoords(); + finalboundaryMesh=MEDCouplingUMesh::MergeUMeshes(boundaryMesh,boundaryMesh1); + DataArrayInt* commonNodes=0; + commonNodes=finalboundaryMesh->zipCoordsTraducer(); + boundaryMesh=finalboundaryMesh; + cout<<"zipcoords"<repr()<decrRef(); + + + if (boundaryMesh!=0) + { + //doing that testMesh becomes second mesh sorted by alphabetical order of name + MEDLoader::WriteUMesh(distfilename.c_str(), boundaryMesh, false); + boundaryMesh->decrRef(); + } + + //cout<<"familyInfo :\n"<getFamilyInfo())<getGroupInfo())<getMesh(idomain)->getName()); + + /*example of adding new family + (_collection->getFamilyInfo())["FaceNotOnCell"]=-500; + vector FaceNotOnCell; + FaceNotOnCell.push_back("FaceNotOnCell"); + (_collection->getGroupInfo())["FaceNotOnCell"]=FaceNotOnCell; + */ + + mfm->setFamilyInfo(_collection->getFamilyInfo()); + mfm->setGroupInfo(_collection->getGroupInfo()); + + //cvwat08 + //without filter mfm->setFamilyFieldArr(-1,(_collection->getFaceFamilyIds())[idomain]); + + string cle=cle1ToStr("faceFamily_toArray",idomain); + if (_collection->getMapDataArrayInt().find(cle)!=_collection->getMapDataArrayInt().end()) + { + DataArrayInt* fam=_collection->getMapDataArrayInt().find(cle)->second; + DataArrayInt* famFilter=0; + if (filter!=0) + { + int* index=filter->getPointer(); + int nbTuples=filter->getNbOfElems(); + //not the good one...buildPartOfMySelf do not exist for DataArray + //Filter=fam->renumberAndReduce(index, filter->getNbOfElems()); + famFilter=DataArrayInt::New(); + famFilter->alloc(nbTuples,1); + int* pfamFilter=famFilter->getPointer(); + int* pfam=fam->getPointer(); + for (int i=0; isetFamilyFieldArr(-1,fam); + famFilter->decrRef(); + } + //cout<<"proc "<setFamilyFieldArr(-1,fam); + // if (famFilter!=0) famFilter->decrRef(); + } + + /*example visualisation of filter + if (_collection->getMapDataArrayInt().find(cle)!=_collection->getMapDataArrayInt().end()) + { + DataArrayInt* fam=_collection->getMapDataArrayInt().find(cle)->second; + string cle2=cle1ToStr("filterNotFaceOnCell",idomain); + if (_collection->getMapDataArrayInt().find(cle2)!=_collection->getMapDataArrayInt().end()) + { + DataArrayInt* filter=_collection->getMapDataArrayInt().find(cle2)->second; + int* index=filter->getPointer(); + int* pfam=fam->getPointer(); + for (int i=0; igetNbOfElems(); i++) pfam[index[i]]=-500; + } + mfm->setFamilyFieldArr(-1,fam); + //mfm->setFamilyFieldArr(-1,_collection->getMapDataArrayInt().find(cle)->second); + } + */ + + cle=cle1ToStr("cellFamily_toArray",idomain); + if (_collection->getMapDataArrayInt().find(cle)!=_collection->getMapDataArrayInt().end()) + mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(cle)->second); + + mfm->write(distfilename.c_str(),0); + cle="/inewFieldDouble="+intToStr(idomain)+"/"; + + map::iterator it; + int nbfFieldFound=0; + for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++) + { + string desc=(*it).first; + size_t found=desc.find(cle); + if (found==string::npos) continue; + if (MyGlobals::_verbose>20) cout<<"proc "<setName(fieldName.c_str()); + field->setMesh(mfm->getLevel0Mesh(false)); + DataArrayDouble* da=(*it).second; + + //get information for components etc.. + vector r1; + r1=selectTagsInVectorOfString(MyGlobals::_generalInformations,"fieldName="+fieldName); + r1=selectTagsInVectorOfString(r1,"typeField="+intToStr(typeField)); + r1=selectTagsInVectorOfString(r1,"DT="+intToStr(DT)); + r1=selectTagsInVectorOfString(r1,"IT="+intToStr(IT)); + //not saved in file? field->setDescription(extractFromDescription(r1[0], "fieldDescription=").c_str()); + int nbc=strToInt(extractFromDescription(r1[0], "nbComponents=")); + //double time=strToDouble(extractFromDescription(r1[0], "time=")); + if (nbc==da->getNumberOfComponents()) + { + for (int i=0; isetInfoOnComponent(i,extractFromDescription(r1[0], "componentInfo"+intToStr(i)+"=").c_str()); + } + else + { + cerr<<"Problem On field "<getNumberOfComponents()<setArray(da); + field->setTime(time,DT,IT); + field->checkCoherency(); + try + { + MEDLoader::WriteField(distfilename.c_str(),field,false); + //if entityName=="MED_NODE_ELEMENT" + //AN INTERP_KERNEL::EXCEPTION HAS BEEN THROWN : Not implemented other profile fitting from already written mesh for fields than on NODES and on CELLS.********** + //modification MEDLoader.cxx done + } + catch(INTERP_KERNEL::Exception& e) + { + //cout trying rewrite all data, only one field defined + string tmp,newName=distfilename; + tmp+="_"+fieldName+"_"+intToStr(nbfFieldFound)+".med"; + newName.replace(newName.find(".med"),4,tmp); + cout<<"WARNING : writeMedFile : new file name with only one field :"<decrRef(); + +} + +/* +void writeFieldNodeCellTryingToFitExistingMesh(const char *fileName, const ParaMEDMEM::MEDCouplingFieldDouble *f) +{ + med_int numdt,numo; + med_float dt; + int nbComp=f->getNumberOfComponents(); + med_idt fid=appendFieldSimpleAtt(fileName,f,numdt,numo,dt); + std::list split; + prepareCellFieldDoubleForWriting(f,thisMeshCellIdsPerType,split); + const double *pt=f->getArray()->getConstPointer(); + int number=0; + for(std::list::const_iterator iter=split.begin();iter!=split.end();iter++) + { + INTERP_KERNEL::AutoPtr nommaa=MEDLoaderBase::buildEmptyString(MED_NAME_SIZE); + MEDLoaderBase::safeStrCpy(f->getMesh()->getName(),MED_NAME_SIZE,nommaa,MEDLoader::_TOO_LONG_STR); + INTERP_KERNEL::AutoPtr profileName=MEDLoaderBase::buildEmptyString(MED_NAME_SIZE); + std::ostringstream oss; oss << "Pfl" << f->getName() << "_" << number++; + MEDLoaderBase::safeStrCpy(oss.str().c_str(),MED_NAME_SIZE,profileName,MEDLoader::_TOO_LONG_STR); + const std::vector& ids=(*iter).getCellIdPerType(); + int *profile=new int [ids.size()]; + std::transform(ids.begin(),ids.end(),profile,std::bind2nd(std::plus(),1)); + MEDprofileWr(fid,profileName,ids.size(),profile); + delete [] profile; + MEDfieldValueWithProfileWr(fid,f->getName(),numdt,numo,dt,MED_NODE_CELL,typmai3[(int)(*iter).getType()],MED_COMPACT_PFLMODE,profileName, + MED_NO_LOCALIZATION,MED_FULL_INTERLACE,MED_ALL_CONSTITUENT,(*iter).getNbOfTuple(),(const unsigned char*)pt); + pt+=(*iter).getNbOfTuple()*nbComp; + } + MEDfileClose(fid); +}*/ + diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionDriver.hxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionDriver.hxx new file mode 100644 index 000000000..53ec045db --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionDriver.hxx @@ -0,0 +1,64 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef MEDPARTITIONER_MESHCOLLECTIONDRIVER_HXX_ +#define MEDPARTITIONER_MESHCOLLECTIONDRIVER_HXX_ + +#include "MEDPARTITIONER.hxx" + +#include +#include + +namespace MEDPARTITIONER +{ + class MESHCollection; + class ParaDomainSelector; + + class MEDPARTITIONER_EXPORT MESHCollectionDriver + { + public: + + MESHCollectionDriver(MESHCollection*); + virtual ~MESHCollectionDriver(){} + + virtual int read(const char*, ParaDomainSelector* sel=0)=0; + int readSeq(const char*,const char*); + + virtual void write(const char*, ParaDomainSelector* sel=0)=0; + + protected: + + void readSubdomain(std::vector& cellglobal, + std::vector& faceglobal, + std::vector& nodeglobal, int idomain); + void readSubdomain(int idomain); + void writeMedFile(int idomain, const std::string& distfilename); + + + MESHCollection* _collection; + + //to Globals + //std::vector _filename; + //std::vector _meshname; + + }; + +} + + +#endif /*MESHCOLLECTIONDRIVER_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedAsciiDriver.H b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedAsciiDriver.H new file mode 100644 index 000000000..c253b586f --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedAsciiDriver.H @@ -0,0 +1,66 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef MEDSPLITTER_MESHCOLLECTIONMEDASCIIDRIVER_H +#define MEDSPLITTER_MESHCOLLECTIONMEDASCIIDRIVER_H + + +#include +#include +#include +#include + +/*!reads a distributed field + * + * \param fields vector of fields (one field per subdomain) + * \param fieldname name of the field + * \param itnumber number of iteration + * \param ordernumber internal number inside the iteration + * */ +template +void MESHCollectionMedAsciiDriver::_readFields(vector* >& fields,char* fieldname, int itnumber, int ordernumber) +{ + for (int i=0; i<_collection->getMesh().size(); i++) + { + char filename[256]; + strcpy(filename,_filename[i].c_str()); + cout << "maillage : " << filename << " champ : " << fieldname << endl; + // MEDMEM::FIELD* field = new MEDMEM::FIELD(MEDMEM::MED_DRIVER,filename,fieldname,itnumber,ordernumber); + fields.push_back (new MEDMEM::FIELD(MEDMEM::MED_DRIVER,filename,fieldname,itnumber,ordernumber)); + } +} + + +/*!writes a distributed field + * + * \param fields vector of fields (one field per subdomain) + * \param fieldname name of the field + * */ +template +void MESHCollectionMedAsciiDriver::_writeFields(vector* >& fields,char* fieldname) +{ + for (int i=0; i<_collection->getMesh().size(); i++) + { + char filename[256]; + strcpy(filename,_filename[i].c_str()); + int driverid = fields[i]->addDriver(MEDMEM::MED_DRIVER, filename, fieldname); + fields[i]->write(driverid); + } +} + +#endif diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedAsciiDriver.cxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedAsciiDriver.cxx new file mode 100644 index 000000000..515b1f8ec --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedAsciiDriver.cxx @@ -0,0 +1,182 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "MEDCouplingUMesh.hxx" +#include "MEDLoader.hxx" + +//MEDPARTITIONER includes +#include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_MESHCollectionDriver.hxx" +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" +#include "MEDPARTITIONER_utils.hxx" + +using namespace MEDPARTITIONER; +using namespace std; + +//template inclusion +//#include "MEDPARTITIONER_MESHCollectionMedAsciiDriver.H" + + +MESHCollectionMedAsciiDriver::MESHCollectionMedAsciiDriver(MESHCollection* collection):MESHCollectionDriver(collection) +{ +} + +/*!reads a MED File v>=2.3 + * and mounts the corresponding meshes in memory + * the connect zones are created from the joints + * + *\param filename ascii file containing the list of MED v2.3 files + * */ + +int MESHCollectionMedAsciiDriver::read(const char* filename, ParaDomainSelector* domainSelector) +{ + //distributed meshes + vector cellglobal; + vector nodeglobal; + vector faceglobal; + int nbdomain; + + // reading ascii master file + try + { + ifstream asciiinput(filename); + if (!asciiinput) throw INTERP_KERNEL::Exception(LOCALIZED("Master ASCII File does not exist")); + char charbuffer[512]; + asciiinput.getline(charbuffer,512); + + while (charbuffer[0]=='#') + { + asciiinput.getline(charbuffer,512); + } + + //reading number of domains + nbdomain=atoi(charbuffer); + //cout << "nb domain "<>nbdomain; + MyGlobals::_fileNames.resize(nbdomain); + MyGlobals::_meshNames.resize(nbdomain); + (_collection->getMesh()).resize(nbdomain); + cellglobal.resize(nbdomain); + nodeglobal.resize(nbdomain); + faceglobal.resize(nbdomain); + + if (nbdomain == 0) throw INTERP_KERNEL::Exception(LOCALIZED("Empty ASCII master file")); + for (int i=0; i> mesh >> idomain >> MyGlobals::_meshNames[i] >> host >> MyGlobals::_fileNames[i]; + + //Setting the name of the global mesh (which should be is the same for all the subdomains) + if (i==0) + _collection->setName(mesh); + + if (idomain!=i+1) + { + throw INTERP_KERNEL::Exception(LOCALIZED("domain must be written from 1 to N in ASCII file descriptor")); + } + if ( !domainSelector || domainSelector->isMyDomain(i)) + readSubdomain(cellglobal,faceglobal,nodeglobal, i); + + }//loop on domains + }//of try + catch(...) + { + throw INTERP_KERNEL::Exception(LOCALIZED("I/O error reading parallel MED file")); + } + + //creation of topology from mesh and connect zones + ParallelTopology* aPT = new ParallelTopology + ((_collection->getMesh()), (_collection->getCZ()), cellglobal, nodeglobal, faceglobal); + _collection->setTopology(aPT); + + for (int i=0; igetMesh().size(); + vector filenames; + filenames.resize(nbdomains); + + //loop on the domains + for (int idomain=0; idomainisMyDomain( idomain ) ) + { + if ( !_collection->getMesh()[idomain]->getNumberOfCells()==0 ) continue;//empty domain + MEDLoader::WriteUMesh(distfilename.c_str(),(_collection->getMesh())[idomain],true); + // writeSubdomain(idomain, nbdomains, distfilename.c_str(), domainSelector); + } + } + + // write master file + if ( !domainSelector || domainSelector->rank() == 0 ) + { + ofstream file(filename); + file << "#MED Fichier V 2.3"<<" "<getMesh().size()<<" "<getName() <<" "<< idomain+1 << " " + << (_collection->getMesh())[idomain]->getName() << " localhost " + << filenames[idomain] << " "< *>& filenames, char* fieldname, +// int itnumber, int ordernumber) +// { +// _readFields(filenames,fieldname,itnumber,ordernumber); +// } +// void readFields(vector *>& filenames, char* fieldname, +// int itnumber, int ordernumber) +// { +// _readFields(filenames,fieldname,itnumber,ordernumber); +// } + +// void writeFields(vector *>& filenames, char* fieldname) +// { +// _writeFields( filenames, fieldname); +// } + +// void writeFields(vector *>& filenames, char* fieldname) +// { +// _writeFields( filenames, fieldname); +// } + + + private : +// template void _readFields(vector *>& filenames, char* fieldname, +// int itnumber, int ordernumber); + +// template +// void _writeFields(vector *>& filenames, char* fieldname); + + + + std::string _master_filename; + }; + +} + + +#endif /*MESHCOLLECTIONDRIVER_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.H b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.H new file mode 100644 index 000000000..54b9c18d5 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.H @@ -0,0 +1,122 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef MEDSPLITTER_MESHCOLLECTIONMEDXMLDRIVER_H +#define MEDSPLITTER_MESHCOLLECTIONMEDXMLDRIVER_H + + +#include +#include +#include +#include + + + +/*!reads a distributed field + * + * \param fields vector of fields (one field per subdomain) + * \param fieldname name of the field + * \param itnumber number of iteration + * \param ordernumber internal number inside the iteration + * */ +template +void MESHCollectionMedXMLDriver::_readFields(vector* >& fields,char* fieldname, int itnumber, int ordernumber) +{ + for (int i=0; i<_collection->getMesh().size(); i++) + { + char filename[256]; + strcpy(filename,_filename[i].c_str()); + cout << "maillage : " << filename << " champ : " << fieldname << endl; + // MEDMEM::FIELD* field = new MEDMEM::FIELD(MEDMEM::MED_DRIVER,filename,fieldname,itnumber,ordernumber); + fields.push_back (new MEDMEM::FIELD(MEDMEM::MED_DRIVER,filename,fieldname,itnumber,ordernumber)); + } +} + + +/*!writes a distributed field + * + * \param fields vector of fields (one field per subdomain) + * \param fieldname name of the field + * */ +template +void MESHCollectionMedXMLDriver::_writeFields(vector* >& fields,const char* fieldname) +{ + xmlDocPtr master_doc=xmlParseFile(_master_filename.c_str()); + + if (!master_doc) + throw MEDEXCEPTION("MEDSPLITTER writeFields - Master File does not exist"); + + //number of domains + + xmlXPathContextPtr xpathCtx = xmlXPathNewContext(master_doc); + xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression(BAD_CAST "//mapping/mesh", xpathCtx); + //assuming there is only one mesh in the XML file + xmlNodePtr mesh_node= xpathObj->nodesetval->nodeTab[0]; + + //adds the field to the master file if necessary + bool exist_field =false; + xpathObj = xmlXPathEvalExpression(BAD_CAST "//mapping/mesh/field", xpathCtx); + //assuming there is only one mesh in the XML file + int field_nr = xpathObj->nodesetval->nodeNr; + for (int i=0; inodesetval->nodeTab[i]->properties->children->content, fieldname)==0) + exist_field = true; + } + + xmlNodePtr field_node; + if (!exist_field) + { + field_node = xmlNewChild(mesh_node, 0, BAD_CAST "field",0); + xmlNewProp(field_node,BAD_CAST "name",BAD_CAST fieldname); + } + + + for (int i=0; i<_collection->getMesh().size(); i++) + { + char filename[256]; + strcpy(filename,_filename[i].c_str()); + int driverid = fields[i]->addDriver(MEDMEM::MED_DRIVER, filename, fieldname); + fields[i]->write(driverid); + + //adds the partition to the master file if the field had not been + //added already + if (!exist_field) + { + xmlNodePtr chunk_node= xmlNewChild(field_node,0,BAD_CAST "chunk",0); + char id[8]; + sprintf(id,"%d",i+1); + xmlNewProp(chunk_node,BAD_CAST "subdomain",BAD_CAST id); + //xmlNewProp(chunk_node,BAD_CAST "name", BAD_CAST fieldname); + //xmlNodePtr fieldname_node= + xmlNewChild(chunk_node,0,BAD_CAST "name", BAD_CAST fieldname); + } + } + //Writing file + xmlKeepBlanksDefault(0); + xmlSaveFormatFileEnc(_master_filename.c_str(), master_doc, "UTF-8", 1); + + //Clean up + xmlXPathFreeContext(xpathCtx); + xmlFreeDoc(master_doc); + //xmlCleanupParser(); + +} + +#endif diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.cxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.cxx new file mode 100644 index 000000000..49c2ac35d --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.cxx @@ -0,0 +1,335 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +//MEDSPLITTER includes +#include "MEDCouplingUMesh.hxx" +#include "MEDLoader.hxx" +#include "MEDFileMesh.hxx" +#include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_MESHCollectionDriver.hxx" +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" +#include "MEDPARTITIONER_utils.hxx" + +using namespace MEDPARTITIONER; +using namespace std; + +//template inclusion +//#include "MEDPARTITIONER_MESHCollectionMedXMLDriver.H" + +/*!\class MESHCollectionMedXMLDriver + * + *\brief Driver for MED 3.2 files having XML master files + * + * Driver for reading and writing distributed files + * for which the master file is written in an XML format compliant with + * the MED 3.2 specification. + * The reading and writing of the meshes and fields are apart : + * the meshes must always be written/read before the fields. Reading/Writing fields + * is optional and is done field after field. API for reading/writing fields + * is written with a template so that MEDMEM::FIELD and MEDMEM::FIELD + * can be conveniently handled. +*/ + +MESHCollectionMedXMLDriver::MESHCollectionMedXMLDriver(MESHCollection* collection):MESHCollectionDriver(collection) +{ +} + +/*!reads a MED File XML Master File v>=2.3 + * and mounts the corresponding meshes in memory + * the connect zones are created from the joints + * + *\param filename XML file containing the list of MED v2.3 files + * */ + +int MESHCollectionMedXMLDriver::read(const char* filename, ParaDomainSelector* domainSelector) +{ + + //distributed meshes //cvwat02 + /*cvw + vector cellglobal; + vector nodeglobal; + vector faceglobal;*/ + + int nbdomain; + + _master_filename=filename; + + // reading ascii master file + try{ + // Setting up the XML tree corresponding to filename + xmlDocPtr master_doc=xmlParseFile(filename); + + if (!master_doc) + throw INTERP_KERNEL::Exception(LOCALIZED("XML Master File does not exist or is not compliant with XML scheme")); + + //////////////////// + //number of domains + //////////////////// + xmlXPathContextPtr xpathCtx = xmlXPathNewContext(master_doc); + xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression(BAD_CAST "//splitting/subdomain", xpathCtx); + if (xpathObj==0 || xpathObj->nodesetval->nodeNr ==0) + throw INTERP_KERNEL::Exception(LOCALIZED("XML Master File does not contain /MED/splitting/subdomain node")); + + /* as subdomain has only one property which is "number" + * it suffices to take the content of its first child */ + const char* mystring = (const char*)xpathObj->nodesetval->nodeTab[0]->properties->children->content; + sscanf(mystring, "%d", &nbdomain); + + ////////////////// + //mesh name + ////////////////// + xmlXPathFreeObject(xpathObj); + xpathObj = xmlXPathEvalExpression(BAD_CAST "//content/mesh", xpathCtx); + if (xpathObj==0 || xpathObj->nodesetval->nodeNr ==0) + throw INTERP_KERNEL::Exception(LOCALIZED("XML Master File does not contain /MED/content/mesh node")); + _collection->setName( (const char*)xpathObj->nodesetval->nodeTab[0]->properties->children->content); + + //cout << "nb domain " << nbdomain << endl; + MyGlobals::_fileNames.resize(nbdomain); + MyGlobals::_meshNames.resize(nbdomain); + (_collection->getMesh()).resize(nbdomain); + (_collection->getFaceMesh()).resize(nbdomain); + (_collection->getCellFamilyIds()).resize(nbdomain); + (_collection->getFaceFamilyIds()).resize(nbdomain); + + /*cvw + cellglobal.resize(nbdomain); + nodeglobal.resize(nbdomain); + faceglobal.resize(nbdomain); + */ + + // retrieving the node which contains the file names + const char filechar[]="//files/subfile"; + xmlXPathFreeObject(xpathObj); + xpathObj = xmlXPathEvalExpression(BAD_CAST filechar, xpathCtx); + if (xpathObj==0 || xpathObj->nodesetval->nodeNr ==0) + throw INTERP_KERNEL::Exception(LOCALIZED("XML Master File does not contain /MED/files/subfile nodes")); + int nbfiles = xpathObj->nodesetval ->nodeNr; + + for (int i=0; inodesetval ==0) + throw INTERP_KERNEL::Exception(LOCALIZED("Error retrieving a file name from subfile of XML Master File")); + MyGlobals::_fileNames[i]=(const char*)xpathObjfilename->nodesetval->nodeTab[0]->children->content; + + //////////////////////////////// + //reading the local mesh names + //////////////////////////////// + ostringstream mesh_search_string; + mesh_search_string<<"//mapping/mesh/chunk[@subdomain=\""<nodesetval ==0) + throw INTERP_KERNEL::Exception(LOCALIZED("Error retrieving mesh name from chunk of XML Master File")); + MyGlobals::_meshNames[i]=(const char*)xpathMeshObj->nodesetval->nodeTab[0]->children->content; + + if ( !domainSelector || domainSelector->isMyDomain(i)) + readSubdomain(i); //cvwat03 + //cvw readSubdomain(cellglobal, faceglobal, nodeglobal, i); //cvwat03 + xmlXPathFreeObject(xpathObjfilename); + xmlXPathFreeObject(xpathMeshObj); + }//loop on domains + + // LIBXML cleanup + xmlXPathFreeObject(xpathObj); + xmlXPathFreeContext(xpathCtx); + xmlFreeDoc(master_doc); + + }//of try + catch(...) + { + throw INTERP_KERNEL::Exception(LOCALIZED("I/O error reading parallel MED file")); + } + + + ParallelTopology* aPT = new ParallelTopology(_collection->getMesh()); + //creation of topology from mesh and connect zones + if ( _collection->isParallelMode() ) + { + //to know nb of cells on each proc to compute global cell ids from locally global + domainSelector->gatherNbOf(_collection->getMesh()); + } + /*cvw + ParallelTopology* aPT = new ParallelTopology + ((_collection->getMesh()), (_collection->getCZ()), cellglobal, nodeglobal, faceglobal); + */ + _collection->setTopology(aPT); + _collection->setDomainNames(_collection->getName()); + /*cvw + for (int i=0; itm_year + ,time_asc->tm_mon+1 + ,time_asc->tm_mday); + + node = xmlNewChild(root_node,0, BAD_CAST "description",0); + + xmlNewProp(node, BAD_CAST "what", BAD_CAST _collection->getDescription().c_str()); + xmlNewProp(node, BAD_CAST "when", BAD_CAST date); + + //Content tag + node =xmlNewChild(root_node,0, BAD_CAST "content",0); + node2 = xmlNewChild(node, 0, BAD_CAST "mesh",0); + xmlNewProp(node2, BAD_CAST "name", BAD_CAST _collection->getName().c_str()); + + //Splitting tag + node=xmlNewChild(root_node,0,BAD_CAST "splitting",0); + node2=xmlNewChild(node,0,BAD_CAST "subdomain",0); + sprintf(buff, "%d", _collection->getMesh().size()); + xmlNewProp(node2, BAD_CAST "number", BAD_CAST buff); + node2=xmlNewChild(node,0,BAD_CAST "global_numbering",0); + xmlNewProp(node2, BAD_CAST "present", BAD_CAST "yes"); + + //Files tag + xmlNodePtr file_node=xmlNewChild(root_node,0,BAD_CAST "files",0); + + //Mapping tag + node = xmlNewChild(root_node,0,BAD_CAST "mapping",0); + xmlNodePtr mesh_node = xmlNewChild(node, 0, BAD_CAST "mesh",0); + xmlNewProp(mesh_node, BAD_CAST "name", BAD_CAST _collection->getName().c_str()); + + int nbdomains= _collection->getMesh().size(); + //vector filenames; + //filenames.resize(nbdomains); + + //loop on the domains + string finalMeshName=extractFromDescription(MyGlobals::_generalInformations[0], "finalMeshName="); + for (int idomain=nbdomains-1; idomain>=0;idomain--) + { + string distfilename; + ostringstream suffix; + suffix<isMyDomain( idomain ) ) + { + if ( (_collection->getMesh())[idomain]->getNumberOfCells()==0 ) continue; //empty domain + if (MyGlobals::_verbose>1) + cout<<"proc "<rank()<<" : writeMedFile "<getMesh())[idomain]->getNumberOfCells()<<" cells" + << " "<<(_collection->getFaceMesh())[idomain]->getNumberOfCells()<<" faces" + << " "<<(_collection->getMesh())[idomain]->getNumberOfNodes()<<" nodes"<rank()==0) + { + //updating the ascii description file + node = xmlNewChild(file_node, 0, BAD_CAST "subfile",0); + sprintf (buff,"%d",idomain+1); + xmlNewProp(node, BAD_CAST "id", BAD_CAST buff); + xmlNewChild(node,0,BAD_CAST "name",BAD_CAST distfilename.c_str()); + xmlNewChild(node,0,BAD_CAST "machine",BAD_CAST "localhost"); + + node = xmlNewChild(mesh_node,0, BAD_CAST "chunk",0); + xmlNewProp(node, BAD_CAST "subdomain", BAD_CAST buff); + xmlNewChild(node,0,BAD_CAST "name", BAD_CAST finalMeshName.c_str()); + //xmlNewChild(node,0,BAD_CAST "name", BAD_CAST (_collection->getMesh())[idomain]->getName()); + } + } + + //create the ascii description file + if (domainSelector->rank()==0) + { + string myfile(filename); + myfile.append(".xml"); + _master_filename=myfile; + if ( !domainSelector || domainSelector->rank() == 0 ) + xmlSaveFormatFileEnc(myfile.c_str(), master_doc, "UTF-8", 1); + //xmlFreeDoc(master_doc); + //xmlCleanupParser(); + } + xmlFreeDoc(master_doc); + xmlCleanupParser(); +} diff --git a/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx new file mode 100644 index 000000000..8a946b60c --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx @@ -0,0 +1,76 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef MESHCOLLECTIONMEDXMLDRIVER_HXX_ +#define MESHCOLLECTIONMEDXMLDRIVER_HXX_ + +#include "MEDPARTITIONER_MESHCollectionDriver.hxx" + +namespace MEDPARTITIONER +{ + class MESHCollection; + + class MESHCollectionMedXMLDriver:public MESHCollectionDriver + { + public: + + MESHCollectionMedXMLDriver(MESHCollection*); + virtual ~MESHCollectionMedXMLDriver(){} + + + int read(const char*, ParaDomainSelector* sel=0); + + void write(const char*, ParaDomainSelector* sel=0); + +// void readFields(vector *>& filenames, char* fieldname, +// int itnumber, int ordernumber) +// { +// _readFields(filenames,fieldname,itnumber,ordernumber); +// } +// void readFields(vector *>& filenames, char* fieldname, +// int itnumber, int ordernumber) +// { +// _readFields(filenames,fieldname,itnumber,ordernumber); +// } + +// void writeFields(vector *>& filenames, char* fieldname) +// { +// _writeFields( filenames, fieldname); +// } +// void writeFields(vector *>& filenames, char* fieldname) +// { +// _writeFields( filenames, fieldname); +// } + + + + private : + + // template void _readFields(vector *>& filenames, char* fieldname, +// int itnumber, int ordernumber); + +// template +// void _writeFields(vector *>& filenames, char* fieldname); + + std::string _master_filename; + }; + +} + + +#endif /*MESHCOLLECTIONDRIVER_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx b/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx new file mode 100644 index 000000000..84731f495 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_METISGraph.cxx @@ -0,0 +1,232 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifdef ENABLE_PARMETIS +#include +#endif +extern "C" { +#include +} +#include "MEDPARTITIONER_METISGraph.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" +#include "MEDPARTITIONER_utils.hxx" +#include "InterpKernelException.hxx" + +#include + +using namespace MEDPARTITIONER; + +METISGraph::METISGraph():Graph() +{ +} + +METISGraph::METISGraph(MEDPARTITIONER::MEDSKYLINEARRAY* graph, int* edgeweight) + :Graph(graph,edgeweight) +{ +} + +METISGraph::~METISGraph() +{ +} + +void METISGraph::partGraph(int ndomain, //cvwat10 + const std::string& options_string, + ParaDomainSelector* parallelizer) +{ + using std::vector; + vector ran,vx,va; //for randomize + + if (MyGlobals::_verbose>10) cout<<"proc "<getNumberOf(); + //graph + int * xadj=const_cast(_graph->getIndex()); + int * adjncy=const_cast(_graph->getValue()); + //constraints + int * vwgt=_cellweight; + int * adjwgt=_edgeweight; + int wgtflag=(_edgeweight!=0)?1:0+(_cellweight!=0)?2:0; + //base 0 or 1 + int base=0; + //ndomain + int nparts=ndomain; + //options + /* + (0=default_option,option,random_seed) see defs.h + #define PMV3_OPTION_DBGLVL 1 + #define PMV3_OPTION_SEED 2 + #define PMV3_OPTION_IPART 3 + #define PMV3_OPTION_PSR 3 + seems no changes int options[4]={1,0,33,0}; //test for a random seed of 33 + */ + int options[4]={0,0,0,0}; + // output parameters + int edgecut; + int* partition=new int[n]; + + //if (MyGlobals::_verbose>10) cout<<"proc "<getProcVtxdist(); + MPI_Comm comm=MPI_COMM_WORLD; + try + { + if (MyGlobals::_verbose>200) + { + cout<<"proc "<0) + { + cout<<"\nproc "<1000) + { + cout<<"\n -cell "<imaxx) imaxx=ilg; + } + cout<<"\nproc "<10) + cout<<"proc "< index(n+1); + vector value(n); + index[0]=0; + if (ran.size()>0 && MyGlobals::_atomize==0) //there is randomize + { + if (MyGlobals::_is0verbose>100) cout<<"randomize"< + +namespace MEDPARTITIONER { + class MEDSKYLINEARRAY; + class MEDPARTITIONER_EXPORT METISGraph:public Graph + { + public: + METISGraph(); + METISGraph(MEDPARTITIONER::MEDSKYLINEARRAY*, int* edgeweight=0); + virtual ~METISGraph(); + void partGraph(int ndomain, const std::string& options_string="", ParaDomainSelector* sel=0); + //private: + // const MEDMEM::MEDSKYLINEARRAY* m_graph; + }; +} +#endif /*METISGRAPH_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.cxx b/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.cxx new file mode 100644 index 000000000..199d545e2 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_ParaDomainSelector.cxx @@ -0,0 +1,700 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : MEDSPLITTER_ParaDomainSelector.cxx +// Created : Wed Jun 24 12:39:59 2009 +// Author : Edward AGAPOV (eap) + +#include "MEDCouplingUMesh.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" +#include "MEDPARTITIONER_UserGraph.hxx" +#include "MEDPARTITIONER_utils.hxx" + +#include +#include + +#ifdef HAVE_MPI2 +#include +#endif + +#ifndef WIN32 +#include +#endif + +//using namespace MED_EN; +using namespace std; +using namespace MEDPARTITIONER; + +//================================================================================ +/*! + * \brief Constructor. Find out my rank and world size + */ +//================================================================================ + +ParaDomainSelector::ParaDomainSelector(bool mesure_memory) + :_rank(0),_world_size(1), _nb_result_domains(-1), _mesure_memory(mesure_memory), + _init_time(0.0), _init_memory(0), _max_memory(0) +{ +#ifdef HAVE_MPI2 + MPI_Comm_size(MPI_COMM_WORLD,&_world_size) ; + MPI_Comm_rank(MPI_COMM_WORLD,&_rank) ; + _init_time = MPI_Wtime(); +#endif + evaluateMemory(); +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +ParaDomainSelector::~ParaDomainSelector() +{ +} + +//================================================================================ +/*! + * \brief Return true if is running on different hosts + */ +//================================================================================ + +bool ParaDomainSelector::isOnDifferentHosts() const +{ + evaluateMemory(); + if ( _world_size < 2 ) return false; + +#ifdef HAVE_MPI2 + char name_here[ MPI_MAX_PROCESSOR_NAME+1 ], name_there[ MPI_MAX_PROCESSOR_NAME+1 ]; + int size; + MPI_Get_processor_name( name_here, &size); + + int next_proc = (rank() + 1) % nbProcs(); + int prev_proc = (rank() - 1 + nbProcs()) % nbProcs(); + int tag = 1111111; + + MPI_Status status; + MPI_Sendrecv((void*)&name_here[0], MPI_MAX_PROCESSOR_NAME, MPI_CHAR, next_proc, tag, + (void*)&name_there[0], MPI_MAX_PROCESSOR_NAME, MPI_CHAR, prev_proc, tag, + MPI_COMM_WORLD, &status); + + //cout<<"proc "<& domain_meshes) +{ + evaluateMemory(); + // get nb of elems of each domain mesh + int nb_domains=domain_meshes.size(); + //cout<<"proc "< nb_elems(nb_domains*2, 0); //NumberOfCells & NumberOfNodes + for (int i=0; igetNumberOfCells(); + nb_elems[i*2+1] = domain_meshes[i]->getNumberOfNodes(); + } + // receive nb of elems from other procs + vector all_nb_elems( nb_domains*2 ); +#ifdef HAVE_MPI2 + MPI_Allreduce((void*)&nb_elems[0], (void*)&all_nb_elems[0], nb_domains*2, + MPI_INT, MPI_SUM, MPI_COMM_WORLD); +#endif + int total_nb_cells=0, total_nb_nodes=0; + for (int i=0; i10) + cout<<"totalNbCells "<& cell_shift_by_domain=_cell_shift_by_domain; + vector& node_shift_by_domain=_node_shift_by_domain; + vector& face_shift_by_domain=_face_shift_by_domain; + + vector< int > ordered_nbs_cell, ordered_nbs_node, domain_order(nb_domains); + ordered_nbs_cell.push_back(0); + ordered_nbs_node.push_back(0); + for (int iproc=0; iproc300) + { + cout<<"proc "<200) + { + cout<<"proc "< ParaDomainSelector::gatherGraph(const Graph* graph) const +{ + Graph* glob_graph = 0; + + evaluateMemory(); +#ifdef HAVE_MPI2 + + // --------------- + // Gather indices + // --------------- + + vector index_size_of_proc( nbProcs() ); // index sizes - 1 + for ( int i = 1; i < _nb_vert_of_procs.size(); ++i ) + index_size_of_proc[i-1] = _nb_vert_of_procs[ i ] - _nb_vert_of_procs[ i-1 ]; + + int index_size = 1 + _cell_shift_by_domain.back(); + int* graph_index = new int[ index_size ]; + const int* index = graph->getGraph()->getIndex(); + int* proc_index_displacement = (int*) & _nb_vert_of_procs[0]; + + MPI_Allgatherv((void*) (index+1), // send local index except first 0 (or 1) + index_size_of_proc[_rank], // index size on this proc + MPI_INT, + (void*) graph_index, // receive indices + & index_size_of_proc[0], // index size on each proc + proc_index_displacement, // displacement of each proc index + MPI_INT, + MPI_COMM_WORLD); + graph_index[0] = index[0]; // it is not overwritten thanks to proc_index_displacement[0]==1 + + // get sizes of graph values on each proc by the got indices of graphs + vector< int > value_size_of_proc( nbProcs() ), proc_value_displacement(1,0); + for ( int i = 0; i < nbProcs(); ++i ) + { + if ( index_size_of_proc[i] > 0 ) + value_size_of_proc[i] = graph_index[ proc_index_displacement[ i+1 ]-1 ] - graph_index[0]; + else + value_size_of_proc[i] = 0; + proc_value_displacement.push_back( proc_value_displacement.back() + value_size_of_proc[i] ); + } + + // update graph_index + for ( int i = 1; i < nbProcs(); ++i ) + { + int shift = graph_index[ proc_index_displacement[i]-1 ]-graph_index[0]; + for ( int j = proc_index_displacement[i]; j < proc_index_displacement[i+1]; ++j ) + graph_index[ j ] += shift; + } + + // -------------- + // Gather values + // -------------- + + int value_size = graph_index[ index_size-1 ] - graph_index[ 0 ]; + int* graph_value = new int[ value_size ]; + const int* value = graph->getGraph()->getValue(); + + MPI_Allgatherv((void*) value, // send local value + value_size_of_proc[_rank], // value size on this proc + MPI_INT, + (void*) graph_value, // receive values + & value_size_of_proc[0], // value size on each proc + & proc_value_displacement[0], // displacement of each proc value + MPI_INT, + MPI_COMM_WORLD); + + // ----------------- + // Gather partition + // ----------------- + + int * partition = new int[ _cell_shift_by_domain.back() ]; + const int* part = graph->getPart(); + + MPI_Allgatherv((void*) part, // send local partition + index_size_of_proc[_rank], // index size on this proc + MPI_INT, + (void*)(partition-1), // -1 compensates proc_index_displacement[0]==1 + & index_size_of_proc[0], // index size on each proc + proc_index_displacement, // displacement of each proc index + MPI_INT, + MPI_COMM_WORLD); + + // ----------- + // Make graph + // ----------- + +// MEDPARTITIONER::MEDSKYLINEARRAY* array = +// new MEDPARTITIONER::MEDSKYLINEARRAY( index_size-1, value_size, graph_index, graph_value, true ); + +// glob_graph = new UserGraph( array, partition, index_size-1 ); + + evaluateMemory(); + + delete [] partition; + +#endif // HAVE_MPI2 + + return auto_ptr( glob_graph ); +} + + +//================================================================================ +/*! + * \brief Set nb of cell/cell pairs in a joint between domains + */ +//================================================================================ + +void ParaDomainSelector::setNbCellPairs( int nb_cell_pairs, int dist_domain, int loc_domain ) +{ + // This method is needed for further computing global numbers of faces in joint. + // Store if both domains are on this proc else on one of procs only + if ( isMyDomain( dist_domain ) || dist_domain < loc_domain ) + { + if ( _nb_cell_pairs_by_joint.empty() ) + _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0); + + int joint_id = jointId( loc_domain, dist_domain ); + _nb_cell_pairs_by_joint[ joint_id ] = nb_cell_pairs; + } + evaluateMemory(); +} + +//================================================================================ +/*! + * \brief Return nb of cell/cell pairs in a joint between domains on different procs + */ +//================================================================================ + +int ParaDomainSelector::getNbCellPairs( int dist_domain, int loc_domain ) const +{ + evaluateMemory(); + + int joint_id = jointId( loc_domain, dist_domain ); + return _nb_cell_pairs_by_joint[ joint_id ]; +} + +//================================================================================ +/*! + * \brief Gather size of each joint + */ +//================================================================================ + +void ParaDomainSelector::gatherNbCellPairs() +{ + if ( _nb_cell_pairs_by_joint.empty() ) + _nb_cell_pairs_by_joint.resize( _nb_result_domains*(_nb_result_domains+1), 0); + evaluateMemory(); + + vector< int > send_buf = _nb_cell_pairs_by_joint; +#ifdef HAVE_MPI2 + MPI_Allreduce((void*)&send_buf[0], + (void*)&_nb_cell_pairs_by_joint[0], + _nb_cell_pairs_by_joint.size(), + MPI_INT, MPI_SUM, MPI_COMM_WORLD); +#endif + // check that the set nbs of cell pairs are correct, + // namely that each joint is treated on one proc only + for ( int j = 0; j < _nb_cell_pairs_by_joint.size(); ++j ) + if ( _nb_cell_pairs_by_joint[j] != send_buf[j] && send_buf[j]>0 ) + throw INTERP_KERNEL::Exception(LOCALIZED("invalid nb of cell pairs")); +} + +//================================================================================ +/*! + * \brief Return the first global id of sub-entity for the joint + */ +//================================================================================ + +int ParaDomainSelector::getFisrtGlobalIdOfSubentity( int loc_domain, int dist_domain ) const +{ + // total_nb_faces includes faces existing before creation of joint faces + // (got in gatherNbOf( MED_FACE )). + evaluateMemory(); + + int total_nb_faces = _face_shift_by_domain.empty() ? 0 : _face_shift_by_domain.back(); + int id = total_nb_faces + 1; + + if ( _nb_cell_pairs_by_joint.empty() ) + throw INTERP_KERNEL::Exception(LOCALIZED("gatherNbCellPairs() must be called before")); + int joint_id = jointId( loc_domain, dist_domain ); + for ( int j = 0; j < joint_id; ++j ) + id += _nb_cell_pairs_by_joint[ j ]; + + return id; +} + +//================================================================================ +/*! + * \brief Send-receive local ids of joint faces + */ +//================================================================================ + +int* ParaDomainSelector::exchangeSubentityIds( int loc_domain, int dist_domain, + const vector& loc_ids_here ) const +{ + int* loc_ids_dist = new int[ loc_ids_here.size()]; + int dest = getProcessorID( dist_domain ); + int tag = 2002 + jointId( loc_domain, dist_domain ); +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Sendrecv((void*)&loc_ids_here[0], loc_ids_here.size(), MPI_INT, dest, tag, + (void*) loc_ids_dist, loc_ids_here.size(), MPI_INT, dest, tag, + MPI_COMM_WORLD, &status); +#endif + evaluateMemory(); + + return loc_ids_dist; +} + +//================================================================================ +/*! + * \brief Return identifier for a joint + */ +//================================================================================ + +int ParaDomainSelector::jointId( int local_domain, int distant_domain ) const +{ + evaluateMemory(); + if (_nb_result_domains < 0) + throw INTERP_KERNEL::Exception(LOCALIZED("setNbDomains() must be called before")); + + if ( local_domain < distant_domain ) + swap( local_domain, distant_domain ); + return local_domain * _nb_result_domains + distant_domain; +} + + +//================================================================================ +/*! + * \brief Return time passed from construction in seconds + */ +//================================================================================ + +double ParaDomainSelector::getPassedTime() const +{ +#ifdef HAVE_MPI2 + return MPI_Wtime() - _init_time; +#else + return 0.0; +#endif +} + +//================================================================================ +/*! + * \brief Evaluate current memory usage and return the maximal one in KB + */ +//================================================================================ + +int ParaDomainSelector::evaluateMemory() const +{ + if ( _mesure_memory ) + { + int used_memory = 0; +#ifndef WIN32 + struct sysinfo si; + int err = sysinfo( &si ); + if ( !err ) + used_memory = + (( si.totalram - si.freeram + si.totalswap - si.freeswap ) * si.mem_unit ) / 1024; +#endif + if ( used_memory > _max_memory ) + ((ParaDomainSelector*) this)->_max_memory = used_memory; + + if ( !_init_memory ) + ((ParaDomainSelector*) this)->_init_memory = used_memory; + } + return _max_memory - _init_memory; +} + +/*! +Sends content of \a mesh to processor \a target. To be used with \a recvMesh method. +\param mesh mesh to be sent +\param target processor id of the target +*/ + +void ParaDomainSelector::sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const +{ + if (MyGlobals::_verbose>600) cout<<"proc "<<_rank<<" : sendMesh '"< tinyInfoLocal; + vector tinyInfoLocalS; + vector tinyInfoLocalD; + //Getting tiny info of local mesh to allow the distant proc to initialize and allocate + //the transmitted mesh. + mesh.getTinySerializationInformation(tinyInfoLocalD,tinyInfoLocal,tinyInfoLocalS); + //cout<<"sendMesh getTinySerializationInformation "<0) //no sends if empty + { + ParaMEDMEM::DataArrayInt *v1Local=0; + ParaMEDMEM::DataArrayDouble *v2Local=0; + //serialization of local mesh to send data to distant proc. + mesh.serialize(v1Local,v2Local); + int nbLocalElems=0; + int* ptLocal=0; + if(v1Local) //cvw if empty getNbOfElems() is 1! + { + nbLocalElems=v1Local->getNbOfElems(); //cvw if empty be 1! + ptLocal=v1Local->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Send(ptLocal, nbLocalElems, MPI_INT, target, 1111, MPI_COMM_WORLD); +#endif + int nbLocalElems2=0; + double *ptLocal2=0; + if(v2Local) //cvw if empty be 0! + { + nbLocalElems2=v2Local->getNbOfElems(); + ptLocal2=v2Local->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Send(ptLocal2, nbLocalElems2, MPI_DOUBLE, target, 1110, MPI_COMM_WORLD); +#endif + if(v1Local) v1Local->decrRef(); + if(v2Local) v2Local->decrRef(); + } + else + { + //cout<<"sendMesh empty Mesh cvw3344 "< tinyInfoDistant; + vector tinyInfoLocalS; + vector tinyInfoDistantD(1); + //Getting tiny info of local mesh to allow the distant proc to initialize and allocate + //the transmitted mesh. +#ifdef HAVE_MPI2 + MPI_Status status; + int tinyVecSize; + MPI_Recv(&tinyVecSize, 1, MPI_INT, source, 1113, MPI_COMM_WORLD, &status); + tinyInfoDistant.resize(tinyVecSize); +#endif + std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0); + +#ifdef HAVE_MPI2 + MPI_Recv(&tinyInfoDistant[0], tinyVecSize, MPI_INT,source,1112,MPI_COMM_WORLD, &status); +#endif + //there was tinyInfoLocal.push_back(mesh.getNumberOfCells()); + int NumberOfCells=tinyInfoDistant[tinyVecSize-1]; + //cout<<"recvMesh NumberOfCells "<0) + { + ParaMEDMEM::DataArrayInt *v1Distant=ParaMEDMEM::DataArrayInt::New(); + ParaMEDMEM::DataArrayDouble *v2Distant=ParaMEDMEM::DataArrayDouble::New(); + //Building the right instance of copy of distant mesh. + ParaMEDMEM::MEDCouplingPointSet *distant_mesh_tmp= + ParaMEDMEM::MEDCouplingPointSet::BuildInstanceFromMeshType( + (ParaMEDMEM::MEDCouplingMeshType) tinyInfoDistant[0]); + std::vector unusedTinyDistantSts; + mesh=dynamic_cast (distant_mesh_tmp); + + mesh->resizeForUnserialization(tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); + int nbDistElem=0; + int *ptDist=0; + if(v1Distant) + { + nbDistElem=v1Distant->getNbOfElems(); + ptDist=v1Distant->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Recv(ptDist, nbDistElem, MPI_INT, source,1111, MPI_COMM_WORLD, &status); +#endif + double *ptDist2=0; + nbDistElem=0; + if(v2Distant) + { + nbDistElem=v2Distant->getNbOfElems(); + ptDist2=v2Distant->getPointer(); + } +#ifdef HAVE_MPI2 + MPI_Recv(ptDist2, nbDistElem, MPI_DOUBLE,source, 1110, MPI_COMM_WORLD, &status); +#endif + //finish unserialization + mesh->unserialization(tinyInfoDistantD,tinyInfoDistant,v1Distant,v2Distant,unusedTinyDistantSts); + if(v1Distant) v1Distant->decrRef(); + if(v2Distant) v2Distant->decrRef(); + } + else + { + mesh=createEmptyMEDCouplingUMesh(); + } + if (MyGlobals::_verbose>600) cout<<"proc "<<_rank<<" : recvMesh '"<getName()<<"' size "<getNumberOfCells()<<" from "< +#include + +namespace ParaMEDMEM +{ + class MEDCouplingUMesh; +} + + +namespace MEDPARTITIONER +{ + class Graph; + class JointExchangeData; + +/*! + * \brief Communication helper in parallel mode + */ +class MEDPARTITIONER_EXPORT ParaDomainSelector +{ + +public: + + ParaDomainSelector(bool mesure_memory=false); + ~ParaDomainSelector(); + + // return processor rank + int rank() const { return _rank; } + // return number of processors + int nbProcs() const { return _world_size; } + // Return true if is running on different hosts + bool isOnDifferentHosts() const; + // Return true if the domain with domainIndex is to be loaded on this proc + bool isMyDomain(int domainIndex) const; + // Return processor id where the domain with domainIndex resides + int getProcessorID(int domainIndex) const; + //Set nb of required domains. (Used to sort joints via jointId()) + void setNbDomains(int nb) { _nb_result_domains = nb; } + // Return identifier for a joint + int jointId( int local_domain, int distant_domain ) const; + + int getNbTotalCells() { return _cell_shift_by_domain.back(); } + int getNbTotalNodes() { return _node_shift_by_domain.back(); }; + int getNbTotalFaces() { return _face_shift_by_domain.back(); }; + + // Return domain order + //int getDomainOrder(int domainIndex, int nb_domains) const; + + // Collect nb of entities on procs + void gatherNbOf(const std::vector& domain_meshes); + + // Return distribution of the graph vertices among the processors + int* getProcVtxdist() const; + + // Return nb of nodes on processors with lower rank + int getProcNodeShift() const; + // Return nb of cells in domains with lower index + int getDomainCellShift(int domainIndex) const; + // Return nb of nodes in domains with lower index + int getDomainNodeShift(int domainIndex) const; + + // Gather graphs from all processors into one + std::auto_ptr gatherGraph(const Graph* graph) const; + + // Set nb of cell/cell pairs in a joint between domains + void setNbCellPairs( int nb_cell_pairs, int dist_domain, int loc_domain ); + // Gather size of each proc/proc joint + void gatherNbCellPairs(); + // Return nb of cell/cell pairs in a joint between domains on different procs + int getNbCellPairs( int dist_domain, int loc_domain ) const; + + // Send-receive joint data + // void exchangeJoint( JointExchangeData* joint ) const; + + // Return the first global id of sub-entity for the joint + int getFisrtGlobalIdOfSubentity( int loc_domain, int dist_domain ) const; + // Send-receive local ids of joint faces + int* exchangeSubentityIds( int loc_domain, int dist_domain, + const std::vector& loc_ids_here ) const; + // Return time passed from construction in seconds + double getPassedTime() const; + + // Evaluate current memory usage and return the maximal one in KB + int evaluateMemory() const; + + void sendMesh(const ParaMEDMEM::MEDCouplingUMesh& mesh, int target) const; + void recvMesh(ParaMEDMEM::MEDCouplingUMesh*& mesh, int source) const; + +private: + + int _rank, _world_size; // my rank and nb of processors + int _nb_result_domains; // required nb of domains + + //int _total_nb_faces; // nb of faces in the whole mesh without proc/proc joint faces + + std::vector< int > _nb_cell_pairs_by_joint; + std::vector< int > _nb_vert_of_procs; // graph vertices + std::vector< int > _cell_shift_by_domain; + std::vector< int > _node_shift_by_domain; + std::vector< int > _face_shift_by_domain; + + double _init_time; + bool _mesure_memory; + int _init_memory, _max_memory; +}; +} + +#endif diff --git a/src/MEDPartitioner/MEDPARTITIONER_ParallelTopology.cxx b/src/MEDPartitioner/MEDPARTITIONER_ParallelTopology.cxx new file mode 100644 index 000000000..42055958f --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_ParallelTopology.cxx @@ -0,0 +1,574 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#include +#include +#include +#include + +#include "InterpKernelHashMap.hxx" +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_Topology.hxx" +#include "MEDPARTITIONER_Graph.hxx" +#include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_ConnectZone.hxx" + +#include "MEDCouplingUMesh.hxx" +#include "MEDPARTITIONER_utils.hxx" + +#ifdef HAVE_MPI2 +#include +#endif + +using namespace MEDPARTITIONER; +using namespace std; + +//empty constructor +ParallelTopology::ParallelTopology():_nb_domain(0),_mesh_dimension(0) +{} + +//!constructing topology according to mesh collection without global numerotation (use setGlobalNumerotation later) +ParallelTopology::ParallelTopology(const vector& meshes) +{ + _nb_domain=meshes.size(); + _nb_cells.resize(_nb_domain); + _nb_nodes.resize(_nb_domain); + // _nb_faces.resize(_nb_domain); + + if (MyGlobals::_is0verbose>100) cout<<"new ParallelTopology\n"; + _loc_to_glob.resize(0); //precaution, need gatherNbOf() setGlobalNumerotation() + _node_loc_to_glob.resize(0); //precaution, need gatherNbOf() setGlobalNumerotation() + //_face_loc_to_glob.resize(_nb_domain); + _mesh_dimension = -1; + bool parallel_mode = false; + for (int idomain=0; !parallel_mode && idomain<_nb_domain; idomain++) + parallel_mode = (!meshes[idomain]); + + if (MyGlobals::_is0verbose>20 && !parallel_mode) cout<<"WARNING : ParallelTopology contructor without parallel_mode"<getMeshDimension(); + } + else + { + if (_mesh_dimension!=meshes[idomain]->getMeshDimension()) + throw INTERP_KERNEL::Exception(LOCALIZED("meshes dimensions incompatible")); + } + _nb_cells[idomain]=meshes[idomain]->getNumberOfCells(); + _nb_nodes[idomain]=meshes[idomain]->getNumberOfNodes(); + //note: in parallel mode _nb_cells and _nb_nodes are not complete now, needs gatherNbOf() + } +} + +//!constructing _loc_to_glob etc by default, needs gatherNbOf() done +void ParallelTopology::setGlobalNumerotationDefault(ParaDomainSelector* domainSelector) +{ + if (MyGlobals::_is0verbose>100) cout<<"setGlobalNumerotationDefault on "<<_nb_domain<<" domains\n"; + if (_loc_to_glob.size()!=0) throw INTERP_KERNEL::Exception(LOCALIZED("a global numerotation is done yet")); + _loc_to_glob.resize(_nb_domain); + _node_loc_to_glob.resize(_nb_domain); + + //warning because _nb_cells[idomain] is 0 if not my domain(s) + //we set loc_to_glob etc.. only for my domain(s) + if (MyGlobals::_is0verbose>500) cout<<"(c)idomain|ilocalCell|iglobalCell"<getDomainCellShift(idomain); + for (int i=0; i<_nb_cells[idomain]; i++) + { + int global=domainCellShift+i ; + _glob_to_loc.insert(make_pair(global,make_pair(idomain,i))); + _loc_to_glob[idomain][i]=global; + if (MyGlobals::_verbose>500) cout<<"c"<500) MPI_Barrier(MPI_COMM_WORLD); + if (MyGlobals::_is0verbose>500) cout<500) cout<<"(n)idomain|ilocalNode|iglobalNode"<getDomainNodeShift(idomain); + for (int i=0; i<_nb_nodes[idomain]; i++) + { + int global=domainNodeShift+i ; + _node_glob_to_loc.insert(make_pair(global,make_pair(idomain,i))); + _node_loc_to_glob[idomain][i]=global; + if (MyGlobals::_verbose>500) cout<<"n"<500) MPI_Barrier(MPI_COMM_WORLD); + if (MyGlobals::_is0verbose>500) cout<getNbTotalCells(); + _nb_total_nodes=domainSelector->getNbTotalNodes(); + _nb_total_faces=domainSelector->getNbTotalFaces(); + if (MyGlobals::_is0verbose>200) cout<<"globalNumerotation default done meshDimension "<<_mesh_dimension<<" nbTotalCells "<<_nb_total_cells<<" nbTotalNodes "<<_nb_total_nodes<& meshes, + const vector& cz, + vector& cellglobal, + vector& nodeglobal, + vector& faceglobal) +{ + _nb_domain=meshes.size(); + int index_global=0; + int index_node_global=0; + int index_face_global=0; + + _nb_cells.resize(_nb_domain); + _nb_nodes.resize(_nb_domain); + // _nb_faces.resize(_nb_domain); + + _loc_to_glob.resize(_nb_domain); + _node_loc_to_glob.resize(_nb_domain); + // _face_loc_to_glob.resize(_nb_domain); + + bool parallel_mode = false; + for (int idomain=0; !parallel_mode && idomain<_nb_domain; idomain++) + parallel_mode = (!meshes[idomain]); + + for (int idomain=0; idomain<_nb_domain; idomain++) + { + if ( !meshes[idomain] ) continue; + _mesh_dimension = meshes[idomain]->getMeshDimension(); + + //creating cell maps + _nb_cells[idomain]=meshes[idomain]->getNumberOfCells(); + // cout << "Nb cells (domain "< ("<getNumberOfNodes()); + for (int i=0; igetNumberOfNodes(); i++) + { + _node_glob_to_loc.insert(make_pair(i,make_pair(0,i))); + _node_loc_to_glob[0][i]=i; + } + _nb_total_nodes=meshes[idomain]->getNumberOfNodes(); + _nb_nodes[0]=_nb_total_nodes; + return; + } + + //creating node maps + _nb_nodes[idomain]=meshes[idomain]->getNumberOfNodes(); + INTERP_KERNEL::HashMap > local2distant; + _node_loc_to_glob[idomain].resize(_nb_nodes[idomain]); + for (int icz=0; iczgetLocalDomainNumber() == idomain && + cz[icz]->getLocalDomainNumber()>cz[icz]->getDistantDomainNumber()) + { + int nb_node= cz[icz]->getNodeNumber(); + const int* node_corresp=cz[icz]->getNodeCorrespValue(); + int distant_ip = cz[icz]->getDistantDomainNumber(); + for (int i=0; i< nb_node; i++) + { + int local= node_corresp[i*2]; + int distant = node_corresp[i*2+1]; + local2distant.insert(make_pair(local, make_pair(distant_ip,distant))); + } + } + } + // setting mappings for all nodes + if (nodeglobal[idomain]==0) + { + for (int inode=0; inode<_nb_nodes[idomain]; inode++) + { + if (local2distant.find(inode)==local2distant.end()) + { + index_node_global++; + _node_glob_to_loc.insert(make_pair(index_node_global,make_pair(idomain,inode))); + //_node_loc_to_glob[make_pair(idomain,inode+1)]=index_node_global; + _node_loc_to_glob[idomain][inode]=index_node_global; + } + else + { + int ip = (local2distant.find(inode)->second).first; + int distant = (local2distant.find(inode)->second).second; + int global_number=_loc_to_glob[ip][distant]; + _node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode))); + _node_loc_to_glob[idomain][inode]=global_number; + } + } + } + //using former node numbering + else + { + for (int inode=0; inode<_nb_nodes[idomain]; inode++) + { + int global_number=nodeglobal[idomain][inode]; + _node_glob_to_loc.insert(make_pair(global_number,make_pair(idomain,inode))); + _node_loc_to_glob[idomain][inode]=global_number; + } + } + } + + _nb_total_cells=index_global; + _nb_total_nodes=index_node_global; + _nb_total_faces=index_face_global; +} + + +//!constructing ParallelTopology from an old topology and a graph +ParallelTopology::ParallelTopology(Graph* graph, Topology* oldTopology, int nb_domain, int mesh_dimension) +{ + + _nb_domain=nb_domain; + //cvw !!whatisit! _nb_cells=graph->nbVertices(); + _mesh_dimension=mesh_dimension; + + if (MyGlobals::_verbose>200) + cout<<"proc "<nbDomain()<<" newNbDomain "<<_nb_domain<getPart(); //all cells for this proc (may be more domains) + _nb_total_cells=graph->nbVertices(); //all cells for this proc (may be more domains) + if (MyGlobals::_verbose>300) + cout<<"proc "<nbDomain(); iold++) + { + int ioldNbCell=oldTopology->getCellNumber(iold); + //cout<<"proc "< globalids(ioldNbCell); + oldTopology->getCellList(iold, &globalids[0]); //unique global numerotation + for (int icell=0; icell300) + for (int idomain=0; idomain<_nb_domain; idomain++) + cout<<"proc "<nbVertices(); + _mesh_dimension=mesh_dimension; + + cout<<"proc "<nbDomain()<<" new nbDomain "<<_nb_domain<getPart(); + _nb_total_cells=graph->nbVertices(); + + cout<<"proc "< local_node = _node_glob_to_loc.find(node_list[i])->second; + ip[i]=local_node.first; + local[i]=local_node.second; + } +} + +/*!Converts a list of global node numbers on domain ip + * to a distributed array with local cell numbers. + * + * If a node in the list is represented on several domains, + * only the value with domain ip is returned + * + * */ +void ParallelTopology::convertGlobalNodeList(const int* node_list, int nbnode, int* local, int ip) +{ + if (_node_glob_to_loc.empty()) + throw INTERP_KERNEL::Exception(LOCALIZED("Node mapping has not yet been built")); + + for (int i=0; i< nbnode; i++) + { + typedef INTERP_KERNEL::HashMultiMap >::iterator mmiter; + pair range=_node_glob_to_loc.equal_range(node_list[i]); + for (mmiter it=range.first; it !=range.second; it++) + { + int ipfound=(it->second).first; + if (ipfound==ip) + local[i]=(it->second).second; + } + } +} + +/*!Converts a list of global node numbers + * to a distributed array with local cell numbers. + * + * If a node in the list is represented on several domains, + * all the values are put in the array + * */ +void ParallelTopology::convertGlobalNodeListWithTwins(const int* node_list, int nbnode, int*& local, int*& ip,int*& full_array, int& size) +{ + if (_node_glob_to_loc.empty()) + throw INTERP_KERNEL::Exception(LOCALIZED("Node mapping has not yet been built")); + + size=0; + for (int i=0; i< nbnode; i++) + { + int count= _node_glob_to_loc.count(node_list[i]); + // if (count > 1) + // cout << "noeud " << node_list[i]<< " doublon d'ordre " << count< >::iterator mmiter; + pair range=_node_glob_to_loc.equal_range(node_list[i]); + for (mmiter it=range.first; it !=range.second; it++) + { + ip[index]=(it->second).first; + local[index]=(it->second).second; + full_array [index]=node_list[i]; + index++; + } + + } +} + +/*!Converts a list of global face numbers + * to a distributed array with local face numbers. + * + * If a face in the list is represented on several domains, + * all the values are put in the array + * */ +void ParallelTopology::convertGlobalFaceListWithTwins(const int* face_list, int nbface, int*& local, int*& ip, int*& full_array,int& size) +{ + size=0; + for (int i=0; i< nbface; i++) + { + //int count = _face_glob_to_loc.count(face_list[i]); + //if (count >1) MESSAGE_MED("face en doublon "< >::iterator mmiter; + pair range=_face_glob_to_loc.equal_range(face_list[i]); + for (mmiter it=range.first; it !=range.second; it++) + { + ip[index]=(it->second).first; + local[index]=(it->second).second; + full_array[index]=face_list[i]; + index++; + } + + } +} + +//!converts a list of global cell numbers +//!to a distributed array with local cell numbers +void ParallelTopology::convertGlobalCellList(const int* cell_list, int nbcell, int* local, int* ip) +{ + for (int i=0; i >::const_iterator iter = _glob_to_loc.find(cell_list[i]); + INTERP_KERNEL::HashMap >::const_iterator iter = _glob_to_loc.find(cell_list[i]); + if (iter == _glob_to_loc.end()) + { + cerr<<"proc "<second).first; //no domain + local[i]=(iter->second).second; //no local cell + //cout<<"proc "< >::const_iterator iter = _face_glob_to_loc.find(face_list[i]); + if (iter == _face_glob_to_loc.end()) + { + throw INTERP_KERNEL::Exception(LOCALIZED("ParallelTopology::convertGlobalFaceList : Face not found")); + } + ip[i]=(iter->second).first; + local[i]=(iter->second).second; + // cout << " in convertGlobalFAceList face global "< ("< >::iterator mmiter; + pair range=_face_glob_to_loc.equal_range(face_list[i]); + for (mmiter it=range.first; it !=range.second; it++) + { + int ipfound=(it->second).first; + if (ipfound==ip) + local[i]=(it->second).second; + + } + } +} + +//replacing a table of global numbering with a table with local numberings +// type_connectivity contains global connectivity for each type in input +// type_connectivity contains local connectivity for each type in output +void ParallelTopology::convertToLocal2ndVersion(int* nodes, int nbnodes, int idomain) +{ + for (int inode=0; inode struct hash< std::pair > + { + size_t operator()( const std::pair& x ) const + { + return hash< int >()( x.first*1000000+x.second ); + } + }; +}*/ + +namespace MEDPARTITIONER { + + class Graph; + class MESHCollection; + class MEDPARTITIONER_FaceModel; + + class ParallelTopology:public Topology + { + + public: + + ParallelTopology(); + + ParallelTopology(const std::vector&); + ParallelTopology(const std::vector&, + const std::vector&, + std::vector&, + std::vector&, + std::vector&); + + ParallelTopology(Graph* graph, Topology* oldTopology, int nbdomain, int mesh_dimension); + + ~ParallelTopology(); + + void setGlobalNumerotationDefault(ParaDomainSelector* domainSelector); + + //!converts a list of global cell numbers + //!to a distributed array with local cell numbers + void convertGlobalNodeList(const int*, int,int*,int*); + void convertGlobalNodeList(const int*, int,int*,int); + void convertGlobalNodeListWithTwins(const int* face_list, int nbnode, int*& local, int*& ip, int*& full_array, int& size); + + //!converts a list of global node numbers + //!to a distributed array with local cell numbers + void convertGlobalCellList(const int*, int , int*, int *); + + //!converts a list of global face numbers + //!to a distributed array with local face numbers + void convertGlobalFaceList(const int*, int , int*, int *); + void convertGlobalFaceList(const int*, int , int*, int); + void convertGlobalFaceListWithTwins(const int* face_list, int nbface, int*& local, int*& ip, int*& full_array,int& size); + +// //!converting node global numberings to local numberings + void convertToLocal2ndVersion(int* nodes, int nbnodes, int idomain); + + //!converting node local numbering to global + inline int convertNodeToGlobal(int ip,int icell) const + { + //return _node_loc_to_glob.find(make_pair(ip,icell))->second; + return _node_loc_to_glob[ip][icell]; + } + + //!converting face local numbering to global + inline int convertFaceToGlobal(int ip,int iface) const + { + // if (_face_loc_to_glob.find(make_pair(ip,icell))==_face_loc_to_glob.end()) + // return -1; + // else + //return _face_loc_to_glob.find(make_pair(ip,icell))->second; + return _face_loc_to_glob[ip][iface]; + } + + //converting cell global numbering to local + inline int convertCellToGlobal(int ip,int icell) const + { + // if (_loc_to_glob.find(make_pair(ip,icell))==_loc_to_glob.end()) + // return -1; + // else + //return _loc_to_glob.find(make_pair(ip,icell))->second; + return _loc_to_glob[ip][icell]; + } + + inline void convertNodeToGlobal(int ip, const int* local, int n, int* global)const + { + for (int i=0; i keys; + for (INTERP_KERNEL::HashMultiMap >::const_iterator iter= _node_glob_to_loc.begin(); + iter!=_node_glob_to_loc.end(); + iter++) { + keys.insert(iter->first); + } + return keys.size(); + } + + //!retrieving list of nodes in global numbers + inline void getNodeList(int idomain, int* list) const + { + for (int i=0; i<_nb_nodes[idomain]; i++) + list[i]=_node_loc_to_glob[idomain][i]; + } + + //!< retrieving cell numbers after fusing in parallel mode + std::vector & getFusedCellNumbers(int idomain) + { + return _cell_loc_to_glob_fuse[idomain]; + } + + const std::vector & getFusedCellNumbers(int idomain) const + { + return _cell_loc_to_glob_fuse[idomain]; + } + + //!< retrieving face numbers after fusing in parallel mode + std::vector & getFusedFaceNumbers(int idomain) + { + return _face_loc_to_glob_fuse[idomain]; + } + const std::vector & getFusedFaceNumbers(int idomain) const + { + return _face_loc_to_glob_fuse[idomain]; + } + + + //!retrieving number of nodes + inline int getCellNumber(int idomain) const + { + return _nb_cells[idomain]; + } + + inline int getCellDomainNumber(int global) const + { + return (_glob_to_loc.find(global)->second).first; + } + + //!retrieving list of nodes in global numbers + inline void getCellList(int idomain, int* list) const + { + for (int i=0; i<_nb_cells[idomain];i++) + list[i]=_loc_to_glob[idomain][i]; + } + + inline int getFaceNumber(int idomain) const + { + return _nb_faces[idomain]; + } + + inline int getFaceNumber() const + { + if (_face_glob_to_loc.empty()) return 0; + std::set keys; + for (INTERP_KERNEL::HashMultiMap >::const_iterator iter= _face_glob_to_loc.begin(); + iter!=_face_glob_to_loc.end(); + iter++) { + keys.insert(iter->first); + } + return keys.size(); + } + + + //!retrieving list of faces in global numbers + inline void getFaceList(int idomain, int* list) const + { + for (int i=0; i<_nb_faces[idomain];i++) + list[i]=_face_loc_to_glob[idomain][i]; + } + + //! converting a global cell number to a local representation (domain + local number) + inline std::pair convertGlobalCell(int iglobal) const + { + return _glob_to_loc.find(iglobal)->second; + } + + inline int convertGlobalFace(int iglobal, int idomain) + { + typedef INTERP_KERNEL::HashMultiMap >::const_iterator MMiter; + std::pair eq = _face_glob_to_loc.equal_range(iglobal); + for (MMiter it=eq.first; it != eq.second; it++) + if (it->second.first == idomain) return it->second.second; + return -1; + } + + inline int convertGlobalNode(int iglobal, int idomain) + { + typedef INTERP_KERNEL::HashMultiMap >::const_iterator MMiter; + std::pair eq = _node_glob_to_loc.equal_range(iglobal); + for (MMiter it=eq.first; it != eq.second; it++) + { + if (it->second.first == idomain) return it->second.second; + } + return -1; + } + + //!adding a face to the topology + inline void appendFace(int idomain, int ilocal, int iglobal) + { + _face_loc_to_glob[idomain].push_back(iglobal); + _face_glob_to_loc.insert(std::make_pair(iglobal,std::make_pair(idomain,ilocal))); + } + + //return max global face number + int getMaxGlobalFace() const; + + private: + bool hasCellWithNodes( const MESHCollection&, int dom, const std::set& nodes ); + + private: + //!mapping global -> local + typedef INTERP_KERNEL::HashMultiMap > TGlob2DomainLoc; + + TGlob2DomainLoc _glob_to_loc; + + std::vector > _loc_to_glob; + + INTERP_KERNEL::HashMultiMap > _node_glob_to_loc; + + //!mapping local -> global + std::vector > _node_loc_to_glob; + + // global numbers in parallel mode + std::vector > _cell_loc_to_glob_fuse; // glob nums after fusing + std::vector > _face_loc_to_glob_fuse; // glob nums after fusing + + //!mapping global -> local + typedef INTERP_KERNEL::HashMultiMap > TGlob2LocsMap; + TGlob2LocsMap _face_glob_to_loc; + + //!mapping local -> global + std::vector > _face_loc_to_glob; + std::vector _nb_cells; + std::vector _nb_nodes; + std::vector _nb_faces; + int _nb_total_cells; + int _nb_total_nodes; + int _nb_total_faces; + int _nb_domain; + int _mesh_dimension; + }; + + +} +#endif /*PARALLELTOPOLOGY_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_SCOTCHGraph.cxx b/src/MEDPartitioner/MEDPARTITIONER_SCOTCHGraph.cxx new file mode 100644 index 000000000..eba5ef4a1 --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_SCOTCHGraph.cxx @@ -0,0 +1,105 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#include +extern "C" { +#define restrict +#include "scotch.h" +} +#include "MEDPARTITIONER_Graph.hxx" +#include "MEDPARTITIONER_SCOTCHGraph.hxx" + +using namespace MEDPARTITIONER; + +SCOTCHGraph::SCOTCHGraph():Graph() +{ +} + +SCOTCHGraph::SCOTCHGraph(MEDPARTITIONER::MEDSKYLINEARRAY* graph, int* edgeweight):Graph(graph,edgeweight) +{ +} + +SCOTCHGraph::~SCOTCHGraph() +{ +} + +void SCOTCHGraph::partGraph(int ndomain, const std::string& options_string, ParaDomainSelector* sel) +{ + // number of graph vertices + int n = _graph->getNumberOf(); + + //graph + int * xadj=const_cast(_graph->getIndex()); + int * adjncy = const_cast(_graph->getValue()); + + //ndomain + int nparts = ndomain; + + // output parameters + int* partition = new int[n+1]; + + SCOTCH_Graph scotch_graph; + + SCOTCH_graphInit(&scotch_graph); + + + SCOTCH_graphBuild(&scotch_graph, + 1, //premier indice 1 + n, // nb of graph nodes + xadj, + 0, + _cellweight, //graph vertices loads + 0, + xadj[n], // number of edges + adjncy, + _edgeweight); + + SCOTCH_Strat scotch_strategy; + SCOTCH_stratInit(&scotch_strategy); + + //!user-defined options for the strategy + if (options_string!="") + SCOTCH_stratGraphMap(&scotch_strategy,options_string.c_str()); + + + if (nparts>1) + SCOTCH_graphPart(&scotch_graph,nparts,&scotch_strategy,partition); + else + // partition for 1 subdomain + for (int i=0; i index(n+1); + std::vector value(n); + index[0]=0; + for (int i=0; i + +namespace MEDPARTITIONER { + class MEDSKYLINEARRAY; + class MEDPARTITIONER_EXPORT SCOTCHGraph:public Graph + { + public: + SCOTCHGraph(); + SCOTCHGraph(MEDPARTITIONER::MEDSKYLINEARRAY*, int* edgeweight=0); + virtual ~SCOTCHGraph(); + void partGraph(int ndomain, const std::string& options_string="", ParaDomainSelector* sel=0); + }; +} +#endif /*SCOTCHGRAPH_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_SkyLineArray.cxx b/src/MEDPartitioner/MEDPARTITIONER_SkyLineArray.cxx new file mode 100644 index 000000000..3f0cf14fd --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_SkyLineArray.cxx @@ -0,0 +1,55 @@ +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#include "MEDPARTITIONER_SkyLineArray.hxx" +#include + + +using namespace MEDPARTITIONER; + +MEDSKYLINEARRAY::MEDSKYLINEARRAY() +{ + //MESSAGE_MED("Constructeur MEDSKYLINEARRAY sans parametre"); +} + +MEDSKYLINEARRAY::MEDSKYLINEARRAY(const MEDSKYLINEARRAY &myArray) +{ + _index=myArray._index; + _value=myArray._value; +} + +MEDSKYLINEARRAY::~MEDSKYLINEARRAY() +{ + //MESSAGE_MED("Destructeur ~MEDSKYLINEARRAY"); + + //if (_index != NULL) delete [] _index; + //if (_value != NULL) delete [] _value; +} + + +MEDSKYLINEARRAY::MEDSKYLINEARRAY(const std::vector& index, const std::vector& value) +{ + _value=value; + _index=index; +} + + diff --git a/src/MEDPartitioner/MEDPARTITIONER_SkyLineArray.hxx b/src/MEDPartitioner/MEDPARTITIONER_SkyLineArray.hxx new file mode 100644 index 000000000..d56fd94bb --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_SkyLineArray.hxx @@ -0,0 +1,81 @@ +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +# ifndef __MEDPARTITIONER_MEDSKYLINEARRAY_H__ +# define __MEDPARTITIONER_MEDSKYLINEARRAY_H__ + +#include "MEDPARTITIONER.hxx" + +#include + +namespace MEDPARTITIONER { +class MEDPARTITIONER_EXPORT MEDSKYLINEARRAY +{ +private : + std::vector _index; + std::vector _value; + +public : + // Attention, avec ce constructeur, il n'est possible de remplir le MEDSKYLINEARRAY + MEDSKYLINEARRAY(); + + // Constructeur par recopie + MEDSKYLINEARRAY( const MEDSKYLINEARRAY &myArray ); + + // Avec ce constructeur la m�moire pour le tableau de valeur et le + // tableau d'index est r�serv�e. Il suffit d'effectuer les s�quences + // d'appels suivantes pour initialiser le MEDSKYLINEARRAY + // 1) setIndex(index) puis fois setI(i,&listValeurN�I) avec i dans 1..count + // rem : listValeurN�I est dupliqu�e + // 2) appeler fois setIJ(i,j,valeur) avec i dans 1..count et avec j dans 1..count + MEDSKYLINEARRAY( const std::vector& index, const std::vector& value ); + + ~MEDSKYLINEARRAY(); + //void setMEDSKYLINEARRAY( const int count, const int length, int* index , int* value ) ; + + inline int getNumberOf() const; + inline int getLength() const; + inline const int* getIndex() const; + inline const int* getValue() const; +}; + +// --------------------------------------- +// Methodes Inline +// --------------------------------------- +inline int MEDSKYLINEARRAY::getNumberOf() const +{ + return _index.size()-1; +} +inline int MEDSKYLINEARRAY::getLength() const +{ + return _value.size() ; +} +inline const int* MEDSKYLINEARRAY::getIndex() const +{ + return (const int*)(&_index[0]) ; +} +inline const int* MEDSKYLINEARRAY::getValue() const +{ + return (const int*)(&_value[0]) ; +} +} +# endif diff --git a/src/MEDPartitioner/MEDPARTITIONER_Topology.hxx b/src/MEDPartitioner/MEDPARTITIONER_Topology.hxx new file mode 100644 index 000000000..d0821e7af --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_Topology.hxx @@ -0,0 +1,171 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef MEDPARTITIONER_TOPOLOGY_HXX_ +#define MEDPARTITIONER_TOPOLOGY_HXX_ + +//#include "boost/shared_ptr.hpp" + +#include +#include + +namespace MEDPARTITIONER +{ + class CONNECTZONE; + class MEDSKYLINEARRAY; +} +namespace ParaMEDMEM { + class MEDCouplingUMesh; +} + +namespace MEDPARTITIONER { + + class Graph; + class MESHCollection; + class MEDPARTITIONER_FaceModel; + + // typedef std::map > TGeom2Faces; + // typedef std::vector< TGeom2Faces > TGeom2FacesByDomain; + + class Topology + { + public: + Topology(){} + Topology(std::vector, std::vector){} + + virtual ~Topology(){} + + //!converts a list of global cell numbers + //!to a distributed array with local cell numbers + virtual void convertGlobalNodeList(const int* list, int nb, int* local, int*ip)=0; + virtual void convertGlobalNodeList(const int* list, int nb, int* local, int ip)=0; + //!converts a list of global node numbers + //!to a distributed array with local cell numbers + virtual void convertGlobalCellList(const int*list , int nb, int* local, int*ip)=0; + + //!converts a list of global face numbers + //!to a distributed array with local face numbers + virtual void convertGlobalFaceList(const int*list , int nb, int* local, int*ip)=0; + virtual void convertGlobalFaceList(const int*list , int nb, int* local, int ip)=0; + virtual void convertGlobalFaceListWithTwins(const int* face_list, int nbface, int*& local, int*& ip, int*& full_array, int& size)=0; + virtual void convertGlobalNodeListWithTwins(const int* face_list, int nbnode, int*& local, int*& ip, int*& full_array, int& size)=0; + + + //number of doamins + virtual int nbDomain() const =0; + + //number of cells + virtual int nbCells() const=0; + + //number of nodes + virtual int nbNodes() const=0; + + //number of cells on a specific domain + virtual int nbCells(int idomain) const=0; + + // ////creating node mapping + // virtual void createNodeMapping(std::map& type_connectivity, + // std::map& present_type_numbers, + // std::vector& polygon_conn, +// std::vector& polygon_conn_index, +// std::vector& polyhedron_conn, +// std::vector& polyhedron_conn_index, +// std::vector& polyhedron_face_index, +// int domain)=0; + + ////creating face mapping + // virtual void createFaceMapping(std::map& type_connectivity, + // std::map& present_type_numbers,int domain)=0; + // + // virtual void createFaceMapping(const MESHCollection&,const MESHCollection&)=0; + + // //converting node global numberings to local numberings + // virtual void convertToLocal(std::map& type_connectivity, + // std::map& present_type_numbers, + // int idomain, + // MED_EN::medEntityMesh entity)=0; + //converting node global numberings to local numberings + virtual void convertToLocal2ndVersion(int*,int,int)=0; + + virtual int convertNodeToGlobal(int ip,int icell)const=0; + virtual int convertFaceToGlobal(int ip,int icell)const=0; + virtual int convertCellToGlobal(int ip,int icell)const=0; + + virtual void convertNodeToGlobal(int ip,const int* local, int n, int* global)const=0 ; + virtual void convertCellToGlobal(int ip,const int* local, int n, int* global)const=0 ; + virtual void convertFaceToGlobal(int ip,const int* local, int n, int* global)const=0 ; + + //retrieving number of nodes + virtual int getNodeNumber(int idomain) const =0; + virtual int getNodeNumber() const=0; + //retrieving list of nodes + virtual void getNodeList(int idomain, int* list) const =0; + + virtual std::vector & getFusedCellNumbers(int idomain) = 0; + virtual const std::vector & getFusedCellNumbers(int idomain) const = 0; + + virtual std::vector & getFusedFaceNumbers(int idomain) = 0; + virtual const std::vector & getFusedFaceNumbers(int idomain) const = 0; + + //retrieving number of nodes + virtual int getCellNumber(int idomain) const =0; + + //retrieving list of nodes + virtual void getCellList(int idomain, int* list) const =0; + + //retrieving number of faces + virtual int getFaceNumber(int idomain) const =0; + virtual int getFaceNumber()const =0; + + //retrieving list of nodes + virtual void getFaceList(int idomain, int* list) const =0; + + //adding a face to the mapping + virtual void appendFace(int idomain, int ilocal, int iglobal)=0; + + //return max global face number + virtual int getMaxGlobalFace()const=0; + + //return next free global face number + //virtual int nextGlobalFace(int start_num) const=0; + + //!converting a global cell number to a local representation + virtual std::pair convertGlobalCell(int iglobal) const =0; + + //converting a global face number to a local representation + virtual int convertGlobalFace(int iglobal, int idomain)=0; + + //converting a global node number to a local representation + virtual int convertGlobalNode(int iglobal, int idomain)=0; + + //! computing arrays with node/node correspondencies + // virtual void computeNodeNodeCorrespondencies(int nbdomain, std::vector&) const =0; + + //! computing arrays with cell/cell correspondencies + // virtual void computeCellCellCorrespondencies(int nbdomain, std::vector&, const Graph*) const =0; + + + //!recreating a face mapping from scratch + // virtual void recreateFaceMapping(const TGeom2FacesByDomain& )=0; + + //!recreating cell and node mapping after send-reveive and fusion of domain meshes + // virtual void recreateMappingAfterFusion(const std::vector& ) = 0; + }; +} + +#endif /*TOPOLOGY_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_UserGraph.cxx b/src/MEDPartitioner/MEDPARTITIONER_UserGraph.cxx new file mode 100644 index 000000000..948ec64cb --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_UserGraph.cxx @@ -0,0 +1,54 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#include "MEDPARTITIONER_Graph.hxx" +#include "MEDPARTITIONER_UserGraph.hxx" +#include +#include +using namespace MEDPARTITIONER; + +/*! constructor that allows to specify the desired partition + * \param partition as table giving the domain number for each cell + * (domain numbers range from 0 to ndomain-1 + * \param n number of cells in the mesh + */ +UserGraph::UserGraph(MEDPARTITIONER::MEDSKYLINEARRAY* array, const int* partition, int n):Graph(array,0) +{ + + std::vector index(n+1),value(n); + + index[0]=0; + for (int i=0; i +namespace MEDPARTITIONER +{ + class MEDSKYLINEARRAY; + class ParaDomainSelector; + class MEDPARTITIONER_EXPORT UserGraph : public Graph + { + public: + UserGraph(MEDPARTITIONER::MEDSKYLINEARRAY*, const int*, int); + virtual ~UserGraph(); + void partGraph(int, const std::string& options=std::string(""), ParaDomainSelector* sel=0); + }; +} + +#endif /*MEDPARTITIONER_USERGRAPH_HXX_*/ diff --git a/src/MEDPartitioner/MEDPARTITIONER_utils.cxx b/src/MEDPartitioner/MEDPARTITIONER_utils.cxx new file mode 100644 index 000000000..19843caba --- /dev/null +++ b/src/MEDPartitioner/MEDPARTITIONER_utils.cxx @@ -0,0 +1,1627 @@ +// Copyright (C) 2007-2010 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#include "MEDLoader.hxx" +#include "MEDLoaderBase.hxx" +#include "MEDFileUtilities.hxx" +#include "CellModel.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" +#include "MEDPARTITIONER_utils.hxx" +#include "InterpKernelException.hxx" +#include "MEDCouplingAutoRefCountObjectPtr.hxx" +#include "InterpKernelAutoPtr.hxx" + +#ifdef HAVE_MPI2 +#include +#endif + +#include +#include +#include +#include +#include + +using namespace std; +using namespace MEDPARTITIONER; + +int MEDPARTITIONER::MyGlobals::_verbose=0; +int MEDPARTITIONER::MyGlobals::_is0verbose=0; +int MEDPARTITIONER::MyGlobals::_rank=-1; +int MEDPARTITIONER::MyGlobals::_world_size=-1; +int MEDPARTITIONER::MyGlobals::_randomize=0; +int MEDPARTITIONER::MyGlobals::_atomize=0; +int MEDPARTITIONER::MyGlobals::_creates_boundary_faces=0; +vector MEDPARTITIONER::MyGlobals::_fileNames; +vector MEDPARTITIONER::MyGlobals::_meshNames; +vector MEDPARTITIONER::MyGlobals::_fieldDescriptions; +vector MEDPARTITIONER::MyGlobals::_generalInformations; + +string MEDPARTITIONER::trim(string& s, const string& drop) +{ + string r=s.erase(s.find_last_not_of(drop)+1); + return r.erase(0,r.find_first_not_of(drop)); +} + +string MEDPARTITIONER::intToStr(int i) +{ + ostringstream oss; + oss<>res; + return res; +} + +string MEDPARTITIONER::doubleToStr(double i) +{ + ostringstream oss; + oss<>res; + return res; +} + +bool MEDPARTITIONER::testArg(const char *arg, const char *argExpected, string& argValue) +{ + argValue=""; + int i; + for (i=0; i MEDPARTITIONER::createRandomSize(int size) +{ + vector res(size); + for (int i=0; i& ran, vector& vx, vector& va) +//randomize a xadj and adjncy, renumbering vertices belong rand. +//work only on one processor!!!! +{ + if (MyGlobals::_world_size>1) + { + cerr<<"randomizeAdj only on one proc!"< invran(size); + for (int i=0; i r=createRandomSize(size); + vector vx,va; + randomizeAdj(&xadj[0],&adjncy[0],r,vx,va); + for (int i=0; i& vec) +{ + if (vec.size()==0) return string(" NONE\n"); + ostringstream oss; + for (vector::const_iterator i=vec.begin(); i!=vec.end(); ++i) + oss<<" -> '"<<*i<<"'"<& vec, string sep) +{ + if (vec.size()==0) return string(" NONE\n"); + ostringstream oss; + for (vector::const_iterator i=vec.begin(); i!=vec.end(); ++i) + oss<& mymap) +{ + if (mymap.size()==0) return string(" NONE\n"); + ostringstream oss; + for (map::const_iterator i=mymap.begin(); i!=mymap.end(); ++i) + oss<<" -> ["<<(*i).first<<"]="<<(*i).second< >& mymap) +{ + if (mymap.size()==0) return string(" NONE\n"); + ostringstream oss; + for (map< string,vector >::const_iterator i=mymap.begin(); i!=mymap.end(); ++i) + oss<<" -> ["<<(*i).first<<"]="<& vec, string sep) +{ + if (vec.size()==0) return string(" NONE\n"); + ostringstream oss; + for (vector::const_iterator i=vec.begin(); i!=vec.end(); ++i) + { + oss<<" ->"; + oss< vec2=deserializeToVectorOfString(*i); + for (vector::const_iterator j=vec2.begin(); j!=vec2.end(); ++j) + oss<.push_back("toto") on serialized_FromVectorOfString_string +{ + ostringstream oss; + oss<& vec) +//a vector of string gives a string +{ + ostringstream oss; + for (vector::const_iterator i=vec.begin(); i!=vec.end(); ++i) + oss< MEDPARTITIONER::deserializeToVectorOfString(const string& str) +//a string gives a vector of string +{ + //vector res=new vector; + vector res; + size_t pos=0; + size_t posmax=str.size(); + if (posmax==0) return res; //empty vector + size_t length; + while (pos < posmax-6) //setw(5)+" " + { + istringstream iss(str.substr(pos,5)); + iss>>length; + //cout<<"length "< vec=deserializeToVectorOfString(fromStr); + vector res; + for (int i=0; i MEDPARTITIONER::vectorizeFromMapOfStringInt(const map& mymap) +//elements first and second of map give one elements in result vector of string +//converting formatted the int second as firsts characters ending at first slash +{ + vector res; + for (map::const_iterator i=mymap.begin(); i!=mymap.end(); ++i) + { + ostringstream oss; + oss<<(*i).second<<"/"<<(*i).first; + res.push_back(oss.str()); + } + return res; +} + +map MEDPARTITIONER::devectorizeToMapOfStringInt(const vector& vec) +//if existing identicals (first,second) in vector no problem, else Exception +{ + map res; + for (vector::const_iterator i=vec.begin(); i!=vec.end(); ++i) + { + size_t pos=0; + size_t posmax=(*i).size(); + size_t found=(*i).find('/'); //first slash + if ((found==string::npos) || (found<1)) + throw INTERP_KERNEL::Exception(LOCALIZED("Error aIntNumber/anyString is expected")); + int second; + istringstream iss((*i).substr(pos,found)); + iss>>second; + string first=(*i).substr(pos+found+1,posmax-found); + map::iterator it=res.find(first); + if (it!=res.end()) + if ((*it).second!=second) + throw INTERP_KERNEL::Exception(LOCALIZED("Error not the same map value")); + res[first]=second; + } + return res; +} + +vector MEDPARTITIONER::vectorizeFromMapOfStringVectorOfString(const map< string,vector >& mymap) +//elements first and second of map give one elements in result vector of string +//adding key map and length of second vector as first string in each serialized vector +//one serialized vector per key map +{ + vector res; + for (map< string,vector >::const_iterator i=mymap.begin(); i!=mymap.end(); ++i) + { + vector vs=(*i).second; //a vector of string; + ostringstream oss; + oss<<"Keymap/"<<(*i).first<<"/"<<(*i).second.size(); + vs.insert(vs.begin(), oss.str()); + res.push_back(serializeFromVectorOfString(vs)); + } + return res; +} + +map< string,vector > MEDPARTITIONER::devectorizeToMapOfStringVectorOfString(const vector& vec) +//if existing identicals keymap in vector no problem +//duplicates in second vector +{ + map< string,vector > res; + for (vector::const_iterator i=vec.begin(); i!=vec.end(); ++i) + { + vector vs=deserializeToVectorOfString(*i); + + string enTete=vs[0]; + size_t posmax=enTete.size(); + size_t foundKey=enTete.find("Keymap/"); + size_t foundSizeVector=enTete.find_last_of('/'); + if ((foundKey==string::npos) || (foundKey!=0) || ((foundKey+7)>=foundSizeVector)) + throw INTERP_KERNEL::Exception(LOCALIZED("Error Keymap/anyString/aIntNumber is expected")); + int sizeVector; + istringstream iss(enTete.substr(foundSizeVector+1,posmax-foundSizeVector)); + iss>>sizeVector; + string keymap=enTete.substr(foundKey+7,foundSizeVector-foundKey-7); + //cout< MEDPARTITIONER::recvVectorOfString(const int source) +//non conseille, interblocages, utiliser sendAndReceive +{ + throw INTERP_KERNEL::Exception(LOCALIZED("use sendAndReceiveVectorOfString please.")); + int recSize=0; + int tag=111000; + MPI_Status status; + MPI_Recv(&recSize, 1, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + string recData(recSize,'x'); + MPI_Recv((void*)recData.data(), recSize, MPI_CHAR, 1, tag+100, MPI_COMM_WORLD, &status); + //cout<<"proc "< MEDPARTITIONER::sendAndReceiveVectorOfString(const vector& vec, const int source, const int target) +//not optimized but suffisant +//return empty vector if i am not target +{ + int rank=MyGlobals::_rank; + + /*for test + ostringstream oss; + oss<<"sendAndReceive from "< res; + return res; //empty one for other proc +} + +vector MEDPARTITIONER::allgathervVectorOfString(const vector& vec) +//strings NO need all same size!!!! +{ + int world_size=MyGlobals::_world_size; + string str=serializeFromVectorOfString(vec); + + /*for test + int rank=MyGlobals::_rank; + ostringstream oss; + oss<<"allgatherv from "< indexes(world_size); + int size=str.length(); + MPI_Allgather(&size, 1, MPI_INT, + &indexes[0], 1, MPI_INT, MPI_COMM_WORLD); + + /*{ + ostringstream oss; + for (int i=0; i disp(1,0); + for (int i=0; i deserial=deserializeToVectorOfString(recData); + if (MyGlobals::_verbose>1000) + { + cout<<"proc "<& vec, const int source, const int target) +//TODO + +string MEDPARTITIONER::cle1ToStr(string s, int inew) +{ + ostringstream oss; + oss<>inew; +} + +string MEDPARTITIONER::cle2ToStr(string s, int inew, int iold) +{ + ostringstream oss; + oss<>inew>>iold; +} + +string MEDPARTITIONER::extractFromDescription(string description,string tag) +{ + size_t found=description.find(tag); + if ((found==string::npos) || (found<1)) + { + cerr<<"ERROR : not found '"<alloc(v.size()/nbComponents,nbComponents); + std::copy(v.begin(),v.end(),p->getPointer()); + return p; +} + +ParaMEDMEM::DataArrayDouble* MEDPARTITIONER::createDataArrayDoubleFromVector(vector& v) +{ + ParaMEDMEM::DataArrayDouble* p=DataArrayDouble::New(); + p->alloc(v.size(),1); + std::copy(v.begin(),v.end(),p->getPointer()); + return p; +} + +vector MEDPARTITIONER::browseFieldDouble(const MEDCouplingFieldDouble* fd) +//quick almost human readable information on a field double +/* example done by fd->simpleRepr() : +FieldDouble with name : "VectorFieldOnCells" +Description of field is : "" +FieldDouble space discretization is : P0 +FieldDouble time discretization is : One time label. Time is defined by iteration=0 order=1 and time=2. +Time unit is : "" +FieldDouble nature of field is : NoNature +FieldDouble default array has 3 components and 30000 tuples. +FieldDouble default array has following info on components : "vx" "vy" "vz" +Mesh support information : +__________________________ +Unstructured mesh with name : "testMesh" +Description of mesh : "" +Time attached to the mesh [unit] : 0 [] +Iteration : -1 Order : -1 +Mesh dimension : 3 +Space dimension : 3 +Info attached on space dimension : "" "" "" +Number of nodes : 33201 +Number of cells : 30000 +Cell types present : NORM_HEXA8 +*/ +{ + vector res; + //res.push_back("fieldName="); res.back()+=fd->getName(); + //not saved in file? res.push_back("fieldDescription="); res.back()+=fd->getDescription(); + //ret << "FieldDouble space discretization is : " << _type->getStringRepr() << "\n"; + //ret << "FieldDouble time discretization is : " << _time_discr->getStringRepr() << "\n"; + //ret << "FieldDouble nature of field is : " << MEDCouplingNatureOfField::getRepr(_nature) << "\n"; + if (fd->getArray()) + { + int nb=fd->getArray()->getNumberOfComponents(); + res.push_back("nbComponents="); res.back()+=intToStr(nb); + //ret << "FieldDouble default array has " << nbOfCompo << " components and " << getArray()->getNumberOfTuples() << " tuples.\n"; + //ret << "FieldDouble default array has following info on components : "; + for (int i=0; igetInfoOnComponent(i) << "\" "; + res.push_back("componentInfo"); + res.back()+=intToStr(i)+"="+fd->getArray()->getInfoOnComponent(i); + } + } + else + { + res.push_back("nbComponents=0"); //unknown + } + return res; +} + +vector MEDPARTITIONER::browseAllFields(const string& myfile) +//quick almost human readable information on all fields in a .med file +{ + vector res; + vector meshNames=MEDLoader::GetMeshNames(myfile.c_str()); + + for (int i=0; i fieldNames= + MEDLoader::GetAllFieldNamesOnMesh(myfile.c_str(),meshNames[i].c_str()); + for (int j = 0; j < fieldNames.size(); j++) + { + vector< ParaMEDMEM::TypeOfField > typeFields= + MEDLoader::GetTypesOfField(myfile.c_str(), meshNames[i].c_str(), fieldNames[j].c_str()); + for (int k = 0; k < typeFields.size(); k++) + { + vector< pair< int, int > > its= + MEDLoader::GetFieldIterations(typeFields[k], myfile.c_str(), meshNames[i].c_str(), fieldNames[j].c_str()); + if (MyGlobals::_is0verbose>100) cout<<"fieldName "< resi; + resi.push_back("fileName="); resi.back()+=myfile; + resi.push_back("meshName="); resi.back()+=meshNames[i]; + resi.push_back("fieldName="); resi.back()+=fieldNames[j]; + resi.push_back("typeField="); resi.back()+=intToStr((int)typeFields[k]); + resi.push_back("DT="); resi.back()+=intToStr((int)its[m].first); + resi.push_back("IT="); resi.back()+=intToStr((int)its[m].second); + res.push_back(serializeFromVectorOfString(resi)); + } + } + } + } + return res; +} + +/* +vector MEDPARTITIONER::browseAllFieldsOnMesh(const string& myfile, const string& mymesh, int idomain) +//quick almost human readable information on all fields on a mesh in a .med file using MEDFILEBROWSER +{ + vector res; + vector meshNames; + meshNames.push_back(mymesh); + MEDMEM::MEDFILEBROWSER myMed(myfile); + for (int i=0; i MEDPARTITIONER::GetInfosOfField(const char *fileName, const char *meshName, int idomain) +{ +const int lggeom=10; +const med_geometry_type GEOMTYPE[lggeom]={ //MED_N_CELL_FIXED_GEO] = { + //MED_POINT1, + //MED_SEG2, + //MED_SEG3, + //MED_SEG4, + //MED_TRIA3, + //MED_QUAD4, + //MED_TRIA6, + //MED_TRIA7, + //MED_QUAD8, + //MED_QUAD9, + MED_TETRA4, + MED_PYRA5, + MED_PENTA6, + MED_HEXA8, + MED_OCTA12, + MED_TETRA10, + MED_PYRA13, + MED_PENTA15, + MED_HEXA20, + MED_HEXA27, + //MED_POLYGON, + //MED_POLYHEDRON +}; + +const char * const GEOMTYPENAME[lggeom]={ + //"MED_POINT1", + //"MED_SEG2", + //"MED_SEG3", + //"MED_SEG4", + //"MED_TRIA3", + //"MED_QUAD4", + //"MED_TRIA6", + //"MED_TRIA7", + //"MED_QUAD8", + //"MED_QUAD9", + "MED_TETRA4", + "MED_PYRA5", + "MED_PENTA6", + "MED_HEXA8", + "MED_OCTA12", + "MED_TETRA10", + "MED_PYRA13", + "MED_PENTA15", + "MED_HEXA20", + "MED_HEXA27", + //"MED_POLYGONE", + //"MED_POLYEDRE", +}; + + +const int lgentity=3; +const med_entity_type ENTITYTYPE[lgentity]={ //MED_N_ENTITY_TYPES+2]={ + //MED_UNDEF_ENTITY_TYPE, + MED_CELL, + //MED_DESCENDING_FACE, + //MED_DESCENDING_EDGE, + MED_NODE, + MED_NODE_ELEMENT, + //MED_STRUCT_ELEMENT, + //MED_UNDEF_ENTITY_TYPE +}; + +const char * const ENTITYTYPENAME[lgentity]={ //MED_N_ENTITY_TYPES+2]={ + //"MED_UNDEF_ENTITY_TYPE", + "MED_CELL", + //"MED_FACE", + //"MED_ARETE", + "MED_NODE", + "MED_NODE_ELEMENT", + //"MED_STRUCT_ELEMENT", + //"MED_UNDEF_ENTITY_TYPE" +}; + + vector res; + med_idt fid=MEDfileOpen(fileName,MED_ACC_RDONLY); + med_int nbFields=MEDnField(fid); + if (MyGlobals::_verbose>20) cout<<"on filename "< comp=new char[ncomp*MED_SNAME_SIZE+1]; + INTERP_KERNEL::AutoPtr unit=new char[ncomp*MED_SNAME_SIZE+1]; + INTERP_KERNEL::AutoPtr dt_unit=new char[MED_LNAME_SIZE+1]; + med_int nbPdt; + MEDfieldInfo(fid,i,nomcha,maa_ass,&localmesh,&typcha,comp,unit,dt_unit,&nbPdt); + std::string curFieldName=MEDLoaderBase::buildStringFromFortran(nomcha,MED_NAME_SIZE+1); + std::string curMeshName=MEDLoaderBase::buildStringFromFortran(maa_ass,MED_NAME_SIZE+1); + for (int k=1; k<=nbPdt; k++) + { + MEDfieldComputingStepInfo(fid,nomcha,k,&numdt,&numo,&dt); + if (MyGlobals::_verbose>20) + cout<<"on filename "<0) + { + if (MyGlobals::_verbose>20) + cout<<"on filename "< resi; + resi.push_back("idomain="); resi.back()+=intToStr(idomain); + resi.push_back("fileName="); resi.back()+=fileName; + resi.push_back("meshName="); resi.back()+=curMeshName; + resi.push_back("fieldName="); resi.back()+=curFieldName; + resi.push_back("typeField="); resi.back()+=intToStr((int)ON_NODES); + resi.push_back("typeData="); resi.back()+=intToStr((int)typcha); //6 for double? + resi.push_back("nbComponent="); resi.back()+=intToStr((int)ncomp); + resi.push_back("DT="); resi.back()+=intToStr((int)numdt); + resi.push_back("IT="); resi.back()+=intToStr((int)numo); + resi.push_back("time="); resi.back()+=doubleToStr(dt); + resi.push_back("entity="); resi.back()+=intToStr((int)enttype); + resi.push_back("entityName="); resi.back()+=ENTITYTYPENAME[ie]; + resi.push_back("nbOfVal="); resi.back()+=intToStr((int)nbOfVal); + resi.push_back("profilName="); resi.back()+=pflname; + resi.push_back("profileSize="); resi.back()+=intToStr((int)profilesize); + resi.push_back("nbPtGauss="); resi.back()+=intToStr((int)nbi); + res.push_back(serializeFromVectorOfString(resi)); + } + break; //on nodes no need to scute all geomtype + } + else + { + med_geometry_type mygeomtype=GEOMTYPE[j]; + med_int nbOfVal=MEDfieldnValueWithProfile(fid,nomcha,numdt,numo,enttype,mygeomtype,profileit, + MED_COMPACT_PFLMODE,pflname,&profilesize,locname,&nbi); + if (nbOfVal>0) + { + if (MyGlobals::_verbose>20) + cout<<"on filename "< resi; + resi.push_back("idomain="); resi.back()+=intToStr(idomain); + resi.push_back("fileName="); resi.back()+=fileName; + resi.push_back("meshName="); resi.back()+=curMeshName; + resi.push_back("fieldName="); resi.back()+=curFieldName; + resi.push_back("typeField="); resi.back()+=intToStr((int)typeField); + resi.push_back("typeData="); resi.back()+=intToStr((int)typcha); //6 for double? + resi.push_back("nbComponent="); resi.back()+=intToStr((int)ncomp); + resi.push_back("DT="); resi.back()+=intToStr((int)numdt); + resi.push_back("IT="); resi.back()+=intToStr((int)numo); + resi.push_back("time="); resi.back()+=doubleToStr(dt); + resi.push_back("entity="); resi.back()+=intToStr((int)enttype); + resi.push_back("entityName="); resi.back()+=ENTITYTYPENAME[ie]; + resi.push_back("geomType="); resi.back()+=intToStr((int)GEOMTYPE[j]); + resi.push_back("geomTypeName="); resi.back()+=GEOMTYPENAME[j]; + resi.push_back("nbOfVal="); resi.back()+=intToStr((int)nbOfVal); + resi.push_back("profilName="); resi.back()+=pflname; + resi.push_back("profileSize="); resi.back()+=intToStr((int)profilesize); + resi.push_back("nbPtGauss="); resi.back()+=intToStr((int)nbi); + if (typeField==-1) + { + cout<<"WARNING : unknown typeField for entity type "<10) cout<<"detected fields:\n"< MEDPARTITIONER::browseAllFieldsOnMesh(const string& myfile, const string& mymesh, int idomain) +//quick almost human readable information on all fields on a mesh in a .med file +{ + vector res=GetInfosOfField(myfile.c_str(),mymesh.c_str(),idomain); + return res; + + /*obsolete do no work on GetTypesOfField ON_GAUSS_NE + vector res; + vector meshNames; + meshNames.push_back(mymesh); + + for (int i=0; i fieldNames= + MEDLoader::GetAllFieldNamesOnMesh(myfile.c_str(),meshNames[i].c_str()); + for (int j=0; j typeFields= + MEDLoader::GetTypesOfField(myfile.c_str(), meshNames[i].c_str(), fieldNames[j].c_str()); + //if (MyGlobals::_is0verbose>100) cout<<"fieldName "< > its; + its=MEDLoader::GetFieldIterations(typeFields[k], myfile.c_str(), meshNames[i].c_str(), fieldNames[j].c_str()); + //if (typeFields[k]==ON_GAUSS_NE) its.push_back(make_pair(5,6)); + if (MyGlobals::_is0verbose>100) cout<<"fieldName "< resi; + resi.push_back("idomain="); resi.back()+=intToStr(idomain); + resi.push_back("fileName="); resi.back()+=myfile; + resi.push_back("meshName="); resi.back()+=meshNames[i]; + resi.push_back("fieldName="); resi.back()+=fieldNames[j]; + resi.push_back("typeField="); resi.back()+=intToStr((int)typeFields[k]); + resi.push_back("DT="); resi.back()+=intToStr((int)its[m].first); + resi.push_back("IT="); resi.back()+=intToStr((int)its[m].second); + //cout<<"browseAllFieldsOnMesh add "<& vec, int target) +{ + int tag = 111002; + int size=vec.size(); + if (MyGlobals::_verbose>1000) + cout<<"proc "< sendDoubleVec "<(&vec[0]), size, MPI_DOUBLE, target, tag+100, MPI_COMM_WORLD); +#endif +} + +/*! Receives messages from proc \a source to fill vector vec. +To be used with \a sendDoubleVec method. + +\param vec vector that is filled +\param source processor id of the incoming messages + */ +std::vector* MEDPARTITIONER::recvDoubleVec(int source) +{ + int tag = 111002; + int size; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(&size, 1, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + if (MyGlobals::_verbose>1000) + cout<<"proc "<* vec=new std::vector; + vec->resize(size); + MPI_Recv(&vec[0], size, MPI_DOUBLE, source, tag+100, MPI_COMM_WORLD, &status); +#endif + return vec; +} + +void MEDPARTITIONER::recvDoubleVec(std::vector& vec, int source) +{ + int tag = 111002; + int size; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(&size, 1, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + if (MyGlobals::_verbose>1000) + cout<<"proc "<& vec, int target) +{ + int tag = 111003; + int size=vec.size(); + if (MyGlobals::_verbose>1000) + cout<<"proc "< sendIntVec "<(&vec[0]), size,MPI_INT, target, tag+100, MPI_COMM_WORLD); +#endif +} + +/*! Receives messages from proc \a source to fill vector vec. +To be used with \a sendIntVec method. +\param vec vector that is filled +\param source processor id of the incoming messages + */ +std::vector* MEDPARTITIONER::recvIntVec(int source) +{ + int tag = 111003; + int size; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(&size, 1, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + if (MyGlobals::_verbose>1000) + cout<<"proc "<* vec=new std::vector; + vec->resize(size); + MPI_Recv(&vec[0], size, MPI_INT, source, tag+100, MPI_COMM_WORLD, &status); +#endif + return vec; +} + +void MEDPARTITIONER::recvIntVec(std::vector& vec, int source) +{ + int tag = 111003; + int size; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(&size, 1, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + if (MyGlobals::_verbose>1000) + cout<<"proc "<getNbOfElems(); + size[1]=da->getNumberOfTuples(); + size[2]=da->getNumberOfComponents(); + if (MyGlobals::_verbose>1000) + cout<<"proc "< sendDataArrayInt "<getPointer(); + MPI_Send(const_cast(&p[0]), size[0] ,MPI_INT, target, tag+100, MPI_COMM_WORLD); +#endif +} + +/*! Receives messages from proc \a source to fill dataArrayInt da. +To be used with \a sendIntVec method. +\param da dataArrayInt that is filled +\param source processor id of the incoming messages + */ +ParaMEDMEM::DataArrayInt* MEDPARTITIONER::recvDataArrayInt(int source) +//std::vector& vec, int source)const +{ + int tag = 111004; + int size[3]; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(size, 3, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + if (MyGlobals::_verbose>1000) + cout<<"proc "<alloc(size[1],size[2]); + int * p=da->getPointer(); + MPI_Recv(const_cast(&p[0]), size[0], MPI_INT, source, tag+100, MPI_COMM_WORLD, &status); +#endif + return da; +} + +/*! +Sends content of \a dataArrayInt to processor \a target. +To be used with \a recvDataArrayDouble method. +\param da dataArray to be sent +\param target processor id of the target +*/ +void MEDPARTITIONER::sendDataArrayDouble(ParaMEDMEM::DataArrayDouble* da, int target) +{ + if (da==0) throw INTERP_KERNEL::Exception(LOCALIZED("Problem send DataArrayDouble* NULL")); + int tag = 111005; + int size[3]; + size[0]=da->getNbOfElems(); + size[1]=da->getNumberOfTuples(); + size[2]=da->getNumberOfComponents(); + if (MyGlobals::_verbose>1000) + cout<<"proc "< sendDataArrayDouble "<getPointer(); + MPI_Send(const_cast(&p[0]), size[0] ,MPI_DOUBLE, target, tag+100, MPI_COMM_WORLD); +#endif +} + +/*! Receives messages from proc \a source to fill dataArrayDouble da. +To be used with \a sendDoubleVec method. +\param da dataArrayDouble that is filled +\param source processor id of the incoming messages + */ +ParaMEDMEM::DataArrayDouble* MEDPARTITIONER::recvDataArrayDouble(int source) +//std::vector& vec, int source)const +{ + int tag = 111005; + int size[3]; +#ifdef HAVE_MPI2 + MPI_Status status; + MPI_Recv(size, 3, MPI_INT, source, tag, MPI_COMM_WORLD, &status); + if (MyGlobals::_verbose>1000) + cout<<"proc "<alloc(size[1],size[2]); + double * p=da->getPointer(); + MPI_Recv(const_cast(&p[0]), size[0], MPI_DOUBLE, source, tag+100, MPI_COMM_WORLD, &status); +#endif + return da; +} + +ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::createEmptyMEDCouplingUMesh() + //create empty MEDCouplingUMesh* dim 3 +{ + ParaMEDMEM::MEDCouplingUMesh* umesh=ParaMEDMEM::MEDCouplingUMesh::New(); + umesh->setMeshDimension(3); + umesh->allocateCells(0); + umesh->finishInsertingCells(); + ParaMEDMEM::DataArrayDouble *myCoords=ParaMEDMEM::DataArrayDouble::New(); + myCoords->alloc(0,3); + umesh->setCoords(myCoords); + umesh->setName("EMPTY"); + myCoords->decrRef(); + umesh->checkCoherency(); + return umesh; +} + +void MEDPARTITIONER::testVectorOfStringMPI() +{ + int rank=MyGlobals::_rank; + int world_size=MyGlobals::_world_size; + vector myVector; + ostringstream oss; + oss<<"hello from "< v0=deserializeToVectorOfString(s0); + cout<<"v0 is : a vector of size "< res=deserializeToVectorOfString(s0); + if (res.size()!=myVector.size()) + throw INTERP_KERNEL::Exception(LOCALIZED("Problem in (de)serialise VectorOfString incoherent sizes")); + for (int i=0; i res=sendAndReceiveVectorOfString(myVector, i, j); + if ((rank==j) && MyGlobals::_verbose>20) + cout<<"proc "< res=allgathervVectorOfString(myVector); + //sometimes for test + res=allgathervVectorOfString(myVector); + res=allgathervVectorOfString(myVector); + if (rank==0 && MyGlobals::_verbose>20) + cout<<"proc "< myMap; + myMap["one"]=1; + myMap["two"]=22; //a bug + myMap["three"]=3; + myMap["two"]=2; //last speaking override + + if (rank==0) + { + vector v2=vectorizeFromMapOfStringInt(myMap); + /* + cout<<"v2 is : a vector of size "< m3=devectorizeToMapOfStringInt(v2); + if (reprMapOfStringInt(m3)!=reprMapOfStringInt(myMap)) + throw INTERP_KERNEL::Exception(LOCALIZED("Problem in (de)vectorize MapOfStringInt")); + } + + vector v2=allgathervVectorOfString(vectorizeFromMapOfStringInt(myMap)); + if (rank==0 && MyGlobals::_verbose>20) + { + cout<<"v2 is : a vector of size "< m2=devectorizeToMapOfStringInt(v2); + cout<<"m2 is : a map of size "< myVector; + ostringstream oss; + oss<<"hello from "< > m2; + m2["first key"]=myVector; + m2["second key"]=myVector; + vector v2=vectorizeFromMapOfStringVectorOfString(m2); + map< string,vector > m3=devectorizeToMapOfStringVectorOfString(v2); + if (rank==0 && MyGlobals::_verbose>20) + cout<<"m2 is : a MapOfStringVectorOfString of size "< > m4; + m4["1rst key"]=myVector; + m4["2snd key"]=myVector; + vector v4=allgathervVectorOfString(vectorizeFromMapOfStringVectorOfString(m4)); + if (rank==0 && MyGlobals::_verbose>20) + { + map< string,vector > m5=devectorizeToMapOfStringVectorOfString(v4); + map< string,vector > m6=deleteDuplicatesInMapOfStringVectorOfString(m5); + cout<<"m5 is : a map of size "<alloc(nbOfTuples,numberOfComponents); + vector vals; + for (int j=0; jgetPointer()); + if (rank==0) sendDataArrayInt(send, 1); + if (rank==1) recv=recvDataArrayInt(0); + if (rank==1 && MyGlobals::_verbose>20) + { + cout<repr()<repr()<repr()!=recv->repr()) + throw INTERP_KERNEL::Exception(LOCALIZED("Problem in send&recv DataArrayInt")); + } + send->decrRef(); + if (rank==1) recv->decrRef(); + } + //double + { + ParaMEDMEM::DataArrayDouble* send=ParaMEDMEM::DataArrayDouble::New(); + ParaMEDMEM::DataArrayDouble* recv=0; + int nbOfTuples=5; + int numberOfComponents=3; + send->alloc(nbOfTuples,numberOfComponents); + vector vals; + for (int j=0; jgetPointer()); + if (rank==0) sendDataArrayDouble(send, 1); + if (rank==1) recv=recvDataArrayDouble(0); + if (rank==1 && MyGlobals::_verbose>20) + { + cout<repr()<repr()<repr()!=recv->repr()) + throw INTERP_KERNEL::Exception(LOCALIZED("Problem in send&recv DataArrayDouble")); + } + send->decrRef(); + if (rank==1) recv->decrRef(); + } + + if (MyGlobals::_verbose) cout<<"proc "< x, y; + int tag=111111; + MPI_Request requete0, requete1; + MPI_Status statut; + int ok=0; + string res; + if (rang==0) + { + x.resize(taille); + MPI_Ssend_init(&x[0], taille, MPI_INT, 1, tag, MPI_COMM_WORLD , &requete0); + for(int k=0; k cela peut prendre du temps + MPI_Start(&requete0); + //Traitement sequentiel independant de "x" + //... + MPI_Wait(&requete0, &statut); + //Traitement sequentiel impliquant une modification de "x" en memoire + //x=... + } + MPI_Request_free(&requete0); + } + else if (rang == 1) + { + y.resize(taille); + MPI_Recv_init(&y[0], taille, MPI_INT, 0, tag, MPI_COMM_WORLD , &requete1); + for(int k=0; k cela peut prendre du temps + MPI_Start(&requete1); + //Traitement sequentiel independant de "y" + //... + MPI_Wait(&requete1, &statut); + //Traitement sequentiel dependant de "y" + //...=f(y) + int nb=0; + for (int i=0; i9) + { + res="0K"; if (nb!=taille) res="KO"; + cout<1) + cout<<"resultat "<=wsize) next=0; + vector x, y; + tagbefo=111111+befo; + tagnext=111111+rang; + MPI_Request requete0, requete1; + MPI_Status statut1, statut2; + int ok=0; + string res; + //cout<<"ini|"< cela peut prendre du temps + MPI_Start(&requete0); + //Reception du gros message --> cela peut prendre du temps + for (int i=0; i9) + { + res="0K"+intToStr(rang); if (nb!=taille) res="KO"+intToStr(rang); + cout<1) + cout<<"resultat proc "<=rangMax) couleur=MPI_UNDEFINED; + cout<<"coul|"<MyGlobals::_world_size) wsize=MyGlobals::_world_size; + befo=rang-1; if (befo<0) befo=wsize-1; + next=rang+1; if (next>=wsize) next=0; + vector x, y; + tagbefo=111111+befo; + tagnext=111111+rang; + MPI_Request requete0, requete1; + MPI_Status statut1, statut2; + int ok=0; + string res; + + //cout<<"ini|"< cela peut prendre du temps + MPI_Start(&requete0); + //Reception du gros message --> cela peut prendre du temps + for (int i=0; i9) + { + res="0K"+intToStr(rang); if (nb!=taille) res="KO"+intToStr(rang); + cout<1) + cout<<"resultat proc "< +#include +#include + +#ifdef LOCALIZED +#undef LOCALIZED +#endif + +#if defined(_DEBUG_) || defined(_DEBUG) +//# define LOCALIZED(message) #message , __FILE__ , __FUNCTION__ , __LINE__ +# define LOCALIZED(message) #message , __FUNCTION__ , __LINE__ +#else +# define LOCALIZED(message) #message +#endif + +namespace MEDPARTITIONER +{ + using namespace std; + using namespace ParaMEDMEM; + + string trim(string& s,const string& drop); + string intToStr(int i); + string doubleToStr(double i); + int strToInt(string s); + double strToDouble(string s); + bool testArg(const char *arg, const char *argExpected, string& argValue); + vector createRandomSize(int size); + void randomizeAdj(int* xadj, int* adjncy, vector& ran, vector& vx, vector& va); + void testRandomize(); + + string reprVectorOfString(const vector& vec); + string reprVectorOfString(const vector& vec, string sep); + string reprMapOfStringInt(const map& mymap); + string reprMapOfStringVectorOfString(const map< string,vector >& mymap); + string reprFieldDescriptions(const vector& vec, string sep); + + string serializeFromString(const string& s); + string serializeFromVectorOfString(const vector& vec); + vector deserializeToVectorOfString(const string& str); + string eraseTagSerialized(string fromStr, string tag); + + vector vectorizeFromMapOfStringInt(const map& mymap); + map devectorizeToMapOfStringInt(const vector& vec); + + vector vectorizeFromMapOfStringVectorOfString(const map< string,vector >& mymap); + map< string,vector > devectorizeToMapOfStringVectorOfString(const vector& vec); + + vector selectTagsInVectorOfString(const vector& vec, string tag); + vector deleteDuplicatesInVectorOfString(const vector& vec); + map< string,vector > deleteDuplicatesInMapOfStringVectorOfString(const map< string,vector >& mymap); + + string cle1ToStr(string s, int inew); + void cle1ToData(string cle, string& s, int& inew); + + string cle2ToStr(string s, int inew, int iold); + void cle2ToData(string cle, string& s, int& inew, int& iold); + + string extractFromDescription(string description, string tag); + void fieldDescriptionToData(string description, + int& idomain, string& fileName, string& meshName, string& fieldName, + int& typeField, int& DT, int& IT); + void fieldShortDescriptionToData(string description, + string& fieldName, int& typeField, int& entity, int& DT, int& IT); + + ParaMEDMEM::DataArrayInt* createDataArrayIntFromVector(vector& v); + ParaMEDMEM::DataArrayInt* createDataArrayIntFromVector(vector& v, int nbComponents); + ParaMEDMEM::DataArrayDouble* createDataArrayDoubleFromVector(vector& v); + + void sendVectorOfString(const vector& vec, const int target); + vector recvVectorOfString(const int source); + //TODO void sendrecvVectorOfString(const vector& vec, const int source, const int target); + vector sendAndReceiveVectorOfString(const vector& vec, const int source, const int target); + vector allgathervVectorOfString(const vector& vec); + + vector browseFieldDouble(const MEDCouplingFieldDouble* fd); + vector browseAllFields(const string& myfile); + vector browseAllFieldsOnMesh(const string& myfile, const string& mymesh, int idomain); + vector GetInfosOfField(const char *fileName, const char *meshName, int idomain ); + + void sendDoubleVec(const std::vector& vec, int target); + std::vector* recvDoubleVec(int source); + void recvDoubleVec(std::vector& vec, int source); + + void sendIntVec(const std::vector& vec, int target); + std::vector* recvIntVec(int source); + void recvIntVec(std::vector& vec, int source); + + void sendDataArrayInt(ParaMEDMEM::DataArrayInt* da, int target); + ParaMEDMEM::DataArrayInt* recvDataArrayInt(int source); + void sendDataArrayDouble(ParaMEDMEM::DataArrayDouble* da, int target); + ParaMEDMEM::DataArrayDouble* recvDataArrayDouble(int source); + + ParaMEDMEM::MEDCouplingUMesh* createEmptyMEDCouplingUMesh(); + + void testVectorOfStringMPI(); + void testMapOfStringIntMPI(); + void testMapOfStringVectorOfStringMPI(); + void testDataArrayMPI(); + void testPersistantMpi0To1(int taille, int nb); + void testPersistantMpiRing(int taille, int nb); + void testPersistantMpiRingOnCommSplit(int taille, int nb); + + class MyGlobals + { + public : static int _verbose; //0 to 1000 over 200 is debug + public : static int _rank; + public : static int _world_size; + public : static int _randomize; + public : static int _atomize; + public : static int _creates_boundary_faces; + public : static int _is0verbose; //cout if rank 0 and verbose + public : static vector _fileNames; //on [iold] + public : static vector _meshNames; //on [iold] + public : static vector _fieldDescriptions; + //used for descriptions of components of fields for example... + public : static vector _generalInformations; + + }; + + /*int MyGlobals::_verbose=0; + int MyGlobals::_is0verbose=0; + int MyGlobals::_rank=-1; + int MyGlobals::_world_size=-1;*/ +} +#endif /*MEDPARTITIONER_UTILS_HXX_*/ diff --git a/src/MEDPartitioner/Makefile.am b/src/MEDPartitioner/Makefile.am new file mode 100644 index 000000000..c731db984 --- /dev/null +++ b/src/MEDPartitioner/Makefile.am @@ -0,0 +1,129 @@ + # Copyright (C) 2007-2010 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. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# MED : MED files in memory +# +include $(top_srcdir)/adm_local/unix/make_common_starter.am + +# this directory must be recompiled before Test folder + +if CPPUNIT_IS_OK + SUBDIRS=. Test +endif + +lib_LTLIBRARIES= libmedpartitioner.la + +salomeinclude_HEADERS= \ +MEDPARTITIONER_MESHCollection.hxx \ +MEDPARTITIONER_MESHCollectionMedXMLDriver.H \ +MEDPARTITIONER_MESHCollectionMedAsciiDriver.H \ +MEDPARTITIONER_MESHCollectionDriver.hxx \ +MEDPARTITIONER_MESHCollectionMedXMLDriver.hxx \ +MEDPARTITIONER_MESHCollectionMedAsciiDriver.hxx \ +MEDPARTITIONER_ParallelTopology.hxx \ +MEDPARTITIONER_JointFinder.hxx \ +MEDPARTITIONER_Graph.hxx\ +MEDPARTITIONER_UserGraph.hxx\ +MEDPARTITIONER_utils.hxx \ +MEDPARTITIONER.hxx \ +MEDPARTITIONER_ParaDomainSelector.hxx \ +MEDPARTITIONER_ConnectZone.hxx \ +MEDPARTITIONER_SkyLineArray.hxx + +if MED_ENABLE_METIS + salomeinclude_HEADERS+= MEDPARTITIONER_METISGraph.hxx +endif +if MED_ENABLE_SCOTCH + salomeinclude_HEADERS+= MEDPARTITIONER_SCOTCHGraph.hxx +endif + +dist_libmedpartitioner_la_SOURCES= \ +MEDPARTITIONER_utils.cxx \ +MEDPARTITIONER_MESHCollection.cxx \ +MEDPARTITIONER_MESHCollectionDriver.cxx \ +MEDPARTITIONER_MESHCollectionMedXMLDriver.cxx \ +MEDPARTITIONER_MESHCollectionMedAsciiDriver.cxx \ +MEDPARTITIONER_ParallelTopology.cxx \ +MEDPARTITIONER_Graph.cxx\ +MEDPARTITIONER_UserGraph.cxx\ +MEDPARTITIONER_ParaDomainSelector.cxx \ +MEDPARTITIONER_JointFinder.cxx \ +MEDPARTITIONER_SkyLineArray.cxx \ +MEDPARTITIONER_ConnectZone.cxx + +if MED_ENABLE_METIS + dist_libmedpartitioner_la_SOURCES+= MEDPARTITIONER_METISGraph.cxx +endif +if MED_ENABLE_SCOTCH + dist_libmedpartitioner_la_SOURCES+= MEDPARTITIONER_SCOTCHGraph.cxx +endif + +libmedpartitioner_la_CPPFLAGS= $(MPI_INCLUDES) $(MED2_INCLUDES) $(HDF5_INCLUDES) @CXXTMPDPTHFLAGS@ \ + $(BOOST_CPPFLAGS) $(LIBXML_INCLUDES) \ + -I$(srcdir)/../INTERP_KERNEL/Bases -I$(srcdir)/../MEDCoupling \ + -I$(srcdir)/../MEDLoader -I$(srcdir)/../INTERP_KERNEL + +libmedpartitioner_la_LDFLAGS= +#libmedpartitioner_la_LDFLAGS= $(MED2_LIBS) $(HDF5_LIBS) $(STDLIB) $(LIBXML_LIBS) \ +# ../MEDMEM/libmedmem.la ../MEDWrapper/V2_1/Core/libmed_V2_1.la + +if MED_ENABLE_PARMETIS + libmedpartitioner_la_CPPFLAGS+= $(PARMETIS_CPPFLAGS) + libmedpartitioner_la_LDFLAGS+= $(PARMETIS_LIBS) +endif +if MED_ENABLE_METIS + libmedpartitioner_la_CPPFLAGS+= $(METIS_CPPFLAGS) + libmedpartitioner_la_LDFLAGS+= $(METIS_LIBS) +endif +if MED_ENABLE_SCOTCH + libmedpartitioner_la_CPPFLAGS+= $(SCOTCH_CPPFLAGS) + libmedpartitioner_la_LDFLAGS+= $(SCOTCH_LIBS) +endif +if MED_ENABLE_KERNEL + libmedpartitioner_la_CPPFLAGS+= ${KERNEL_CXXFLAGS} + libmedpartitioner_la_LDFLAGS+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace +endif + +#libmedpartitioner_la_LDFLAGS+= $(MED2_LIBS) $(HDF5_LIBS) $(STDLIB) $(LIBXML_LIBS) $(MPI_LIBS) \ +# ../MEDMEM/libmedmem.la ../INTERP_KERNEL/libinterpkernel.la ../MEDCoupling/libmedcoupling.la ../MEDLoader/libmedloader.la +libmedpartitioner_la_LDFLAGS+= $(MED2_LIBS) $(HDF5_LIBS) $(STDLIB) $(LIBXML_LIBS) $(MPI_LIBS) \ + ../INTERP_KERNEL/libinterpkernel.la ../MEDCoupling/libmedcoupling.la ../MEDLoader/libmedloader.la + +# Executables targets +bin_PROGRAMS= medpartitioner + +dist_medpartitioner_SOURCES= medpartitioner.cxx + +medpartitioner_CPPFLAGS= $(libmedpartitioner_la_CPPFLAGS) +medpartitioner_LDADD= $(libmedpartitioner_la_LDFLAGS) -lm $(BOOST_LIBS) libmedpartitioner.la +if MED_ENABLE_KERNEL + medpartitioner_LDADD+= -lSALOMEBasics +endif + +if MPI_IS_OK + bin_PROGRAMS+=medpartitioner_para + dist_medpartitioner_para_SOURCES= medpartitioner_para.cxx + medpartitioner_para_CPPFLAGS= $(medpartitioner_CPPFLAGS) + medpartitioner_para_LDADD= $(medpartitioner_LDADD) +endif + +OBSOLETE_FILES = \ + MEDPARTITIONER_SequentialTopology.cxx \ + test_HighLevelAPI.cxx + +EXTRA_DIST += $(OBSOLETE_FILES) diff --git a/src/MEDPartitioner/Test/MEDPARTITIONERTest.cxx b/src/MEDPartitioner/Test/MEDPARTITIONERTest.cxx new file mode 100644 index 000000000..7b8bde6cc --- /dev/null +++ b/src/MEDPartitioner/Test/MEDPARTITIONERTest.cxx @@ -0,0 +1,1109 @@ +// Copyright (C) 2007-2011 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + + +#include +#include +#include +#include +#include +#include + +//#include "MEDPARTITIONERTest.hxx" +#include + +//#include "MEDMEM_STRING.hxx" + +#include "MEDPARTITIONERTest.hxx" + +#include "CellModel.hxx" +#include "MEDFileMesh.hxx" +#include "MEDLoader.hxx" +#include "MEDLoaderBase.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingExtrudedMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" +#include "MEDCouplingMemArray.hxx" +#include "MEDCouplingMultiFields.hxx" + +#include "MEDPARTITIONER_MESHCollection.hxx" +//#include "MEDPARTITIONER_Topology.hxx" +#include "MEDPARTITIONER_ParallelTopology.hxx" +#include "MEDPARTITIONER_ParaDomainSelector.hxx" + +#include "MEDPARTITIONER_utils.hxx" + +#ifdef HAVE_MPI2 +#include +#endif + +using namespace std; +//using namespace MEDPARTITIONER; +using namespace ParaMEDMEM; +using namespace MEDPARTITIONER; + +/*string intToStr(int i) +{ + ostringstream oss; + oss<_ni=ni; //nb of hexa9 + this->_nj=nj; + this->_nk=nk; + this->_ntot=_ni*_nj*_nk; + string ijk=intToStr(ni)+"x"+intToStr(nj)+"x"+intToStr(nk); + this->_fileName="tmp_testMesh_"+ijk+".med"; + this->_fileNameWithFaces="tmp_testMeshWithFaces_"+ijk+".med"; + string ij=intToStr(ni)+"x"+intToStr(nj); + this->_fileName2="tmp_testMesh_"+ij+".med"; + this->_meshName="testMesh"; +} + +void MEDPARTITIONERTest::setSmallSize() +{ + setSize(2,3,5); //nb of hexa9 +} + +void MEDPARTITIONERTest::setMedianSize() +{ + setSize(20,30,50); //nb of hexa9 +} + +void MEDPARTITIONERTest::setbigSize() +{ + setSize(200,300,500); //nb of hexa9 +} + + +// ============================================================================ +/*! + * Set up the environment called at every CPPUNIT_TEST () + */ +// ============================================================================ +void MEDPARTITIONERTest::setUp() +{ + this->_verbose=0; +} + +// ============================================================================ +/*! + * - delete + */ +// ============================================================================ +void MEDPARTITIONERTest::tearDown() +{ +} + +ParaMEDMEM::MEDCouplingUMesh * MEDPARTITIONERTest::buildCUBE3DMesh() +//only hexa8 +{ + vector conn; + vector coor; + for (int k=0; k<=_nk; k++) + for (int j=0; j<=_nj; j++) + for (int i=0; i<=_ni; i++) + { + coor.push_back(i+.1); + coor.push_back(j+.2); + coor.push_back(k+.3); + } + int ii; + for (int k=0; k<_nk; k++) + for (int j=0; j<_nj; j++) + for (int i=0; i<_ni; i++) + { + ii=i + j*(_ni+1) + k*(_ni+1)*(_nj+1); + conn.push_back(ii); + conn.push_back(ii+1); + ii=ii + _ni + 2 ; + conn.push_back(ii); + conn.push_back(ii-1); + + ii=i + j*(_ni+1) + (k+1)*(_ni+1)*(_nj+1); + conn.push_back(ii); + conn.push_back(ii+1); + ii=ii + _ni + 2 ; + conn.push_back(ii); + conn.push_back(ii-1); + } + + if (false) //(_verbose) + { + cout<<"\nnb coor "<<(_ni+1)*(_nj+1)*(_nk+1)*3<<" "<setMeshDimension(3); + int nbc=conn.size()/8; //nb of cells + int nbv=coor.size()/3; //nb of vertices + mesh->allocateCells(nbc); + for(int i=0; iinsertNextCell(INTERP_KERNEL::NORM_HEXA8,8,onehexa); + } + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(nbv,3); + std::copy(coor.begin(),coor.end(),myCoords->getPointer()); + mesh->setCoords(myCoords); + mesh->setName(_meshName.c_str()); + myCoords->decrRef(); + mesh->checkCoherency(); + return mesh; +} + +ParaMEDMEM::MEDCouplingUMesh * MEDPARTITIONERTest::buildCARRE3DMesh() +//only quad4 in oblique (k=j) +{ + vector conn; + vector coor; + for (int j=0; j<=_nj; j++) + for (int i=0; i<=_ni; i++) + { + int k=j; + coor.push_back(i+.1); + coor.push_back(j+.2); + coor.push_back(k+.3); + } + int ii; + int k=0; + for (int j=0; j<_nj; j++) + for (int i=0; i<_ni; i++) + { + ii=i + j*(_ni+1) + k*(_ni+1)*(_nj+1); + conn.push_back(ii); + conn.push_back(ii+1); + ii=ii + _ni + 2 ; + conn.push_back(ii); + conn.push_back(ii-1); + } + + if (false) //(_verbose) + { + cout<<"\nnb coor "<<(_ni+1)*(_nj+1)*3<<" "<setMeshDimension(2); + int nbc=conn.size()/4; //nb of cells + int nbv=coor.size()/3; //nb of vertices + mesh->allocateCells(nbc); + for(int i=0; iinsertNextCell(INTERP_KERNEL::NORM_QUAD4,4,onequa); + } + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(nbv,3); + std::copy(coor.begin(),coor.end(),myCoords->getPointer()); + mesh->setCoords(myCoords); + mesh->setName(_meshName.c_str()); + myCoords->decrRef(); + mesh->checkCoherency(); + return mesh; +} + +ParaMEDMEM::MEDCouplingUMesh * MEDPARTITIONERTest::buildFACE3DMesh() +//only quad4 on a global face of the CUBE3D (k=0) +{ + vector conn; + vector coor; + for (int j=0; j<=_nj; j++) + for (int i=0; i<=_ni; i++) + { + int k=0; + coor.push_back(i+.1); + coor.push_back(j+.2); + coor.push_back(k+.3); + } + int ii; + int k=0; + for (int j=0; j<_nj; j++) + for (int i=0; i<_ni; i++) + { + ii=i + j*(_ni+1) + k*(_ni+1)*(_nj+1); + conn.push_back(ii); + conn.push_back(ii+1); + ii=ii + _ni + 2 ; + conn.push_back(ii); + conn.push_back(ii-1); + } + + if (false) //(_verbose) + { + cout<<"\nnb coor "<<(_ni+1)*(_nj+1)*3<<" "<setMeshDimension(2); + int nbc=conn.size()/4; //nb of cells + int nbv=coor.size()/3; //nb of vertices + mesh->allocateCells(nbc); + for(int i=0; iinsertNextCell(INTERP_KERNEL::NORM_QUAD4,4,onequa); + } + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(nbv,3); + std::copy(coor.begin(),coor.end(),myCoords->getPointer()); + mesh->setCoords(myCoords); + mesh->setName(_meshName.c_str()); + myCoords->decrRef(); + mesh->checkCoherency(); + return mesh; +} + +MEDCouplingFieldDouble * MEDPARTITIONERTest::buildVecFieldOnCells(string myfileName) +{ + //int ni=2,nj=3,nk=5; //nb of hexa9 + vector field; + for (int k=0; k<_nk; k++) + for (int j=0; j<_nj; j++) + for (int i=0; i<_ni; i++) + { + field.push_back(i+.1); + field.push_back(j+.2); + field.push_back(k+.3); + } + //cvwat + MEDCouplingUMesh *mesh=MEDLoader::ReadUMeshFromFile(myfileName.c_str(),_meshName.c_str(),0); + int nbOfCells=mesh->getNumberOfCells(); + MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME); + f1->setName("VectorFieldOnCells"); + f1->setDescription("DescriptionOfFieldOnCells"); //not saved in file? + f1->setMesh(mesh); + DataArrayDouble *myField=DataArrayDouble::New(); + myField->alloc(nbOfCells,3); + std::copy(field.begin(),field.end(),myField->getPointer()); + f1->setArray(myField); + myField->setInfoOnComponent(0,"vx"); + myField->setInfoOnComponent(1,"vy"); + myField->setInfoOnComponent(2,"vz"); + myField->decrRef(); + f1->setTime(2.,0,1); + f1->checkCoherency(); + mesh->decrRef(); + return f1; +} + +MEDCouplingFieldDouble * MEDPARTITIONERTest::buildVecFieldOnNodes() +{ + //int ni=2,nj=3,nk=5; //nb of hexa9 + vector field; + for (int k=0; k<=_nk; k++) + for (int j=0; j<=_nj; j++) + for (int i=0; i<=_ni; i++) + { + field.push_back(i+.1); + field.push_back(j+.2); + field.push_back(k+.3); + } + + MEDCouplingUMesh *mesh=MEDLoader::ReadUMeshFromFile(_fileName.c_str(),_meshName.c_str(),0); + int nbOfNodes=mesh->getNumberOfNodes(); + MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_NODES,ONE_TIME); + f1->setName("VectorFieldOnNodes"); + f1->setDescription("DescriptionOfFieldOnNodes"); //not saved in file? + f1->setMesh(mesh); + DataArrayDouble *myField=DataArrayDouble::New(); + myField->alloc(nbOfNodes,3); + std::copy(field.begin(),field.end(),myField->getPointer()); + f1->setArray(myField); + myField->setInfoOnComponent(0,"vx"); + myField->setInfoOnComponent(1,"vy"); + myField->setInfoOnComponent(2,"vz"); + myField->decrRef(); + f1->setTime(2.,0,1); + f1->checkCoherency(); + mesh->decrRef(); + return f1; +} + + +void MEDPARTITIONERTest::createTestMeshWithoutField() +{ + { + MEDCouplingUMesh * mesh = buildCUBE3DMesh(); + MEDLoader::WriteUMesh(_fileName.c_str(),mesh,true); + if (_verbose) cout<getName(),0); + if (_verbose) cout<<_fileName<<" reread"<isEqual(mesh_rw,1e-12)); + mesh_rw->decrRef(); + } + mesh->decrRef(); + } + + { + vector meshes; + MEDCouplingUMesh * mesh1 = buildCUBE3DMesh(); + MEDCouplingUMesh * mesh2 = buildFACE3DMesh(); + mesh1->setName("testMesh"); + mesh2->setName("theFaces"); + mesh2->tryToShareSameCoordsPermute(*mesh1, 1e-9); + mesh2->checkCoherency(); + mesh1->checkCoherency(); + meshes.push_back(mesh1); + meshes.push_back(mesh2); + MEDLoader::WriteUMeshes(_fileNameWithFaces.c_str(), meshes, true); + + ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(_fileNameWithFaces.c_str(), mesh1->getName()); + DataArrayInt* FacesFam=DataArrayInt::New(); + FacesFam->alloc(mfm->getSizeAtLevel(-1),1); + FacesFam->fillWithValue(-1); + DataArrayInt* CellsFam=DataArrayInt::New(); + CellsFam->alloc(mfm->getSizeAtLevel(0),1); + CellsFam->fillWithValue(1); + mfm->setFamilyFieldArr(-1,FacesFam); + mfm->setFamilyFieldArr(0,CellsFam); + map theFamilies; + theFamilies["FAMILLE_ZERO"]=0; + theFamilies["FamilyFaces"]=-1; + theFamilies["FamilyCells"]=1; + map > theGroups; + theGroups["GroupFaces"].push_back("FamilyFaces"); + theGroups["GroupCells"].push_back("FamilyCells"); + mfm->setFamilyInfo(theFamilies); + mfm->setGroupInfo(theGroups); + mfm->write(_fileNameWithFaces.c_str(),0); + FacesFam->decrRef(); + CellsFam->decrRef(); + + /*ce truc marche pas! + ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(_fileNameWithFaces.c_str(), mesh1->getName()); + vector ms; + ms.push_back(mesh2); + mfm->setGroupsFromScratch(-1, ms); + mfm->write(_fileNameWithFaces.c_str(),0); + */ + + if (_verbose) cout<getName(),0); + if (_verbose) cout<<_fileNameWithFaces<<" reread"<isEqual(mesh_rw,1e-12)); + mesh_rw->decrRef(); + } + mesh1->decrRef(); + mesh2->decrRef(); + } + + { + MEDCouplingUMesh * mesh = buildCARRE3DMesh(); + MEDLoader::WriteUMesh(_fileName2.c_str(),mesh,true); + if (_verbose) cout<getName(),0); + if (_verbose) cout<<_fileName2<<" reread"<isEqual(mesh_rw,1e-12)); + mesh_rw->decrRef(); + mesh->decrRef(); + } +} + +void MEDPARTITIONERTest::createHugeTestMesh(int ni, int nj, int nk, int nbx, int nby, int nbz, int nbTarget) +{ + setSize(ni,nj,nk); + _nbTargetHuge=nbTarget; + MEDCouplingUMesh * mesh = buildCUBE3DMesh(); + //int nbx=1, nby=1, nbz=2; + std::vector< double > cooDep,cooFin; + mesh->getCoordinatesOfNode(0, cooDep); + mesh->getCoordinatesOfNode(mesh->getNumberOfNodes()-1, cooFin); + //cout<\n \ +\n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n \ + \n$tagSubfile \ + \n \ + \n$tagMesh \ + \n \ +\n"; + + string tagSubfiles, tagSubfile="\ + \n \ + $fileName\n \ + localhost\n \ + \n"; + string tagMeshes, tagMesh="\ + \n \ + \n \ + testMesh\n \ + \n \ + \n"; + + int xyz=1; + string sxyz; + DataArrayDouble* coordsInit=mesh->getCoords()->deepCpy(); + double* ptrInit=coordsInit->getPointer(); + double deltax=cooFin[0]-cooDep[0]; + double deltay=cooFin[1]-cooDep[1]; + double deltaz=cooFin[2]-cooDep[2]; + + double dz=0.; + for (int z=0; zgetCoords(); + //int nbOfComp=coords->getNumberOfComponents(); //be 3D + int nbOfTuple=coords->getNumberOfTuples(); + double* ptr=coords->getPointer(); + double* ptrini=ptrInit; + for (int i=0; idecrRef(); + + tagXml.replace(tagXml.find("$subdomainNumber"),16,sxyz); + tagXml.replace(tagXml.find("$tagSubfile"),11,tagSubfiles); + tagXml.replace(tagXml.find("$tagMesh"),8,tagMeshes); + + string nameFileXml; + _fileNameHugeXml="tmp_testMeshHuge_"+intToStr(_ni)+"x"+intToStr(_nj)+"x"+intToStr(_nk)+"_"+sxyz+".xml"; + std::ofstream f(_fileNameHugeXml.c_str()); + f<setTime(3.,1,1); //time,it,order + f1->applyFunc("x/2."); + MEDLoader::WriteField(name.c_str(),f1,false); + if (_verbose) cout<getMesh()->getName(),0,f1->getName(),0,1); + //DataArrayDouble *res=f2->getArray(); + if (_verbose) cout<isEqual(f2,1e-12,1e-12)); + f2->decrRef(); + } + f1->decrRef(); + } + { + string name=_fileName; + MEDCouplingFieldDouble *f1=buildVecFieldOnCells(name); + name.replace(name.find(".med"),4,"_WithVecFieldOnGaussNe.med"); + MEDCouplingFieldDouble *f3=MEDCouplingFieldDouble::New(ON_GAUSS_NE,ONE_TIME); + f3->setMesh(f1->getMesh()); + //cout<<"\nNumberOfMeshPlacesExpected "<getNumberOfMeshPlacesExpected()<<" " + // /*<getMesh())<<" "*/ + // <getMesh()->getNumberOfNodes()<<" " + // <getMesh()->getNumberOfCells()<setName("MyFieldOnGaussNE"); + f3->setDescription("MyDescriptionNE"); + DataArrayDouble *array=DataArrayDouble::New(); + //int nb=f1->getMesh()->getNumberOfNodes(); + + /*8 pt de gauss by cell + int nb=f3->getMesh()->getNumberOfCells()*8; + array->alloc(nb,2); + double *ptr=array->getPointer(); + for (int i=0; igetMesh()->getNumberOfCells(); + int nb=nbcell*nbptgauss; + int nbcomp=2; + array->alloc(nb,nbcomp); + double *ptr=array->getPointer(); + int ii=0; + for (int i=0; isetInfoOnComponent(0,"vGx"); + array->setInfoOnComponent(1,"vGy"); + f3->setTime(4.,5,6); + f3->setArray(array); + array->decrRef(); + MEDLoader::WriteField(name.c_str(),f3,true); + if (_verbose) cout<checkCoherency(); + f1->decrRef(); + if (_ntot<1000000) //too long + { + MEDCouplingFieldDouble* f4=MEDLoader::ReadField(ON_GAUSS_NE, + name.c_str(), f3->getMesh()->getName(), 0, "MyFieldOnGaussNE", 5, 6); + if (_verbose) cout<<"MyFieldOnGaussNE reread"<decrRef(); + } + f3->decrRef(); + } + { + string name=_fileNameWithFaces; + MEDCouplingFieldDouble *f1=buildVecFieldOnCells(name); + name.replace(name.find(".med"),4,"_WithVecFieldOnCells.med"); + MEDLoader::WriteField(name.c_str(),f1,true); + if (_verbose) cout<getMesh()->getName(),0,f1->getName(),0,1); + if (_verbose) cout<isEqual(f2,1e-12,1e-12)); assertion failed!! + f2->decrRef(); + } + f1->decrRef(); + } +} + +void MEDPARTITIONERTest::createTestMeshWithVecFieldOnNodes() +{ + MEDCouplingFieldDouble *f1=buildVecFieldOnNodes(); + string name=_fileName; + name.replace(name.find(".med"),4,"_WithVecFieldOnNodes.med"); + MEDLoader::WriteField(name.c_str(),f1,true); + if (_verbose) cout<getMesh()->getName(),0,f1->getName(),0,1); + if (_verbose) cout<isEqual(f2,1e-12,1e-12)); assertion failed!! + f2->decrRef(); + } + f1->decrRef(); +} + +void MEDPARTITIONERTest::verifyTestMeshWithVecFieldOnNodes() +{ + string name=_fileName; + name.replace(name.find(".med"),4,"_WithVecFieldOnNodes.med"); + MEDCouplingUMesh * m=MEDLoader::ReadUMeshFromFile(name.c_str(),_meshName.c_str(),0); + const std::set& types=m->getAllTypes(); + if (_verbose) + { + cout<<"\n types in "<::iterator t=types.begin(); t!=types.end(); ++t) cout<<" "<<*t; + for (std::set::iterator t=types.begin(); t!=types.end(); ++t) + { + //INTERP_KERNEL::CellModel essai=INTERP_KERNEL::CellModel::GetCellModel(*t); + cout<<" "<<(INTERP_KERNEL::CellModel::GetCellModel(*t)).getRepr(); + } + cout<decrRef(); + + MEDFileUMesh * mf; + mf->New(_fileName.c_str(),_meshName.c_str(),-1,-1); + vector lev; + lev=mf->getNonEmptyLevels(); + if (_verbose) + { + cout<<" levels in "<::iterator l=lev.begin(); l!=lev.end(); ++l) cout<<" "<<*l; + cout<decrRef(); +} + +void MEDPARTITIONERTest::verifyMedpartitionerOnSmallSizeForMesh() +{ + int res; + string fileName,cmd,execName,sourceName,targetName,input; + execName=getenv("MED_ROOT_DIR"); //.../INSTALL/MED + execName+="/bin/salome/medpartitioner_para"; + fileName=_fileNameWithFaces; + + ParaMEDMEM::MEDFileUMesh* initialMesh=ParaMEDMEM::MEDFileUMesh::New(fileName.c_str(),_meshName.c_str()); + ParaMEDMEM::MEDCouplingUMesh* cellMesh=initialMesh->getLevel0Mesh(false); + ParaMEDMEM::MEDCouplingUMesh* faceMesh=initialMesh->getLevelM1Mesh(false); + + cmd="mpirun -np 5 "+execName+" --ndomains=5 --split-method=metis"; //on same proc + sourceName=fileName; + targetName=fileName; + targetName.replace(targetName.find(".med"),4,"_partitionedTo5_"); + cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+intToStr(_verbose); + if (_verbose) cout<cellMeshes=collection.getMesh(); + CPPUNIT_ASSERT_EQUAL(5, (int) cellMeshes.size()); + int nbcells=0; + for (int i = 0; i < cellMeshes.size(); i++) nbcells+=cellMeshes[i]->getNumberOfCells(); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), nbcells); + + std::vectorfaceMeshes=collection.getFaceMesh(); + CPPUNIT_ASSERT_EQUAL(5, (int) faceMeshes.size()); + int nbfaces=0; + for (int i = 0; i < faceMeshes.size(); i++) nbfaces+=faceMeshes[i]->getNumberOfCells(); + CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), nbfaces); + + //merge split meshes and test equality + cmd="mpirun -np 1 "+execName+" --ndomains=1 --split-method=metis"; //on same proc + sourceName=targetName+".xml"; + targetName=fileName; + targetName.replace(targetName.find(".med"),4,"_remergedFrom5_"); + cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+intToStr(_verbose); + if (_verbose) cout<getLevel0Mesh(false); + ParaMEDMEM::MEDCouplingUMesh* refusedFaceMesh=refusedMesh->getLevelM1Mesh(false); + + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), refusedCellMesh->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), refusedFaceMesh->getNumberOfCells()); + + /*not the good job + ParaMEDMEM::MEDCouplingMesh* mergeCell=cellMesh->mergeMyselfWith(refusedCellMesh); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), mergeCell->getNumberOfCells()); + + ParaMEDMEM::MEDCouplingMesh* mergeFace=faceMesh->mergeMyselfWith(refusedFaceMesh); + CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), mergeFace->getNumberOfCells()); + + CPPUNIT_ASSERT(faceMesh->isEqual(refusedFaceMesh,1e-12)); + */ + + std::vector meshes; + std::vector corr; + meshes.push_back(cellMesh); + refusedCellMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-9); + meshes.push_back(refusedCellMesh); + MEDCouplingUMesh* fusedCell=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), fusedCell->getNumberOfCells()); + + meshes.resize(0); + for (int i = 0; i < corr.size(); i++) corr[i]->decrRef(); + corr.resize(0); + meshes.push_back(faceMesh); + refusedFaceMesh->tryToShareSameCoordsPermute(*faceMesh, 1e-9); + meshes.push_back(refusedFaceMesh); + MEDCouplingUMesh* fusedFace=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr); + CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), fusedFace->getNumberOfCells()); + + for (int i = 0; i < corr.size(); i++) corr[i]->decrRef(); + fusedFace->decrRef(); + refusedFaceMesh->decrRef(); + faceMesh->decrRef(); + fusedCell->decrRef(); + refusedCellMesh->decrRef(); + cellMesh->decrRef(); + //done in ~collection + //for (int i = 0; i < faceMeshes.size(); i++) faceMeshes[i]->decrRef(); + //for (int i = 0; i < cellMeshes.size(); i++) cellMeshes[i]->decrRef(); +} + +void MEDPARTITIONERTest::verifyMedpartitionerOnSmallSizeForFieldOnCells() +{ + int res; + string fileName,cmd,execName,sourceName,targetName,input; + execName=getenv("MED_ROOT_DIR"); //.../INSTALL/MED + execName+="/bin/salome/medpartitioner_para"; + fileName=_fileName; + fileName.replace(fileName.find(".med"),4,"_WithVecFieldOnCells.med"); + + ParaMEDMEM::MEDFileUMesh* initialMesh=ParaMEDMEM::MEDFileUMesh::New(fileName.c_str(),_meshName.c_str()); + ParaMEDMEM::MEDCouplingUMesh* cellMesh=initialMesh->getLevel0Mesh(false); + + cmd="mpirun -np 5 "+execName+" --ndomains=5 --split-method=metis"; //on same proc + sourceName=fileName; + targetName=fileName; + targetName.replace(targetName.find(".med"),4,"_partitionedTo5_"); + cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+intToStr(_verbose); + if (_verbose) cout<getLevel0Mesh(false); + + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), refusedCellMesh->getNumberOfCells()); + + std::vector meshes; + std::vector corr; + meshes.push_back(cellMesh); + refusedCellMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-9); + meshes.push_back(refusedCellMesh); + MEDCouplingUMesh* fusedCell=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), fusedCell->getNumberOfCells()); + + MEDCouplingFieldDouble* field1=MEDLoader::ReadFieldCell(fileName.c_str(),initialMesh->getName(),0,"VectorFieldOnCells",0,1); + MEDCouplingFieldDouble* field2=MEDLoader::ReadFieldCell(refusedName.c_str(),refusedCellMesh->getName(),0,"VectorFieldOnCells",0,1); + + int nbcells=corr[1]->getNumberOfTuples(); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), nbcells); + //use corr to test equality of field + DataArrayDouble* f1=field1->getArray(); + DataArrayDouble* f2=field2->getArray(); + if (_verbose) + { + cout<<"\nf1 : "<reprZip(); + cout<<"\nf2 : "<reprZip(); //field2->advancedRepradvancedRepr(); + for (int i = 0; i < corr.size(); i++) cout<<"\ncorr "<reprZip(); + + } + int nbequal=0; + int nbcomp=field1->getNumberOfComponents(); + double* p1=f1->getPointer(); + double* p2=f2->getPointer(); + int* pc=corr[1]->getPointer(); + for (int i = 0; i < nbcells; i++) + { + int i1=pc[i]*nbcomp; + int i2=i*nbcomp; + for (int j = 0; j < nbcomp; j++) + { + if (p1[i1+j]==p2[i2+j]) nbequal++; + //cout<<" "<decrRef(); + field1->decrRef(); + field2->decrRef(); + fusedCell->decrRef(); + refusedCellMesh->decrRef(); + cellMesh->decrRef(); +} + +void MEDPARTITIONERTest::verifyMedpartitionerOnSmallSizeForFieldOnGaussNe() +{ + int res; + string fileName,cmd,execName,sourceName,targetName,input; + execName=getenv("MED_ROOT_DIR"); //.../INSTALL/MED + execName+="/bin/salome/medpartitioner_para"; + fileName=_fileName; + fileName.replace(fileName.find(".med"),4,"_WithVecFieldOnGaussNe.med"); + + ParaMEDMEM::MEDFileUMesh* initialMesh=ParaMEDMEM::MEDFileUMesh::New(fileName.c_str(),_meshName.c_str()); + ParaMEDMEM::MEDCouplingUMesh* cellMesh=initialMesh->getLevel0Mesh(false); + + cmd="mpirun -np 5 "+execName+" --ndomains=5 --split-method=metis"; //on same proc + sourceName=fileName; + targetName=fileName; + targetName.replace(targetName.find(".med"),4,"_partitionedTo5_"); + cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+intToStr(_verbose); + if (_verbose) cout<getLevel0Mesh(false); + + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), refusedCellMesh->getNumberOfCells()); + + std::vector meshes; + std::vector corr; + meshes.push_back(cellMesh); + refusedCellMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-9); + meshes.push_back(refusedCellMesh); + MEDCouplingUMesh* fusedCell=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), fusedCell->getNumberOfCells()); + + MEDCouplingFieldDouble* field1=MEDLoader::ReadField(ON_GAUSS_NE,fileName.c_str(),initialMesh->getName(),0,"MyFieldOnGaussNE",5,6); + MEDCouplingFieldDouble* field2=MEDLoader::ReadField(ON_GAUSS_NE,refusedName.c_str(),refusedCellMesh->getName(),0,"MyFieldOnGaussNE",5,6); + + int nbcells=corr[1]->getNumberOfTuples(); + CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), nbcells); + //use corr to test equality of field + DataArrayDouble* f1=field1->getArray(); + DataArrayDouble* f2=field2->getArray(); + if (_verbose) + { + cout<<"\nf1 : "<reprZip(); //123.4 for 12th cell,3rd component, 4th gausspoint + cout<<"\nf2 : "<reprZip(); //field2->advancedRepradvancedRepr(); + for (int i = 0; i < corr.size(); i++) cout<<"\ncorr "<reprZip(); + + } + int nbequal=0; + int nbptgauss=8; + int nbcomp=field1->getNumberOfComponents(); + double* p1=f1->getPointer(); + double* p2=f2->getPointer(); + int* pc=corr[1]->getPointer(); + for (int i = 0; i < nbcells; i++) + { + int i1=pc[i]*nbcomp*nbptgauss; + int i2=i*nbcomp*nbptgauss; + for (int j = 0; j < nbcomp*nbptgauss; j++) + { + if (p1[i1+j]==p2[i2+j]) nbequal++; + //cout<<" "<decrRef(); + field1->decrRef(); + field2->decrRef(); + fusedCell->decrRef(); + refusedCellMesh->decrRef(); + cellMesh->decrRef(); +} + +void MEDPARTITIONERTest::createTestMeshes() +{ + createTestMeshWithoutField(); + createTestMeshWithVecFieldOnCells(); + createTestMeshWithVecFieldOnNodes(); +} + + +void MEDPARTITIONERTest::deleteTestMeshes() +{ + string cmd="rm *tmp_testMesh*"; + if (_verbose) cout</dev/null 1>/dev/null"); //no trace + CPPUNIT_ASSERT_EQUAL(0, res); + + execName=getenv("MED_ROOT_DIR"); //.../INSTALL/MED + execName+="/bin/salome/medpartitioner_para"; + + cmd="which "+execName+" 2>/dev/null 1>/dev/null"; //no trace + res=system(cmd.c_str()); + CPPUNIT_ASSERT_EQUAL(0, res); + + cmd="mpirun -np 2 "+execName+" --ndomains=2 --split-method=metis"; //on same proc + sourceName=_fileName; + targetName=_fileName; + targetName.replace(targetName.find(".med"),4,"_partitionedTo2_"); + cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+intToStr(_verbose); + if (_verbose) cout< +namespace po=boost::program_options; +#endif + +#include +#include +#include +#include +#include + +#include "MEDPARTITIONER_Graph.hxx" +#include "MEDPARTITIONER_MESHCollection.hxx" +#include "MEDPARTITIONER_Topology.hxx" + +using namespace std; + +int main(int argc, char** argv) +{ +#ifndef ENABLE_METIS +#ifndef ENABLE_SCOTCH + cout << "Sorry, no one split method is available. Please, compile with METIS or SCOTCH."<(),"name of the input MED file") + ("output-file",po::value(),"name of the resulting file") + ("meshname",po::value(),"name of the input mesh") +#ifdef ENABLE_METIS +#ifdef ENABLE_SCOTCH + ("split-method",po::value(&library)->default_value("metis"),"name of the splitting library (metis,scotch)") +#endif +#endif + ("ndomains",po::value(&ndomains)->default_value(1),"number of subdomains in the output file") + ("plain-master","creates a plain masterfile instead of an XML file") + ("creates-boundary-faces","creates the necessary faces so that faces joints are created in the output files") + ("family-splitting","preserves the family names instead of focusing on the groups") + ("empty-groups","creates empty groups in zones that do not contain a group from the original domain"); + + po::variables_map vm; + po::store(po::parse_command_line(argc,argv,desc),vm); + po::notify(vm); + + if (vm.count("help")) + { + cout<(); + if (!vm.count("input-file") || !vm.count("output-file")) + { + cout << "input-file and output-file names must be specified"<(); + output = vm["output-file"].as(); + + if (vm.count("mesh-only")) + mesh_only=true; + + if (vm.count("distributed")) + is_sequential=false; + + if (is_sequential) + meshname = vm["meshname"].as(); + + if (vm.count("plain-master")) + xml_output_master=false; + + if (vm.count("creates-boundary-faces")) + creates_boundary_faces=true; + + if (vm.count("split-families")) + split_families=true; + + if (vm.count("empty-groups")) + empty_groups=true; + +#else // BOOST_PROGRAM_OPTIONS_LIB + + // Primitive parsing of command-line options + + string desc ("Available options of medpartitioner V1.0:\n" + "\t--help : produces this help message\n" + "\t--mesh-only : do not create the fields contained in the original file(s)\n" + "\t--distributed : specifies that the input file is distributed\n" + "\t--input-file= : name of the input MED file\n" + "\t--output-file= : name of the resulting file\n" + "\t--meshname= : name of the input mesh (not used with --distributed option)\n" + "\t--ndomains= : number of subdomains in the output file, default is 1\n" +#ifdef ENABLE_METIS +#ifdef ENABLE_SCOTCH + "\t--split-method=: name of the splitting library (metis/scotch), default is metis\n" +#endif +#endif + "\t--plain-master : creates a plain masterfile instead of an XML file\n" + "\t--creates-boundary-faces: creates the necessary faces so that faces joints are created in the output files\n" + "\t--family-splitting : preserves the family names instead of focusing on the groups\n" + "\t--empty-groups : creates empty groups in zones that do not contain a group from the original domain" + ); + + if (argc < 4) { + cout << desc.c_str() << endl; + return 1; + } + + for (int i = 1; i < argc; i++) { + if (strlen(argv[i]) < 3) { + cout << desc.c_str() << endl; + return 1; + } + + if (strncmp(argv[i],"--m",3) == 0) { + if (strcmp(argv[i],"--mesh-only") == 0) { + mesh_only = true; + cout << "\tmesh-only = " << mesh_only << endl; // tmp + } + else if (strlen(argv[i]) > 11) { // "--meshname=" + meshname = (argv[i] + 11); + cout << "\tmeshname = " << meshname << endl; // tmp + } + } + else if (strncmp(argv[i],"--d",3) == 0) { + is_sequential = false; + cout << "\tis_sequential = " << is_sequential << endl; // tmp + } + else if (strncmp(argv[i],"--i",3) == 0) { + if (strlen(argv[i]) > 13) { // "--input-file=" + input = (argv[i] + 13); + cout << "\tinput-file = " << input << endl; // tmp + } + } + else if (strncmp(argv[i],"--o",3) == 0) { + if (strlen(argv[i]) > 14) { // "--output-file=" + output = (argv[i] + 14); + cout << "\toutput-file = " << output << endl; // tmp + } + } + else if (strncmp(argv[i],"--s",3) == 0) { + if (strlen(argv[i]) > 15) { // "--split-method=" + library = (argv[i] + 15); + cout << "\tsplit-method = " << library << endl; // tmp + } + } + else if (strncmp(argv[i],"--f",3) == 0) { //"--family-splitting" + split_families=true; + cout << "\tfamily-splitting true" << endl; // tmp + } + else if (strncmp(argv[i],"--n",3) == 0) { + if (strlen(argv[i]) > 11) { // "--ndomains=" + ndomains = atoi(argv[i] + 11); + cout << "\tndomains = " << ndomains << endl; // tmp + } + } + else if (strncmp(argv[i],"--p",3) == 0) { // "--plain-master" + xml_output_master = false; + cout << "\txml_output_master = " << xml_output_master << endl; // tmp + } + else if (strncmp(argv[i],"--c",3) == 0) { // "--creates-boundary-faces" + creates_boundary_faces = true; + cout << "\tcreates_boundary_faces = " << creates_boundary_faces << endl; // tmp + } + else if (strncmp(argv[i],"--e",3) == 0) { // "--empty-groups" + empty_groups = true; + cout << "\tempty_groups = true" << endl; // tmp + } + else { + cout << desc.c_str() << endl; + return 1; + } + } + + if (is_sequential && meshname.empty()) { + cout << "Mesh name must be given for sequential(not distributed) input file." << endl; + cout << desc << endl; + return 1; + } + +#endif // BOOST_PROGRAM_OPTIONS_LIB + + + //testing whether it is possible to write a file at the specified location + string outputtest = output + ".testioms."; + ofstream testfile (outputtest.c_str()); + if (testfile.fail()) + { + cout << "MEDPARTITIONER : output-file directory does not exist or is in read-only access" << endl; + return 1; + }; + //deletes test file + remove(outputtest.c_str()); + + // Beginning of the computation + + // Loading the mesh collection + MEDPARTITIONER::MESHCollection* collection; + cout << "MEDPARTITIONER : reading input files "<createPartition(ndomains,MEDPARTITIONER::Graph::METIS); + else + new_topo = collection->createPartition(ndomains,MEDPARTITIONER::Graph::SCOTCH); + + cout << "MEDPARTITIONER : creating new meshes"< +#include +#include +#include +#include + +#ifdef HAVE_MPI2 +#include +#endif + +#ifdef BOOST_PROGRAM_OPTIONS_LIB +#include +namespace po=boost::program_options; +#endif + +using namespace std; +using namespace MEDPARTITIONER; + +int main(int argc, char** argv) +{ +#ifndef ENABLE_PARMETIS +#ifndef ENABLE_PTSCOTCH + cout << "Sorry, no one split method is available. Please, compile with ParMETIS or PT-SCOTCH."< : name of the input .med file or .xml master file\n" + "\t--output-file= : name of the resulting file (without exension)\n" + "\t--ndomains= : number of subdomains in the output file, default is 1\n" +#ifdef ENABLE_PARMETIS +#ifdef ENABLE_PTSCOTCH + "\t--split-method= : name of the splitting library (metis/scotch), default is metis\n" +#endif +#endif + "\t--creates-boundary-faces : creates boundary faces mesh in the output files\n" + //"\t--family-splitting : preserves the family names instead of focusing on the groups\n" + "\t--dump-cpu-memory : dumps passed CPU time and maximal increase of used memory\n" + //"\t--randomize= : random seed for other partitionning (only on one proc)\n" + //"\t--atomize : do the opposite of a good partitionner (only on one proc)\n" + ); + + if (argc<=1) help=1; + string value; + for (int i = 1; i < argc; i++) + { + if (strlen(argv[i]) < 3) + { + if (MyGlobals::_rank==0) cerr << "bad argument : "<< argv[i] << endl; + MPI_Finalize(); return 1; + } + + if (testArg(argv[i],"--verbose",value)) + { + MyGlobals::_verbose=1; + if (value!="") MyGlobals::_verbose = atoi(value.c_str()); + } + else if (testArg(argv[i],"--help",value)) help=1; + else if (testArg(argv[i],"--test",value)) test=1; + else if (testArg(argv[i],"--input-file",value)) input=value; + else if (testArg(argv[i],"--output-file",value)) output=value; + else if (testArg(argv[i],"--split-method",value)) library=value; + //else if (testArg(argv[i],"--family-splitting",value)) split_family=true; + else if (testArg(argv[i],"--ndomains",value)) ndomains=atoi(value.c_str()); + else if (testArg(argv[i],"--randomize",value)) MyGlobals::_randomize=atoi(value.c_str()); + else if (testArg(argv[i],"--atomize",value)) MyGlobals::_atomize=atoi(value.c_str()); + else if (testArg(argv[i],"--creates-boundary-faces",value)) MyGlobals::_creates_boundary_faces=1; + //else if (testArg(argv[i],"--empty-groups",value)) empty_groups=true; + else if (testArg(argv[i],"--dump-cpu-memory",value)) mesure_memory=true; + else + { + if (MyGlobals::_rank==0) cerr << "unknown argument : "<< argv[i] << endl; + MPI_Finalize(); return 1; + } + } + + + if (MyGlobals::_randomize!=0 && MyGlobals::_world_size!=1) + { + if (MyGlobals::_rank==0) cerr << "randomize only available on 1 proc (mpirun -np 1)" << endl; + MyGlobals::_randomize=0; + } + +#ifdef ENABLE_PARMETIS +#ifndef ENABLE_PTSCOTCH + library = "metis"; +#endif +#else + library = "scotch"; +#endif + + if (help==1) + { + if (MyGlobals::_rank==0) cout<0) cout<<"tests on "<setGlobalNumerotationDefault(collection.getParaDomainSelector()); + //int nbfiles=MyGlobals::_fileMedNames->size(); //nb domains + //to have unique valid fields names/pointers/descriptions for partitionning + collection.prepareFieldDescriptions(); + //int nbfields=collection.getFieldDescriptions().size(); //on all domains + //cout< new_topo; + if (library == "metis") //cvwat06 + new_topo.reset( collection.createPartition(ndomains,MEDPARTITIONER::Graph::METIS)); + else + new_topo.reset( collection.createPartition(ndomains,MEDPARTITIONER::Graph::SCOTCH)); + parallelizer.evaluateMemory(); + + // Creating a new mesh collection from the partitioning + if (MyGlobals::_is0verbose) cout << "Creating new meshes"< finalInformations; + vector r1,r2; + r1=allgathervVectorOfString(MyGlobals::_generalInformations); + //if (MyGlobals::_is0verbose>1000) cout << "generalInformations : \n"<0) cout<<"OK END"<< endl; + MPI_Finalize(); + return 0; + } + catch(const char *mess) + { + cerr<<"proc "<