From eb01444eb94f0ba718529166a66f6c8ee0d86395 Mon Sep 17 00:00:00 2001 From: vbd Date: Tue, 15 Sep 2009 09:59:17 +0000 Subject: [PATCH] creation of a renumbering tool that : - takes as input a MED file (with only one cell type) - as output creates another file with the same topology, but with renumbered cells so that the mesh is more suitable for numerical computation. The renumbering can be computed with reverse cuthill mckee algorithm (boost) or with nested dissection (metis) --- src/RENUMBER/Makefile.am | 92 ++++++ src/RENUMBER/RENUMBER_BOOSTRenumbering.cxx | 35 +++ src/RENUMBER/RENUMBER_BOOSTRenumbering.hxx | 12 + src/RENUMBER/RENUMBER_METISRenumbering.cxx | 15 + src/RENUMBER/RENUMBER_METISRenumbering.hxx | 12 + src/RENUMBER/RENUMBER_Renumbering.cxx | 1 + src/RENUMBER/RENUMBER_Renumbering.hxx | 11 + src/RENUMBER/RenumberingFactory.cxx | 64 ++++ src/RENUMBER/RenumberingFactory.hxx | 12 + src/RENUMBER/renumbering.cxx | 347 +++++++++++++++++++++ src/RENUMBER/testRenumbering.py | 287 +++++++++++++++++ 11 files changed, 888 insertions(+) create mode 100644 src/RENUMBER/Makefile.am create mode 100644 src/RENUMBER/RENUMBER_BOOSTRenumbering.cxx create mode 100644 src/RENUMBER/RENUMBER_BOOSTRenumbering.hxx create mode 100644 src/RENUMBER/RENUMBER_METISRenumbering.cxx create mode 100644 src/RENUMBER/RENUMBER_METISRenumbering.hxx create mode 100644 src/RENUMBER/RENUMBER_Renumbering.cxx create mode 100644 src/RENUMBER/RENUMBER_Renumbering.hxx create mode 100644 src/RENUMBER/RenumberingFactory.cxx create mode 100644 src/RENUMBER/RenumberingFactory.hxx create mode 100644 src/RENUMBER/renumbering.cxx create mode 100755 src/RENUMBER/testRenumbering.py diff --git a/src/RENUMBER/Makefile.am b/src/RENUMBER/Makefile.am new file mode 100644 index 000000000..80d1f7733 --- /dev/null +++ b/src/RENUMBER/Makefile.am @@ -0,0 +1,92 @@ +# Copyright (C) 2007-2008 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 MEDMEM : 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= librenumber.la + +salomeinclude_HEADERS= \ +RENUMBER_Renumbering.hxx \ +RenumberingFactory.hxx + +if MED_ENABLE_METIS + salomeinclude_HEADERS+= RENUMBER_METISRenumbering.hxx +endif +if BOOST_IS_OK + salomeinclude_HEADERS+= RENUMBER_BOOSTRenumbering.hxx +endif + +dist_librenumber_la_SOURCES= \ +RENUMBER_Renumbering.cxx \ +RenumberingFactory.cxx + +if MED_ENABLE_METIS + dist_librenumber_la_SOURCES+= RENUMBER_METISRenumbering.cxx +endif +if BOOST_IS_OK + dist_librenumber_la_SOURCES+= RENUMBER_BOOSTRenumbering.cxx +endif + +librenumber_la_CPPFLAGS= $(MED2_INCLUDES) $(HDF5_INCLUDES) @CXXTMPDPTHFLAGS@ \ + $(BOOST_CPPFLAGS) \ + -I$(srcdir)/../MEDMEM -I$(srcdir)/../MEDWrapper/V2_1/Core + +librenumber_la_LDFLAGS= +#libmedsplitter_la_LDFLAGS= $(MED2_LIBS) $(HDF5_LIBS) $(STDLIB) \ +# ../MEDMEM/libmedmem.la ../MEDWrapper/V2_1/Core/libmed_V2_1.la + +if MED_ENABLE_METIS + librenumber_la_CPPFLAGS+= $(METIS_CPPFLAGS) + librenumber_la_LDFLAGS+= $(METIS_LIBS) +endif +if BOOST_IS_OK + librenumber_la_CPPFLAGS+= $(BOOST_CPPFLAGS) -DENABLE_BOOST + librenumber_la_LDFLAGS+= $(BOOST_LIBS) +endif +if MED_ENABLE_KERNEL + librenumber_la_CPPFLAGS+= ${KERNEL_CXXFLAGS} + librenumber_la_LDFLAGS+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace +endif + +librenumber_la_LDFLAGS+= $(MED2_LIBS) $(HDF5_LIBS) $(STDLIB) \ + ../MEDMEM/libmedmem.la ../MEDWrapper/V2_1/Core/libmed_V2_1.la ../INTERP_KERNEL/libinterpkernel.la + +# Executables targets +bin_PROGRAMS= renumber + +dist_renumber_SOURCES= renumbering.cxx + +renumber_CPPFLAGS= $(librenumber_la_CPPFLAGS) +renumber_LDADD= $(librenumber_la_LDFLAGS) -lm $(BOOST_LIBS) librenumber.la +if MED_ENABLE_KERNEL + renumber_LDADD+= -lSALOMEBasics +endif + +OBSOLETE_FILES = +# MEDSPLITTER_SequentialTopology.cxx \ +# test_HighLevelAPI.cxx + +EXTRA_DIST += $(OBSOLETE_FILES) diff --git a/src/RENUMBER/RENUMBER_BOOSTRenumbering.cxx b/src/RENUMBER/RENUMBER_BOOSTRenumbering.cxx new file mode 100644 index 000000000..f54a7a1d0 --- /dev/null +++ b/src/RENUMBER/RENUMBER_BOOSTRenumbering.cxx @@ -0,0 +1,35 @@ + +#include +#include +#include +#include +#include + +#include "RENUMBER_BOOSTRenumbering.hxx" + +void BOOSTRenumbering::renumber(const int* graph,const int* index_graph,int nb_cell,std::vector& iperm,std::vector& perm) +{ + iperm.resize(nb_cell,0); + perm.resize(nb_cell,0); + + typedef boost::adjacency_list > > Graph; + typedef boost::graph_traits::vertex_descriptor Vertex; + typedef boost::graph_traits::vertices_size_type size_type; + Graph G(nb_cell); + for (int i=0;i::type + index_map = boost::get(boost::vertex_index, G); + boost::cuthill_mckee_ordering(G, iperm.rbegin(), boost::get(boost::vertex_color, G), + boost::make_degree_map(G)); + for (size_type c = 0; c != iperm.size(); ++c) + perm[index_map[iperm[c]]] = c; + for(int i=0;i& iperm,std::vector& perm); +}; + +#endif /*BOOSTRENUMBERING_HXX_*/ diff --git a/src/RENUMBER/RENUMBER_METISRenumbering.cxx b/src/RENUMBER/RENUMBER_METISRenumbering.cxx new file mode 100644 index 000000000..5c8daab9d --- /dev/null +++ b/src/RENUMBER/RENUMBER_METISRenumbering.cxx @@ -0,0 +1,15 @@ +extern "C" +{ +#include "metis.h" +} + +#include "RENUMBER_METISRenumbering.hxx" + +void METISRenumbering::renumber(const int* graph,const int* index_graph,int nb_cell,std::vector& iperm,std::vector& perm) +{ + iperm.resize(nb_cell,0); + perm.resize(nb_cell,0); + int num_flag=1; + int options=0; + METIS_NodeND(&nb_cell,(int*)index_graph,(int*)graph,&num_flag,&options,&iperm[0],&perm[0]); +} diff --git a/src/RENUMBER/RENUMBER_METISRenumbering.hxx b/src/RENUMBER/RENUMBER_METISRenumbering.hxx new file mode 100644 index 000000000..560a5c5a9 --- /dev/null +++ b/src/RENUMBER/RENUMBER_METISRenumbering.hxx @@ -0,0 +1,12 @@ +#ifndef METISRENUMBERING_HXX_ +#define METISRENUMBERING_HXX_ + +#include "RENUMBER_Renumbering.hxx" + +class METISRenumbering:public Renumbering +{ +public: + virtual void renumber(const int* graph,const int* index_graph,int nb_cell,std::vector& iperm,std::vector& perm); +}; + +#endif /*METISRENUMBERING_HXX_*/ diff --git a/src/RENUMBER/RENUMBER_Renumbering.cxx b/src/RENUMBER/RENUMBER_Renumbering.cxx new file mode 100644 index 000000000..8d1c8b69c --- /dev/null +++ b/src/RENUMBER/RENUMBER_Renumbering.cxx @@ -0,0 +1 @@ + diff --git a/src/RENUMBER/RENUMBER_Renumbering.hxx b/src/RENUMBER/RENUMBER_Renumbering.hxx new file mode 100644 index 000000000..c623d340f --- /dev/null +++ b/src/RENUMBER/RENUMBER_Renumbering.hxx @@ -0,0 +1,11 @@ +#ifndef RENUMBERING_HXX_ +#define RENUMBERING_HXX_ +#include + +class Renumbering +{ +public: + virtual void renumber(const int* graphe,const int* index_graphe,int nb_cell,std::vector& iperm,std::vector& perm)=0; +}; + +#endif /*RENUMBERING_HXX_*/ diff --git a/src/RENUMBER/RenumberingFactory.cxx b/src/RENUMBER/RenumberingFactory.cxx new file mode 100644 index 000000000..4dc9f9582 --- /dev/null +++ b/src/RENUMBER/RenumberingFactory.cxx @@ -0,0 +1,64 @@ +#include "RenumberingFactory.hxx" +#include "RENUMBER_Renumbering.hxx" +#ifdef ENABLE_METIS +#include "RENUMBER_METISRenumbering.hxx" +#endif +#ifdef ENABLE_BOOST +#include "RENUMBER_BOOSTRenumbering.hxx" +#endif + +#include + +using namespace std; + +namespace MED_RENUMBER +{ + Renumbering* RenumberingFactory(const string &s) + { +#ifdef ENABLE_METIS +#ifdef ENABLE_BOOST + if (s=="METIS") + { + return new METISRenumbering; + } + else if(s=="BOOST") + { + return new BOOSTRenumbering; + } + else + { + std::cerr << "The method has to be METIS or BOOST" << std::endl; + return 0; + } +#endif +#ifndef ENABLE_BOOST + if (s=="METIS") + { + return new METISRenumbering; + } + else + { + std::cerr << "The method has to be METIS!" << std::endl; + return 0; + } +#endif +#endif +#ifndef ENABLE_METIS +#ifdef ENABLE_BOOST + if (s=="BOOST") + { + return new BOOSTRenumbering; + } + else + { + std::cerr << "The method has to be BOOST!" << std::endl; + return 0; + } +#endif +#ifndef ENABLE_BOOST + std::cerr << "Error, no method compiled" << std::endl; + return 0; +#endif +#endif + } +} diff --git a/src/RENUMBER/RenumberingFactory.hxx b/src/RENUMBER/RenumberingFactory.hxx new file mode 100644 index 000000000..fe0f13f0f --- /dev/null +++ b/src/RENUMBER/RenumberingFactory.hxx @@ -0,0 +1,12 @@ +#ifndef RENUMBERINGFACTORY_HXX_ +#define RENUMBERINGFACTORY_HXX_ + +#include +#include "RENUMBER_Renumbering.hxx" + +namespace MED_RENUMBER +{ + Renumbering* RenumberingFactory(const std::string& s); +} + +#endif /*RENUMBERINGFACTORY_HXX_*/ diff --git a/src/RENUMBER/renumbering.cxx b/src/RENUMBER/renumbering.cxx new file mode 100644 index 000000000..44eac9a2e --- /dev/null +++ b/src/RENUMBER/renumbering.cxx @@ -0,0 +1,347 @@ +#include +#include +#include +#include + +#include "MEDMEM_Family.hxx" +#include "MEDMEM_Mesh.hxx" +#include "MEDMEM_Meshing.hxx" +#include "MEDMEM_MedMeshDriver.hxx" +#include "MEDMEM_Connectivity.hxx" +#include "MEDMEM_Field.hxx" +#include "MEDMEM_DriversDef.hxx" +#include "MEDMEM_Med.hxx" +#include "MEDMEM_MedMeshDriver22.hxx" + +#include "RenumberingFactory.hxx" + +#include +using namespace MEDMEM; +using namespace std; +using namespace MED_EN; +using namespace MED_RENUMBER; + +void computeNeighbour(const MESH* mesh,const medGeometryElement& Type, vector >& neighbour, int& ntot,int& nb_cell) +{ + CONNECTIVITY* conn = (CONNECTIVITY*)mesh->getConnectivityptr(); + conn->calculateFullDescendingConnectivity(MED_CELL); + const int* rev_conn=mesh->getReverseConnectivity(MED_EN::MED_DESCENDING, MED_EN::MED_CELL); + const int* rev_conn_index=mesh->getReverseConnectivityIndex(MED_EN::MED_DESCENDING, MED_EN::MED_CELL); + int nb_face= mesh->getNumberOfElementsWithPoly(MED_FACE,MED_ALL_ELEMENTS); + int nb_edge = mesh->getNumberOfElementsWithPoly(MED_EDGE,MED_ALL_ELEMENTS); + nb_cell= mesh->getNumberOfElementsWithPoly(MED_CELL,Type); + + int nb_constituent; + if(mesh->getMeshDimension()==2) + nb_constituent = nb_edge; + else if (mesh->getMeshDimension()==3) + nb_constituent = nb_face; + else + throw MEDEXCEPTION("Wrong dimension"); + + neighbour.resize(nb_cell,(list)0); + ntot=0; + for(int i=0;i& iperm) +{ + if(Type==MED_POLYHEDRA) + { + int *conn_face_index_init=(int*)mesh.getPolyhedronFacesIndex(); + int *conn_index_init=(int*)mesh.getPolyhedronIndex(MED_FULL_INTERLACE); + int *conn_init=(int*)mesh.getPolyhedronConnectivity(MED_FULL_INTERLACE); + + int *conn_index_renum=new int[nb_cell+1]; + int *conn_face_index_renum=new int[conn_index_init[nb_cell]]; + int *conn_renum=new int[conn_face_index_init[conn_index_init[nb_cell]-1]-1]; + + int i_cell,i_face,i_conn; + int iter_face=0; + int iter_conn=0; + int i2; + conn_index_renum[0]=1; + conn_face_index_renum[0]=1; + for(i_cell=0;i_cellsetNodal(conn_renum,MED_CELL,Type); + delete[] conn_renum; + delete[] conn_index_renum; + } +} + +void changeFamily(MESH* mesh, const medGeometryElement& Type, const vector& perm) +{ + int nb_families=mesh->getNumberOfFamilies(MED_CELL); + for (int i=0;igetFamily(MED_CELL,i+1); + if (!family->isOnAllElements()) + { + int nb_elem=family->getNumberOfElements(Type); + int *number=(int *)family->getNumber(Type); + for(int j=0;j mesh_names,f_names; + nb_mesh=med_struct.getNumberOfMeshes(); + nb_fields=med_struct.getNumberOfFields(); + mesh_names=med_struct.getMeshNames(); + f_names=med_struct.getFieldNames(); + if(nb_mesh!=1) + { + cout << "There are many meshes in the file" << endl; + return -1; + } + if(mesh_names[0].c_str()!=meshname) + { + cout << "Mesh name does not match" << endl; + return -1; + } + vector field_names; + vector iternumber; + vector ordernumber; + vector types; + for (int ifield = 0; ifield < nb_fields; ifield++) + { + deque dtit=med_struct.getFieldIteration(f_names[ifield]); + for (deque::const_iterator iter =dtit.begin(); iter!=dtit.end(); iter++) + { + field_names.push_back(f_names[ifield]); + iternumber.push_back(iter->dt); + ordernumber.push_back(iter->it); + + FIELD_* field = med_struct.getField(f_names[ifield],iter->dt,iter->it); + if (dynamic_cast*>(field)) + types.push_back(1); + else + types.push_back(0); + + } + } + t_read_st=clock(); + + // Reading mesh + MESH myMesh; + myMesh.setName(meshname); + MED_MESH_RDONLY_DRIVER22 *drv22=new MED_MESH_RDONLY_DRIVER22(filename_in,&myMesh); + drv22->desactivateFacesComputation(); + int newDrv=myMesh.addDriver(*drv22); + delete drv22; + myMesh.read(newDrv); + int nb_type=myMesh.getNumberOfTypesWithPoly(MED_CELL); + if (nb_type!=1) + { + cout << "Mesh must have only one type of cell" << endl; + return -1; + } + medGeometryElement *Types = myMesh.getTypesWithPoly(MED_CELL); + medGeometryElement Type=Types[0]; + delete[] Types; + + t_read_mesh=clock(); + MESH* workMesh=new MESH(myMesh); + cout << "Building the graph "; + cout.flush(); + int ntot,nb_cell; + vector > neighbour; + computeNeighbour(workMesh,Type,neighbour,ntot,nb_cell); + int* graph=new int[ntot]; + int* graph_index=new int[nb_cell+1]; + graph_index[0]=1; + int count=0; + for(int i=0;i::const_iterator it=neighbour[i].begin();it!=neighbour[i].end();++it) + { + graph[count]=*it; + ++count; + } + graph_index[i+1]=count+1; + } + + + // Compute permutation + vector iperm,perm; + Renumbering* renumb= RenumberingFactory(type_renum); + renumb->renumber(graph,graph_index,nb_cell,iperm,perm); + delete renumb; + delete workMesh; + t_compute_graph=clock(); + cout << " : " << (t_compute_graph-t_read_mesh)/(double) CLOCKS_PER_SEC << "s" << endl; + cout.flush(); + + // Connectivity + cout << "Computing connectivity"; + cout.flush(); + MESH meshRenum(myMesh); + changeConnectivity(meshRenum,Type,nb_cell,iperm); + t_connectiv=clock(); + cout << " : " << (t_connectiv-t_compute_graph)/(double) CLOCKS_PER_SEC << "s" << endl; + cout.flush(); + + // Familles + cout << "Computing families "; + cout.flush(); + changeFamily(&meshRenum,Type,perm); + int drv3=meshRenum.addDriver(MED_DRIVER,filename_out,meshRenum.getName()); + meshRenum.write(drv3); + t_family=clock(); + cout << " : " << (t_family-t_connectiv)/(double) CLOCKS_PER_SEC << "s" << endl; + cout.flush(); + + // Fields + cout << "Computing fields "; + cout.flush(); + bool exist_type; + for(int ifield=0;ifield myField(MED_DRIVER,filename_in,field_names[ifield],iternumber[ifield],ordernumber[ifield]); + FIELD newField(myField); + const SUPPORT* mySupport=newField.getSupport(); + const medGeometryElement *typesOfSupport = mySupport->getTypes(); + for(int t=0;tgetNumberOfTypes();++t) + { + if(typesOfSupport[t]==Type) + { + exist_type=true; + break; + } + } + if(exist_type) + { + for(int i=0;igetNumberOfElements(Type);++i) + { + for(int j=0;j