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"
31 template<
typename Scalar>
36 template<
typename Scalar>
37 SpaceSharedPtr<Scalar>::SpaceSharedPtr(
const SpaceSharedPtr& other) :
std::tr1::shared_ptr<
Hermes::Hermes2D::Space<Scalar> >(other)
41 template<
typename Scalar>
42 void SpaceSharedPtr<Scalar>::operator=(
const SpaceSharedPtr& other)
44 std::tr1::shared_ptr<Hermes::Hermes2D::Space<Scalar> >::operator=(other);
47 const char* spaceTypeToString(SpaceType spaceType)
53 case HERMES_HCURL_SPACE:
55 case HERMES_HDIV_SPACE:
59 case HERMES_L2_MARKERWISE_CONST_SPACE:
60 return "l2-markerwise";
61 case HERMES_INVALID_SPACE:
66 SpaceType spaceTypeFromString(
const char* spaceTypeString)
68 if (!strcmp(spaceTypeString,
"h1"))
69 return HERMES_H1_SPACE;
70 if (!strcmp(spaceTypeString,
"hcurl"))
71 return HERMES_HCURL_SPACE;
72 if (!strcmp(spaceTypeString,
"hdiv"))
73 return HERMES_HDIV_SPACE;
74 if (!strcmp(spaceTypeString,
"l2"))
75 return HERMES_L2_SPACE;
76 if (!strcmp(spaceTypeString,
"l2-markerwise"))
77 return HERMES_L2_MARKERWISE_CONST_SPACE;
78 if (!strcmp(spaceTypeString,
"invalid"))
79 return HERMES_INVALID_SPACE;
82 unsigned g_space_seq = 0;
84 template<
typename Scalar>
87 this->ndata =
nullptr;
88 this->edata =
nullptr;
89 this->nsize = esize = 0;
91 this->seq = g_space_seq++;
92 this->seq_assigned = -1;
94 this->proj_mat =
nullptr;
95 this->chol_p =
nullptr;
96 this->vertex_functions_count = this->edge_functions_count = this->bubble_functions_count = 0;
98 if (essential_bcs !=
nullptr)
101 for (
unsigned int i = 0; i < (*it)->markers.size(); i++)
102 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)
103 throw Hermes::Exceptions::Exception(
"A boundary condition defined on a non-existent marker %s.", (*it)->markers.at(i).c_str());
105 own_shapeset = (shapeset ==
nullptr);
108 template<
typename Scalar>
114 free_with_check(ndata,
true);
119 free_with_check(edata,
true);
132 Space<std::complex<double> >::Space() : shapeset(nullptr), essential_bcs(nullptr)
137 template<
typename Scalar>
138 Space<Scalar>::Space(MeshSharedPtr mesh, Shapeset* shapeset, EssentialBCs<Scalar>* essential_bcs)
139 : shapeset(shapeset), essential_bcs(essential_bcs), mesh(mesh)
142 throw Hermes::Exceptions::NullException(0);
146 template<
typename Scalar>
153 throw Hermes::Exceptions::Exception(
"Mesh not present in Space.");
157 if (ndata ==
nullptr || edata ==
nullptr || !nsize || !esize)
166 if (edata ==
nullptr)
168 throw Hermes::Exceptions::Exception(
"nullptr edata detected in Space<Scalar>.");
172 if (!this->is_up_to_date())
174 throw Hermes::Exceptions::Exception(
"Space is not up to date.");
181 template<
typename Scalar>
185 free_with_check(this->proj_mat,
true);
186 free_with_check(this->chol_p);
188 if (this->own_shapeset)
189 delete this->shapeset;
192 template<
typename Scalar>
195 unsigned char prev = e->prev_vert(i);
196 if (e->is_triangle())
197 return e->
sons[3]->
vn[e->prev_vert(i)];
198 else if (e->
sons[2] ==
nullptr)
199 return i == 1 ? e->
sons[0]->
vn[2] : i == 3 ? e->
sons[0]->
vn[3] :
nullptr;
200 else if (e->
sons[0] ==
nullptr)
201 return i == 0 ? e->
sons[2]->
vn[1] : i == 2 ? e->
sons[2]->
vn[2] :
nullptr;
202 else return e->
sons[i]->
vn[j];
205 template<
typename Scalar>
208 if ((nsize < mesh->get_max_node_id()) || (ndata ==
nullptr))
213 while (nsize < mesh->get_max_node_id())
214 nsize = nsize * 3 / 2;
215 ndata = realloc_with_check<Space<Scalar>,
NodeData>(ndata, nsize,
this);
216 for (
int i = oldsize; i < nsize; i++)
218 ndata[i].edge_bc_proj =
nullptr;
219 ndata[i].baselist =
nullptr;
220 ndata[i].base =
nullptr;
224 if ((esize < mesh->get_max_element_id()) || (edata ==
nullptr))
229 while (esize < mesh->get_max_element_id())
230 esize = esize * 3 / 2;
231 edata = realloc_with_check<Space<Scalar>,
ElementData>(edata, nsize,
this);
232 for (
int i = oldsize; i < esize; i++)
237 edata[i].changed_in_last_adaptation =
true;
242 template<
typename Scalar>
246 this->vertex_functions_count = this->edge_functions_count = this->bubble_functions_count = 0;
248 this->essential_bcs = space->essential_bcs;
250 if (new_mesh->get_seq() != space->get_mesh()->get_seq())
252 new_mesh->copy(space->get_mesh());
253 this->mesh = new_mesh;
256 this->mesh = space->get_mesh();
258 this->resize_tables();
261 for_all_active_elements(e, this->mesh)
263 this->set_element_order_internal(e->
id, space->get_element_order(e->
id));
266 this->seq = g_space_seq++;
268 for_all_active_elements(e, this->mesh)
270 if (space->edata[e->
id].changed_in_last_adaptation)
271 this->edata[e->
id].changed_in_last_adaptation =
true;
273 this->edata[e->
id].changed_in_last_adaptation =
false;
279 template<
typename Scalar>
282 return this->shapeset;
285 template<
typename Scalar>
292 template<
typename Scalar>
299 template<
typename Scalar>
305 template<
typename Scalar>
308 return (seq_assigned == this->seq) && (mesh_seq == mesh->get_seq());
311 template<
typename Scalar>
315 return essential_bcs;
318 template<
typename Scalar>
321 set_element_order_internal(
id, order, order_v);
324 template<
typename Scalar>
328 if (id < 0 || id >= mesh->get_max_element_id())
329 throw Hermes::Exceptions::Exception(
"Space<Scalar>::set_element_order_internal: Invalid element id.");
334 if (mesh->get_element(
id)->is_quad() && H2D_GET_V_ORDER(order) == 0)
336 order = H2D_MAKE_QUAD_ORDER(order, order_v);
338 order = H2D_MAKE_QUAD_ORDER(order, order);
340 edata[
id].order = order;
344 template<
typename Scalar>
347 int n = spaces.size();
348 for (
int i = 0; i < n; i++)
350 if (spaces[i]->get_essential_bcs() !=
nullptr)
351 spaces[i]->get_essential_bcs()->set_current_time(time);
352 spaces[i]->update_essential_bc_values();
356 template<
typename Scalar>
359 space->get_essential_bcs()->set_current_time(time);
360 space->update_essential_bc_values();
363 template<
typename Scalar>
367 for (
unsigned char i = 0; i < spaces.size(); i++)
368 ndof += spaces[i]->get_num_dofs();
372 template<
typename Scalar>
375 return space->get_num_dofs();
378 template<
typename Scalar>
381 int n = spaces.size();
384 for (
int i = 0; i < n; i++) {
385 ndof += spaces[i]->assign_dofs(ndof);
391 template<
typename Scalar>
394 if (marker == HERMES_ANY)
395 set_uniform_order_internal(order, -1234);
397 set_uniform_order_internal(order, mesh->element_markers_conversion.get_internal_marker(marker).marker);
400 template<
typename Scalar>
404 int quad_order = H2D_MAKE_QUAD_ORDER(order, order);
407 for_all_active_elements(e, mesh)
409 if (marker == HERMES_ANY_INT || e->
marker == marker)
412 if (e->is_triangle())
415 ed->order = quad_order;
421 template<
typename Scalar>
425 for_all_active_elements(e, this->get_mesh())
427 if (e->is_triangle())
428 set_element_order_internal(e->
id, std::max<int>(min_order, get_element_order(e->
id) + order_change));
431 int h_order, v_order;
438 if (H2D_GET_V_ORDER(get_element_order(e->
id)) + order_change < min_order)
441 v_order = H2D_GET_V_ORDER(get_element_order(e->
id)) + order_change;
443 set_element_order_internal(e->
id, H2D_MAKE_QUAD_ORDER(h_order, v_order));
448 template<
typename Scalar>
452 for_all_active_elements(e, this->get_mesh())
454 if (e->is_triangle())
456 this->warn(
"Using quad version of Space<Scalar>::adjust_element_order(), only horizontal orders will be used.");
457 set_element_order_internal(e->
id, std::max<int>(horizontal_min_order, get_element_order(e->
id) + horizontal_order_change));
460 if (get_element_order(e->
id) == -1)
461 set_element_order_internal(e->
id, H2D_MAKE_QUAD_ORDER(horizontal_min_order, vertical_min_order));
463 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)));
467 template<
typename Scalar>
471 std::vector<int> list;
473 for_all_inactive_elements(e, this->mesh)
476 for (
unsigned int i = 0; i < 4; i++)
477 if (e->
sons[i] !=
nullptr &&
478 (!e->
sons[i]->
active || (keep_initial_refinements && e->
sons[i]->
id < this->mesh->ninitial))
481 found =
false;
break;
484 if (found) list.push_back(e->
id);
488 for (
unsigned int i = 0; i < list.size(); i++)
490 unsigned int order = 0, h_order = 0, v_order = 0;
491 unsigned int num_sons = 0;
492 if (this->mesh->get_element_fast(list[i])->bsplit())
495 for (
int sons_i = 0; sons_i < 4; sons_i++)
497 if (this->mesh->get_element_fast(list[i])->sons[sons_i]->active)
499 if (this->mesh->get_element_fast(list[i])->sons[sons_i]->is_triangle())
500 order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[sons_i]->id);
503 h_order +=
H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[sons_i]->id));
504 v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[sons_i]->id));
511 if (this->mesh->get_element_fast(list[i])->hsplit())
514 if (this->mesh->get_element_fast(list[i])->sons[0]->active)
516 if (this->mesh->get_element_fast(list[i])->sons[0]->is_triangle())
517 order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[0]->id);
520 h_order +=
H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[0]->id));
521 v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[0]->id));
524 if (this->mesh->get_element_fast(list[i])->sons[1]->active)
526 if (this->mesh->get_element_fast(list[i])->sons[1]->is_triangle())
527 order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[1]->id);
530 h_order +=
H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[1]->id));
531 v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[1]->id));
538 if (this->mesh->get_element_fast(list[i])->sons[2]->active)
540 if (this->mesh->get_element_fast(list[i])->sons[2]->is_triangle())
541 order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[2]->id);
544 h_order +=
H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[2]->id));
545 v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[2]->id));
548 if (this->mesh->get_element_fast(list[i])->sons[3]->active)
550 if (this->mesh->get_element_fast(list[i])->sons[3]->is_triangle())
551 order += this->get_element_order(this->mesh->get_element_fast(list[i])->sons[3]->id);
554 h_order +=
H2D_GET_H_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[3]->id));
555 v_order += H2D_GET_V_ORDER(this->get_element_order(this->mesh->get_element_fast(list[i])->sons[3]->id));
560 order = (
unsigned int)(order / num_sons);
561 h_order = (
unsigned int)(h_order / num_sons);
562 v_order = (
unsigned int)(v_order / num_sons);
564 if (this->mesh->get_element_fast(list[i])->is_triangle())
565 edata[list[i]].order = order;
567 edata[list[i]].order = H2D_MAKE_QUAD_ORDER(h_order, v_order);
569 if (!only_unrefine_space_data)
570 this->mesh->unrefine_element_id(list[i]);
574 for_all_active_elements(e, this->mesh)
575 this->edata[e->
id].changed_in_last_adaptation =
true;
578 template<
typename Scalar>
582 this->unrefine_all_mesh_elements_internal(keep_initial_refinements,
false);
585 template<
typename Scalar>
588 for (
unsigned char i = 0; i < spaces.size() - 1; i++)
589 spaces[i]->unrefine_all_mesh_elements_internal(keep_initial_refinements,
true);
591 spaces[spaces.size() - 1]->unrefine_all_mesh_elements_internal(keep_initial_refinements,
false);
594 template<
typename Scalar>
599 template<
typename Scalar>
604 template<
typename Scalar>
608 for_all_active_elements(e, coarse_space->get_mesh())
611 int coarse_element_id = e->
id;
614 int current_order = coarse_space->get_element_order(coarse_element_id);
616 if (current_order < 0)
617 throw Hermes::Exceptions::Exception(
"Source space has an uninitialized order (element id = %d)", coarse_element_id);
621 if (e->is_triangle())
624 int current_order_triangular = current_order;
628 new_order = current_order_triangular + this->order_increase;
634 int current_order_quadrilateral_horizontal =
H2D_GET_H_ORDER(current_order);
636 int current_order_quadrilateral_vertical = H2D_GET_V_ORDER(current_order);
639 int new_order_quadrilateral_horizontal = current_order_quadrilateral_horizontal + this->order_increase;
640 int new_order_quadrilateral_vertical = current_order_quadrilateral_vertical + this->order_increase;
641 new_order = H2D_MAKE_QUAD_ORDER(new_order_quadrilateral_horizontal, new_order_quadrilateral_vertical);
645 ref_space->update_orders_recurrent(ref_space->mesh->get_element(coarse_element_id), new_order);
649 template<
typename Scalar>
652 this->coarse_space = coarse_space;
653 this->ref_mesh = ref_mesh;
654 return this->create_ref_space(assign_dofs);
657 template<
typename Scalar>
663 return this->coarse_space;
666 if (
dynamic_cast<L2Space<Scalar>*
>(this->coarse_space.get()) !=
nullptr)
667 ref_space = this->init_construction_l2();
668 if (
dynamic_cast<H1Space<Scalar>*
>(this->coarse_space.get()) !=
nullptr)
669 ref_space = this->init_construction_h1();
671 ref_space = this->init_construction_hcurl();
673 ref_space = this->init_construction_hdiv();
675 if (ref_space ==
nullptr)
676 throw Exceptions::Exception(
"Something went wrong in ReferenceSpaceCreator::create_ref_space().");
679 this->handle_orders(ref_space);
682 this->finish_construction(ref_space);
686 ref_space->assign_dofs();
692 template<
typename Scalar>
695 if (this->coarse_space->own_shapeset)
701 template<
typename Scalar>
704 if (this->coarse_space->own_shapeset)
710 template<
typename Scalar>
711 SpaceSharedPtr<Scalar> Space<Scalar>::ReferenceSpaceCreator::init_construction_hcurl()
713 if (this->coarse_space->own_shapeset)
714 return SpaceSharedPtr<Scalar>(
new HcurlSpace<Scalar>(this->ref_mesh, this->coarse_space->get_essential_bcs(), 1));
716 return SpaceSharedPtr<Scalar>(
new HcurlSpace<Scalar>(this->ref_mesh, this->coarse_space->get_essential_bcs(), 1, this->coarse_space->get_shapeset()));
719 template<
typename Scalar>
720 SpaceSharedPtr<Scalar> Space<Scalar>::ReferenceSpaceCreator::init_construction_hdiv()
722 if (this->coarse_space->own_shapeset)
723 return SpaceSharedPtr<Scalar>(
new HdivSpace<Scalar>(this->ref_mesh, this->coarse_space->get_essential_bcs(), 1));
725 return SpaceSharedPtr<Scalar>(
new HdivSpace<Scalar>(this->ref_mesh, this->coarse_space->get_essential_bcs(), 1, this->coarse_space->get_shapeset()));
728 template<
typename Scalar>
729 void Space<Scalar>::ReferenceSpaceCreator::finish_construction(SpaceSharedPtr<Scalar> ref_space)
731 ref_space->seq = g_space_seq++;
734 for_all_active_elements(e, coarse_space->get_mesh())
736 bool to_set = this->coarse_space->edata[e->id].changed_in_last_adaptation;
738 if (ref_space->mesh->get_element(e->id)->active)
739 ref_space->edata[e->id].changed_in_last_adaptation = to_set;
741 for (
unsigned int i = 0; i < 4; i++)
742 if (ref_space->mesh->get_element(e->id)->sons[i] !=
nullptr)
743 if (ref_space->mesh->get_element(e->id)->sons[i]->active)
744 ref_space->edata[ref_space->mesh->get_element(e->id)->sons[i]->id].changed_in_last_adaptation = to_set;
749 template<
typename Scalar>
750 void Space<Scalar>::update_orders_recurrent(Element* e,
int order)
753 int mo = shapeset->get_max_order();
755 int lower_limit = (
get_type() == HERMES_L2_SPACE ||
get_type() == HERMES_L2_MARKERWISE_CONST_SPACE ||
get_type() == HERMES_HCURL_SPACE) ? 0 : 1;
757 int vo = std::max(lower_limit, std::min(H2D_GET_V_ORDER(order), mo));
758 order = e->is_triangle() ? ho : H2D_MAKE_QUAD_ORDER(ho, vo);
761 edata[e->id].order = order;
763 for (
int i = 0; i < 4; i++)
764 if (e->sons[i] !=
nullptr)
765 update_orders_recurrent(e->sons[i], order);
768 template<
typename Scalar>
772 if (en->
id >=
nsize || edge >= e->get_nvert())
782 template<
typename Scalar>
785 assert(en->
type == HERMES_TYPE_EDGE);
787 int o1 = 1000, o2 = 1000;
788 assert(e[0] !=
nullptr || e[1] !=
nullptr);
792 if (e[0]->is_triangle() || en == e[0]->en[0] || en == e[0]->en[2])
795 o1 = H2D_GET_V_ORDER(
edata[e[0]->
id].order);
800 if (e[1]->is_triangle() || en == e[1]->en[0] || en == e[1]->en[2])
803 o2 = H2D_GET_V_ORDER(
edata[e[1]->
id].order);
806 if (o1 == 0)
return o2 == 1000 ? 0 : o2;
807 if (o2 == 0)
return o1 == 1000 ? 0 : o1;
808 return std::min(o1, o2);
811 template<
typename Scalar>
814 if (this->mesh == mesh)
822 template<
typename Scalar>
826 this->mesh->set_seq(seq);
829 template<
typename Scalar>
834 template<
typename Scalar>
839 template<
typename Scalar>
845 template<
typename Scalar>
848 int num = mesh->get_max_element_id();
849 int* orders = malloc_with_check<Space<Scalar>,
int>(num + 1,
this);
851 for_all_active_elements(e, mesh)
854 if (e->is_triangle() && (H2D_GET_V_ORDER(p) != 0))
858 for_all_active_elements(e, mesh)
860 free_with_check(orders);
863 template<
typename Scalar>
870 if (this->mesh ==
nullptr)
875 if (
edata ==
nullptr)
877 throw Hermes::Exceptions::Exception(
"nullptr edata detected in Space<Scalar>::assign_dofs().");
882 throw Hermes::Exceptions::ValueException(
"first_dof", first_dof, 0);
889 assign_vertex_dofs();
891 assign_bubble_dofs();
906 template<
typename Scalar>
913 for (i = 0; i < mesh->get_max_node_id(); i++)
923 for_all_active_elements(e, mesh)
925 for (
unsigned char i = 0; i < e->get_nvert(); i++)
928 if (essential_bcs !=
nullptr)
929 if (essential_bcs->get_boundary_condition(mesh->boundary_markers_conversion.get_user_marker(e->
en[i]->
marker).marker) !=
nullptr)
939 template<
typename Scalar>
945 template<
typename Scalar>
948 return this->edge_functions_count;
951 template<
typename Scalar>
954 return this->bubble_functions_count;
957 template<
typename Scalar>
963 throw Hermes::Exceptions::Exception(
"Uninitialized element order in get_element_assembly_list(id = #%d).", e->
id);
965 throw Hermes::Exceptions::Exception(
"The space in get_element_assembly_list() is out of date. You need to update it with assign_dofs()"
966 " any time the mesh changes.");
970 for (
unsigned char i = 0; i < e->get_nvert(); i++)
971 get_vertex_assembly_list(e, i, al);
972 for (
unsigned char i = 0; i < e->get_nvert(); i++)
973 get_boundary_assembly_list_internal(e, i, al);
974 get_bubble_assembly_list(e, al);
977 template<
typename Scalar>
982 get_vertex_assembly_list(e, surf_num, al);
983 get_vertex_assembly_list(e, e->
next_vert(surf_num), al);
984 get_boundary_assembly_list_internal(e, surf_num, al);
987 template<
typename Scalar>
991 ElementData* ed = &
edata[e->
id];
995 short* indices = shapeset->get_bubble_indices(ed->order, e->get_mode());
996 for (
int i = 0, dof = ed->bdof; i < ed->n; i++, dof++, indices++)
1000 template<
typename Scalar>
1006 template<
typename Scalar>
1009 unsigned char n = shapeset->get_max_order() + 1 - nv;
1010 mat = new_matrix<double>(n, n);
1011 int component = (
get_type() == HERMES_HDIV_SPACE) ? 1 : 0;
1014 for (
unsigned char i = 0; i < n; i++)
1016 for (
unsigned char j = i; j < n; j++)
1019 double2* pt = quad1d.get_points(o);
1020 int ii = shapeset->get_edge_index(0, 0, i + nv, HERMES_MODE_QUAD);
1021 int ij = shapeset->get_edge_index(0, 0, j + nv, HERMES_MODE_QUAD);
1023 for (
int k = 0; k < quad1d.get_num_points(o); k++)
1025 double val_ii = shapeset->get_fn_value(ii, pt[k][0], -1.0, component, HERMES_MODE_QUAD);
1026 double val_ij = shapeset->get_fn_value(ij, pt[k][0], -1.0, component, HERMES_MODE_QUAD);
1027 val += pt[k][1] * val_ii * val_ij;
1033 p = malloc_with_check<Space<Scalar>,
double>(n,
this);
1037 template<
typename Scalar>
1038 void Space<Scalar>::update_edge_bc(Element* e, SurfPos* surf_pos)
1045 Node* en = e->en[surf_pos->surf_num];
1046 NodeData* nd = &
ndata[en->id];
1047 nd->edge_bc_proj =
nullptr;
1050 if (essential_bcs !=
nullptr)
1052 EssentialBoundaryCondition<Scalar> *bc = this->essential_bcs->get_boundary_condition(this->mesh->boundary_markers_conversion.get_user_marker(en->marker).marker);
1056 surf_pos->marker = en->marker;
1057 nd->edge_bc_proj = get_bc_projection(surf_pos, order, bc);
1060 int i = surf_pos->surf_num, j = e->next_vert(i);
1061 ndata[e->vn[i]->id].vertex_bc_coef = nd->edge_bc_proj + 0;
1062 ndata[e->vn[j]->id].vertex_bc_coef = nd->edge_bc_proj + 1;
1067 template<
typename Scalar>
1071 for_all_active_elements(e, mesh)
1073 for (
unsigned char i = 0; i < e->get_nvert(); i++)
1078 SurfPos surf_pos = { 0, i, e, e->
vn[i]->
id, e->vn[j]->id, 0.0, 0.0, 1.0 };
1079 update_edge_bc(e, &surf_pos);
1085 template<
typename Scalar>
1090 for (
unsigned int i = 0; i < bc_data_base_components.size(); i++)
1091 free_with_check(bc_data_base_components[i]);
1093 bc_data_base_components.clear();
1096 template<
typename Scalar>
1106 for_all_used_elements(e, this->get_mesh())
1110 std::ofstream out(filename);
1112 XMLSpace::space_(out, xmlspace, namespace_info_map,
"UTF-8", parsing_flags);
1117 template<
typename Scalar>
1128 bson_append_string(&bw,
"type", spaceTypeToString(this->
get_type()));
1131 bson_append_int(&bw,
"element_data_count", this->mesh->get_max_element_id());
1134 bson_append_start_array(&bw,
"orders");
1135 for (
int _id = 0, _max = this->get_mesh()->get_max_element_id(); _id < _max; _id++)
1136 bson_append_int(&bw,
"c", this->
edata[_id].order);
1137 bson_append_finish_array(&bw);
1139 bson_append_start_array(&bw,
"bdofs");
1140 for (
int _id = 0, _max = this->get_mesh()->get_max_element_id(); _id < _max; _id++)
1141 bson_append_int(&bw,
"c", this->
edata[_id].bdof);
1142 bson_append_finish_array(&bw);
1144 bson_append_start_array(&bw,
"ns");
1145 for (
int _id = 0, _max = this->get_mesh()->get_max_element_id(); _id < _max; _id++)
1146 bson_append_int(&bw,
"c", this->
edata[_id].n);
1147 bson_append_finish_array(&bw);
1149 bson_append_start_array(&bw,
"changed");
1150 for (
int _id = 0, _max = this->get_mesh()->get_max_element_id(); _id < _max; _id++)
1151 bson_append_bool(&bw,
"c", this->
edata[_id].changed_in_last_adaptation);
1152 bson_append_finish_array(&bw);
1159 fpw = fopen(filename,
"wb");
1160 const char *dataw = (
const char *)bson_data(&bw);
1161 fwrite(dataw, bson_size(&bw), 1, fpw);
1168 template<
typename Scalar>
1173 if (spaceType == HERMES_H1_SPACE)
1178 if (shapeset ==
nullptr)
1181 space->own_shapeset =
true;
1186 throw Hermes::Exceptions::SpaceLoadFailureException(
"Wrong shapeset / Wrong spaceType in Space loading subroutine.");
1188 space->shapeset = shapeset;
1191 space->precalculate_projection_matrix(2, space->proj_mat, space->chol_p);
1194 else if (spaceType == HERMES_HCURL_SPACE)
1199 if (shapeset ==
nullptr)
1202 space->own_shapeset =
true;
1207 throw Hermes::Exceptions::Exception(
"HcurlSpace requires a vector shapeset in Space::load.");
1209 throw Hermes::Exceptions::SpaceLoadFailureException(
"Wrong shapeset / Wrong spaceType in Space loading subroutine.");
1211 space->shapeset = shapeset;
1214 space->precalculate_projection_matrix(0, space->proj_mat, space->chol_p);
1217 else if (spaceType == HERMES_HDIV_SPACE)
1222 if (shapeset ==
nullptr)
1225 space->own_shapeset =
true;
1230 throw Hermes::Exceptions::Exception(
"HdivSpace requires a vector shapeset in Space::load.");
1232 throw Hermes::Exceptions::SpaceLoadFailureException(
"Wrong shapeset / Wrong spaceType in Space loading subroutine.");
1234 space->shapeset = shapeset;
1237 space->precalculate_projection_matrix(0, space->proj_mat, space->chol_p);
1240 else if (spaceType == HERMES_L2_SPACE)
1245 if (shapeset ==
nullptr)
1248 space->own_shapeset =
true;
1252 throw Hermes::Exceptions::SpaceLoadFailureException(
"Wrong shapeset / Wrong spaceType in Space loading subroutine.");
1254 space->shapeset = shapeset;
1258 else if (spaceType == HERMES_L2_MARKERWISE_CONST_SPACE)
1263 Hermes::Mixins::Loggable::Static::warn(
"L2MarkerWiseConstSpace does not need a shapeset when loading.");
1268 throw Exceptions::SpaceLoadFailureException(
"Wrong spaceType in Space loading subroutine.");
1275 template<
typename Scalar>
1285 parsing_flags = xml_schema::flags::dont_validate;
1287 std::auto_ptr<XMLSpace::space> parsed_xml_space(
XMLSpace::space_(filename, parsing_flags));
1289 if (spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) == HERMES_H1_SPACE)
1291 else if (spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) == HERMES_HCURL_SPACE)
1293 else if (spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) == HERMES_HDIV_SPACE)
1295 else if (spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) == HERMES_L2_SPACE)
1297 else if (spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) == HERMES_L2_MARKERWISE_CONST_SPACE)
1300 throw Hermes::Exceptions::IOException(Exceptions::IOException::Read, filename);
1303 space->mesh_seq = space->mesh->get_seq();
1304 space->init(shapeset, 1,
false);
1306 if (essential_bcs !=
nullptr && spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) != HERMES_L2_SPACE && spaceTypeFromString(parsed_xml_space->spaceType().get().c_str()) != HERMES_L2_MARKERWISE_CONST_SPACE)
1310 for (
unsigned int i = 0; i < (*it)->markers.size(); i++)
1311 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())
1312 throw Hermes::Exceptions::Exception(
"A boundary condition defined on a non-existent marker.");
1315 space->resize_tables();
1318 unsigned int elem_data_count = parsed_xml_space->element_data().size();
1319 for (
unsigned int elem_data_i = 0; elem_data_i < elem_data_count; elem_data_i++)
1321 space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].order = parsed_xml_space->element_data().at(elem_data_i).ord();
1322 space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].bdof = parsed_xml_space->element_data().at(elem_data_i).bd();
1323 space->edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].n = parsed_xml_space->element_data().at(elem_data_i).n();
1324 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();
1327 space->seq = g_space_seq++;
1329 space->assign_dofs();
1335 throw Hermes::Exceptions::SpaceLoadFailureException(e.what());
1339 template<
typename Scalar>
1347 parsing_flags = xml_schema::flags::dont_validate;
1349 std::auto_ptr<XMLSpace::space> parsed_xml_space(
XMLSpace::space_(filename, parsing_flags));
1351 if (strcmp(parsed_xml_space->spaceType().get().c_str(), spaceTypeToString(this->
get_type())))
1352 throw Exceptions::Exception(
"Saved Space is not of the same type as the current one in loading.");
1357 unsigned int elem_data_count = parsed_xml_space->element_data().size();
1358 for (
unsigned int elem_data_i = 0; elem_data_i < elem_data_count; elem_data_i++)
1360 this->
edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].order = parsed_xml_space->element_data().at(elem_data_i).ord();
1361 this->
edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].bdof = parsed_xml_space->element_data().at(elem_data_i).bd();
1362 this->
edata[parsed_xml_space->element_data().at(elem_data_i).e_id()].n = parsed_xml_space->element_data().at(elem_data_i).n();
1363 this->
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();
1366 this->
seq = g_space_seq++;
1372 throw Hermes::Exceptions::SpaceLoadFailureException(e.what());
1377 template<
typename Scalar>
1381 fpr = fopen(filename,
"rb");
1384 fseek(fpr, 0, SEEK_END);
1385 int size = ftell(fpr);
1389 char *datar = malloc_with_check<char>(size);
1390 fread(datar, size, 1, fpr);
1394 bson_init_finished_data(&br, datar, 0);
1398 bson_find(&it, &br,
"type");
1401 space->mesh_seq = space->mesh->get_seq();
1402 space->resize_tables();
1405 if (essential_bcs !=
nullptr && space->get_type() != HERMES_L2_SPACE && space->get_type() != HERMES_L2_MARKERWISE_CONST_SPACE)
1409 for (
unsigned int i = 0; i < (*it)->markers.size(); i++)
1410 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())
1411 throw Hermes::Exceptions::Exception(
"A boundary condition defined on a non-existent marker.");
1415 bson_find(&it, &br,
"element_data_count");
1416 if (bson_iterator_int(&it) != mesh->get_max_element_id())
1417 throw Exceptions::Exception(
"Mesh and saved space mixed in Space<Scalar>::load_bson.");
1420 bson_iterator it_coeffs;
1421 bson_find(&it_coeffs, &br,
"orders");
1422 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1423 bson_iterator_init(&it, &sub);
1424 int index_coeff = 0;
1425 while (bson_iterator_next(&it))
1426 space->edata[index_coeff++].order = bson_iterator_int(&it);
1428 bson_find(&it_coeffs, &br,
"bdofs");
1429 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1430 bson_iterator_init(&it, &sub);
1432 while (bson_iterator_next(&it))
1433 space->edata[index_coeff++].bdof = bson_iterator_int(&it);
1435 bson_find(&it_coeffs, &br,
"ns");
1436 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1437 bson_iterator_init(&it, &sub);
1439 while (bson_iterator_next(&it))
1440 space->edata[index_coeff++].n = bson_iterator_int(&it);
1442 bson_find(&it_coeffs, &br,
"changed");
1443 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1444 bson_iterator_init(&it, &sub);
1446 while (bson_iterator_next(&it))
1447 space->edata[index_coeff++].changed_in_last_adaptation = bson_iterator_int(&it);
1450 free_with_check(datar);
1452 space->seq = g_space_seq++;
1454 space->assign_dofs();
1459 template<
typename Scalar>
1460 void Space<Scalar>::load_bson(
const char *filename)
1463 fpr = fopen(filename,
"rb");
1466 fseek(fpr, 0, SEEK_END);
1467 int size = ftell(fpr);
1471 char *datar = malloc_with_check<char>(size);
1472 fread(datar, size, 1, fpr);
1476 bson_init_finished_data(&br, datar, 0);
1480 bson_find(&it, &br,
"type");
1482 if (strcmp(bson_iterator_string(&it), spaceTypeToString(this->
get_type())))
1483 throw Exceptions::Exception(
"Saved Space is not of the same type as the current one in loading.");
1486 bson_find(&it, &br,
"element_data_count");
1487 if (bson_iterator_int(&it) != mesh->get_max_element_id())
1488 throw Exceptions::Exception(
"Current and saved space mixed in Space<Scalar>::load_bson.");
1493 bson_iterator it_coeffs;
1494 bson_find(&it_coeffs, &br,
"orders");
1495 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1496 bson_iterator_init(&it, &sub);
1497 int index_coeff = 0;
1498 while (bson_iterator_next(&it))
1499 this->
edata[index_coeff++].order = bson_iterator_int(&it);
1501 bson_find(&it_coeffs, &br,
"bdofs");
1502 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1503 bson_iterator_init(&it, &sub);
1505 while (bson_iterator_next(&it))
1506 this->
edata[index_coeff++].bdof = bson_iterator_int(&it);
1508 bson_find(&it_coeffs, &br,
"ns");
1509 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1510 bson_iterator_init(&it, &sub);
1512 while (bson_iterator_next(&it))
1513 this->
edata[index_coeff++].n = bson_iterator_int(&it);
1515 bson_find(&it_coeffs, &br,
"changed");
1516 bson_iterator_subobject_init(&it_coeffs, &sub, 0);
1517 bson_iterator_init(&it, &sub);
1519 while (bson_iterator_next(&it))
1520 this->
edata[index_coeff++].changed_in_last_adaptation = bson_iterator_int(&it);
1523 free_with_check(datar);
1525 this->
seq = g_space_seq++;
1533 template<
typename Scalar>
1536 throw Hermes::Exceptions::MethodNotOverridenException(
"SettableSpaces<Scalar>::get_spaces()");
1539 template<
typename Scalar>
1542 return this->get_spaces()[n];
1545 template<
typename Scalar>
1548 std::vector<SpaceSharedPtr<Scalar> > spaces;
1549 spaces.push_back(space);
1550 this->set_spaces(spaces);
1553 template class HERMES_API SettableSpaces < double > ;
1554 template class HERMES_API SettableSpaces < std::complex<double> > ;
1557 template class HERMES_API Space < double > ;
1558 template class HERMES_API Space < std::complex<double> > ;
1560 template class HERMES_API SpaceSharedPtr < double > ;
1561 template class HERMES_API SpaceSharedPtr < std::complex<double> > ;
::xsd::cxx::tree::id< char, ncname > id
C++ type corresponding to the ID XML Schema built-in type.
int first_dof
For equation systems.
H1ShapesetJacobi H1Shapeset
Experimental.
unsigned type
0 = vertex node; 1 = edge node
static const int H2D_UNASSIGNED_DOF
DOF which was not assigned yet.
int ndof
Number of degrees of freedom (dimension of the space).
unsigned char next_vert(unsigned char i) const
Helper functions to obtain the index of the next or previous vertex/edge.
int get_bubble_functions_count()
Returns the total (global) number of bubble functions.
virtual bool isOkay() const
State querying helpers.
void save(const char *filename) const
Saves this space into a file.
::xsd::cxx::tree::flags flags
Parsing and serialization flags.
int vertex_functions_count
For statistics.
int mesh_seq
Tracking changes - mesh.
int get_edge_functions_count()
Returns the total (global) number of edge functions.
const spaceType_optional & spaceType() const
Return a read-only (constant) reference to the attribute container.
Stores one element of a mesh.
int get_element_order(int id) const
Returns element polynomial order.
HdivShapesetLegendre HdivShapeset
This is the default Hdiv shapeset typedef.
void free_bc_data()
Free BC data.
int get_vertex_functions_count()
void set_uniform_order(int order, std::string marker=HERMES_ANY)
void init()
Common code for constructors.
Element * elem[2]
elements sharing the edge node
void get_boundary_assembly_list(Element *e, int surf_num, AsmList< Scalar > *al) const
Obtains an edge assembly list (contains shape functions that are nonzero on the specified edge)...
int get_max_dof() const
Returns the DOF number of the last basis function.
virtual ~Space()
Destructor.
void set_mesh(MeshSharedPtr mesh)
Sets a (new) mesh and calls assign_dofs().
virtual void get_element_assembly_list(Element *e, AsmList< Scalar > *al) const
Obtains an assembly list for the given element.
Generated from space_h2d_xml.xsd.
::xsd::cxx::tree::exception< char > exception
Root of the C++/Tree exception hierarchy.
void add_triplet(int i, int d, Scalar c)
Adds a record for one basis function (shape functions index, basis functions index, coefficient).
std::vector< Scalar * > bc_data_projections
Used for bc projection.
Used to pass the instances of Space around.
HcurlShapesetGradLeg HcurlShapeset
This is the default Hcurl shapeset typedef.
void free()
Freeing the data.
Stores one node of a mesh.
void adjust_element_order(int order_change, int min_order)
const element_data_sequence & element_data() const
Return a read-only (constant) reference to the element sequence.
::xsd::cxx::tree::time< char, simple_type > time
C++ type corresponding to the time XML Schema built-in type.
virtual int get_edge_order(Element *e, int edge) const
Internal. Obtains the order of an edge, according to the minimum rule.
virtual int get_edge_order_internal(Node *en) const
Internal.
static SpaceSharedPtr< Scalar > init_empty_space(SpaceType spaceType, MeshSharedPtr mesh, Shapeset *shapeset)
virtual void copy(SpaceSharedPtr< Scalar > space, MeshSharedPtr new_mesh)
int get_num_dofs() const
Returns the number of basis functions contained in the space.
Determines the position on an element surface (edge in 2D and Face in 3D).
Node * en[H2D_MAX_NUMBER_EDGES]
edge node pointers
#define H2D_GET_H_ORDER(encoded_order)
Macros for combining quad horizontal and vertical encoded_orders.
void distribute_orders(MeshSharedPtr mesh, int *parents)
Sets polynomial orders to elements created by Mesh::regularize() using "parents". ...
virtual SpaceType get_space_type() const =0
virtual void handle_orders(SpaceSharedPtr< Scalar > ref_space)
::xsd::cxx::xml::dom::namespace_infomap< char > namespace_infomap
Namespace serialization information map.
unsigned short cnt
the number of items in the arrays idx, dof and coef
unsigned int seq_assigned
Tracking changes - mark call to assign_dofs().
virtual void resize_tables()
Updates internal node and element tables.
virtual SpaceSharedPtr< Scalar > create_ref_space(bool assign_dofs=true)
Methods that user calls to get the reference space pointer (has to be properly casted if necessary)...
NodeData * ndata
node data table
Should be exactly the same as is the count of enum ShapesetType.
unsigned bnd
1 = boundary node; 0 = inner node
void update_essential_bc_values()
virtual int assign_dofs(int first_dof=0)
Builds basis functions and assigns DOF numbers to them.
EssentialBCs< Scalar > * get_essential_bcs() const
Obtains an boundary conditions.
void unrefine_all_mesh_elements(bool keep_initial_refinements=true)
Class corresponding to the element_data schema type.
int esize
element data table size
ElementData * edata
element data table
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.
unsigned char get_num_components() const
Returns 2 if this is a vector shapeset, 1 otherwise.
int get_seq() const
Internal. Used by DiscreteProblem to detect changes in the space.
virtual void set_element_order_internal(int id, int order, int order_v=-1)
Represents a finite element space over a domain.
virtual void reset_dof_assignment()
Resets assignment of DOF to an unassigned state.
MeshSharedPtr mesh
FE mesh.
Class corresponding to the space schema type.
static SpaceSharedPtr< Scalar > load(const char *filename, MeshSharedPtr mesh, bool validate=false, EssentialBCs< Scalar > *essential_bcs=nullptr, Shapeset *shapeset=nullptr)
Loads a space from a file in XML format.
Element * sons[H2D_MAX_ELEMENT_SONS]
son elements (up to four)
void set_uniform_order_internal(int order, int marker)
ReferenceSpaceCreator(unsigned int order_increase=1)
bool active
0 = active, no sons; 1 = inactive (refined), has sons
bool is_up_to_date() const
Returns true if the space is ready for computation, false otherwise.
virtual void post_assign()
::std::auto_ptr< ::XMLSpace::space > space_(const ::std::string &uri,::xml_schema::flags f=0, const ::xml_schema::properties &p=::xml_schema::properties())
Parse a URI or a local file.
void set_mesh_seq(int seq)
Sets a (new) mesh seq, and mesh_seq.
EssentialBCs< Scalar > * essential_bcs
Boundary conditions.
void set_essential_bcs(EssentialBCs< Scalar > *essential_bcs)
Sets the boundary condition.
unsigned int seq
Tracking changes.
L2ShapesetLegendre L2Shapeset
This is the default shapeset typedef.
virtual SpaceType get_type() const =0
Internal. Return type of this space. See enum SpaceType.
int nsize
node data table size
1D quadrature points on the standard reference domain (-1,1)
virtual void set_element_order(int id, int order, int order_v=-1)
Node * vn[H2D_MAX_NUMBER_VERTICES]
vertex node pointers
void unrefine_all_mesh_elements_internal(bool keep_initial_refinements, bool only_unrefine_space_data)
virtual void update_constraints()