Hermes2D  3.0
mesh_reader_h2d_bson.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_reader_h2d_bson.h"
17 
18 #ifdef WITH_BSON
19 
20 #include "mesh.h"
21 #include "refmap.h"
22 #include "api2d.h"
23 
24 using namespace std;
25 
26 namespace Hermes
27 {
28  namespace Hermes2D
29  {
30  MeshReaderH2DBSON::MeshReaderH2DBSON()
31  {
32  }
33 
34  MeshReaderH2DBSON::~MeshReaderH2DBSON()
35  {
36  }
37 
38  void MeshReaderH2DBSON::load(const char *filename, MeshSharedPtr mesh)
39  {
40  if (!mesh)
41  throw Exceptions::NullException(1);
42 
43  mesh->free();
44 
45  FILE *fpr;
46  fpr = fopen(filename, "rb");
47 
48  // file size:
49  fseek(fpr, 0, SEEK_END);
50  int size = ftell(fpr);
51  rewind(fpr);
52 
53  // allocate memory to contain the whole file:
54  char *datar = malloc_with_check<char>(size);
55  fread(datar, size, 1, fpr);
56  fclose(fpr);
57 
58  bson br;
59  bson_init_finished_data(&br, datar, 0);
60 
61  bson_iterator it;
62  bson sub;
63 
64  // Vertices //
65  bson_find(&it, &br, "vertex-count");
66  int vertices_count = bson_iterator_int(&it);
67 
68  // Initialize mesh.
69  int mesh_size = HashTable::H2D_DEFAULT_HASH_SIZE;
70  while (mesh_size < 8 * vertices_count)
71  mesh_size *= 2;
72  mesh->init(mesh_size);
73 
74  // Create top-level vertex nodes.
75  double* vertex_xes = new double[vertices_count];
76  double* vertex_yes = new double[vertices_count];
77 
78  // Xes
79  bson_iterator it_coeffs;
80  bson_find(&it_coeffs, &br, "vertex-x");
81  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
82  bson_iterator_init(&it, &sub);
83  int index_coeff = 0;
84  while (bson_iterator_next(&it))
85  vertex_xes[index_coeff++] = bson_iterator_double(&it);
86 
87  // Yes
88  bson_find(&it_coeffs, &br, "vertex-y");
89  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
90  bson_iterator_init(&it, &sub);
91  index_coeff = 0;
92  while (bson_iterator_next(&it))
93  vertex_yes[index_coeff++] = bson_iterator_double(&it);
94 
95  // Assign.
96  for (int vertex_i = 0; vertex_i < vertices_count; vertex_i++)
97  {
98  Node* node = mesh->nodes.add();
99  assert(node->id == vertex_i);
100  node->ref = TOP_LEVEL_REF;
101  node->type = HERMES_TYPE_VERTEX;
102  node->bnd = 0;
103  node->p1 = node->p2 = -1;
104  node->next_hash = nullptr;
105  node->x = vertex_xes[vertex_i];
106  node->y = vertex_yes[vertex_i];
107  }
108  mesh->ntopvert = vertices_count;
109 
110  delete[] vertex_xes;
111  delete[] vertex_yes;
112 
113  // Elements //
114  bson_find(&it, &br, "element-count");
115  int element_count = bson_iterator_int(&it);
116  mesh->nbase = mesh->nactive = mesh->ninitial = element_count;
117 
118  int* vertex_0 = new int[element_count];
119  int* vertex_1 = new int[element_count];
120  int* vertex_2 = new int[element_count];
121  int* vertex_3 = new int[element_count];
122  const char** markers = new const char*[element_count];
123 
124  // Vertex_0s
125  bson_find(&it_coeffs, &br, "element-1");
126  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
127  bson_iterator_init(&it, &sub);
128  index_coeff = 0;
129  while (bson_iterator_next(&it))
130  vertex_0[index_coeff++] = bson_iterator_int(&it);
131  // Vertex_1s
132  bson_find(&it_coeffs, &br, "element-2");
133  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
134  bson_iterator_init(&it, &sub);
135  index_coeff = 0;
136  while (bson_iterator_next(&it))
137  vertex_1[index_coeff++] = bson_iterator_int(&it);
138  // Vertex_2s
139  bson_find(&it_coeffs, &br, "element-3");
140  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
141  bson_iterator_init(&it, &sub);
142  index_coeff = 0;
143  while (bson_iterator_next(&it))
144  vertex_2[index_coeff++] = bson_iterator_int(&it);
145  // Vertex_3s
146  bson_find(&it_coeffs, &br, "element-4");
147  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
148  bson_iterator_init(&it, &sub);
149  index_coeff = 0;
150  while (bson_iterator_next(&it))
151  vertex_3[index_coeff++] = bson_iterator_int(&it);
152  // Markers
153  bson_find(&it_coeffs, &br, "element-marker");
154  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
155  bson_iterator_init(&it, &sub);
156  index_coeff = 0;
157  while (bson_iterator_next(&it))
158  markers[index_coeff++] = bson_iterator_string(&it);
159 
160  Element* e;
161  for (int element_i = 0; element_i < element_count; element_i++)
162  {
163  mesh->element_markers_conversion.insert_marker(markers[element_i]);
164  if (vertex_3[element_i] != -1)
165  e = mesh->create_quad(mesh->element_markers_conversion.get_internal_marker(markers[element_i]).marker,
166  &mesh->nodes[vertex_0[element_i]],
167  &mesh->nodes[vertex_1[element_i]],
168  &mesh->nodes[vertex_2[element_i]],
169  &mesh->nodes[vertex_3[element_i]],
170  nullptr);
171  else
172  e = mesh->create_triangle(mesh->element_markers_conversion.get_internal_marker(markers[element_i]).marker,
173  &mesh->nodes[vertex_0[element_i]],
174  &mesh->nodes[vertex_1[element_i]],
175  &mesh->nodes[vertex_2[element_i]],
176  nullptr);
177  }
178 
179  delete[] vertex_0;
180  delete[] vertex_1;
181  delete[] vertex_2;
182  delete[] vertex_3;
183  delete[] markers;
184 
185  // Boundaries //
186  bson_find(&it, &br, "edge-count");
187  int edges_count = bson_iterator_int(&it);
188 
189  int* vertex_0_edge = new int[edges_count];
190  int* vertex_1_edge = new int[edges_count];
191  const char** edge_markers = new const char*[edges_count];
192 
193  // Vertex_0s
194  bson_find(&it_coeffs, &br, "edge-1");
195  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
196  bson_iterator_init(&it, &sub);
197  index_coeff = 0;
198  while (bson_iterator_next(&it))
199  vertex_0_edge[index_coeff++] = bson_iterator_int(&it);
200  // Vertex_1s
201  bson_find(&it_coeffs, &br, "edge-2");
202  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
203  bson_iterator_init(&it, &sub);
204  index_coeff = 0;
205  while (bson_iterator_next(&it))
206  vertex_1_edge[index_coeff++] = bson_iterator_int(&it);
207  // Markers
208  bson_find(&it_coeffs, &br, "edge-marker");
209  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
210  bson_iterator_init(&it, &sub);
211  index_coeff = 0;
212  while (bson_iterator_next(&it))
213  edge_markers[index_coeff++] = bson_iterator_string(&it);
214 
215  Node* en;
216  for (unsigned int edge_i = 0; edge_i < edges_count; edge_i++)
217  {
218  en = mesh->peek_edge_node(vertex_0_edge[edge_i], vertex_1_edge[edge_i]);
219  if (en == nullptr)
220  throw Hermes::Exceptions::MeshLoadFailureException("Boundary data #%d: edge %d-%d does not exist.", edge_i, vertex_0_edge[edge_i], vertex_1_edge[edge_i]);
221 
222  std::string edge_marker = edge_markers[edge_i];
223 
224  // This functions check if the user-supplied marker on this element has been
225  // already used, and if not, inserts it in the appropriate structure.
226  mesh->boundary_markers_conversion.insert_marker(edge_marker);
227 
228  en->marker = mesh->boundary_markers_conversion.get_internal_marker(edge_marker).marker;
229  }
230 
231  Node* node;
232  for_all_edge_nodes(node, mesh)
233  {
234  if (node->ref < 2)
235  {
236  mesh->nodes[node->p1].bnd = 1;
237  mesh->nodes[node->p2].bnd = 1;
238  node->bnd = 1;
239  }
240  }
241 
242  // check that all boundary edges have a marker assigned
243  for_all_edge_nodes(en, mesh)
244  if (en->ref < 2 && en->marker == 0)
245  this->warn("Boundary edge node does not have a boundary marker.");
246 
247  delete[] vertex_0_edge;
248  delete[] vertex_1_edge;
249  delete[] edge_markers;
250 
251  // Curves //
252 
253  // Save arcs.
254  // - count.
255  bson_find(&it, &br, "arc-count");
256  int arc_count = bson_iterator_int(&it);
257 
258  int* p1s = new int[arc_count];
259  int* p2s = new int[arc_count];
260  double* angles = new double[arc_count];
261 
262  // p1s
263  bson_find(&it_coeffs, &br, "arcs-p1");
264  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
265  bson_iterator_init(&it, &sub);
266  index_coeff = 0;
267  while (bson_iterator_next(&it))
268  p1s[index_coeff++] = bson_iterator_int(&it);
269  // p2s
270  bson_find(&it_coeffs, &br, "arcs-p2");
271  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
272  bson_iterator_init(&it, &sub);
273  index_coeff = 0;
274  while (bson_iterator_next(&it))
275  p2s[index_coeff++] = bson_iterator_int(&it);
276  // angles
277  bson_find(&it_coeffs, &br, "arcs-angle");
278  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
279  bson_iterator_init(&it, &sub);
280  index_coeff = 0;
281  while (bson_iterator_next(&it))
282  angles[index_coeff++] = bson_iterator_double(&it);
283 
284  // Arcs //
285  for (unsigned int curves_i = 0; curves_i < arc_count; curves_i++)
286  {
287  // load the control points, knot vector, etc.
288  Node* en;
289  int p1, p2;
290  double angle;
291 
292  Curve* curve;
293 
294  // read the end point indices
295  p1 = p1s[curves_i];
296  p2 = p2s[curves_i];
297  angle = angles[curves_i];
298 
299  curve = MeshUtil::load_arc(mesh, curves_i, &en, p1, p2, angle);
300 
301  // assign the arc to the elements sharing the edge node
302  MeshUtil::assign_curve(en, curve, p1, p2);
303  }
304 
305  // update refmap coeffs of curvilinear elements
306  for_all_used_elements(e, mesh)
307  {
308  if (e->cm != nullptr)
309  e->cm->update_refmap_coeffs(e);
310  RefMap::set_element_iro_cache(e);
311  }
312 
313  delete[] p1s;
314  delete[] p2s;
315  delete[] angles;
316 
317  // Refinements.
318  bson_find(&it, &br, "refinement-count");
319  int refinement_count = bson_iterator_int(&it);
320  int* refined_elements = new int[refinement_count];
321  int* refinement_types = new int[refinement_count];
322 
323  // refinement-id
324  bson_find(&it_coeffs, &br, "refinement-id");
325  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
326  bson_iterator_init(&it, &sub);
327  index_coeff = 0;
328  while (bson_iterator_next(&it))
329  refined_elements[index_coeff++] = bson_iterator_int(&it);
330  // refinement-id
331  bson_find(&it_coeffs, &br, "refinement-type");
332  bson_iterator_subobject_init(&it_coeffs, &sub, 0);
333  bson_iterator_init(&it, &sub);
334  index_coeff = 0;
335  while (bson_iterator_next(&it))
336  refinement_types[index_coeff++] = bson_iterator_int(&it);
337 
338  bson_destroy(&br);
339  free_with_check(datar);
340 
341  // perform initial refinements
342  for (unsigned int refinement_i = 0; refinement_i < refinement_count; refinement_i++)
343  {
344  int element_id = refined_elements[refinement_i];
345  int refinement_type = refinement_types[refinement_i];
346  if (refinement_type == -1)
347  mesh->unrefine_element_id(element_id);
348  else
349  mesh->refine_element_id(element_id, refinement_type);
350  }
351 
352  delete[] refined_elements;
353  delete[] refinement_types;
354 
355  if (HermesCommonApi.get_integral_param_value(checkMeshesOnLoad))
356  mesh->initial_single_check();
357  }
358 
359  void MeshReaderH2DBSON::save(const char *filename, MeshSharedPtr mesh)
360  {
361  // Utility pointer.
362  Element* e;
363 
364  // Init bson
365  bson bw;
366  bson_init(&bw);
367 
368  // Save vertices
369  // - count.
370  bson_append_int(&bw, "vertex-count", mesh->ntopvert);
371  // - xes
372  bson_append_start_array(&bw, "vertex-x");
373  for (int i = 0; i < mesh->ntopvert; i++)
374  bson_append_double(&bw, "c", mesh->nodes[i].x);
375  bson_append_finish_array(&bw);
376  // - yes
377  bson_append_start_array(&bw, "vertex-y");
378  for (int i = 0; i < mesh->ntopvert; i++)
379  bson_append_double(&bw, "c", mesh->nodes[i].y);
380  bson_append_finish_array(&bw);
381 
382  // Save elements
383  // - count.
384  bson_append_int(&bw, "element-count", mesh->get_num_base_elements());
385  // - first vertex ids
386  bson_append_start_array(&bw, "element-1");
387  for (int i = 0; i < mesh->get_num_base_elements(); i++)
388  {
389  e = mesh->get_element_fast(i);
390  bson_append_int(&bw, "c", e->vn[0]->id);
391  }
392  bson_append_finish_array(&bw);
393  // - second vertex ids
394  bson_append_start_array(&bw, "element-2");
395  for (int i = 0; i < mesh->get_num_base_elements(); i++)
396  {
397  e = mesh->get_element_fast(i);
398  bson_append_int(&bw, "c", e->vn[1]->id);
399  }
400  bson_append_finish_array(&bw);
401  // - third vertex ids
402  bson_append_start_array(&bw, "element-3");
403  for (int i = 0; i < mesh->get_num_base_elements(); i++)
404  {
405  e = mesh->get_element_fast(i);
406  bson_append_int(&bw, "c", e->vn[2]->id);
407  }
408  bson_append_finish_array(&bw);
409  // - first vertex ids
410  bson_append_start_array(&bw, "element-4");
411  for (int i = 0; i < mesh->get_num_base_elements(); i++)
412  {
413  e = mesh->get_element_fast(i);
414  if (e->is_triangle())
415  bson_append_int(&bw, "c", -1);
416  else
417  bson_append_int(&bw, "c", e->vn[3]->id);
418  }
419  bson_append_finish_array(&bw);
420  // - markers
421  bson_append_start_array(&bw, "element-marker");
422  for (int i = 0; i < mesh->get_num_base_elements(); i++)
423  {
424  e = mesh->get_element_fast(i);
425  bson_append_string(&bw, "c", mesh->get_element_markers_conversion().get_user_marker(e->marker).marker.c_str());
426  }
427  bson_append_finish_array(&bw);
428 
429  // Save boundary markers
430  // - count.
431  int edge_count = 0;
432  for_all_base_elements(e, mesh)
433  for (unsigned i = 0; i < e->get_nvert(); i++)
434  if (MeshUtil::get_base_edge_node(e, i)->marker)
435  edge_count++;
436  bson_append_int(&bw, "edge-count", edge_count);
437  // - ids 1
438  bson_append_start_array(&bw, "edge-1");
439  for_all_base_elements(e, mesh)
440  for (unsigned i = 0; i < e->get_nvert(); i++)
441  if (MeshUtil::get_base_edge_node(e, i)->marker)
442  bson_append_int(&bw, "c", e->vn[i]->id);
443  bson_append_finish_array(&bw);
444  // - ids 2
445  bson_append_start_array(&bw, "edge-2");
446  for_all_base_elements(e, mesh)
447  for (unsigned i = 0; i < e->get_nvert(); i++)
448  if (MeshUtil::get_base_edge_node(e, i)->marker)
449  bson_append_int(&bw, "c", e->vn[e->next_vert(i)]->id);
450  bson_append_finish_array(&bw);
451  // - markers
452  bson_append_start_array(&bw, "edge-marker");
453  for_all_base_elements(e, mesh)
454  for (unsigned i = 0; i < e->get_nvert(); i++)
455  if (MeshUtil::get_base_edge_node(e, i)->marker)
456  bson_append_string(&bw, "c", mesh->boundary_markers_conversion.get_user_marker(MeshUtil::get_base_edge_node(e, i)->marker).marker.c_str());
457  bson_append_finish_array(&bw);
458 
459  // Save arcs.
460  // - count.
461  int arc_count = 0;
462  for_all_base_elements(e, mesh)
463  {
464  if (e->is_curved())
465  for (unsigned i = 0; i < e->get_nvert(); i++)
466  if (e->cm->curves[i] != nullptr)
467  if (e->cm->curves[i]->type == ArcType)
468  arc_count++;
469  }
470  bson_append_int(&bw, "arc-count", arc_count);
471  // - p1.
472  bson_append_start_array(&bw, "arcs-p1");
473  for_all_base_elements(e, mesh)
474  {
475  if (e->is_curved())
476  for (unsigned i = 0; i < e->get_nvert(); i++)
477  if (e->cm->curves[i] != nullptr)
478  {
479  if (e->cm->curves[i]->type != ArcType)
480  continue;
481  bson_append_int(&bw, "c", e->vn[i]->id);
482  }
483  }
484  bson_append_finish_array(&bw);
485  // - p2.
486  bson_append_start_array(&bw, "arcs-p2");
487  for_all_base_elements(e, mesh)
488  {
489  if (e->is_curved())
490  for (unsigned i = 0; i < e->get_nvert(); i++)
491  if (e->cm->curves[i] != nullptr)
492  {
493  if (e->cm->curves[i]->type != ArcType)
494  continue;
495  bson_append_int(&bw, "c", e->vn[e->next_vert(i)]->id);
496  }
497  }
498  bson_append_finish_array(&bw);
499  // - angles.
500  bson_append_start_array(&bw, "arcs-angle");
501  for_all_base_elements(e, mesh)
502  {
503  if (e->is_curved())
504  for (unsigned i = 0; i < e->get_nvert(); i++)
505  if (e->cm->curves[i] != nullptr)
506  {
507  if (e->cm->curves[i]->type != ArcType)
508  continue;
509  bson_append_int(&bw, "c", ((Arc*)e->cm->curves[i])->angle);
510  }
511  }
512  bson_append_finish_array(&bw);
513 
514  // Save general NURBS.
515  // - count.
516  for_all_base_elements(e, mesh)
517  if (e->is_curved())
518  for (unsigned i = 0; i < e->get_nvert(); i++)
519  if (e->cm->curves[i] != nullptr)
520  if (e->cm->curves[i]->type != ArcType)
521  throw Exceptions::Exception("BSON mesh loader can not operate with general NURBS so far.");
522 
523  // Save refinements
524  // - count.
525  for (unsigned int refinement_i = 0; refinement_i < mesh->refinements.size(); refinement_i++)
526  bson_append_int(&bw, "refinement-count", mesh->refinements.size());
527  // - ids
528  bson_append_start_array(&bw, "refinement-id");
529  for (unsigned int refinement_i = 0; refinement_i < mesh->refinements.size(); refinement_i++)
530  bson_append_int(&bw, "c", mesh->refinements[refinement_i].first);
531  bson_append_finish_array(&bw);
532  // - type
533  bson_append_start_array(&bw, "refinement-type");
534  for (unsigned int refinement_i = 0; refinement_i < mesh->refinements.size(); refinement_i++)
535  bson_append_int(&bw, "c", mesh->refinements[refinement_i].second);
536  bson_append_finish_array(&bw);
537 
538  // Done.
539  bson_finish(&bw);
540 
541  // Write to disk.
542  FILE *fpw;
543  fpw = fopen(filename, "wb");
544  const char *dataw = (const char *)bson_data(&bw);
545  fwrite(dataw, bson_size(&bw), 1, fpw);
546  fclose(fpw);
547 
548  bson_destroy(&bw);
549  }
550 
551  void MeshReaderH2DBSON::save(const char *filename, std::vector<MeshSharedPtr> meshes)
552  {
553  // For mapping of physical coordinates onto top vertices.
554  std::map<std::pair<double, double>, unsigned int> points_to_vertices;
555  // For mapping of vertex pairs onto boundary edges.
556  std::map<std::pair<unsigned int, unsigned int>, unsigned int> vertices_to_boundaries;
557  // For mapping of vertex pairs onto curves.
558  std::map<std::pair<unsigned int, unsigned int>, bool> vertices_to_curves;
559 
560  std::vector<element_BSON> elements;
561  std::vector<edge_BSON> edges;
562  std::vector<vertex_BSON> vertices;
563  std::vector<arc_BSON> arcs;
564  std::vector<subdomain_BSON> subdomains;
565 
566  bool* baseElementsSaved = new bool[meshes[0]->get_num_base_elements()];
567  memset(baseElementsSaved, 0, sizeof(bool)* meshes[0]->get_num_base_elements());
568 
569 #pragma region Save to structures
570 
571  for (unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
572  {
573  bool hasAllElements = (meshes[meshes_i]->get_num_used_base_elements() == meshes[meshes_i]->get_num_base_elements());
574 
575  subdomain_BSON subdomain;
576 
577  // Mapping of top vertices of subdomains to the global mesh.
578  std::map<unsigned int, unsigned int> vertices_to_vertices;
579 
580  // Utility pointer.
581  Element* e;
582 
583  // save vertices
584  for (int i = 0; i < meshes[meshes_i]->ntopvert; i++)
585  {
586  // Look for the coordinates of this vertex.
587  // If found, then insert the pair <this vertex number, the found vertex number> into vertices_to_vertices dictionary.
588  // If not, insert.
589  std::map<std::pair<double, double>, unsigned int>::iterator it = points_to_vertices.find(std::pair<double, double>(meshes[meshes_i]->nodes[i].x, meshes[meshes_i]->nodes[i].y));
590  if (it != points_to_vertices.end())
591  vertices_to_vertices.insert(std::pair<unsigned int, unsigned int>(i, it->second));
592  else
593  {
594  int new_i = points_to_vertices.size();
595  vertices_to_vertices.insert(std::pair<unsigned int, unsigned int>(i, new_i));
596  points_to_vertices.insert(std::pair<std::pair<double, double>, unsigned int>(std::pair<double, double>(meshes[meshes_i]->nodes[i].x, meshes[meshes_i]->nodes[i].y), points_to_vertices.size()));
597 
598  vertices.push_back(vertex_BSON(meshes[meshes_i]->nodes[i].x, meshes[meshes_i]->nodes[i].y, new_i));
599  }
600  if (!hasAllElements)
601  subdomain.vertices.push_back(vertices_to_vertices.find(i)->second);
602  }
603 
604  // save elements
605  for (int i = 0; i < meshes[meshes_i]->get_num_base_elements(); i++)
606  {
607  e = &(meshes[meshes_i]->elements[i]);
608  if (!e->used)
609  continue;
610  if (!baseElementsSaved[e->id])
611  {
612  if (e->nvert == 3)
613  elements.push_back(element_BSON(vertices_to_vertices.find(e->vn[0]->id)->second, vertices_to_vertices.find(e->vn[1]->id)->second, vertices_to_vertices.find(e->vn[2]->id)->second, -1, meshes[meshes_i]->get_element_markers_conversion().get_user_marker(e->marker).marker, e->id));
614  else
615  elements.push_back(element_BSON(vertices_to_vertices.find(e->vn[0]->id)->second, vertices_to_vertices.find(e->vn[1]->id)->second, vertices_to_vertices.find(e->vn[2]->id)->second, vertices_to_vertices.find(e->vn[3]->id)->second, meshes[meshes_i]->get_element_markers_conversion().get_user_marker(e->marker).marker, e->id));
616  baseElementsSaved[e->id] = true;
617  }
618 
619  if (!hasAllElements)
620  subdomain.elements.push_back(e->id);
621  }
622 
623  // save boundary edge markers
624  bool has_inner_edges = false;
625  for_all_base_elements(e, meshes[meshes_i])
626  {
627  for (unsigned i = 0; i < e->get_nvert(); i++)
628  {
629  if (MeshUtil::get_base_edge_node(e, i)->bnd)
630  {
631  if (vertices_to_boundaries.find(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second))) == vertices_to_boundaries.end())
632  {
633  unsigned int edge_i = edges.size();
634  vertices_to_boundaries.insert(std::pair<std::pair<unsigned int, unsigned int>, unsigned int>(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second)), edge_i));
635  edges.push_back(edge_BSON(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second, meshes[meshes_i]->boundary_markers_conversion.get_user_marker(MeshUtil::get_base_edge_node(e, i)->marker).marker, edge_i));
636  }
637  if (!hasAllElements)
638  subdomain.boundary_edges.push_back(vertices_to_boundaries.find(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second)))->second);
639  }
640  else
641  has_inner_edges = true;
642  }
643  }
644 
645  if (has_inner_edges)
646  {
647  for_all_base_elements(e, meshes[meshes_i])
648  for (unsigned i = 0; i < e->get_nvert(); i++)
649  {
650  if (!MeshUtil::get_base_edge_node(e, i)->bnd)
651  {
652  if (vertices_to_boundaries.find(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second))) == vertices_to_boundaries.end())
653  {
654  unsigned int edge_i = edges.size();
655  vertices_to_boundaries.insert(std::pair<std::pair<unsigned int, unsigned int>, unsigned int>(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second)), edge_i));
656  edges.push_back(edge_BSON(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second, meshes[meshes_i]->boundary_markers_conversion.get_user_marker(MeshUtil::get_base_edge_node(e, i)->marker).marker, edge_i));
657  }
658  if (!hasAllElements)
659  subdomain.inner_edges.push_back(vertices_to_boundaries.find(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second)))->second);
660  }
661  }
662  }
663 
664  // save curved edges
665  for_all_base_elements(e, meshes[meshes_i])
666  {
667  if (e->is_curved())
668  for (unsigned i = 0; i < e->get_nvert(); i++)
669  if (e->cm->curves[i] != nullptr)
670  if (vertices_to_curves.find(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second))) == vertices_to_curves.end())
671  {
672  if (e->cm->curves[i]->type == ArcType)
673  arcs.push_back(arc_BSON(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second, ((Arc*)e->cm->curves[i])->angle));
674  else
675  throw Exceptions::Exception("BSON mesh loader can not operate with general NURBS so far.");
676 
677  vertices_to_curves.insert(std::pair<std::pair<unsigned int, unsigned int>, bool>(std::pair<unsigned int, unsigned int>(std::min(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second), std::max(vertices_to_vertices.find(e->vn[i]->id)->second, vertices_to_vertices.find(e->vn[e->next_vert(i)]->id)->second)), true));
678  }
679  }
680 
681  // save refinements
682  for (unsigned int refinement_i = 0; refinement_i < meshes[meshes_i]->refinements.size(); refinement_i++)
683  subdomain.refinements.push_back(refinement_BSON(meshes[meshes_i]->refinements[refinement_i].first, meshes[meshes_i]->refinements[refinement_i].second));
684 
685  subdomains.push_back(subdomain);
686  }
687 
688  delete[] baseElementsSaved;
689 
690  std::sort(elements.begin(), elements.end(), elementCompare);
691 #pragma endregion
692 
693  // Init bson
694  bson bw;
695  bson_init(&bw);
696 
697  // Save domain
698  // - vertices.
699  bson_append_start_array(&bw, "vertices");
700  for (int i = 0; i < vertices.size(); i++)
701  vertices[i].save_to_BSON(bw);
702  bson_append_finish_array(&bw);
703  // - elements.
704  bson_append_start_array(&bw, "elements");
705  for (int i = 0; i < elements.size(); i++)
706  elements[i].save_to_BSON(bw);
707  bson_append_finish_array(&bw);
708  // - edges
709  bson_append_start_array(&bw, "edges");
710  for (int i = 0; i < edges.size(); i++)
711  edges[i].save_to_BSON(bw);
712  bson_append_finish_array(&bw);
713  // - curves
714  bson_append_start_array(&bw, "arcs");
715  for (int i = 0; i < arcs.size(); i++)
716  arcs[i].save_to_BSON(bw);
717  bson_append_finish_array(&bw);
718 
719  // Save subdomains
720  bson_append_start_array(&bw, "subdomains");
721  for (int i = 0; i < subdomains.size(); i++)
722  subdomains[i].save_to_BSON(bw);
723  bson_append_finish_array(&bw);
724 
725  // Done.
726  bson_finish(&bw);
727 
728  // Write to disk.
729  FILE *fpw;
730  fpw = fopen(filename, "wb");
731  const char *dataw = (const char *)bson_data(&bw);
732  fwrite(dataw, bson_size(&bw), 1, fpw);
733  fclose(fpw);
734 
735  bson_destroy(&bw);
736  }
737 
738  void MeshReaderH2DBSON::load(const char *filename, std::vector<MeshSharedPtr> meshes)
739  {
740  FILE *fpr;
741  fpr = fopen(filename, "rb");
742 
743  // file size:
744  fseek(fpr, 0, SEEK_END);
745  int size = ftell(fpr);
746  rewind(fpr);
747 
748  // allocate memory to contain the whole file:
749  char *datar = malloc_with_check<char>(size);
750  fread(datar, size, 1, fpr);
751  fclose(fpr);
752 
753  bson br;
754  bson_init_finished_data(&br, datar, 0);
755 
756  for (unsigned int meshes_i = 0; meshes_i < meshes.size(); meshes_i++)
757  meshes.at(meshes_i)->free();
758 
759  MeshSharedPtr global_mesh(new Mesh);
760 
761  std::map<int, int> vertex_is;
762  std::map<int, int> element_is;
763  std::map<int, int> edge_is;
764 
765  std::vector<element_BSON> elements;
766  std::vector<edge_BSON> edges;
767  std::vector<vertex_BSON> vertices;
768  std::vector<arc_BSON> arcs;
769  std::vector<subdomain_BSON> subdomains;
770 
771  // load
772  this->load_domain(br, global_mesh, vertex_is, element_is, edge_is, elements, edges, vertices, arcs, subdomains);
773 
774  int max_vertex_i = -1;
775  for (std::map<int, int>::iterator it = vertex_is.begin(); it != vertex_is.end(); ++it)
776  if (it->first > max_vertex_i)
777  max_vertex_i = it->first;
778  int max_element_i = -1;
779  for (std::map<int, int>::iterator it = element_is.begin(); it != element_is.end(); ++it)
780  if (it->first > max_element_i)
781  max_element_i = it->first;
782  int max_edge_i = -1;
783  for (std::map<int, int>::iterator it = edge_is.begin(); it != edge_is.end(); ++it)
784  if (it->first > max_edge_i)
785  max_edge_i = it->first;
786 
787  // Subdomains //
788  unsigned int subdomains_count = subdomains.size();
789  if (subdomains_count != meshes.size())
790  throw Hermes::Exceptions::MeshLoadFailureException("Number of subdomains( = %u) does not equal the number of provided meshes in the vector( = %u).", subdomains_count, meshes.size());
791 
792  for (unsigned int subdomains_i = 0; subdomains_i < subdomains_count; subdomains_i++)
793  {
794  for (int element_i = 0; element_i < elements.size(); element_i++)
795  {
796  element_BSON element = elements.at(element_i);
797 
798  meshes[subdomains_i]->element_markers_conversion.insert_marker(element.marker);
799  }
800  for (unsigned int edge_i = 0; edge_i < edges.size(); edge_i++)
801  {
802  edge_BSON edge = edges.at(edge_i);
803 
804  meshes[subdomains_i]->boundary_markers_conversion.insert_marker(edge.marker);
805  }
806  }
807 
808  for (unsigned int subdomains_i = 0; subdomains_i < subdomains_count; subdomains_i++)
809  {
810  unsigned int vertex_number_count = subdomains.at(subdomains_i).vertices.empty() ? 0 : subdomains.at(subdomains_i).vertices.size();
811  unsigned int element_number_count = subdomains.at(subdomains_i).elements.empty() ? 0 : subdomains.at(subdomains_i).elements.size();
812  unsigned int boundary_edge_number_count = subdomains.at(subdomains_i).boundary_edges.empty() ? 0 : subdomains.at(subdomains_i).boundary_edges.size();
813  unsigned int inner_edge_number_count = subdomains.at(subdomains_i).inner_edges.empty() ? 0 : subdomains.at(subdomains_i).inner_edges.size();
814 
815  // copy the whole mesh if the subdomain is the whole mesh.
816  if (element_number_count == 0 || element_number_count == elements.size())
817  {
818  meshes[subdomains_i]->copy(global_mesh);
819  // refinements.
820  if (!subdomains.at(subdomains_i).refinements.empty() && subdomains.at(subdomains_i).refinements.size() > 0)
821  {
822  // perform initial refinements
823  for (unsigned int i = 0; i < subdomains.at(subdomains_i).refinements.size(); i++)
824  {
825  int element_id = subdomains.at(subdomains_i).refinements.at(i).id;
826  int refinement_type = subdomains.at(subdomains_i).refinements.at(i).type;
827  if (refinement_type == -1)
828  meshes[subdomains_i]->unrefine_element_id(element_id);
829  else
830  meshes[subdomains_i]->refine_element_id(element_id, refinement_type);
831  }
832  }
833  }
834  else
835  {
836  // Vertex numbers //
837  // create a mapping order-in-the-whole-domain <-> order-in-this-subdomain.
838  std::map<unsigned int, unsigned int> vertex_vertex_numbers;
839 
840  // Initialize mesh.
841  int size = HashTable::H2D_DEFAULT_HASH_SIZE;
842  while (size < 8 * vertex_number_count)
843  size *= 2;
844  meshes[subdomains_i]->init(size);
845 
846  // Create top-level vertex nodes.
847  if (vertex_number_count == 0)
848  vertex_number_count = vertices.size();
849  for (unsigned int vertex_numbers_i = 0; vertex_numbers_i < vertex_number_count; vertex_numbers_i++)
850  {
851  unsigned int vertex_number;
852  if (vertex_number_count == vertices.size())
853  vertex_number = vertex_is[vertex_numbers_i];
854  else
855  {
856  vertex_number = subdomains.at(subdomains_i).vertices.at(vertex_numbers_i);
857  if (vertex_number > max_vertex_i)
858  throw Exceptions::MeshLoadFailureException("Wrong vertex number:%u in subdomain %u.", vertex_number, subdomains_i);
859  }
860 
861  vertex_vertex_numbers.insert(std::pair<unsigned int, unsigned int>(vertex_number, vertex_numbers_i));
862  Node* node = meshes[subdomains_i]->nodes.add();
863  assert(node->id == vertex_numbers_i);
864  node->ref = TOP_LEVEL_REF;
865  node->type = HERMES_TYPE_VERTEX;
866  node->bnd = 0;
867  node->p1 = node->p2 = -1;
868  node->next_hash = nullptr;
869 
870  // assignment.
871  node->x = vertices[vertex_number].x;
872  node->y = vertices[vertex_number].y;
873  }
874  meshes[subdomains_i]->ntopvert = vertex_number_count;
875 
876  // Element numbers //
877  unsigned int element_count = elements.size();
878  meshes[subdomains_i]->nbase = element_count;
879  meshes[subdomains_i]->nactive = meshes[subdomains_i]->ninitial = element_number_count;
880 
881  Element* e;
882  int* elements_existing = new int[element_count];
883  for (int i = 0; i < element_count; i++)
884  elements_existing[i] = -1;
885  for (int element_number_i = 0; element_number_i < element_number_count; element_number_i++)
886  {
887  int elementI = subdomains.at(subdomains_i).elements.at(element_number_i);
888  if (elementI > max_element_i)
889  throw Exceptions::MeshLoadFailureException("Wrong element number:%i in subdomain %u.", elementI, subdomains_i);
890 
891  elements_existing[element_is[subdomains.at(subdomains_i).elements.at(element_number_i)]] = elementI;
892  }
893  for (int element_i = 0; element_i < element_count; element_i++)
894  {
895  bool found = false;
896  if (element_number_count == 0)
897  found = true;
898  else
899  found = elements_existing[element_i] != -1;
900 
901  if (!found)
902  {
903  meshes[subdomains_i]->elements.skip_slot()->cm = nullptr;
904  continue;
905  }
906 
907  element_BSON element(-1, -1, -1, -1, "", -1);
908  for (int searched_element_i = 0; searched_element_i < element_count; searched_element_i++)
909  {
910  element = elements.at(searched_element_i);
911  if (element.i == elements_existing[element_i])
912  break;
913  }
914  if (element.i == -1)
915  throw Exceptions::MeshLoadFailureException("Element number wrong in the mesh file.");
916 
917  if (element.v4 != -1)
918  e = meshes[subdomains_i]->create_quad(meshes[subdomains_i]->element_markers_conversion.get_internal_marker(element.marker).marker,
919  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v1)->second],
920  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v2)->second],
921  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v3)->second],
922  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v4)->second],
923  nullptr, element_i);
924  else
925  e = meshes[subdomains_i]->create_triangle(meshes[subdomains_i]->element_markers_conversion.get_internal_marker(element.marker).marker,
926  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v1)->second],
927  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v2)->second],
928  &meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(element.v3)->second],
929  nullptr, element_i);
930  }
931 
932  // Boundary Edge numbers //
933  if (boundary_edge_number_count == 0)
934  boundary_edge_number_count = edges.size();
935 
936  for (int boundary_edge_number_i = 0; boundary_edge_number_i < boundary_edge_number_count; boundary_edge_number_i++)
937  {
938  edge_BSON edge(-1, -1, "", -1);
939  int domain_edge_count = edges.size();
940  for (unsigned int to_find_i = 0; to_find_i < domain_edge_count; to_find_i++)
941  {
942  if (boundary_edge_number_count != domain_edge_count)
943  {
944  if (edges.at(to_find_i).i == subdomains.at(subdomains_i).boundary_edges.at(boundary_edge_number_i))
945  {
946  edge = edges.at(to_find_i);
947  break;
948  }
949  }
950  else
951  {
952  if (edges.at(to_find_i).i == edge_is[boundary_edge_number_i])
953  {
954  edge = edges.at(to_find_i);
955  break;
956  }
957  }
958  }
959 
960  if (edge.i == -1)
961  throw Exceptions::MeshLoadFailureException("Wrong boundary-edge number:%i in subdomain %u.", subdomains.at(subdomains_i).boundary_edges.at(boundary_edge_number_i), subdomains_i);
962 
963  Node* en = meshes[subdomains_i]->peek_edge_node(vertex_vertex_numbers.find(edge.v1)->second, vertex_vertex_numbers.find(edge.v2)->second);
964  if (en == nullptr)
965  throw Hermes::Exceptions::MeshLoadFailureException("Boundary data error (edge %i does not exist).", boundary_edge_number_i);
966 
967  en->marker = meshes[subdomains_i]->boundary_markers_conversion.get_internal_marker(edge.marker).marker;
968 
969  meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(edge.v1)->second].bnd = 1;
970  meshes[subdomains_i]->nodes[vertex_vertex_numbers.find(edge.v2)->second].bnd = 1;
971  en->bnd = 1;
972  }
973 
974  // Inner Edge numbers //
975  for (int inner_edge_number_i = 0; inner_edge_number_i < inner_edge_number_count; inner_edge_number_i++)
976  {
977  edge_BSON edge(-1, -1, "", -1);
978 
979  for (unsigned int to_find_i = 0; to_find_i < edges.size(); to_find_i++)
980  {
981  if (edges.at(to_find_i).i == subdomains.at(subdomains_i).inner_edges.at(inner_edge_number_i))
982  {
983  edge = edges.at(to_find_i);
984  break;
985  }
986  }
987 
988  if (edge.i == -1)
989  throw Exceptions::MeshLoadFailureException("Wrong inner-edge number:%i in subdomain %u.", subdomains.at(subdomains_i).boundary_edges.at(inner_edge_number_i), subdomains_i);
990 
991  Node* en = meshes[subdomains_i]->peek_edge_node(vertex_vertex_numbers.find(edge.v1)->second, vertex_vertex_numbers.find(edge.v2)->second);
992  if (en == nullptr)
993  throw Hermes::Exceptions::MeshLoadFailureException("Inner data error (edge %i does not exist).", inner_edge_number_i);
994 
995  en->marker = meshes[subdomains_i]->boundary_markers_conversion.get_internal_marker(edge.marker).marker;
996  en->bnd = 0;
997  }
998 
999  // Curves //
1000  // Arcs & NURBSs //
1001  unsigned int arc_count = arcs.empty() ? 0 : arcs.size();
1002 
1003  for (unsigned int curves_i = 0; curves_i < arc_count; curves_i++)
1004  {
1005  // load the control points, knot vector, etc.
1006  Node* en;
1007  int p1, p2;
1008 
1009  // first do arcs, then NURBSs.
1010  Curve* curve;
1011  if (vertex_vertex_numbers.find(arcs.at(curves_i).p1) == vertex_vertex_numbers.end() ||
1012  vertex_vertex_numbers.find(arcs.at(curves_i).p2) == vertex_vertex_numbers.end())
1013  continue;
1014  else
1015  {
1016  // read the end point indices
1017  p1 = vertex_vertex_numbers.find(arcs.at(curves_i).p1)->second;
1018  p2 = vertex_vertex_numbers.find(arcs.at(curves_i).p2)->second;
1019 
1020  curve = MeshUtil::load_arc(meshes[subdomains_i], curves_i, &en, p1, p2, arcs.at(curves_i).angle, true);
1021  if (curve == nullptr)
1022  continue;
1023  }
1024 
1025  // assign the arc to the elements sharing the edge node
1026  MeshUtil::assign_curve(en, curve, p1, p2);
1027  }
1028 
1029  // update refmap coeffs of curvilinear elements
1030  for_all_used_elements(e, meshes[subdomains_i])
1031  {
1032  if (e->cm != nullptr)
1033  e->cm->update_refmap_coeffs(e);
1034  RefMap::set_element_iro_cache(e);
1035  }
1036 
1037  // refinements.
1038  if (!subdomains.at(subdomains_i).refinements.empty() && subdomains.at(subdomains_i).refinements.size() > 0)
1039  {
1040  // perform initial refinements
1041  for (unsigned int i = 0; i < subdomains.at(subdomains_i).refinements.size(); i++)
1042  {
1043  int element_id = subdomains.at(subdomains_i).refinements.at(i).id;
1044  int refinement_type = subdomains.at(subdomains_i).refinements.at(i).type;
1045  if (refinement_type == -1)
1046  meshes[subdomains_i]->unrefine_element_id(element_id);
1047  else
1048  meshes[subdomains_i]->refine_element_id(element_id, refinement_type);
1049  }
1050  }
1051 
1052  delete[] elements_existing;
1053  }
1054  meshes[subdomains_i]->seq = g_mesh_seq++;
1055  if (HermesCommonApi.get_integral_param_value(checkMeshesOnLoad))
1056  meshes[subdomains_i]->initial_single_check();
1057  }
1058  bson_destroy(&br);
1059  free_with_check(datar);
1060  }
1061 
1062  void MeshReaderH2DBSON::load_domain(bson& br, MeshSharedPtr mesh, std::map<int, int>& vertex_is, std::map<int, int>& element_is, std::map<int, int>& edge_is,
1063  std::vector<element_BSON>& elements, std::vector<edge_BSON>& edges, std::vector<vertex_BSON>& vertices, std::vector<arc_BSON>& arcs, std::vector<subdomain_BSON>& subdomains)
1064  {
1065  bson_iterator it, it_sub;
1066  bson b_sub, b_sub_sub;
1067  bson_find(&it, &br, "vertices");
1068  bson_iterator_subobject_init(&it, &b_sub, 0);
1069  bson_iterator_init(&it_sub, &b_sub);
1070  while (bson_iterator_next(&it_sub))
1071  {
1072  vertex_BSON vertex;
1073  bson_iterator_subobject_init(&it_sub, &b_sub_sub, 0);
1074  vertex.load_from_BSON(b_sub_sub);
1075  vertices.push_back(vertex);
1076  }
1077 
1078  bson_find(&it, &br, "elements");
1079  bson_iterator_subobject_init(&it, &b_sub, 0);
1080  bson_iterator_init(&it_sub, &b_sub);
1081  while (bson_iterator_next(&it_sub))
1082  {
1083  element_BSON element;
1084  bson_iterator_subobject_init(&it_sub, &b_sub_sub, 0);
1085  element.load_from_BSON(b_sub_sub);
1086  elements.push_back(element);
1087  }
1088 
1089  bson_find(&it, &br, "edges");
1090  bson_iterator_subobject_init(&it, &b_sub, 0);
1091  bson_iterator_init(&it_sub, &b_sub);
1092  while (bson_iterator_next(&it_sub))
1093  {
1094  edge_BSON edge;
1095  bson_iterator_subobject_init(&it_sub, &b_sub_sub, 0);
1096  edge.load_from_BSON(b_sub_sub);
1097  edges.push_back(edge);
1098  }
1099 
1100  bson_find(&it, &br, "arcs");
1101  bson_iterator_subobject_init(&it, &b_sub, 0);
1102  bson_iterator_init(&it_sub, &b_sub);
1103  while (bson_iterator_next(&it_sub))
1104  {
1105  arc_BSON arc;
1106  bson_iterator_subobject_init(&it_sub, &b_sub_sub, 0);
1107  arc.load_from_BSON(b_sub_sub);
1108  arcs.push_back(arc);
1109  }
1110 
1111  bson_find(&it, &br, "subdomains");
1112  bson_iterator_subobject_init(&it, &b_sub, 0);
1113  bson_iterator_init(&it_sub, &b_sub);
1114  while (bson_iterator_next(&it_sub))
1115  {
1116  subdomain_BSON subdomain;
1117  bson_iterator_subobject_init(&it_sub, &b_sub_sub, 0);
1118  subdomain.load_from_BSON(b_sub_sub);
1119  subdomains.push_back(subdomain);
1120  }
1121 
1122  // Vertices //
1123  int vertices_count = vertices.size();
1124 
1125  // Initialize mesh.
1126  int size = HashTable::H2D_DEFAULT_HASH_SIZE;
1127  while (size < 8 * vertices_count)
1128  size *= 2;
1129  mesh->init(size);
1130 
1131  // Create top-level vertex nodes.
1132  for (int vertex_i = 0; vertex_i < vertices_count; vertex_i++)
1133  {
1134  Node* node = mesh->nodes.add();
1135  assert(node->id == vertex_i);
1136  node->ref = TOP_LEVEL_REF;
1137  node->type = HERMES_TYPE_VERTEX;
1138  node->bnd = 0;
1139  node->p1 = node->p2 = -1;
1140  node->next_hash = nullptr;
1141 
1142  if (vertices[vertex_i].i > H2D_MAX_NODE_ID - 1)
1143  throw Exceptions::MeshLoadFailureException("The index 'i' of vertex in the mesh file must be lower than %i.", H2D_MAX_NODE_ID);
1144 
1145  // insert
1146  vertex_is.insert(std::pair<int, int>(vertices[vertex_i].i, vertex_i));
1147 
1148  // assignment.
1149  node->x = vertices[vertex_i].x;
1150  node->y = vertices[vertex_i].y;
1151  }
1152  mesh->ntopvert = vertices_count;
1153 
1154  // Elements //
1155  unsigned int element_count = elements.size();
1156  mesh->nbase = mesh->nactive = mesh->ninitial = element_count;
1157 
1158  Element* e;
1159  for (int element_i = 0; element_i < element_count; element_i++)
1160  {
1161  element_BSON element = elements[element_i];
1162 
1163  // insert.
1164  if (element.i > H2D_MAX_NODE_ID - 1)
1165  throw Exceptions::MeshLoadFailureException("The index 'i' of element in the mesh file must be lower than %i.", H2D_MAX_NODE_ID);
1166 
1167  element_is.insert(std::pair<int, int>(element.i, element_i));
1168 
1169  mesh->element_markers_conversion.insert_marker(element.marker);
1170 
1171  if (element.v4 != -1)
1172  e = mesh->create_quad(mesh->element_markers_conversion.get_internal_marker(element.marker).marker,
1173  &mesh->nodes[element.v1],
1174  &mesh->nodes[element.v2],
1175  &mesh->nodes[element.v3],
1176  &mesh->nodes[element.v4],
1177  nullptr);
1178  else
1179  e = mesh->create_triangle(mesh->element_markers_conversion.get_internal_marker(element.marker).marker,
1180  &mesh->nodes[element.v1],
1181  &mesh->nodes[element.v2],
1182  &mesh->nodes[element.v3],
1183  nullptr);
1184  }
1185 
1186  // Boundaries //
1187  unsigned int edges_count = edges.size();
1188 
1189  Node* en;
1190  for (unsigned int edge_i = 0; edge_i < edges_count; edge_i++)
1191  {
1192  int v1 = edges[edge_i].v1;
1193  int v2 = edges[edge_i].v2;
1194 
1195  // insert
1196  if (edges[edge_i].i > H2D_MAX_NODE_ID - 1)
1197  throw Exceptions::MeshLoadFailureException("The index 'i' of edge in the mesh file must be lower than %i.", H2D_MAX_NODE_ID);
1198 
1199  edge_is.insert(std::pair<int, int>(edge_i, edges[edge_i].i));
1200 
1201  en = mesh->peek_edge_node(v1, v2);
1202  if (en == nullptr)
1203  throw Hermes::Exceptions::MeshLoadFailureException("Boundary data #%d: edge %d-%d does not exist.", edge_i, v1, v2);
1204 
1205  std::string edge_marker = edges[edge_i].marker;
1206 
1207  // Trim whitespaces.
1208  unsigned int begin = edge_marker.find_first_not_of(" \t\n");
1209  unsigned int end = edge_marker.find_last_not_of(" \t\n");
1210  edge_marker.erase(end + 1, edge_marker.length());
1211  edge_marker.erase(0, begin);
1212 
1213  // This functions check if the user-supplied marker on this element has been
1214  // already used, and if not, inserts it in the appropriate structure.
1215  mesh->boundary_markers_conversion.insert_marker(edge_marker);
1216  int marker = mesh->boundary_markers_conversion.get_internal_marker(edge_marker).marker;
1217 
1218  en->marker = marker;
1219  }
1220 
1221  Node* node;
1222  for_all_edge_nodes(node, mesh)
1223  {
1224  if (node->ref < 2)
1225  {
1226  mesh->nodes[node->p1].bnd = 1;
1227  mesh->nodes[node->p2].bnd = 1;
1228  node->bnd = 1;
1229  }
1230  }
1231 
1232  // check that all boundary edges have a marker assigned
1233  for_all_edge_nodes(en, mesh)
1234  if (en->ref < 2 && en->marker == 0)
1235  this->warn("Boundary edge node does not have a boundary marker.");
1236 
1237  // Curves //
1238  // Arcs & NURBSs //
1239  unsigned int arc_count = arcs.size();
1240 
1241  for (unsigned int curves_i = 0; curves_i < arc_count; curves_i++)
1242  {
1243  // load the control points, knot vector, etc.
1244  Node* en;
1245  int p1, p2;
1246 
1247  // first do arcs, then NURBSs.
1248  Curve* curve;
1249  // read the end point indices
1250  p1 = arcs[curves_i].p1;
1251  p2 = arcs[curves_i].p2;
1252 
1253  curve = MeshUtil::load_arc(mesh, curves_i, &en, p1, p2, arcs[curves_i].angle);
1254 
1255  // assign the arc to the elements sharing the edge node
1256  MeshUtil::assign_curve(en, curve, p1, p2);
1257  }
1258 
1259  // update refmap coeffs of curvilinear elements
1260  for_all_used_elements(e, mesh)
1261  {
1262  if (e->cm != nullptr)
1263  e->cm->update_refmap_coeffs(e);
1264  RefMap::set_element_iro_cache(e);
1265  }
1266  }
1267  }
1268 }
1269 
1270 #endif
Definition: adapt.h:24
::xsd::cxx::tree::type type
C++ type corresponding to the anyType XML Schema built-in type.
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.