Hermes2D  3.0
mesh_reader_h2d.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 "mesh.h"
17 #include "mesh_reader_h2d.h"
18 #include "refmap.h"
19 
20 namespace Hermes
21 {
22  namespace Hermes2D
23  {
24  extern unsigned g_mesh_seq;
25 
26  MeshReaderH2D::MeshReaderH2D()
27  {
28  }
29 
30  MeshReaderH2D::~MeshReaderH2D()
31  {
32  }
33 
34  Arc* MeshReaderH2D::load_arc(MeshSharedPtr mesh, MeshData *m, int id, Node** en, int &p1, int &p2, Arc* arc)
35  {
36  // read the end point indices
37  p1 = m->curv_first[id];
38  p2 = m->curv_second[id];
39 
40  *en = mesh->peek_edge_node(p1, p2);
41  if (*en == nullptr)
42  throw Hermes::Exceptions::MeshLoadFailureException("Curve #%d: edge %d-%d does not exist.", id, p1, p2);
43 
44  // get the number of control points
45  std::vector<double> vCP;
46 
47  // edge endpoints are also control points, with weight 1.0
48  arc->pt[0][0] = mesh->nodes[p1].x;
49  arc->pt[0][1] = mesh->nodes[p1].y;
50  arc->pt[0][2] = 1.0;
51  arc->pt[2][0] = mesh->nodes[p2].x;
52  arc->pt[2][1] = mesh->nodes[p2].y;
53  arc->pt[2][2] = 1.0;
54 
55  // read the arc angle
56  arc->angle = m->curv_third[id];
57  double a = (180.0 - arc->angle) / 180.0 * M_PI;
58 
59  // generate one control point
60  double x = 1.0 / tan(a * 0.5);
61  arc->pt[1][0] = 0.5*((arc->pt[2][0] + arc->pt[0][0]) + (arc->pt[2][1] - arc->pt[0][1]) * x);
62  arc->pt[1][1] = 0.5*((arc->pt[2][1] + arc->pt[0][1]) - (arc->pt[2][0] - arc->pt[0][0]) * x);
63  arc->pt[1][2] = cos((M_PI - a) * 0.5);
64 
65  // get the number of knot vector points
66  std::vector<double> vKnots;
67 
68  return arc;
69  }
70 
71  Nurbs* MeshReaderH2D::load_nurbs(MeshSharedPtr mesh, MeshData *m, int id, Node** en, int &p1, int &p2, Nurbs* nurbs)
72  {
73  // read the end point indices
74  p1 = m->curv_first[id];
75  p2 = m->curv_second[id];
76 
77  *en = mesh->peek_edge_node(p1, p2);
78  if (*en == nullptr)
79  throw Hermes::Exceptions::MeshLoadFailureException("Curve #%d: edge %d-%d does not exist.", id, p1, p2);
80 
81  nurbs->degree = m->curv_third[id];
82 
83  // get the number of control points
84  std::vector<double> vCP;
85 
86  int inner = 1, outer;
87  for (unsigned int i = 0; i < m->vars_[m->curv_inner_pts[id]].size(); ++i)
88  {
89  std::istringstream istr(m->vars_[m->curv_inner_pts[id]][i]);
90 
91  double dummy_dbl;
92  if (!(istr >> dummy_dbl))
93  vCP.push_back(atof(m->vars_[m->vars_[m->curv_inner_pts[id]][i]][0].c_str()));
94  else
95  vCP.push_back(atof(m->vars_[m->curv_inner_pts[id]][i].c_str()));
96  }
97  inner = vCP.size() / 3;
98 
99  nurbs->np = inner + 2;
100 
101  // edge endpoints are also control points, with weight 1.0
102  nurbs->pt = new double3[nurbs->np];
103  nurbs->pt[0][0] = mesh->nodes[p1].x;
104  nurbs->pt[0][1] = mesh->nodes[p1].y;
105  nurbs->pt[0][2] = 1.0;
106  nurbs->pt[inner + 1][0] = mesh->nodes[p2].x;
107  nurbs->pt[inner + 1][1] = mesh->nodes[p2].y;
108  nurbs->pt[inner + 1][2] = 1.0;
109 
110  // read inner control points
111  for (int i = 0; i < inner; i++)
112  {
113  for (int j = 0; j < 3; ++j)
114  nurbs->pt[i + 1][j] = vCP[3 * i + j];
115  }
116 
117  // get the number of knot vector points
118  std::vector<double> vKnots;
119 
120  inner = 0;
121  for (unsigned int i = 0; i < m->vars_[m->curv_knots[id]].size(); ++i)
122  {
123  std::istringstream istr(m->vars_[m->curv_knots[id]][i]);
124  double dummy_dbl;
125 
126  if (!(istr >> dummy_dbl))
127  vKnots.push_back(atof(m->vars_[m->vars_[m->curv_knots[id]][i]][0].c_str()));
128  else
129  vKnots.push_back(atof(m->vars_[m->curv_knots[id]][i].c_str()));
130  }
131  inner = vKnots.size();
132 
133  nurbs->nk = nurbs->degree + nurbs->np + 1;
134  outer = nurbs->nk - inner;
135  if ((outer & 1) == 1)
136  throw Hermes::Exceptions::MeshLoadFailureException("Curve #%d: incorrect number of knot points.", id);
137 
138  // knot vector is completed by 0.0 on the left and by 1.0 on the right
139  nurbs->kv = new double[nurbs->nk];
140 
141  for (int i = 0; i < outer / 2; i++)
142  nurbs->kv[i] = 0.0;
143 
144  if (inner) {
145  for (int i = outer / 2; i < inner + outer / 2; i++) {
146  nurbs->kv[i] = vKnots[i - (outer / 2)];
147  }
148  }
149 
150  for (int i = outer / 2 + inner; i < nurbs->nk; i++)
151  nurbs->kv[i] = 1.0;
152 
153  return nurbs;
154  }
155 
156  Curve* MeshReaderH2D::load_curve(MeshSharedPtr mesh, MeshData *m, int id, Node** en, int &p1, int &p2)
157  {
158  if (m->curv_nurbs[id])
159  return load_nurbs(mesh, m, id, en, p1, p2, new Nurbs);
160  else
161  return load_arc(mesh, m, id, en, p1, p2, new Arc);
162  }
163 
164  void MeshReaderH2D::load(const char *filename, MeshSharedPtr mesh)
165  {
166  // Check if file exists
167  std::ifstream s(filename);
168  if (!s.good())
169  throw Hermes::Exceptions::MeshLoadFailureException("Mesh file not found.");
170  s.close();
171 
172  int i, j, n;
173  Node* en;
174  bool debug = false;
175 
176  mesh->free();
177 
178  std::string fName(filename);
179  MeshData m(fName);
180  m.parse_mesh();
181 
183 
184  n = m.n_vert;
185  if (n < 0) throw Hermes::Exceptions::MeshLoadFailureException("File %s: 'vertices' must be a list.", filename);
186  if (n < 2) throw Hermes::Exceptions::MeshLoadFailureException("File %s: invalid number of vertices.", filename);
187 
188  // create a hash table large enough
190  while (size < 8 * n) size *= 2;
191  mesh->init(size);
192 
193  // create top-level vertex nodes
194  for (i = 0; i < n; i++)
195  {
196  Node* node = mesh->nodes.add();
197  assert(node->id == i);
198  node->ref = TOP_LEVEL_REF;
199  node->type = HERMES_TYPE_VERTEX;
200  node->bnd = 0;
201  node->p1 = node->p2 = -1;
202  node->next_hash = nullptr;
203  node->x = m.x_vertex[i];
204  node->y = m.y_vertex[i];
205  }
206  mesh->ntopvert = n;
207 
209 
210  n = m.n_el;
211  if (n < 0) throw Hermes::Exceptions::MeshLoadFailureException("File %s: 'elements' must be a list.", filename);
212  if (n < 1) throw Hermes::Exceptions::MeshLoadFailureException("File %s: no elements defined.", filename);
213 
214  // create elements
215  mesh->nactive = 0;
216  for (i = 0; i < n; i++)
217  {
218  // read and check vertex indices
219  int nv;
220  if (m.en4[i] == -1) nv = 4;
221  else nv = 5;
222 
223  int* idx = new int[nv - 1];
224  std::string el_marker;
225  if (!nv) {
226  mesh->elements.skip_slot()->cm = nullptr;
227  continue;
228  }
229 
230  if (nv < 4 || nv > 5)
231  {
232  delete[] idx;
233  throw Hermes::Exceptions::MeshLoadFailureException("File %s: element #%d: wrong number of vertex indices.", filename, i);
234  }
235 
236  if (nv == 4) {
237  idx[0] = m.en1[i];
238  idx[1] = m.en2[i];
239  idx[2] = m.en3[i];
240 
241  el_marker = m.e_mtl[i];
242  }
243  else {
244  idx[0] = m.en1[i];
245  idx[1] = m.en2[i];
246  idx[2] = m.en3[i];
247  idx[3] = m.en4[i];
248 
249  el_marker = m.e_mtl[i];
250  }
251  for (j = 0; j < nv - 1; j++)
252  if (idx[j] < 0 || idx[j] >= mesh->ntopvert)
253  {
254  delete[] idx;
255  throw Hermes::Exceptions::MeshLoadFailureException("File %s: error creating element #%d: vertex #%d does not exist.", filename, i, idx[j]);
256  }
257 
258  Node *v0 = &mesh->nodes[idx[0]], *v1 = &mesh->nodes[idx[1]], *v2 = &mesh->nodes[idx[2]];
259 
260  int marker;
261 
262  // This functions check if the user-supplied marker on this element has been
263  // already used, and if not, inserts it in the appropriate structure.
264  mesh->element_markers_conversion.insert_marker(el_marker);
265  marker = mesh->element_markers_conversion.get_internal_marker(el_marker).marker;
266 
267  if (nv == 4) {
268  Mesh::check_triangle(i, v0, v1, v2);
269  mesh->create_triangle(marker, v0, v1, v2, nullptr);
270  }
271  else {
272  Node *v3 = &mesh->nodes[idx[3]];
273  Mesh::check_quad(i, v0, v1, v2, v3);
274  mesh->create_quad(marker, v0, v1, v2, v3, nullptr);
275  }
276 
277  mesh->nactive++;
278 
279  delete[] idx;
280  }
281  mesh->nbase = n;
282 
284  if (m.n_bdy > 0)
285  {
286  n = m.n_bdy;
287 
288  // read boundary data
289  for (i = 0; i < n; i++)
290  {
291  int v1, v2, marker;
292  v1 = m.bdy_first[i];
293  v2 = m.bdy_second[i];
294 
295  en = mesh->peek_edge_node(v1, v2);
296  if (en == nullptr)
297  throw Hermes::Exceptions::MeshLoadFailureException("File %s: boundary data #%d: edge %d-%d does not exist", filename, i, v1, v2);
298 
299  std::string bnd_marker;
300  bnd_marker = m.bdy_type[i];
301 
302  // This functions check if the user-supplied marker on this element has been
303  // already used, and if not, inserts it in the appropriate structure.
304  mesh->boundary_markers_conversion.insert_marker(bnd_marker);
305  marker = mesh->boundary_markers_conversion.get_internal_marker(bnd_marker).marker;
306 
307  en->marker = marker;
308  }
309  }
310 
311  Node* node;
312  for_all_edge_nodes(node, mesh)
313  {
314  if (node->ref < 2)
315  {
316  mesh->nodes[node->p1].bnd = 1;
317  mesh->nodes[node->p2].bnd = 1;
318  node->bnd = 1;
319  }
320  }
321 
323  if (m.n_curv > 0)
324  {
325  n = m.n_curv;
326  if (n < 0) throw Hermes::Exceptions::MeshLoadFailureException("File %s: 'curves' must be a list.", filename);
327 
328  // load curved edges
329  for (i = 0; i < n; i++)
330  {
331  // load the control points, knot vector, etc.
332  Node* en;
333  int p1, p2;
334 
335  Curve* curve = load_curve(mesh, &m, i, &en, p1, p2);
336 
337  // assign the arc to the elements sharing the edge node
338  MeshUtil::assign_curve(en, curve, p1, p2);
339  }
340  }
341 
342  // update refmap coeffs of curvilinear elements
343  Element* e;
344  for_all_used_elements(e, mesh)
345  {
346  if (e->cm != nullptr)
347  e->cm->update_refmap_coeffs(e);
348  this->ref_map.set_element_iro_cache(e);
349  }
350 
352  if (m.n_ref > 0)
353  {
354  n = m.n_ref;
355  if (n < 0) throw Hermes::Exceptions::MeshLoadFailureException("File %s: 'refinements' must be a list.", filename);
356 
357  // perform initial refinements
358  for (i = 0; i < n; i++)
359  {
360  int id, ref;
361  id = m.ref_elt[i];
362  ref = m.ref_type[i];
363  mesh->refine_element_id(id, ref);
364  }
365  }
366  mesh->ninitial = mesh->elements.get_num_items();
367 
368  mesh->seq = g_mesh_seq++;
369  if (HermesCommonApi.get_integral_param_value(checkMeshesOnLoad))
370  mesh->initial_single_check();
371  }
372 
373  void MeshReaderH2D::save_refinements(MeshSharedPtr mesh, FILE* f, Element* e, int id, bool& first)
374  {
375  if (e->active) return;
376  fprintf(f, first ? "refinements =\n{\n" : ",\n"); first = false;
377  if (e->bsplit())
378  {
379  fprintf(f, " { %d, 0 }", id);
380  int sid = mesh->seq; mesh->seq += 4;
381  for (int i = 0; i < 4; i++)
382  save_refinements(mesh, f, e->sons[i], sid + i, first);
383  }
384  else if (e->hsplit())
385  {
386  fprintf(f, " { %d, 1 }", id);
387  int sid = mesh->seq; mesh->seq += 2;
388  save_refinements(mesh, f, e->sons[0], sid, first);
389  save_refinements(mesh, f, e->sons[1], sid + 1, first);
390  }
391  else
392  {
393  fprintf(f, " { %d, 2 }", id);
394  int sid = mesh->seq; mesh->seq += 2;
395  save_refinements(mesh, f, e->sons[2], sid, first);
396  save_refinements(mesh, f, e->sons[3], sid + 1, first);
397  }
398  }
399 
400  void MeshReaderH2D::save_curve(MeshSharedPtr mesh, FILE* f, int p1, int p2, Curve* curve)
401  {
402  if (curve->type == ArcType)
403  {
404  fprintf(f, " [ %d, %d, %.16g ]", p1, p2, ((Arc*)curve)->angle);
405  }
406  else
407  {
408  int inner = ((Nurbs*)curve)->np - 2;
409  int outer = ((Nurbs*)curve)->nk - inner;
410  fprintf(f, " ú %d, %d, %d, ú ", p1, p2, ((Nurbs*)curve)->degree);
411  for (int i = 1; i < ((Nurbs*)curve)->np - 1; i++)
412  fprintf(f, "ú %.16g, %.16g, %.16g ]%s ",
413  ((Nurbs*)curve)->pt[i][0], ((Nurbs*)curve)->pt[i][1], ((Nurbs*)curve)->pt[i][2],
414  i < ((Nurbs*)curve)->np - 2 ? "," : "");
415 
416  fprintf(f, "],[ ");
417  int max = ((Nurbs*)curve)->nk - (((Nurbs*)curve)->degree + 1);
418  for (int i = ((Nurbs*)curve)->degree + 1; i < max; i++)
419  fprintf(f, "%.16g%s", ((Nurbs*)curve)->kv[i], i < max - 1 ? "," : "");
420  fprintf(f, "] ]");
421  }
422  }
423 
424  void MeshReaderH2D::save(const char* filename, MeshSharedPtr mesh)
425  {
426  int i, mrk;
427  Element* e;
428 
429  // open output file
430  FILE* f = fopen(filename, "w");
431  if (f == nullptr) throw Hermes::Exceptions::MeshLoadFailureException("Could not create mesh file.");
432  //fprintf(f, "# hermes2d saved mesh\n\n");
433 
434  // save vertices
435  fprintf(f, "vertices =\n[\n");
436  for (i = 0; i < mesh->ntopvert; i++)
437  fprintf(f, " [ %.16g, %.16g ]%s\n", mesh->nodes[i].x, mesh->nodes[i].y, (i < mesh->ntopvert - 1 ? "," : ""));
438 
439  // save elements
440  fprintf(f, "]\n\nelements =\n[");
441  bool first = true;
442  for (i = 0; i < mesh->get_num_base_elements(); i++)
443  {
444  const char* nl = first ? "\n" : ",\n"; first = false;
445  e = mesh->get_element_fast(i);
446  if (!e->used)
447  fprintf(f, "%s [ ]", nl);
448  else if (e->is_triangle())
449  fprintf(f, "%s [ %d, %d, %d, \"%s\" ]", nl, e->vn[0]->id, e->vn[1]->id, e->vn[2]->id, mesh->get_element_markers_conversion().get_user_marker(e->marker).marker.c_str());
450  else
451  fprintf(f, "%s [ %d, %d, %d, %d, \"%s\" ]", nl, e->vn[0]->id, e->vn[1]->id, e->vn[2]->id, e->vn[3]->id, mesh->get_element_markers_conversion().get_user_marker(e->marker).marker.c_str());
452  }
453 
454  // save boundary markers
455  fprintf(f, "\n]\n\nboundaries =\n[");
456  first = true;
457  for_all_base_elements(e, mesh)
458  {
459  for (unsigned char i = 0; i < e->get_nvert(); i++)
460  if ((mrk = MeshUtil::get_base_edge_node(e, i)->marker))
461  {
462  const char* nl = first ? "\n" : ",\n"; first = false;
463  fprintf(f, "%s [ %d, %d, \"%s\" ]", nl, e->vn[i]->id, e->vn[e->next_vert(i)]->id, mesh->boundary_markers_conversion.get_user_marker(mrk).marker.c_str());
464  }
465  }
466  fprintf(f, "\n]\n\n");
467 
468  // save curved edges
469  first = true;
470  for_all_base_elements(e, mesh)
471  {
472  if (e->is_curved())
473  {
474  for (unsigned char i = 0; i < e->get_nvert(); i++)
475  if (e->cm->curves[i] != nullptr)
476  {
477  fprintf(f, first ? "curves =\n[\n" : ",\n"); first = false;
478  save_curve(mesh, f, e->vn[i]->id, e->vn[e->next_vert(i)]->id, e->cm->curves[i]);
479  }
480  if (!first) fprintf(f, "\n]\n\n");
481  }
482  }
483  // save refinements
484  unsigned temp = mesh->seq;
485  mesh->seq = mesh->nbase;
486  first = true;
487  for_all_base_elements(e, mesh)
488  save_refinements(mesh, f, e, e->id, first);
489  if (!first) fprintf(f, "\n]\n\n");
490 
491  mesh->seq = temp;
492  fclose(f);
493  }
494  }
495 }
::xsd::cxx::tree::id< char, ncname > id
C++ type corresponding to the ID XML Schema built-in type.
std::vector< std::string > e_mtl
Element markers – single word strings.
Definition: mesh_data.h:82
double x
vertex node coordinates
Definition: element.h:63
void parse_mesh(void)
This function parses a given input mesh file line by line and extracts the necessary information into...
Definition: mesh_data.cpp:152
unsigned type
0 = vertex node; 1 = edge node
Definition: element.h:52
int id
node id number
Definition: element.h:48
Definition: adapt.h:24
std::vector< int > ref_elt
List of elements to be refined.
Definition: mesh_data.h:107
Node * next_hash
next node in hash synonym list
Definition: element.h:77
Stores one element of a mesh.
Definition: element.h:107
int n_ref
Number of elements with specified refinements.
Definition: mesh_data.h:65
RefMap ref_map
Reference mapping for detecting the inverse reference mapping order.
Definition: mesh_reader.h:43
int n_vert
Number of vertices.
Definition: mesh_data.h:57
CurvMap * cm
curved mapping, nullptr if not curvilinear
Definition: element.h:188
int n_el
Number of elements.
Definition: mesh_data.h:59
std::vector< std::string > bdy_type
Boundary name.
Definition: mesh_data.h:89
Class to stored 2d mesh parameters. The MeshData class organizes all the necessary data structures re...
Definition: mesh_data.h:40
Stores one node of a mesh.
Definition: element.h:45
std::vector< int > en4
Nodes with local node number 4. Only for quadrilateral elements. For triangular elements it is set to...
Definition: mesh_data.h:79
int n_curv
Number of curved edges (including NURBS curves)
Definition: mesh_data.h:63
std::vector< int > en3
Nodes with local node number 3.
Definition: mesh_data.h:77
std::vector< int > en2
Nodes with local node number 2.
Definition: mesh_data.h:75
std::vector< int > ref_type
List of element refinement type.
Definition: mesh_data.h:109
int n_bdy
Number of boundary edges.
Definition: mesh_data.h:61
std::vector< int > bdy_first
First node of a boundary edge.
Definition: mesh_data.h:85
int p1
parent id numbers
Definition: element.h:75
unsigned bnd
1 = boundary node; 0 = inner node
Definition: element.h:54
std::vector< int > en1
Nodes with local node number 1.
Definition: mesh_data.h:73
std::vector< double > y_vertex
y-coordinate of the vertices
Definition: mesh_data.h:70
std::vector< double > x_vertex
x-coordinate of the vertices
Definition: mesh_data.h:68
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.
std::vector< int > bdy_second
Second node of a boundary edge.
Definition: mesh_data.h:87
unsigned ref
the number of elements using the node
Definition: element.h:50
void update_refmap_coeffs(Element *e)
Definition: curved.cpp:758
virtual void load(const char *filename, MeshSharedPtr mesh)
Element * sons[H2D_MAX_ELEMENT_SONS]
son elements (up to four)
Definition: element.h:140
bool active
0 = active, no sons; 1 = inactive (refined), has sons
Definition: element.h:114
int marker
edge marker
Definition: element.h:68
static Node * get_base_edge_node(Element *base, int edge)
For internal use.
Definition: mesh_util.cpp:95
static const int H2D_DEFAULT_HASH_SIZE
32K entries
Definition: hash.h:49
static void assign_curve(Node *en, Curve *curve, int p1, int p2)
For internal use.
Definition: mesh_util.cpp:23