Hermes2D  2.0
space.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.org/licenses/>.
15 
16 #include "space_h1.h"
17 #include "space_l2.h"
18 #include "shapeset_hc_all.h"
19 #include "shapeset_hd_all.h"
20 #include "shapeset_h1_all.h"
21 #include "shapeset_l2_all.h"
22 #include "space_hcurl.h"
23 #include "space_hdiv.h"
24 #include "space_h2d_xml.h"
25 #include "api2d.h"
26 #include <iostream>
27 
28 namespace Hermes
29 {
30  namespace Hermes2D
31  {
32  unsigned g_space_seq = 0;
33 
34  template<>
35  void Space<double>::init()
36  {
37  this->default_tri_order = -1;
38  this->default_quad_order = -1;
39  this->ndata = NULL;
40  this->edata = NULL;
41  this->nsize = esize = 0;
42  this->ndata_allocated = 0;
43  this->mesh_seq = -1;
44  this->seq = g_space_seq;
45  this->was_assigned = -1;
46  this->ndof = 0;
47  this->proj_mat = NULL;
48  this->chol_p = NULL;
49  this->vertex_functions_count = this->edge_functions_count = this->bubble_functions_count = 0;
50 
51  if(essential_bcs != NULL)
52  for(Hermes::vector<EssentialBoundaryCondition<double>*>::const_iterator it = essential_bcs->begin(); it != essential_bcs->end(); it++)
53  for(unsigned int i = 0; i < (*it)->markers.size(); i++)
54  {
55  if(mesh->boundary_markers_conversion.conversion_table_inverse.find((*it)->markers.at(i)) == mesh->boundary_markers_conversion.conversion_table_inverse.end() && (*it)->markers.at(i) != HERMES_ANY)
56  throw Hermes::Exceptions::Exception("A boundary condition defined on a non-existent marker %s.", (*it)->markers.at(i).c_str());
57  }
58 
59  own_shapeset = (shapeset == NULL);
60  }
61 
62  template<>
63  void Space<std::complex<double> >::init()
64  {
65  this->default_tri_order = -1;
66  this->default_quad_order = -1;
67  this->ndata = NULL;
68  this->edata = NULL;
69  this->nsize = esize = 0;
70  this->ndata_allocated = 0;
71  this->mesh_seq = -1;
72  this->seq = g_space_seq;
73  this->was_assigned = -1;
74  this->ndof = 0;
75  this->proj_mat = NULL;
76  this->chol_p = NULL;
77  this->vertex_functions_count = this->edge_functions_count = this->bubble_functions_count = 0;
78 
79  if(essential_bcs != NULL)
80  for(Hermes::vector<EssentialBoundaryCondition<std::complex<double> >*>::const_iterator it = essential_bcs->begin(); it != essential_bcs->end(); it++)
81  for(unsigned int i = 0; i < (*it)->markers.size(); i++)
82  {
83  if(mesh->boundary_markers_conversion.conversion_table_inverse.find((*it)->markers.at(i)) == mesh->boundary_markers_conversion.conversion_table_inverse.end() && (*it)->markers.at(i) != HERMES_ANY)
84  throw Hermes::Exceptions::Exception("A boundary condition defined on a non-existent marker %s.", (*it)->markers.at(i).c_str());
85  }
86 
87  own_shapeset = (shapeset == NULL);
88  }
89 
90  template<>
91  void Space<double>::free()
92  {
93  free_bc_data();
94  if(nsize) { ::free(ndata); nsize = 0; ndata = NULL; }
95  if(esize) { ::free(edata); edata = 0; edata = NULL; }
96  this->seq = -1;
97  }
98 
99  template<>
100  void Space<std::complex<double> >::free()
101  {
102  free_bc_data();
103  if(nsize) { ::free(ndata); nsize = 0; ndata = NULL; }
104  if(esize) { ::free(edata); edata = 0; edata = NULL; }
105  this->seq = -1;
106  }
107 
108  template<>
109  Space<double>::Space() : shapeset(NULL), essential_bcs(NULL), mesh(NULL)
110  {
111  this->init();
112  }
113 
114  template<>
115  Space<std::complex<double> >::Space() : shapeset(NULL), essential_bcs(NULL), mesh(NULL)
116  {
117  this->init();
118  }
119 
120  template<>
121  Space<double>::Space(const Mesh* mesh, Shapeset* shapeset, EssentialBCs<double>* essential_bcs)
122  : shapeset(shapeset), essential_bcs(essential_bcs), mesh(mesh)
123  {
124  if(mesh == NULL)
125  throw Hermes::Exceptions::NullException(0);
126  this->init();
127  }
128 
129  template<>
130  Space<std::complex<double> >::Space(const Mesh* mesh, Shapeset* shapeset, EssentialBCs<std::complex<double> >* essential_bcs)
131  : shapeset(shapeset), essential_bcs(essential_bcs), mesh(mesh)
132  {
133  if(mesh == NULL)
134  throw Hermes::Exceptions::NullException(0);
135  this->init();
136  }
137 
138  template<typename Scalar>
139  bool Space<Scalar>::isOkay() const
140  {
141  if(ndata == NULL || edata == NULL || !nsize || !esize)
142  return false;
143  if(seq < 0)
144  return false;
145  if(this->mesh == NULL)
146  return false;
147 
148  this->mesh->check();
149 
150  if(edata == NULL)
151  {
152  throw Hermes::Exceptions::Exception("NULL edata detected in Space<Scalar>::get_element_order().");
153  return false;
154  }
155 
156  if(!this->is_up_to_date())
157  {
158  throw Hermes::Exceptions::Exception("Space is not up to date.");
159  return false;
160  }
161  else
162  return true;
163  }
164 
165  template<>
166  Space<double>::~Space()
167  {
168  free();
169 
170  if(this->proj_mat != NULL)
171  delete [] this->proj_mat;
172  if(this->chol_p != NULL)
173  delete [] this->chol_p;
174 
175  }
176 
177  template<>
178  Space<std::complex<double> >::~Space()
179  {
180  free();
181 
182  if(this->proj_mat != NULL)
183  delete [] this->proj_mat;
184  if(this->chol_p != NULL)
185  delete [] this->chol_p;
186 
187  }
188 
189  template<typename Scalar>
190  Node* Space<Scalar>::get_mid_edge_vertex_node(Element* e, int i, int j)
191  {
192  if(e->is_triangle())
193  return e->sons[3]->vn[e->prev_vert(i)];
194  else if(e->sons[2] == NULL)
195  return i == 1 ? e->sons[0]->vn[2] : i == 3 ? e->sons[0]->vn[3] : NULL;
196  else if(e->sons[0] == NULL)
197  return i == 0 ? e->sons[2]->vn[1] : i == 2 ? e->sons[2]->vn[2] : NULL;
198  else return e->sons[i]->vn[j];
199  }
200 
201  template<typename Scalar>
202  void Space<Scalar>::resize_tables()
203  {
204  if((nsize < mesh->get_max_node_id()) || (ndata == NULL))
205  {
206  //HACK: definition of allocated size and the result number of elements
207  nsize = mesh->get_max_node_id();
208  if((nsize > ndata_allocated) || (ndata == NULL))
209  {
210  int prev_allocated = ndata_allocated;
211  if(ndata_allocated == 0)
212  ndata_allocated = 1024;
213  while (ndata_allocated < nsize)
214  ndata_allocated = ndata_allocated * 3 / 2;
215  ndata = (NodeData*)realloc(ndata, ndata_allocated * sizeof(NodeData));
216  for(int i = prev_allocated; i < ndata_allocated; i++)
217  ndata[i].edge_bc_proj = NULL;
218  }
219  }
220 
221  if((esize < mesh->get_max_element_id()) || (edata == NULL))
222  {
223  int oldsize = esize;
224  if(!esize)
225  esize = 1024;
226  while (esize < mesh->get_max_element_id())
227  esize = esize * 3 / 2;
228  edata = (ElementData*) realloc(edata, sizeof(ElementData) * esize);
229  for (int i = oldsize; i < esize; i++)
230  edata[i].order = -1;
231  for (int i = oldsize; i < esize; i++)
232  edata[i].changed_in_last_adaptation = true;
233  }
234  }
235 
236  template<typename Scalar>
237  void Space<Scalar>::copy(const Space<Scalar>* space, Mesh* new_mesh)
238  {
239  this->free();
240 
241  this->vertex_functions_count = this->edge_functions_count = this->bubble_functions_count = 0;
242 
243  this->essential_bcs = space->essential_bcs;
244  this->shapeset = space->shapeset->clone();
245 
246  new_mesh->copy(space->get_mesh());
247  this->mesh = new_mesh;
248 
249  this->resize_tables();
250 
251  Element* e;
252  for_all_active_elements(e, this->mesh)
253  {
254  this->set_element_order_internal(e->id, space->get_element_order(e->id));
255  }
256 
257  this->seq = g_space_seq++;
258 
259  for_all_active_elements(e, this->mesh)
260  {
261  if(space->edata[e->id].changed_in_last_adaptation)
262  this->edata[e->id].changed_in_last_adaptation = true;
263  else
264  this->edata[e->id].changed_in_last_adaptation = false;
265  }
266  }
267 
268  template<typename Scalar>
269  Shapeset* Space<Scalar>::get_shapeset() const
270  {
271  return this->shapeset;
272  }
273 
274  template<typename Scalar>
275  int Space<Scalar>::get_num_dofs() const
276  {
277  check();
278  return ndof;
279  }
280 
281  template<typename Scalar>
282  int Space<Scalar>::get_max_dof() const
283  {
284  check();
285  return next_dof - stride;
286  }
287 
288  template<typename Scalar>
289  Mesh* Space<Scalar>::get_mesh() const
290  {
291  return const_cast<Mesh*>(mesh);
292  }
293 
294  template<typename Scalar>
295  bool Space<Scalar>::is_up_to_date() const
296  {
297  return was_assigned == this->seq && mesh_seq == (int) mesh->get_seq();
298  }
299 
300  template<typename Scalar>
301  EssentialBCs<Scalar>* Space<Scalar>::get_essential_bcs() const
302  {
303  check();
304  return essential_bcs;
305  }
306 
307  template<typename Scalar>
308  void Space<Scalar>::set_element_order(int id, int order)
309  {
310  set_element_order_internal(id, order);
311  }
312 
313  template<typename Scalar>
314  void Space<Scalar>::set_element_order_internal(int id, int order)
315  {
316  if(id < 0 || id >= mesh->get_max_element_id())
317  throw Hermes::Exceptions::Exception("Space<Scalar>::set_element_order_internal: Invalid element id.");
318 
319  resize_tables();
320 
321  if(mesh->get_element(id)->is_quad() && get_type() != HERMES_L2_SPACE && H2D_GET_V_ORDER(order) == 0)
322  order = H2D_MAKE_QUAD_ORDER(order, order);
323 
324  edata[id].order = order;
325  seq = g_space_seq++;
326  }
327 
328  template<typename Scalar>
329  void Space<Scalar>::update_essential_bc_values(Hermes::vector<Space<Scalar>*> spaces, double time)
330  {
331  int n = spaces.size();
332  for (int i = 0; i < n; i++)
333  {
334  if(spaces[i]->get_essential_bcs() != NULL)
335  spaces[i]->get_essential_bcs()->set_current_time(time);
336  spaces[i]->update_essential_bc_values();
337  }
338  }
339 
340  template<typename Scalar>
341  void Space<Scalar>::update_essential_bc_values(Space<Scalar>*s, double time)
342  {
343  s->get_essential_bcs()->set_current_time(time);
344  s->update_essential_bc_values();
345  }
346 
347  template<typename Scalar>
348  int Space<Scalar>::get_num_dofs(Hermes::vector<const Space<Scalar>*> spaces)
349  {
350  int ndof = 0;
351  for (unsigned int i = 0; i<spaces.size(); i++)
352  ndof += spaces[i]->get_num_dofs();
353  return ndof;
354  }
355 
356  template<typename Scalar>
357  int Space<Scalar>::get_num_dofs(Hermes::vector<Space<Scalar>*> spaces)
358  {
359  int ndof = 0;
360  for (unsigned int i = 0; i<spaces.size(); i++)
361  ndof += spaces[i]->get_num_dofs();
362  return ndof;
363  }
364 
365  template<typename Scalar>
366  int Space<Scalar>::get_num_dofs(const Space<Scalar>* space)
367  {
368  return space->get_num_dofs();
369  }
370 
371  template<typename Scalar>
372  int Space<Scalar>::get_num_dofs(Space<Scalar>* space)
373  {
374  return space->get_num_dofs();
375  }
376 
377  template<typename Scalar>
378  int Space<Scalar>::assign_dofs(Hermes::vector<Space<Scalar>*> spaces)
379  {
380  int n = spaces.size();
381 
382  int ndof = 0;
383  for (int i = 0; i < n; i++) {
384  ndof += spaces[i]->assign_dofs(ndof);
385  }
386 
387  return ndof;
388  }
389 
390  template<typename Scalar>
391  int Space<Scalar>::get_element_order(int id) const
392  {
393  if(id >= esize)
394  {
395  this->warn("Element index %d in Space<Scalar>::get_element_order() while maximum is %d.", id, esize);
396  throw Hermes::Exceptions::Exception("Wrong element index in Space<Scalar>::get_element_order().");
397  }
398 
399  return edata[id].order;
400  }
401 
402  template<typename Scalar>
403  void Space<Scalar>::set_uniform_order(int order, std::string marker)
404  {
405  if(marker == HERMES_ANY)
406  set_uniform_order_internal(order, -1234);
407  else
408  set_uniform_order_internal(order, mesh->element_markers_conversion.get_internal_marker(marker).marker);
409  }
410 
411  template<typename Scalar>
412  void Space<Scalar>::set_uniform_order_internal(int order, int marker)
413  {
414  resize_tables();
415  int quad_order = H2D_MAKE_QUAD_ORDER(order, order);
416 
417  Element* e;
418  for_all_active_elements(e, mesh)
419  {
420  if(marker == HERMES_ANY_INT || e->marker == marker)
421  {
422  ElementData* ed = &edata[e->id];
423  if(e->is_triangle())
424  ed->order = order;
425  else
426  ed->order = quad_order;
427  }
428  }
429  seq = g_space_seq++;
430  }
431 
432  template<typename Scalar>
433  void Space<Scalar>::set_element_orders(int* elem_orders_)
434  {
435  resize_tables();
436 
437  Element* e;
438  int counter = 0;
439  for_all_elements(e, mesh)
440  {
441  assert(elem_orders_[counter] >= 0 && elem_orders_[counter] <= shapeset->get_max_order());
442  ElementData* ed = &edata[e->id];
443  if(e->is_triangle())
444  ed->order = elem_orders_[counter];
445  else
446  ed->order = H2D_MAKE_QUAD_ORDER(elem_orders_[counter], elem_orders_[counter]);
447  counter++;
448  }
449  }
450 
451  template<typename Scalar>
452  void Space<Scalar>::adjust_element_order(int order_change, int min_order)
453  {
454  Element* e;
455  for_all_active_elements(e, this->get_mesh())
456  {
457  if(e->is_triangle())
458  set_element_order_internal(e->id, std::max<int>(min_order, get_element_order(e->id) + order_change));
459  else
460  {
461  int h_order, v_order;
462 
463  if(H2D_GET_H_ORDER(get_element_order(e->id)) + order_change < min_order)
464  h_order = min_order;
465  else
466  h_order = H2D_GET_H_ORDER(get_element_order(e->id)) + order_change;
467 
468  if(H2D_GET_V_ORDER(get_element_order(e->id)) + order_change < min_order)
469  v_order = min_order;
470  else
471  v_order = H2D_GET_V_ORDER(get_element_order(e->id)) + order_change;
472 
473  set_element_order_internal(e->id, H2D_MAKE_QUAD_ORDER(h_order, v_order));
474  }
475  }
476  }
477 
478  template<typename Scalar>
479  void Space<Scalar>::adjust_element_order(int horizontal_order_change, int vertical_order_change, unsigned int horizontal_min_order, unsigned int vertical_min_order)
480  {
481  Element* e;
482  for_all_active_elements(e, this->get_mesh())
483  {
484  if(e->is_triangle())
485  {
486  this->warn("Using quad version of Space<Scalar>::adjust_element_order(), only horizontal orders will be used.");
487  set_element_order_internal(e->id, std::max<int>(horizontal_min_order, get_element_order(e->id) + horizontal_order_change));
488  }
489  else
490  if(get_element_order(e->id) == -1)
491  set_element_order_internal(e->id, H2D_MAKE_QUAD_ORDER(horizontal_min_order, vertical_min_order));
492  else
493  set_element_order_internal(e->id, H2D_MAKE_QUAD_ORDER(std::max<int>(H2D_GET_H_ORDER(get_element_order(e->id)) + horizontal_order_change, horizontal_min_order), std::max<int>(H2D_GET_V_ORDER(get_element_order(e->id)) + vertical_order_change, vertical_min_order)));
494  }
495  }
496 
497  template<typename Scalar>
498  void Space<Scalar>::update_element_orders_after_refinement()
499  {
500  Element* e;
501  for_all_active_elements(e, this->mesh)
502  {
503  if(this->get_element_order(e->id) < 0)
504  this->set_element_order_internal(e->id, this->get_element_order(e->parent->id));
505  }
506  }
507 
508  template<typename Scalar>
509  void Space<Scalar>::unrefine_all_mesh_elements(bool keep_initial_refinements)
510  {
511  check();
512  // find inactive elements with active sons
513  Hermes::vector<int> list;
514  Element* e;
515  for_all_inactive_elements(e, this->mesh)
516  {
517  bool found = true;
518  for (unsigned int i = 0; i < 4; i++)
519  if(e->sons[i] != NULL &&
520  (!e->sons[i]->active || (keep_initial_refinements && e->sons[i]->id < this->mesh->ninitial))
521  )
522  { found = false; break; }
523 
524  if(found) list.push_back(e->id);
525  }
526 
527  // unrefine the found elements
528  for (unsigned int i = 0; i < list.size(); i++)
529  {
530  unsigned int order = 0, h_order = 0, v_order = 0;
531  unsigned int num_sons = 0;
532  if(this->mesh->get_element_fast(list[i])->bsplit())
533  {
534  num_sons = 4;
535  for (int sons_i = 0; sons_i < 4; sons_i++)
536  {
537  if(this->mesh->get_element_fast(list[i])->sons[sons_i]->active)
538  {
539  if(this->mesh->get_element_fast(list[i])->sons[sons_i]->is_triangle())
540  order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[sons_i]->id);
541  else
542  {
543  h_order += H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[sons_i]->id));
544  v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[sons_i]->id));
545  }
546  }
547  }
548  }
549  else
550  {
551  if(this->mesh->get_element_fast(list[i])->hsplit())
552  {
553  num_sons = 2;
554  if(this->mesh->get_element_fast(list[i])->sons[0]->active)
555  {
556  if(this->mesh->get_element_fast(list[i])->sons[0]->is_triangle())
557  order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[0]->id);
558  else
559  {
560  h_order += H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[0]->id));
561  v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[0]->id));
562  }
563  }
564  if(this->mesh->get_element_fast(list[i])->sons[1]->active)
565  {
566  if(this->mesh->get_element_fast(list[i])->sons[1]->is_triangle())
567  order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[1]->id);
568  else
569  {
570  h_order += H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[1]->id));
571  v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[1]->id));
572  }
573  }
574  }
575  else
576  {
577  num_sons = 2;
578  if(this->mesh->get_element_fast(list[i])->sons[2]->active)
579  {
580  if(this->mesh->get_element_fast(list[i])->sons[2]->is_triangle())
581  order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[2]->id);
582  else
583  {
584  h_order += H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[2]->id));
585  v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[2]->id));
586  }
587  }
588  if(this->mesh->get_element_fast(list[i])->sons[3]->active)
589  {
590  if(this->mesh->get_element_fast(list[i])->sons[3]->is_triangle())
591  order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[3]->id);
592  else
593  {
594  h_order += H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[3]->id));
595  v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[3]->id));
596  }
597  }
598  }
599  }
600  order = (unsigned int)(order / num_sons);
601  h_order = (unsigned int)(h_order / num_sons);
602  v_order = (unsigned int)(v_order / num_sons);
603 
604  if(this->mesh->get_element_fast(list[i])->is_triangle())
605  edata[list[i]].order = order;
606  else
607  edata[list[i]].order = H2D_MAKE_QUAD_ORDER(h_order, v_order);
608  const_cast<Mesh*>(this->mesh)->unrefine_element_id(list[i]);
609  }
610 
611  // Recalculate all integrals, do not use previous adaptivity step.
612  for_all_active_elements(e, this->mesh)
613  this->edata[e->id].changed_in_last_adaptation = true;
614 
615  }
616 
617  template<typename Scalar>
618  Space<Scalar>::ReferenceSpaceCreator::ReferenceSpaceCreator(const Space<Scalar>* coarse_space, const Mesh* ref_mesh, unsigned int order_increase) : coarse_space(coarse_space), ref_mesh(ref_mesh), order_increase(order_increase)
619  {
620  }
621 
622  template<typename Scalar>
623  void Space<Scalar>::ReferenceSpaceCreator::handle_orders(Space<Scalar>* ref_space)
624  {
625  Element* e;
626  for_all_active_elements(e, coarse_space->get_mesh())
627  {
628  // This is the element id on the COARSE mesh. One can use it in the logic of increasing polynomial order (selectively).
629  int coarse_element_id = e->id;
630 
631  // Get the current order (may be triangular or quadrilateral - in the latter case, it is both horizontal and vertical order encoded in one number).
632  int current_order = coarse_space->get_element_order(coarse_element_id);
633  // Sanity check.
634  if(current_order < 0)
635  throw Hermes::Exceptions::Exception("Source space has an uninitialized order (element id = %d)", coarse_element_id);
636 
637  // New order calculation.
638  int new_order;
639  if(e->is_triangle())
640  {
641  // The triangular order is just the current order.
642  int current_order_triangular = current_order;
646  new_order = current_order_triangular + this->order_increase;
647  }
648  else
649  {
650  // We have to get the proper parts of the encoded order.
651  // - horizontal.
652  int current_order_quadrilateral_horizontal = H2D_GET_H_ORDER(current_order);
653  // - vertical.
654  int current_order_quadrilateral_vertical = H2D_GET_V_ORDER(current_order);
655 
656  // And now we have to create the new encoded order.
657  int new_order_quadrilateral_horizontal = current_order_quadrilateral_horizontal + this->order_increase;
658  int new_order_quadrilateral_vertical = current_order_quadrilateral_vertical + this->order_increase;
659  new_order = H2D_MAKE_QUAD_ORDER(new_order_quadrilateral_horizontal, new_order_quadrilateral_vertical);
660  }
661 
662  // And now call this method that does the magic for us (sets the new_order to .
663  ref_space->update_orders_recurrent(ref_space->mesh->get_element(coarse_element_id), new_order);
664  }
665  }
666 
667  template<typename Scalar>
668  Space<Scalar>* Space<Scalar>::ReferenceSpaceCreator::create_ref_space(bool assign_dofs)
669  {
671  Space<Scalar>* ref_space = NULL;
672  if(dynamic_cast<const L2Space<Scalar>*>(this->coarse_space) != NULL)
673  ref_space = this->init_construction_l2();
674  if(dynamic_cast<const H1Space<Scalar>*>(this->coarse_space) != NULL)
675  ref_space = this->init_construction_h1();
676  if(dynamic_cast<const HcurlSpace<Scalar>*>(this->coarse_space) != NULL)
677  ref_space = this->init_construction_hcurl();
678  if(dynamic_cast<const HdivSpace<Scalar>*>(this->coarse_space) != NULL)
679  ref_space = this->init_construction_hdiv();
680 
681  if(ref_space == NULL)
682  throw Exceptions::Exception("Something went wrong in ReferenceSpaceCreator::create_ref_space().");
683 
685  this->handle_orders(ref_space);
686 
688  this->finish_construction(ref_space);
689 
690  // Assign dofs?
691  if(assign_dofs)
692  ref_space->assign_dofs();
693 
695  return ref_space;
696  }
697 
698  template<typename Scalar>
699  L2Space<Scalar>* Space<Scalar>::ReferenceSpaceCreator::init_construction_l2()
700  {
701  L2Space<Scalar>* ref_space;
702  if(this->coarse_space->own_shapeset)
703  ref_space = new L2Space<Scalar>(this->ref_mesh, 0);
704  else
705  ref_space = new L2Space<Scalar>(this->ref_mesh, 0, this->coarse_space->get_shapeset());
706  return ref_space;
707  }
708 
709  template<typename Scalar>
710  H1Space<Scalar>* Space<Scalar>::ReferenceSpaceCreator::init_construction_h1()
711  {
712  H1Space<Scalar>* ref_space;
713  if(this->coarse_space->own_shapeset)
714  ref_space = new H1Space<Scalar>(this->ref_mesh, this->coarse_space->get_essential_bcs(), 1);
715  else
716  ref_space = new H1Space<Scalar>(this->ref_mesh, this->coarse_space->get_essential_bcs(), 1, this->coarse_space->get_shapeset());
717  return ref_space;
718  }
719 
720  template<typename Scalar>
721  HcurlSpace<Scalar>* Space<Scalar>::ReferenceSpaceCreator::init_construction_hcurl()
722  {
723  HcurlSpace<Scalar>* ref_space;
724  if(this->coarse_space->own_shapeset)
725  ref_space = new HcurlSpace<Scalar>(this->ref_mesh, this->coarse_space->essential_bcs, 1);
726  else
727  ref_space = new HcurlSpace<Scalar>(this->ref_mesh, this->coarse_space->essential_bcs, 1, this->coarse_space->shapeset);
728  return ref_space;
729  }
730 
731  template<typename Scalar>
732  HdivSpace<Scalar>* Space<Scalar>::ReferenceSpaceCreator::init_construction_hdiv()
733  {
734  HdivSpace<Scalar>* ref_space;
735  if(this->coarse_space->own_shapeset)
736  ref_space = new HdivSpace<Scalar>(this->ref_mesh, this->coarse_space->essential_bcs, 1);
737  else
738  ref_space = new HdivSpace<Scalar>(this->ref_mesh, this->coarse_space->essential_bcs, 1, this->coarse_space->shapeset);
739  return ref_space;
740  }
741 
742  template<typename Scalar>
743  void Space<Scalar>::ReferenceSpaceCreator::finish_construction(Space<Scalar>* ref_space)
744  {
745  ref_space->seq = g_space_seq++;
746 
747  Element *e;
748  for_all_active_elements(e, coarse_space->get_mesh())
749  {
750  if(this->coarse_space->edata[e->id].changed_in_last_adaptation)
751  {
752  if(ref_space->mesh->get_element(e->id)->active)
753  ref_space->edata[e->id].changed_in_last_adaptation = true;
754  else
755  for(unsigned int i = 0; i < 4; i++)
756  if(ref_space->mesh->get_element(e->id)->sons[i] != NULL)
757  if(ref_space->mesh->get_element(e->id)->sons[i]->active)
758  ref_space->edata[ref_space->mesh->get_element(e->id)->sons[i]->id].changed_in_last_adaptation = true;
759  }
760  }
761  }
762 
763  template<typename Scalar>
764  void Space<Scalar>::update_orders_recurrent(Element* e, int order)
765  {
766  // Adjust wrt. max and min possible orders.
767  int mo = shapeset->get_max_order();
768  int lower_limit = (get_type() == HERMES_L2_SPACE || get_type() == HERMES_HCURL_SPACE) ? 0 : 1; // L2 and Hcurl may use zero orders.
769  int ho = std::max(lower_limit, std::min(H2D_GET_H_ORDER(order), mo));
770  int vo = std::max(lower_limit, std::min(H2D_GET_V_ORDER(order), mo));
771  order = e->is_triangle() ? ho : H2D_MAKE_QUAD_ORDER(ho, vo);
772 
773  if(e->active)
774  edata[e->id].order = order;
775  else
776  for (int i = 0; i < 4; i++)
777  if(e->sons[i] != NULL)
778  update_orders_recurrent(e->sons[i], order);
779  }
780 
781  template<typename Scalar>
782  int Space<Scalar>::get_edge_order(Element* e, int edge) const
783  {
784  Node* en = e->en[edge];
785  if(en->id >= nsize || edge >= (int)e->get_nvert()) return 0;
786 
787  if(ndata[en->id].n == -1)
788  return get_edge_order_internal(ndata[en->id].base); // constrained node
789  else
790  return get_edge_order_internal(en);
791  }
792 
793  template<typename Scalar>
794  int Space<Scalar>::get_edge_order_internal(Node* en) const
795  {
796  assert(en->type == HERMES_TYPE_EDGE);
797  Element** e = en->elem;
798  int o1 = 1000, o2 = 1000;
799  assert(e[0] != NULL || e[1] != NULL);
800 
801  if(e[0] != NULL)
802  {
803  if(e[0]->is_triangle() || en == e[0]->en[0] || en == e[0]->en[2])
804  o1 = H2D_GET_H_ORDER(edata[e[0]->id].order);
805  else
806  o1 = H2D_GET_V_ORDER(edata[e[0]->id].order);
807  }
808 
809  if(e[1] != NULL)
810  {
811  if(e[1]->is_triangle() || en == e[1]->en[0] || en == e[1]->en[2])
812  o2 = H2D_GET_H_ORDER(edata[e[1]->id].order);
813  else
814  o2 = H2D_GET_V_ORDER(edata[e[1]->id].order);
815  }
816 
817  if(o1 == 0) return o2 == 1000 ? 0 : o2;
818  if(o2 == 0) return o1 == 1000 ? 0 : o1;
819  return std::min(o1, o2);
820  }
821 
822  template<typename Scalar>
823  void Space<Scalar>::set_mesh(Mesh* mesh)
824  {
825  if(this->mesh == mesh)
826  return;
827  free();
828  this->mesh = mesh;
829  this->mesh_seq = mesh->get_seq();
830  seq = g_space_seq++;
831  }
832 
833  template<typename Scalar>
834  void Space<Scalar>::set_mesh_seq(int seq)
835  {
836  this->mesh_seq = seq;
837  const_cast<Mesh*>(this->mesh)->set_seq(seq);
838  }
839 
840  template<typename Scalar>
841  void Space<Scalar>::update_constraints()
842  {
843  }
844 
845  template<typename Scalar>
846  void Space<Scalar>::post_assign()
847  {
848  }
849 
850  template<typename Scalar>
851  int Space<Scalar>::get_seq() const
852  {
853  return seq;
854  }
855 
856  template<typename Scalar>
857  void Space<Scalar>::distribute_orders(Mesh* mesh, int* parents)
858  {
859  int num = mesh->get_max_element_id();
860  int* orders = new int[num + 1];
861  Element* e;
862  for_all_active_elements(e, mesh)
863  {
864  int p = get_element_order(parents[e->id]);
865  if(e->is_triangle() && (H2D_GET_V_ORDER(p) != 0))
866  p = std::max(H2D_GET_H_ORDER(p), H2D_GET_V_ORDER(p));
867  orders[e->id] = p;
868  }
869  for_all_active_elements(e, mesh)
870  set_element_order_internal(e->id, orders[e->id]);
871  delete [] orders;
872  }
873 
874  template<typename Scalar>
875  int Space<Scalar>::assign_dofs(int first_dof, int stride)
876  {
877  if(ndata == NULL || edata == NULL || !nsize || !esize)
878  return false;
879  if(seq < 0)
880  return false;
881  if(this->mesh == NULL)
882  return false;
883 
884  this->mesh->check();
885 
886  if(edata == NULL)
887  {
888  throw Hermes::Exceptions::Exception("NULL edata detected in Space<Scalar>::get_element_order().");
889  return false;
890  }
891 
892  if(first_dof < 0)
893  throw Hermes::Exceptions::ValueException("first_dof", first_dof, 0);
894  if(stride < 1)
895  throw Hermes::Exceptions::ValueException("stride", stride, 1);
896 
897  resize_tables();
898 
899  Element* e;
900 
901  //check validity of orders
902  for_all_active_elements(e, mesh)
903  {
904  if(e->id >= esize || edata[e->id].order < 0)
905  {
906  printf("e->id = %d\n", e->id);
907  printf("esize = %d\n", esize);
908  printf("edata[%d].order = %d\n", e->id, edata[e->id].order);
909  throw
910  Hermes::Exceptions::Exception("Uninitialized element order in Space::assign_dofs().");
911  }
912  this->edata[e->id].changed_in_last_adaptation = true;
913  }
914 
915  this->first_dof = next_dof = first_dof;
916  this->stride = stride;
917 
918  reset_dof_assignment();
919  assign_vertex_dofs();
920  assign_edge_dofs();
921  assign_bubble_dofs();
922 
923  free_bc_data();
924  update_essential_bc_values();
925  update_constraints();
926  post_assign();
927 
928  mesh_seq = mesh->get_seq();
929  was_assigned = this->seq;
930  this->ndof = (next_dof - first_dof) / stride;
931 
932  return this->ndof;
933  check();
934  }
935 
936  template<typename Scalar>
937  void Space<Scalar>::reset_dof_assignment()
938  {
939  // First assume that all vertex nodes are part of a natural BC. the member NodeData::n
940  // is misused for this purpose, since it stores nothing at this point. Also assume
941  // that all DOFs are unassigned.
942  int i, j;
943  for (i = 0; i < mesh->get_max_node_id(); i++)
944  {
945  ndata[i].n = 1; // Natural boundary condition. The point is that it is not (0 == Dirichlet).
946  ndata[i].dof = H2D_UNASSIGNED_DOF;
947  }
948 
949  // next go through all boundary edge nodes constituting an essential BC and mark their
950  // neighboring vertex nodes also as essential
951  Element* e;
952  for_all_active_elements(e, mesh)
953  {
954  for (unsigned int i = 0; i < e->get_nvert(); i++)
955  {
956  if(e->en[i]->bnd)
957  if(essential_bcs != NULL)
958  if(essential_bcs->get_boundary_condition(mesh->boundary_markers_conversion.get_user_marker(e->en[i]->marker).marker) != NULL)
959  {
960  j = e->next_vert(i);
961  ndata[e->vn[i]->id].n = 0;
962  ndata[e->vn[j]->id].n = 0;
963  }
964  }
965  }
966  }
967 
968  template<typename Scalar>
969  int Space<Scalar>::get_vertex_functions_count()
970  {
971  return this->vertex_functions_count;
972  }
973 
974  template<typename Scalar>
975  int Space<Scalar>::get_edge_functions_count()
976  {
977  return this->edge_functions_count;
978  }
979 
980  template<typename Scalar>
981  int Space<Scalar>::get_bubble_functions_count()
982  {
983  return this->bubble_functions_count;
984  }
985 
986  template<typename Scalar>
987  void Space<Scalar>::get_element_assembly_list(Element* e, AsmList<Scalar>* al, unsigned int first_dof) const
988  {
989  this->check();
990  // some checks
991  if(e->id >= esize || edata[e->id].order < 0)
992  throw Hermes::Exceptions::Exception("Uninitialized element order in get_element_assembly_list(id = #%d).", e->id);
993  if(!is_up_to_date())
994  throw Hermes::Exceptions::Exception("The space in get_element_assembly_list() is out of date. You need to update it with assign_dofs()"
995  " any time the mesh changes.");
996 
997  // add vertex, edge and bubble functions to the assembly list
998  al->cnt = 0;
999  for (unsigned int i = 0; i < e->get_nvert(); i++)
1000  get_vertex_assembly_list(e, i, al);
1001  for (unsigned int i = 0; i < e->get_nvert(); i++)
1002  get_boundary_assembly_list_internal(e, i, al);
1003  get_bubble_assembly_list(e, al);
1004  for(unsigned int i = 0; i < al->cnt; i++)
1005  if(al->dof[i] >= 0)
1006  al->dof[i] += first_dof;
1007  }
1008 
1009  template<typename Scalar>
1010  void Space<Scalar>::get_boundary_assembly_list(Element* e, int surf_num, AsmList<Scalar>* al, unsigned int first_dof) const
1011  {
1012  this->check();
1013  al->cnt = 0;
1014  get_vertex_assembly_list(e, surf_num, al);
1015  get_vertex_assembly_list(e, e->next_vert(surf_num), al);
1016  get_boundary_assembly_list_internal(e, surf_num, al);
1017  for(unsigned int i = 0; i < al->cnt; i++)
1018  if(al->dof[i] >= 0)
1019  al->dof[i] += first_dof;
1020  }
1021 
1022  template<typename Scalar>
1023  void Space<Scalar>::get_bubble_assembly_list(Element* e, AsmList<Scalar>* al) const
1024  {
1025  this->check();
1026  ElementData* ed = &edata[e->id];
1027 
1028  if(!ed->n) return;
1029 
1030  int* indices = shapeset->get_bubble_indices(ed->order, e->get_mode());
1031  for (int i = 0, dof = ed->bdof; i < ed->n; i++, dof += stride, indices++)
1032  al->add_triplet(*indices, dof, 1.0);
1033  }
1034 
1035  template<typename Scalar>
1036  void Space<Scalar>::set_essential_bcs(EssentialBCs<Scalar>* essential_bcs)
1037  {
1038  this->essential_bcs = essential_bcs;
1039  }
1040 
1041  template<typename Scalar>
1042  void Space<Scalar>::precalculate_projection_matrix(int nv, double**& mat, double*& p)
1043  {
1044  int n = shapeset->get_max_order() + 1 - nv;
1045  mat = new_matrix<double>(n, n);
1046  int component = (get_type() == HERMES_HDIV_SPACE) ? 1 : 0;
1047 
1048  Quad1DStd quad1d;
1049  for (int i = 0; i < n; i++)
1050  {
1051  for (int j = i; j < n; j++)
1052  {
1053  int o = i + j + 4;
1054  double2* pt = quad1d.get_points(o);
1055  int ii = shapeset->get_edge_index(0, 0, i + nv, HERMES_MODE_QUAD);
1056  int ij = shapeset->get_edge_index(0, 0, j + nv, HERMES_MODE_QUAD);
1057  double val = 0.0;
1058  for (int k = 0; k < quad1d.get_num_points(o); k++)
1059  {
1060  val += pt[k][1] * shapeset->get_fn_value(ii, pt[k][0], -1.0, component, HERMES_MODE_QUAD)
1061  * shapeset->get_fn_value(ij, pt[k][0], -1.0, component, HERMES_MODE_QUAD);
1062  }
1063  mat[i][j] = val;
1064  }
1065  }
1066 
1067  p = new double[n];
1068  choldc(mat, n, p);
1069  }
1070 
1071  template<typename Scalar>
1072  void Space<Scalar>::update_edge_bc(Element* e, SurfPos* surf_pos)
1073  {
1074  if(e->active)
1075  {
1076  Node* en = e->en[surf_pos->surf_num];
1077  NodeData* nd = &ndata[en->id];
1078  nd->edge_bc_proj = NULL;
1079 
1080  if(nd->dof != H2D_UNASSIGNED_DOF && en->bnd)
1081  if(essential_bcs != NULL)
1082  {
1083  EssentialBoundaryCondition<Scalar> *bc = this->essential_bcs->get_boundary_condition(this->mesh->boundary_markers_conversion.get_user_marker(en->marker).marker);
1084  if(bc != NULL)
1085  {
1086  int order = get_edge_order_internal(en);
1087  surf_pos->marker = en->marker;
1088  nd->edge_bc_proj = get_bc_projection(surf_pos, order, bc);
1089  bc_data.push_back(nd->edge_bc_proj);
1090 
1091  int i = surf_pos->surf_num, j = e->next_vert(i);
1092  ndata[e->vn[i]->id].vertex_bc_coef = nd->edge_bc_proj + 0;
1093  ndata[e->vn[j]->id].vertex_bc_coef = nd->edge_bc_proj + 1;
1094  }
1095  }
1096  }
1097  else
1098  {
1099  int son1, son2;
1100  if(mesh->get_edge_sons(e, surf_pos->surf_num, son1, son2) == 2)
1101  {
1102  double mid = (surf_pos->lo + surf_pos->hi) * 0.5, tmp = surf_pos->hi;
1103  surf_pos->hi = mid;
1104  update_edge_bc(e->sons[son1], surf_pos);
1105  surf_pos->lo = mid; surf_pos->hi = tmp;
1106  update_edge_bc(e->sons[son2], surf_pos);
1107  }
1108  else
1109  update_edge_bc(e->sons[son1], surf_pos);
1110  }
1111  }
1112 
1113  template<typename Scalar>
1114  void Space<Scalar>::update_essential_bc_values()
1115  {
1116  Element* e;
1117  for_all_base_elements(e, mesh)
1118  {
1119  for (unsigned int i = 0; i < e->get_nvert(); i++)
1120  {
1121  int j = e->next_vert(i);
1122  if(e->vn[i]->bnd && e->vn[j]->bnd)
1123  {
1124  SurfPos surf_pos = {0, i, e, e->vn[i]->id, e->vn[j]->id, 0.0, 0.0, 1.0};
1125  update_edge_bc(e, &surf_pos);
1126  }
1127  }
1128  }
1129  }
1130 
1131  template<typename Scalar>
1132  void Space<Scalar>::free_bc_data()
1133  {
1134  for (unsigned int i = 0; i < bc_data.size(); i++)
1135  delete (bc_data[i]);
1136  bc_data.clear();
1137  }
1138 
1139  template<typename Scalar>
1140  bool Space<Scalar>::save(const char *filename) const
1141  {
1142  this->check();
1143  XMLSpace::space xmlspace;
1144 
1145  switch(this->get_type())
1146  {
1147  case HERMES_H1_SPACE:
1148  xmlspace.spaceType().set("h1");
1149  break;
1150  case HERMES_HCURL_SPACE:
1151  xmlspace.spaceType().set("hcurl");
1152  break;
1153  case HERMES_HDIV_SPACE:
1154  xmlspace.spaceType().set("hdiv");
1155  break;
1156  case HERMES_L2_SPACE:
1157  xmlspace.spaceType().set("l2");
1158  break;
1159  default:
1160  return false;
1161  }
1162 
1163  // Utility pointer.
1164  Element *e;
1165  for_all_elements(e, this->get_mesh())
1166  xmlspace.element_data().push_back(XMLSpace::space::element_data_type(e->id, this->edata[e->id].order, this->edata[e->id].bdof, this->edata[e->id].n, this->edata[e->id].changed_in_last_adaptation));
1167 
1168  std::string space_schema_location(Hermes2DApi.get_text_param_value(xmlSchemasDirPath));
1169  space_schema_location.append("/space_h2d_xml.xsd");
1170  ::xml_schema::namespace_info namespace_info_space("XMLSpace", space_schema_location);
1171 
1172  ::xml_schema::namespace_infomap namespace_info_map;
1173  namespace_info_map.insert(std::pair<std::basic_string<char>, xml_schema::namespace_info>("space", namespace_info_space));
1174 
1175  std::ofstream out(filename);
1176 
1177  ::xml_schema::flags parsing_flags = ::xml_schema::flags::dont_pretty_print;
1178 
1179  XMLSpace::space_(out, xmlspace, namespace_info_map, "UTF-8", parsing_flags);
1180  out.close();
1181 
1182  return true;
1183  }
1184 
1185  template<typename Scalar>
1186  Space<Scalar>* Space<Scalar>::load(const char *filename, Mesh* mesh, bool validate, EssentialBCs<Scalar>* essential_bcs, Shapeset* shapeset)
1187  {
1188  try
1189  {
1190  Space<Scalar>* space;
1191 
1192  ::xml_schema::flags parsing_flags = 0;
1193 
1194  if(!validate)
1195  parsing_flags = xml_schema::flags::dont_validate;
1196 
1197  std::auto_ptr<XMLSpace::space> parsed_xml_space (XMLSpace::space_(filename, parsing_flags));
1198 
1199  if(!strcmp(parsed_xml_space->spaceType().get().c_str(),"h1"))
1200  {
1201  space = new H1Space<Scalar>();
1202  space->mesh = mesh;
1203 
1204  if(shapeset == NULL)
1205  {
1206  space->shapeset = new H1Shapeset;
1207  space->own_shapeset = true;
1208  }
1209  else
1210  {
1211  if(shapeset->get_space_type() != HERMES_H1_SPACE)
1212  throw Hermes::Exceptions::SpaceLoadFailureException("Wrong shapeset / Wrong spaceType in the Solution XML file %s in Space::load.", filename);
1213  else
1214  space->shapeset = shapeset;
1215  }
1216 
1217  space->precalculate_projection_matrix(2, space->proj_mat, space->chol_p);
1218  }
1219  else if (!strcmp(parsed_xml_space->spaceType().get().c_str(),"hcurl"))
1220  {
1221  space = new HcurlSpace<Scalar>();
1222  space->mesh = mesh;
1223 
1224  if(shapeset == NULL)
1225  {
1226  space->shapeset = new HcurlShapeset;
1227  space->own_shapeset = true;
1228  }
1229  else
1230  {
1231  if(shapeset->get_num_components() < 2)
1232  throw Hermes::Exceptions::Exception("HcurlSpace requires a vector shapeset in Space::load.");
1233  if(shapeset->get_space_type() != HERMES_HCURL_SPACE)
1234  throw Hermes::Exceptions::SpaceLoadFailureException("Wrong shapeset / Wrong spaceType in the Solution XML file %s in Space::load.", filename);
1235  else
1236  space->shapeset = shapeset;
1237  }
1238 
1239  space->precalculate_projection_matrix(0, space->proj_mat, space->chol_p);
1240  }
1241  else if(!!strcmp(parsed_xml_space->spaceType().get().c_str(),"hdiv"))
1242  {
1243  space = new HdivSpace<Scalar>();
1244  space->mesh = mesh;
1245 
1246  if(shapeset == NULL)
1247  {
1248  space->shapeset = new HdivShapeset;
1249  space->own_shapeset = true;
1250  }
1251  else
1252  {
1253  if(shapeset->get_num_components() < 2)
1254  throw Hermes::Exceptions::Exception("HdivSpace requires a vector shapeset in Space::load.");
1255  if(shapeset->get_space_type() != HERMES_HDIV_SPACE)
1256  throw Hermes::Exceptions::SpaceLoadFailureException("Wrong shapeset / Wrong spaceType in the Solution XML file %s in Space::load.", filename);
1257  else
1258  space->shapeset = shapeset;
1259  }
1260 
1261  space->precalculate_projection_matrix(0, space->proj_mat, space->chol_p);
1262  }
1263  else if(strcmp(parsed_xml_space->spaceType().get().c_str(),"l2"))
1264  {
1265  space = new L2Space<Scalar>();
1266  space->mesh = mesh;
1267 
1268  if(shapeset == NULL)
1269  {
1270  space->shapeset = new L2Shapeset;
1271  space->own_shapeset = true;
1272  }
1273  {
1274  if(shapeset->get_space_type() != HERMES_L2_SPACE)
1275  throw Hermes::Exceptions::SpaceLoadFailureException("Wrong shapeset / Wrong spaceType in the Solution XML file %s in Space::load.", filename);
1276  else
1277  space->shapeset = shapeset;
1278  }
1279 
1280  static_cast<L2Space<Scalar>*>(space)->ldata = NULL;
1281  static_cast<L2Space<Scalar>*>(space)->lsize = 0;
1282  }
1283  else
1284  {
1285  throw Exceptions::SpaceLoadFailureException("Wrong spaceType in the Solution XML file %s in Space::load.", filename);
1286  return NULL;
1287  }
1288 
1289  space->essential_bcs = essential_bcs;
1290  space->mesh_seq == space->mesh->get_seq();
1291 
1292  // L2 space does not have any (strong) essential BCs.
1293  if(essential_bcs != NULL && parsed_xml_space->spaceType().get().c_str() != "l2")
1294  for(typename Hermes::vector<EssentialBoundaryCondition<Scalar>*>::const_iterator it = essential_bcs->begin(); it != essential_bcs->end(); it++)
1295  for(unsigned int i = 0; i < (*it)->markers.size(); i++)
1296  if(space->get_mesh()->boundary_markers_conversion.conversion_table_inverse.find((*it)->markers.at(i)) == space->get_mesh()->boundary_markers_conversion.conversion_table_inverse.end())
1297  throw Hermes::Exceptions::Exception("A boundary condition defined on a non-existent marker.");
1298 
1299  space->resize_tables();
1300 
1301  // Element data //
1302  unsigned int elem_data_count = parsed_xml_space->element_data().size();
1303  for (unsigned int elem_data_i = 0; elem_data_i < elem_data_count; elem_data_i++)
1304  {
1305  space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].order = parsed_xml_space->element_data().at(elem_data_i).ord();
1306  space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].bdof = parsed_xml_space->element_data().at(elem_data_i).bd();
1307  space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].n = parsed_xml_space->element_data().at(elem_data_i).n();
1308  space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].changed_in_last_adaptation = parsed_xml_space->element_data().at(elem_data_i).chgd();
1309  }
1310 
1311  space->seq = g_space_seq++;
1312 
1313  space->assign_dofs();
1314 
1315  return space;
1316  }
1317  catch (const xml_schema::exception& e)
1318  {
1319  throw Hermes::Exceptions::SpaceLoadFailureException(e.what());
1320  }
1321  }
1322 
1323  template class HERMES_API Space<double>;
1324  template class HERMES_API Space<std::complex<double> >;
1325  }
1326 }