src/RENUMBER/renumbering.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include <string>
  21 #include <stdlib.h>
  22 #include <iostream>
  23 #include <fstream>
  24
  25 #include "MEDMEM_Family.hxx"
  26 #include "MEDMEM_Mesh.hxx"
  27 #include "MEDMEM_Meshing.hxx"
  28 #include "MEDMEM_MedMeshDriver.hxx"
  29 #include "MEDMEM_Connectivity.hxx"
  30 #include "MEDMEM_Field.hxx"
  31 #include "MEDMEM_DriversDef.hxx"
  32 #include "MEDMEM_MedFileBrowser.hxx"
  33 #include "MEDMEM_MedMeshDriver.hxx"
  34
  35 #include "RenumberingFactory.hxx"
  36
  37 #include <time.h>
  38 using namespace MEDMEM;
  39 using namespace std;
  40 using namespace MED_EN;
  41 using namespace MED_RENUMBER;
  42
  43 void computeNeighbour(const MESH* mesh,const medGeometryElement& Type, vector<list<int> >& neighbour, int& ntot,int& nb_cell)
  44 {
  45   CONNECTIVITY* conn = (CONNECTIVITY*)mesh->getConnectivityptr();
  46   conn->calculateFullDescendingConnectivity(MED_CELL);
  47   const int* rev_conn=mesh->getReverseConnectivity(MED_EN::MED_DESCENDING, MED_EN::MED_CELL);
  48   const int* rev_conn_index=mesh->getReverseConnectivityIndex(MED_EN::MED_DESCENDING, MED_EN::MED_CELL);
  49   int nb_face= mesh->getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
  50   int nb_edge = mesh->getNumberOfElements(MED_EDGE,MED_ALL_ELEMENTS);
  51   nb_cell= mesh->getNumberOfElements(MED_CELL,Type);
  52
  53   int nb_constituent;
  54   if(mesh->getMeshDimension()==2)
  55     nb_constituent = nb_edge;
  56   else if (mesh->getMeshDimension()==3)
  57     nb_constituent = nb_face;
  58   else
  59     throw MEDEXCEPTION("Wrong dimension");
  60
  61   neighbour.resize(nb_cell,(list<int>)0);
  62   ntot=0;
  63   for(int i=0;i<nb_constituent;++i)
  64     {
  65       for(int j=rev_conn_index[i]-1;j<rev_conn_index[i+1]-1;++j)
  66         {
  67           for(int k=j+1;k<rev_conn_index[i+1]-1;++k)
  68             {
  69               if(rev_conn[j]!=0 && rev_conn[k]!=0)
  70                 {
  71                   ntot+=2;
  72                   neighbour[rev_conn[j]-1].push_back(rev_conn[k]);
  73                   neighbour[rev_conn[k]-1].push_back(rev_conn[j]);
  74                 }
  75             }
  76         }
  77     }
  78 }
  79
  80 void changeConnectivity(MESH& mesh, const medGeometryElement& Type, const int& nb_cell, const vector<int>& iperm)
  81 {
  82   /*if(Type==MED_POLYHEDRA)
  83     {
  84       int *conn_face_index_init=(int*)mesh.getPolyhedronFacesIndex();
  85       int *conn_index_init=(int*)mesh.getPolyhedronIndex(MED_FULL_INTERLACE);
  86       int *conn_init=(int*)mesh.getPolyhedronConnectivity(MED_FULL_INTERLACE);
  87
  88       int *conn_index_renum=new int[nb_cell+1];
  89       int *conn_face_index_renum=new int[conn_index_init[nb_cell]];
  90       int *conn_renum=new int[conn_face_index_init[conn_index_init[nb_cell]-1]-1];
  91
  92       int i_cell,i_face,i_conn;
  93       int iter_face=0;
  94       int iter_conn=0;
  95       int i2;
  96       conn_index_renum[0]=1;
  97       conn_face_index_renum[0]=1;
  98       for(i_cell=0;i_cell<nb_cell;++i_cell)
  99         {
 100           i2=iperm[i_cell]-1;
 101           for(i_face=conn_index_init[i2]-1;i_face<conn_index_init[i2+1]-1;++i_face)
 102             {
 103               for(i_conn=conn_face_index_init[i_face]-1;i_conn<conn_face_index_init[i_face+1]-1;++i_conn)
 104                 {
 105                   conn_renum[iter_conn]=conn_init[i_conn];
 106                   ++iter_conn;
 107                 }
 108               conn_face_index_renum[iter_face+1]=iter_conn+1;
 109               ++iter_face;
 110             }
 111           conn_index_renum[i_cell+1]=iter_face+1;
 112         }
 113       memcpy(conn_face_index_init,conn_face_index_renum,sizeof(int)*conn_index_init[nb_cell]);
 114       memcpy(conn_index_init,conn_index_renum,sizeof(int)*(nb_cell+1));
 115       memcpy(conn_init,conn_renum, sizeof(int)*(conn_face_index_init[conn_index_init[nb_cell]-1]-1));
 116
 117       delete[] conn_index_renum;
 118       delete[] conn_face_index_renum;
 119       delete[] conn_renum;
 120     }
 121   else if (Type==MED_POLYGON)
 122     {
 123       int *conn_init=(int*)mesh.getPolygonsConnectivity(MED_FULL_INTERLACE,MED_CELL);
 124       int *conn_index_init=(int*)mesh.getPolygonsConnectivityIndex(MED_FULL_INTERLACE,MED_CELL);
 125       int *conn_index_renum=new int[nb_cell+1];
 126       int *conn_renum=new int[conn_index_init[nb_cell]-1];
 127
 128       int iter=0;
 129       int i2;
 130       conn_index_renum[0]=1;
 131       for(int i=0;i<nb_cell;++i)
 132         {
 133           i2=iperm[i]-1;
 134           for(int k=conn_index_init[i2];k<conn_index_init[i2+1];++k)
 135             {
 136               conn_renum[iter]=conn_init[k-1];
 137               ++iter;
 138             }
 139           conn_index_renum[i+1]=iter+1;
 140         }
 141       memcpy(conn_index_init,conn_index_renum,sizeof(int)*(nb_cell+1));
 142       memcpy(conn_init,conn_renum, sizeof(int)*(conn_index_init[nb_cell]-1));
 143
 144       delete[] conn_renum;
 145       delete[] conn_index_renum;
 146     }
 147     else*/
 148     {
 149       const int *conn_init=mesh.getConnectivity(MED_NODAL,MED_CELL,Type);
 150       const int *conn_index_init=mesh.getConnectivityIndex(MED_NODAL,MED_CELL);
 151       int *conn_renum=new int[conn_index_init[nb_cell]-1];
 152       int *conn_index_renum=new int[nb_cell+1];
 153
 154       int iter=0;
 155       int i2;
 156       conn_index_renum[0]=1;
 157       for(int i=0;i<nb_cell;++i)
 158         {
 159           i2=iperm[i]-1;
 160           for(int k=conn_index_init[i2];k<conn_index_init[i2+1];++k)
 161             {
 162               conn_renum[iter]=conn_init[k-1];
 163               ++iter;
 164             }
 165           conn_index_renum[i+1]=iter+1;
 166         }
 167
 168       CONNECTIVITY* myConnectivity=(CONNECTIVITY*)mesh.getConnectivityptr();
 169       myConnectivity->setNodal(conn_renum,MED_CELL,Type,conn_index_renum);
 170       delete[] conn_renum;
 171       delete[] conn_index_renum;
 172     }
 173 }
 174
 175 void changeFamily(MESH* mesh, const medGeometryElement& Type, const vector<int>& perm)
 176 {
 177   int nb_families=mesh->getNumberOfFamilies(MED_CELL);
 178   for (int i=0;i<nb_families;++i)
 179     {
 180       const FAMILY* family=mesh->getFamily(MED_CELL,i+1);
 181       if (!family->isOnAllElements())
 182         {
 183           int nb_elem=family->getNumberOfElements(Type);
 184           int *number=(int *)family->getNumber(Type);
 185           for(int j=0;j<nb_elem;++j)
 186             number[j]=perm[number[j]-1];
 187         }
 188     }
 189 }
 190
 191 int main (int argc, char** argv)
 192 {
 193   double t_begin,t_read_st,t_read_mesh,t_compute_graph,t_connectiv,t_family,t_field;
 194   t_begin=clock();
 195   if (argc <5)
 196     {
 197       cerr << "Usage : " << argv[0]
 198            << " filename_in meshname method[BOOST/METIS] filename_out" << endl << endl;
 199       return -1;
 200     }
 201   string filename_in = argv[1];
 202   string meshname = argv[2];
 203   string type_renum = argv[3];
 204   string filename_out = argv[4];
 205
 206   if(type_renum!="METIS" && type_renum!="BOOST")
 207     {
 208       cout << "The method has to be METIS or BOOST!" << endl;
 209       exit(-1);
 210     }
 211
 212   string s="rm "+filename_out;
 213   system(s.c_str());
 214
 215   // Reading file structure
 216   const MEDFILEBROWSER med_struct(filename_in);
 217   int nb_mesh, nb_fields;
 218   vector<string> mesh_names,f_names;
 219   nb_mesh=med_struct.getNumberOfMeshes();
 220   nb_fields=med_struct.getNumberOfFields();
 221   mesh_names=med_struct.getMeshNames();
 222   f_names=med_struct.getFieldNames();
 223   if(nb_mesh!=1)
 224     {
 225       cout << "There are many meshes in the file" << endl;
 226       return -1;
 227     }
 228   if(mesh_names[0].c_str()!=meshname)
 229     {
 230       cout << "Mesh name does not match" << endl;
 231       return -1;
 232     }
 233   vector<string> field_names;
 234   vector<int> iternumber;
 235   vector<int> ordernumber;
 236   vector<int> types;
 237   int nb_fields_tot=0;
 238   for (int ifield = 0; ifield < nb_fields; ifield++)
 239     {
 240       vector<DT_IT_> dtit=med_struct.getFieldIteration(f_names[ifield]);
 241       for (vector<DT_IT_>::const_iterator iter =dtit.begin(); iter!=dtit.end(); iter++)
 242         {
 243           field_names.push_back(f_names[ifield]);
 244           iternumber.push_back(iter->dt);
 245           ordernumber.push_back(iter->it);
 246           ++nb_fields_tot;
 247           if(med_struct.getFieldType(f_names[ifield])==MED_EN::MED_REEL64)
 248             types.push_back(1);
 249           else
 250             types.push_back(0);
 251
 252         }
 253     }
 254   t_read_st=clock();
 255
 256   // Reading mesh
 257   MESH myMesh;
 258   myMesh.setName(meshname);
 259   MED_MESH_RDONLY_DRIVER *drv22=new MED_MESH_RDONLY_DRIVER(filename_in,&myMesh);
 260   drv22->desactivateFacesComputation();
 261   int newDrv=myMesh.addDriver(*drv22);
 262   delete drv22;
 263   myMesh.read(newDrv);
 264   int nb_type=myMesh.getNumberOfTypes(MED_CELL);
 265   if (nb_type!=1)
 266     {
 267       cout << "Mesh must have only one type of cell" << endl;
 268       return -1;
 269     }
 270   const medGeometryElement *Types = myMesh.getTypes(MED_CELL);
 271   medGeometryElement Type=Types[0];
 272
 273   t_read_mesh=clock();
 274   MESH* workMesh=new MESH(myMesh);
 275   cout << "Building the graph    ";
 276   cout.flush();
 277   int ntot,nb_cell;
 278   vector<list<int> > neighbour;
 279   computeNeighbour(workMesh,Type,neighbour,ntot,nb_cell);
 280   int* graph=new int[ntot];
 281   int* graph_index=new int[nb_cell+1];
 282   graph_index[0]=1;
 283   int count=0;
 284   for(int i=0;i<nb_cell;++i)
 285     {
 286       for (list<int>::const_iterator it=neighbour[i].begin();it!=neighbour[i].end();++it)
 287         {
 288           graph[count]=*it;
 289           ++count;
 290         }
 291       graph_index[i+1]=count+1;
 292     }
 293
 294
 295   // Compute permutation
 296   vector<int> iperm,perm;
 297   Renumbering* renumb= RenumberingFactory(type_renum);
 298   renumb->renumber(graph,graph_index,nb_cell,iperm,perm);
 299   delete renumb;
 300   delete workMesh;
 301   t_compute_graph=clock();
 302   cout << " : " << (t_compute_graph-t_read_mesh)/(double) CLOCKS_PER_SEC << "s" << endl;
 303   cout.flush();
 304
 305   // Connectivity
 306   cout << "Computing connectivity";
 307   cout.flush();
 308   MESH meshRenum(myMesh);
 309   changeConnectivity(meshRenum,Type,nb_cell,iperm);
 310   t_connectiv=clock();
 311   cout << " : " << (t_connectiv-t_compute_graph)/(double) CLOCKS_PER_SEC << "s" << endl;
 312   cout.flush();
 313
 314   // Familles
 315   cout << "Computing families    ";
 316   cout.flush();
 317   changeFamily(&meshRenum,Type,perm);
 318   int drv3=meshRenum.addDriver(MED_DRIVER,filename_out,meshRenum.getName());
 319   meshRenum.write(drv3);
 320   t_family=clock();
 321   cout << " : " << (t_family-t_connectiv)/(double) CLOCKS_PER_SEC << "s" << endl;
 322   cout.flush();
 323
 324   // Fields
 325   cout << "Computing fields      ";
 326   cout.flush();
 327   bool exist_type;
 328   for(int ifield=0;ifield<nb_fields_tot;++ifield)
 329     {
 330       exist_type=false;
 331       FIELD<double> myField(MED_DRIVER,filename_in,field_names[ifield],iternumber[ifield],ordernumber[ifield]);
 332       FIELD<double> newField(myField);
 333       const SUPPORT* mySupport=newField.getSupport();
 334       const medGeometryElement *typesOfSupport = mySupport->getTypes();
 335       for(int t=0;t<mySupport->getNumberOfTypes();++t)
 336         {
 337           if(typesOfSupport[t]==Type)
 338             {
 339               exist_type=true;
 340               break;
 341             }
 342         }
 343       if(exist_type)
 344         {
 345           for(int i=0;i<mySupport->getNumberOfElements(Type);++i)
 346             {
 347               for(int j=0;j<newField.getNumberOfComponents();++j)
 348                 {
 349                   newField.setValueIJ(i+1,j+1,myField.getValueIJ(iperm[i],j+1));
 350                 }
 351             }
 352         }
 353       int drv=newField.addDriver(MED_DRIVER,filename_out,field_names[ifield]);
 354       newField.write(drv);
 355     }
 356   t_field=clock();
 357   cout << " : " << (t_field-t_family)/(double) CLOCKS_PER_SEC << "s" << endl;
 358   cout.flush();
 359
 360   delete[] graph_index;
 361   delete[] graph;
 362
 363   return 0;
 364 }