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