Back to index

salome-med  6.5.0
renumbering.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include <string>
00021 #include <stdlib.h>
00022 #include <iostream>
00023 #include <fstream>
00024 
00025 #include "MEDMEM_Family.hxx"
00026 #include "MEDMEM_Mesh.hxx"
00027 #include "MEDMEM_Meshing.hxx"
00028 #include "MEDMEM_MedMeshDriver.hxx"
00029 #include "MEDMEM_Connectivity.hxx"
00030 #include "MEDMEM_Field.hxx"
00031 #include "MEDMEM_DriversDef.hxx"
00032 #include "MEDMEM_MedFileBrowser.hxx"
00033 #include "MEDMEM_MedMeshDriver.hxx"
00034 
00035 #include "RenumberingFactory.hxx"
00036 
00037 #include <time.h>
00038 using namespace MEDMEM;
00039 using namespace std;
00040 using namespace MED_EN;
00041 using namespace MED_RENUMBER;
00042 
00043 void computeNeighbour(const MESH* mesh,const medGeometryElement& Type, vector<list<int> >& neighbour, int& ntot,int& nb_cell)
00044 {
00045   CONNECTIVITY* conn = (CONNECTIVITY*)mesh->getConnectivityptr();
00046   conn->calculateFullDescendingConnectivity(MED_CELL);
00047   const int* rev_conn=mesh->getReverseConnectivity(MED_EN::MED_DESCENDING, MED_EN::MED_CELL);
00048   const int* rev_conn_index=mesh->getReverseConnectivityIndex(MED_EN::MED_DESCENDING, MED_EN::MED_CELL);
00049   int nb_face= mesh->getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
00050   int nb_edge = mesh->getNumberOfElements(MED_EDGE,MED_ALL_ELEMENTS);
00051   nb_cell= mesh->getNumberOfElements(MED_CELL,Type);
00052 
00053   int nb_constituent;
00054   if(mesh->getMeshDimension()==2)
00055     nb_constituent = nb_edge;
00056   else if (mesh->getMeshDimension()==3)
00057     nb_constituent = nb_face;
00058   else
00059     throw MEDEXCEPTION("Wrong dimension");
00060 
00061   neighbour.resize(nb_cell,(list<int>)0);
00062   ntot=0;
00063   for(int i=0;i<nb_constituent;++i)
00064     {
00065       for(int j=rev_conn_index[i]-1;j<rev_conn_index[i+1]-1;++j)
00066         {
00067           for(int k=j+1;k<rev_conn_index[i+1]-1;++k)
00068             {
00069               if(rev_conn[j]!=0 && rev_conn[k]!=0)
00070                 {
00071                   ntot+=2;
00072                   neighbour[rev_conn[j]-1].push_back(rev_conn[k]);
00073                   neighbour[rev_conn[k]-1].push_back(rev_conn[j]);
00074                 }
00075             }
00076         }
00077     }
00078 }
00079 
00080 void changeConnectivity(MESH& mesh, const medGeometryElement& Type, const int& nb_cell, const vector<int>& iperm)
00081 {
00082   /*if(Type==MED_POLYHEDRA)
00083     {
00084       int *conn_face_index_init=(int*)mesh.getPolyhedronFacesIndex();
00085       int *conn_index_init=(int*)mesh.getPolyhedronIndex(MED_FULL_INTERLACE);
00086       int *conn_init=(int*)mesh.getPolyhedronConnectivity(MED_FULL_INTERLACE);
00087 
00088       int *conn_index_renum=new int[nb_cell+1];
00089       int *conn_face_index_renum=new int[conn_index_init[nb_cell]];
00090       int *conn_renum=new int[conn_face_index_init[conn_index_init[nb_cell]-1]-1];
00091 
00092       int i_cell,i_face,i_conn;
00093       int iter_face=0;
00094       int iter_conn=0;
00095       int i2;
00096       conn_index_renum[0]=1;
00097       conn_face_index_renum[0]=1;
00098       for(i_cell=0;i_cell<nb_cell;++i_cell)
00099         {
00100           i2=iperm[i_cell]-1;
00101           for(i_face=conn_index_init[i2]-1;i_face<conn_index_init[i2+1]-1;++i_face)
00102             {
00103               for(i_conn=conn_face_index_init[i_face]-1;i_conn<conn_face_index_init[i_face+1]-1;++i_conn)
00104                 {
00105                   conn_renum[iter_conn]=conn_init[i_conn];
00106                   ++iter_conn;
00107                 }
00108               conn_face_index_renum[iter_face+1]=iter_conn+1;
00109               ++iter_face;
00110             }
00111           conn_index_renum[i_cell+1]=iter_face+1;
00112         }
00113       memcpy(conn_face_index_init,conn_face_index_renum,sizeof(int)*conn_index_init[nb_cell]);
00114       memcpy(conn_index_init,conn_index_renum,sizeof(int)*(nb_cell+1));
00115       memcpy(conn_init,conn_renum, sizeof(int)*(conn_face_index_init[conn_index_init[nb_cell]-1]-1));
00116 
00117       delete[] conn_index_renum;
00118       delete[] conn_face_index_renum;
00119       delete[] conn_renum;
00120     }
00121   else if (Type==MED_POLYGON)
00122     {
00123       int *conn_init=(int*)mesh.getPolygonsConnectivity(MED_FULL_INTERLACE,MED_CELL);
00124       int *conn_index_init=(int*)mesh.getPolygonsConnectivityIndex(MED_FULL_INTERLACE,MED_CELL);
00125       int *conn_index_renum=new int[nb_cell+1];
00126       int *conn_renum=new int[conn_index_init[nb_cell]-1];
00127 
00128       int iter=0;
00129       int i2;
00130       conn_index_renum[0]=1;
00131       for(int i=0;i<nb_cell;++i)
00132         {
00133           i2=iperm[i]-1;
00134           for(int k=conn_index_init[i2];k<conn_index_init[i2+1];++k)
00135             {
00136               conn_renum[iter]=conn_init[k-1];
00137               ++iter;
00138             }
00139           conn_index_renum[i+1]=iter+1;
00140         }
00141       memcpy(conn_index_init,conn_index_renum,sizeof(int)*(nb_cell+1));
00142       memcpy(conn_init,conn_renum, sizeof(int)*(conn_index_init[nb_cell]-1));
00143 
00144       delete[] conn_renum;
00145       delete[] conn_index_renum;
00146     }
00147     else*/
00148     {
00149       const int *conn_init=mesh.getConnectivity(MED_NODAL,MED_CELL,Type);
00150       const int *conn_index_init=mesh.getConnectivityIndex(MED_NODAL,MED_CELL);
00151       int *conn_renum=new int[conn_index_init[nb_cell]-1];
00152       int *conn_index_renum=new int[nb_cell+1];
00153 
00154       int iter=0;
00155       int i2;
00156       conn_index_renum[0]=1;
00157       for(int i=0;i<nb_cell;++i)
00158         {
00159           i2=iperm[i]-1;
00160           for(int k=conn_index_init[i2];k<conn_index_init[i2+1];++k)
00161             {
00162               conn_renum[iter]=conn_init[k-1];
00163               ++iter;
00164             }
00165           conn_index_renum[i+1]=iter+1;
00166         }
00167 
00168       CONNECTIVITY* myConnectivity=(CONNECTIVITY*)mesh.getConnectivityptr();
00169       myConnectivity->setNodal(conn_renum,MED_CELL,Type,conn_index_renum);
00170       delete[] conn_renum;
00171       delete[] conn_index_renum;
00172     }
00173 }
00174 
00175 void changeFamily(MESH* mesh, const medGeometryElement& Type, const vector<int>& perm)
00176 {
00177   int nb_families=mesh->getNumberOfFamilies(MED_CELL);
00178   for (int i=0;i<nb_families;++i)
00179     {
00180       const FAMILY* family=mesh->getFamily(MED_CELL,i+1);
00181       if (!family->isOnAllElements())
00182         {
00183           int nb_elem=family->getNumberOfElements(Type);
00184           int *number=(int *)family->getNumber(Type);
00185           for(int j=0;j<nb_elem;++j)
00186             number[j]=perm[number[j]-1];
00187         }
00188     }
00189 }
00190 
00191 int main (int argc, char** argv)
00192 {
00193   double t_begin,t_read_st,t_read_mesh,t_compute_graph,t_connectiv,t_family,t_field;
00194   t_begin=clock();
00195   if (argc <5)
00196     {
00197       cerr << "Usage : " << argv[0] 
00198            << " filename_in meshname method[BOOST/METIS] filename_out" << endl << endl;
00199       return -1;
00200     }
00201   string filename_in = argv[1];
00202   string meshname = argv[2];
00203   string type_renum = argv[3];
00204   string filename_out = argv[4];
00205 
00206   if(type_renum!="METIS" && type_renum!="BOOST")
00207     {
00208       cout << "The method has to be METIS or BOOST!" << endl;
00209       exit(-1);
00210     }
00211 
00212   string s="rm "+filename_out;
00213   system(s.c_str());
00214 
00215   // Reading file structure
00216   const MEDFILEBROWSER med_struct(filename_in);
00217   int nb_mesh, nb_fields;
00218   vector<string> mesh_names,f_names;
00219   nb_mesh=med_struct.getNumberOfMeshes();
00220   nb_fields=med_struct.getNumberOfFields();
00221   mesh_names=med_struct.getMeshNames();
00222   f_names=med_struct.getFieldNames();
00223   if(nb_mesh!=1)
00224     {
00225       cout << "There are many meshes in the file" << endl;
00226       return -1;
00227     }
00228   if(mesh_names[0].c_str()!=meshname)
00229     {
00230       cout << "Mesh name does not match" << endl;
00231       return -1;
00232     }
00233   vector<string> field_names;
00234   vector<int> iternumber;
00235   vector<int> ordernumber;
00236   vector<int> types;
00237   int nb_fields_tot=0;
00238   for (int ifield = 0; ifield < nb_fields; ifield++)
00239     {
00240       vector<DT_IT_> dtit=med_struct.getFieldIteration(f_names[ifield]);
00241       for (vector<DT_IT_>::const_iterator iter =dtit.begin(); iter!=dtit.end(); iter++)
00242         {
00243           field_names.push_back(f_names[ifield]);
00244           iternumber.push_back(iter->dt);
00245           ordernumber.push_back(iter->it);
00246           ++nb_fields_tot;
00247           if(med_struct.getFieldType(f_names[ifield])==MED_EN::MED_REEL64)
00248             types.push_back(1);
00249           else
00250             types.push_back(0);
00251 
00252         }
00253     }
00254   t_read_st=clock();
00255 
00256   // Reading mesh
00257   MESH myMesh;
00258   myMesh.setName(meshname);
00259   MED_MESH_RDONLY_DRIVER *drv22=new MED_MESH_RDONLY_DRIVER(filename_in,&myMesh);
00260   drv22->desactivateFacesComputation();
00261   int newDrv=myMesh.addDriver(*drv22);
00262   delete drv22;
00263   myMesh.read(newDrv);
00264   int nb_type=myMesh.getNumberOfTypes(MED_CELL);
00265   if (nb_type!=1)
00266     {
00267       cout << "Mesh must have only one type of cell" << endl;
00268       return -1;
00269     }
00270   const medGeometryElement *Types = myMesh.getTypes(MED_CELL);
00271   medGeometryElement Type=Types[0];
00272 
00273   t_read_mesh=clock();
00274   MESH* workMesh=new MESH(myMesh);
00275   cout << "Building the graph    ";
00276   cout.flush();
00277   int ntot,nb_cell;
00278   vector<list<int> > neighbour;
00279   computeNeighbour(workMesh,Type,neighbour,ntot,nb_cell);
00280   int* graph=new int[ntot];
00281   int* graph_index=new int[nb_cell+1];
00282   graph_index[0]=1;
00283   int count=0;
00284   for(int i=0;i<nb_cell;++i)
00285     {
00286       for (list<int>::const_iterator it=neighbour[i].begin();it!=neighbour[i].end();++it)
00287         {
00288           graph[count]=*it;
00289           ++count;
00290         }
00291       graph_index[i+1]=count+1;
00292     }
00293 
00294 
00295   // Compute permutation
00296   vector<int> iperm,perm;
00297   Renumbering* renumb= RenumberingFactory(type_renum);
00298   renumb->renumber(graph,graph_index,nb_cell,iperm,perm);
00299   delete renumb;
00300   delete workMesh;
00301   t_compute_graph=clock();
00302   cout << " : " << (t_compute_graph-t_read_mesh)/(double) CLOCKS_PER_SEC << "s" << endl;
00303   cout.flush();
00304 
00305   // Connectivity
00306   cout << "Computing connectivity";
00307   cout.flush();
00308   MESH meshRenum(myMesh);
00309   changeConnectivity(meshRenum,Type,nb_cell,iperm);
00310   t_connectiv=clock();
00311   cout << " : " << (t_connectiv-t_compute_graph)/(double) CLOCKS_PER_SEC << "s" << endl;
00312   cout.flush();
00313 
00314   // Familles
00315   cout << "Computing families    ";
00316   cout.flush();
00317   changeFamily(&meshRenum,Type,perm);
00318   int drv3=meshRenum.addDriver(MED_DRIVER,filename_out,meshRenum.getName());
00319   meshRenum.write(drv3);
00320   t_family=clock();
00321   cout << " : " << (t_family-t_connectiv)/(double) CLOCKS_PER_SEC << "s" << endl;
00322   cout.flush();
00323 
00324   // Fields
00325   cout << "Computing fields      ";
00326   cout.flush();
00327   bool exist_type;
00328   for(int ifield=0;ifield<nb_fields_tot;++ifield)
00329     {
00330       exist_type=false;
00331       FIELD<double> myField(MED_DRIVER,filename_in,field_names[ifield],iternumber[ifield],ordernumber[ifield]);
00332       FIELD<double> newField(myField);
00333       const SUPPORT* mySupport=newField.getSupport();
00334       const medGeometryElement *typesOfSupport = mySupport->getTypes();
00335       for(int t=0;t<mySupport->getNumberOfTypes();++t)
00336         {
00337           if(typesOfSupport[t]==Type)
00338             {
00339               exist_type=true;
00340               break;
00341             }
00342         }
00343       if(exist_type)
00344         {
00345           for(int i=0;i<mySupport->getNumberOfElements(Type);++i)
00346             {
00347               for(int j=0;j<newField.getNumberOfComponents();++j)
00348                 {
00349                   newField.setValueIJ(i+1,j+1,myField.getValueIJ(iperm[i],j+1));
00350                 }
00351             }
00352         }
00353       int drv=newField.addDriver(MED_DRIVER,filename_out,field_names[ifield]);
00354       newField.write(drv);
00355     }
00356   t_field=clock();
00357   cout << " : " << (t_field-t_family)/(double) CLOCKS_PER_SEC << "s" << endl;
00358   cout.flush();
00359 
00360   delete[] graph_index;
00361   delete[] graph;
00362 
00363   return 0;
00364 }