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