1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
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"
35 #include "RenumberingFactory.hxx"
38 using namespace MEDMEM;
40 using namespace MED_EN;
41 using namespace MED_RENUMBER;
43 void computeNeighbour(const MESH* mesh,const medGeometryElement& Type, vector<list<int> >& neighbour, int& ntot,int& nb_cell)
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);
54 if(mesh->getMeshDimension()==2)
55 nb_constituent = nb_edge;
56 else if (mesh->getMeshDimension()==3)
57 nb_constituent = nb_face;
59 throw MEDEXCEPTION("Wrong dimension");
61 neighbour.resize(nb_cell,(list<int>)0);
63 for(int i=0;i<nb_constituent;++i)
65 for(int j=rev_conn_index[i]-1;j<rev_conn_index[i+1]-1;++j)
67 for(int k=j+1;k<rev_conn_index[i+1]-1;++k)
69 if(rev_conn[j]!=0 && rev_conn[k]!=0)
72 neighbour[rev_conn[j]-1].push_back(rev_conn[k]);
73 neighbour[rev_conn[k]-1].push_back(rev_conn[j]);
80 void changeConnectivity(MESH& mesh, const medGeometryElement& Type, const int& nb_cell, const vector<int>& iperm)
82 /*if(Type==MED_POLYHEDRA)
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);
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];
92 int i_cell,i_face,i_conn;
96 conn_index_renum[0]=1;
97 conn_face_index_renum[0]=1;
98 for(i_cell=0;i_cell<nb_cell;++i_cell)
101 for(i_face=conn_index_init[i2]-1;i_face<conn_index_init[i2+1]-1;++i_face)
103 for(i_conn=conn_face_index_init[i_face]-1;i_conn<conn_face_index_init[i_face+1]-1;++i_conn)
105 conn_renum[iter_conn]=conn_init[i_conn];
108 conn_face_index_renum[iter_face+1]=iter_conn+1;
111 conn_index_renum[i_cell+1]=iter_face+1;
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));
117 delete[] conn_index_renum;
118 delete[] conn_face_index_renum;
121 else if (Type==MED_POLYGON)
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];
130 conn_index_renum[0]=1;
131 for(int i=0;i<nb_cell;++i)
134 for(int k=conn_index_init[i2];k<conn_index_init[i2+1];++k)
136 conn_renum[iter]=conn_init[k-1];
139 conn_index_renum[i+1]=iter+1;
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));
145 delete[] conn_index_renum;
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];
156 conn_index_renum[0]=1;
157 for(int i=0;i<nb_cell;++i)
160 for(int k=conn_index_init[i2];k<conn_index_init[i2+1];++k)
162 conn_renum[iter]=conn_init[k-1];
165 conn_index_renum[i+1]=iter+1;
168 CONNECTIVITY* myConnectivity=(CONNECTIVITY*)mesh.getConnectivityptr();
169 myConnectivity->setNodal(conn_renum,MED_CELL,Type,conn_index_renum);
171 delete[] conn_index_renum;
175 void changeFamily(MESH* mesh, const medGeometryElement& Type, const vector<int>& perm)
177 int nb_families=mesh->getNumberOfFamilies(MED_CELL);
178 for (int i=0;i<nb_families;++i)
180 const FAMILY* family=mesh->getFamily(MED_CELL,i+1);
181 if (!family->isOnAllElements())
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];
191 int main (int argc, char** argv)
193 double t_begin,t_read_st,t_read_mesh,t_compute_graph,t_connectiv,t_family,t_field;
197 cerr << "Usage : " << argv[0]
198 << " filename_in meshname method[BOOST/METIS] filename_out" << endl << endl;
201 string filename_in = argv[1];
202 string meshname = argv[2];
203 string type_renum = argv[3];
204 string filename_out = argv[4];
206 if(type_renum!="METIS" && type_renum!="BOOST")
208 cout << "The method has to be METIS or BOOST!" << endl;
212 string s="rm "+filename_out;
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();
225 cout << "There are many meshes in the file" << endl;
228 if(mesh_names[0].c_str()!=meshname)
230 cout << "Mesh name does not match" << endl;
233 vector<string> field_names;
234 vector<int> iternumber;
235 vector<int> ordernumber;
238 for (int ifield = 0; ifield < nb_fields; ifield++)
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++)
243 field_names.push_back(f_names[ifield]);
244 iternumber.push_back(iter->dt);
245 ordernumber.push_back(iter->it);
247 if(med_struct.getFieldType(f_names[ifield])==MED_EN::MED_REEL64)
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);
264 int nb_type=myMesh.getNumberOfTypes(MED_CELL);
267 cout << "Mesh must have only one type of cell" << endl;
270 const medGeometryElement *Types = myMesh.getTypes(MED_CELL);
271 medGeometryElement Type=Types[0];
274 MESH* workMesh=new MESH(myMesh);
275 cout << "Building the graph ";
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];
284 for(int i=0;i<nb_cell;++i)
286 for (list<int>::const_iterator it=neighbour[i].begin();it!=neighbour[i].end();++it)
291 graph_index[i+1]=count+1;
295 // Compute permutation
296 vector<int> iperm,perm;
297 Renumbering* renumb= RenumberingFactory(type_renum);
298 renumb->renumber(graph,graph_index,nb_cell,iperm,perm);
301 t_compute_graph=clock();
302 cout << " : " << (t_compute_graph-t_read_mesh)/(double) CLOCKS_PER_SEC << "s" << endl;
306 cout << "Computing connectivity";
308 MESH meshRenum(myMesh);
309 changeConnectivity(meshRenum,Type,nb_cell,iperm);
311 cout << " : " << (t_connectiv-t_compute_graph)/(double) CLOCKS_PER_SEC << "s" << endl;
315 cout << "Computing families ";
317 changeFamily(&meshRenum,Type,perm);
318 int drv3=meshRenum.addDriver(MED_DRIVER,filename_out,meshRenum.getName());
319 meshRenum.write(drv3);
321 cout << " : " << (t_family-t_connectiv)/(double) CLOCKS_PER_SEC << "s" << endl;
325 cout << "Computing fields ";
328 for(int ifield=0;ifield<nb_fields_tot;++ifield)
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)
337 if(typesOfSupport[t]==Type)
345 for(int i=0;i<mySupport->getNumberOfElements(Type);++i)
347 for(int j=0;j<newField.getNumberOfComponents();++j)
349 newField.setValueIJ(i+1,j+1,myField.getValueIJ(iperm[i],j+1));
353 int drv=newField.addDriver(MED_DRIVER,filename_out,field_names[ifield]);
357 cout << " : " << (t_field-t_family)/(double) CLOCKS_PER_SEC << "s" << endl;
360 delete[] graph_index;