Hermes2D  3.0
mesh_reader_exodusii.cpp
1 // This file is part of Hermes2D
2 //
3 // Hermes2D is free software; you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation; either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D 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
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D; if not, see <http://www.gnu.prg/licenses/>.
15 
16 #include "config.h"
17 #ifdef WITH_EXODUSII
18 
19 #include "mesh_reader_exodusii.h"
20 #include "mesh.h"
21 #include <exodusII.h>
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  MeshReaderExodusII::MeshReaderExodusII()
27  {
28  }
29 
30  MeshReaderExodusII::~MeshReaderExodusII()
31  {
32  }
33 
34  struct Vertex
35  {
36  double x, y;
37  };
38 
39  // needed for comparison of vertices inside std::map
40  struct VCompare
41  {
42  bool operator()(Vertex a, Vertex b) const
43  {
44  if (a.x < b.x) return true;
45  else if (a.x > b.x) return false;
46  else
47  {
48  if (a.y < b.y) return true;
49  return false;
50  }
51  }
52  };
53 
54  void MeshReaderExodusII::load(const char *file_name, MeshSharedPtr mesh)
55  {
56  int err;
57  // use float or double
58  int cpu_ws = sizeof(double);
59  // store variables as doubles
60  int io_ws = 8;
61  float version;
62  int exoid = ex_open(file_name, EX_READ, &cpu_ws, &io_ws, &version);
63 
64  // read initialization parameters
65  int n_dims, n_nodes, n_elems, n_eblocks, n_nodesets, n_sidesets;
66  char title[MAX_LINE_LENGTH + 1];
67  err = ex_get_init(exoid, title, &n_dims, &n_nodes, &n_elems, &n_eblocks, &n_nodesets, &n_sidesets);
68  if (n_dims != 2)
69  throw Hermes::Exceptions::Exception("File '%s' does not contain 2D mesh", file_name);
70 
71  // load coordinates
72  double *x = new double[n_nodes];
73  double *y = new double[n_nodes];
74  err = ex_get_coord(exoid, x, y, nullptr);
75 
76  // remove duplicate vertices and build renumbering map
77  // map for eliminating duplicities
78  std::map<Vertex, int, VCompare> vtx_list;
79  // reindexing map
80  std::map<int, int> vmap;
81  // vertex array
82  std::vector<Vertex> vtx_arr;
83  int vid = 0;
84  for (int i = 0; i < n_nodes; i++)
85  {
86  int k;
87  Vertex v = { x[i], y[i] };
88  if (vtx_list.count(v) == 0)
89  {
90  vtx_arr.push_back(v);
91  k = vid++;
92  vtx_list[v] = k;
93  }
94  else
95  k = vtx_list[v];
96 
97  vmap[i + 1] = k;
98  }
99  delete[] x;
100  delete[] y;
101 
102  int n_vtx = vtx_arr.size();
103  double2 *vtx = new double2[n_vtx];
104  for (int i = 0; i < n_vtx; i++)
105  {
106  vtx[i][0] = vtx_arr[i].x;
107  vtx[i][1] = vtx_arr[i].y;
108  }
109 
110  // number of triangles
111  int n_tri = 0;
112  // number of quads
113  int n_quad = 0;
114 
115  // get info about element blocks
116  int *eid_blocks = new int[n_eblocks];
117  err = ex_get_elem_blk_ids(exoid, eid_blocks);
118  // go over all element blocks
119  for (int i = 0; i < n_eblocks; i++)
120  {
121  int id = eid_blocks[i];
122 
123  // get block info
124  char elem_type[MAX_STR_LENGTH + 1];
125  int n_elems_in_blk, n_elem_nodes, n_attrs;
126  err = ex_get_elem_block(exoid, id, elem_type, &n_elems_in_blk, &n_elem_nodes, &n_attrs);
127 
128  if (n_elem_nodes == 3) n_tri += n_elems_in_blk;
129  else if (n_elem_nodes == 4) n_quad += n_elems_in_blk;
130  else
131  {
132  delete[] vtx;
133  throw Hermes::Exceptions::Exception("Unknown type of element");
134  }
135  }
136  // triangles
137  int3 *tri = n_tri > 0 ? new int3[n_tri] : nullptr;
138  std::string *tri_markers = n_tri > 0 ? new std::string[n_tri] : nullptr;
139  // quads
140  int4 *quad = n_quad > 0 ? new int4[n_quad] : nullptr;
141  std::string *quad_markers = n_quad > 0 ? new std::string[n_quad] : nullptr;
142 
143  // total number of elements
144  int n_els = n_tri + n_quad;
145  // elements
146  int **els = n_els > 0 ? new int *[n_els] : nullptr;
147  // number of vertices for each element
148  int *el_nv = n_els > 0 ? new int[n_els] : nullptr;
149 
150  int it = 0, iq = 0, iel = 0;
151  for (int i = 0; i < n_eblocks; i++)
152  {
153  int id = eid_blocks[i];
154 
155  // get block info
156  char elem_type[MAX_STR_LENGTH + 1];
157  int n_elems_in_blk, n_elem_nodes, n_attrs;
158  err = ex_get_elem_block(exoid, id, elem_type, &n_elems_in_blk, &n_elem_nodes, &n_attrs);
159 
160  // read connectivity array
161  int *connect = new int[n_elem_nodes * n_elems_in_blk];
162  err = ex_get_elem_conn(exoid, id, connect);
163 
164  // Update the mesh' internal array element_markers_conversion.
165  std::ostringstream string_stream;
166  string_stream << id;
167  std::string el_marker = string_stream.str();
168 
169  // This functions check if the user-supplied marker on this element has been
170  // already used, and if not, inserts it in the appropriate structure.
171  mesh->element_markers_conversion.insert_marker(el_marker);
172 
173  int ic = 0;
174  for (int j = 0; j < n_elems_in_blk; j++)
175  {
176  el_nv[iel] = n_elem_nodes;
177  if (n_elem_nodes == 3)
178  {
179  tri[it][0] = vmap[connect[ic++]];
180  tri[it][1] = vmap[connect[ic++]];
181  tri[it][2] = vmap[connect[ic++]];
182  tri_markers[it] = el_marker;
183  els[iel] = tri[it];
184  it++;
185  }
186  else if (n_elem_nodes == 4)
187  {
188  quad[iq][0] = vmap[connect[ic++]];
189  quad[iq][1] = vmap[connect[ic++]];
190  quad[iq][2] = vmap[connect[ic++]];
191  quad[iq][3] = vmap[connect[ic++]];
192  quad_markers[iq] = el_marker;
193  els[iel] = quad[iq];
194  iq++;
195  }
196  else
197  {
198  delete[] vtx;
199  throw Hermes::Exceptions::Exception("Unknown type of element");
200  }
201  iel++;
202  }
203  delete[] connect;
204  }
205  delete[] eid_blocks;
206 
207  // query number of side sets
208  int *sid_blocks = new int[n_sidesets];
209  err = ex_get_side_set_ids(exoid, sid_blocks);
210 
211  // go over the sidesets
212  // number of markers
213  int n_mark = 0;
214  for (int i = 0; i < n_sidesets; i++)
215  {
216  int sid = sid_blocks[i];
217  int n_sides_in_set, n_df_in_set;
218  err = ex_get_side_set_param(exoid, sid, &n_sides_in_set, &n_df_in_set);
219  n_mark += n_sides_in_set;
220  }
221  int2 *marks = new int2[n_mark];
222  std::string *bnd_markers = new std::string[n_mark];
223 
224  int im = 0;
225  for (int i = 0; i < n_sidesets; i++)
226  {
227  int sid = sid_blocks[i];
228  int n_sides_in_set, n_df_in_set;
229  err = ex_get_side_set_param(exoid, sid, &n_sides_in_set, &n_df_in_set);
230  int num_elem_in_set = n_sides_in_set;
231 
232  int *elem_list = new int[num_elem_in_set];
233  int *side_list = new int[n_sides_in_set];
234  err = ex_get_side_set(exoid, sid, elem_list, side_list);
235 
236  // Update the mesh' internal array boundary_markers_conversion.
237  std::ostringstream string_stream;
238  string_stream << sid;
239  std::string bnd_marker = string_stream.str();
240 
241  // This functions check if the user-supplied marker on this element has been
242  // already used, and if not, inserts it in the appropriate structure.
243  mesh->boundary_markers_conversion.insert_marker(bnd_marker);
244 
245  for (int j = 0; j < num_elem_in_set; j++)
246  {
247  // # of vertices of the element
248  int nv = el_nv[side_list[j] - 1];
249  int vt = side_list[j] - 1;
250  marks[im][0] = els[elem_list[j] - 1][vt];
251  marks[im][1] = els[elem_list[j] - 1][(vt + 1) % nv];
252  bnd_markers[im] = bnd_marker;
253  im++;
254  }
255 
256  delete[] elem_list;
257  delete[] side_list;
258  }
259  delete[] sid_blocks;
260 
261  // we are done
262  err = ex_close(exoid);
263 
264  mesh->create(n_vtx, vtx, n_tri, tri, tri_markers, n_quad, quad, quad_markers, n_mark, marks, bnd_markers);
265 
266  // clean-up
267  delete[] marks;
268  delete[] tri;
269  delete[] quad;
270  delete[] tri_markers;
271  delete[] quad_markers;
272  delete[] bnd_markers;
273  delete[] vtx;
274  delete[] el_nv;
275  delete[] els;
276  }
277  }
278 }
279 #endif
::xsd::cxx::tree::id< char, ncname > id
C++ type corresponding to the ID XML Schema built-in type.
Definition: adapt.h:24
virtual void load(const char *file_name, MeshSharedPtr mesh)
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.