17 #include "space_hcurl.h"
20 #include "shapeset_hc_all.h"
21 #include "essential_boundary_conditions.h"
27 template<
typename Scalar>
28 HcurlSpace<Scalar>::HcurlSpace() : Space<Scalar>()
32 template<
typename Scalar>
38 this->own_shapeset =
true;
40 if(this->shapeset->
get_num_components() < 2)
throw Hermes::Exceptions::Exception(
"HcurlSpace requires a vector shapeset.");
42 this->precalculate_projection_matrix(0, this->proj_mat, this->chol_p);
45 if(p_init < 0)
throw Hermes::Exceptions::Exception(
"P_INIT must be >= 0 in an Hcurl space.");
46 else this->set_uniform_order_internal(p_init, HERMES_ANY_INT);
52 template<
typename Scalar>
54 :
Space<Scalar>(mesh, shapeset, essential_bcs)
56 init(shapeset, p_init);
59 template<
typename Scalar>
60 HcurlSpace<Scalar>::HcurlSpace(
const Mesh* mesh,
int p_init, Shapeset* shapeset)
61 : Space<Scalar>(mesh, shapeset, NULL)
63 init(shapeset, p_init);
66 template<
typename Scalar>
67 HcurlSpace<Scalar>::~HcurlSpace()
69 if(this->own_shapeset)
70 delete this->shapeset;
73 template<
typename Scalar>
78 this->precalculate_projection_matrix(0, this->proj_mat, this->chol_p);
81 template<
typename Scalar>
86 this->shapeset = shapeset;
87 this->own_shapeset =
false;
90 throw Hermes::Exceptions::Exception(
"Wrong shapeset type in HcurlSpace<Scalar>::set_shapeset()");
93 template<
typename Scalar>
97 this->edge_functions_count = 0;
98 for_all_edge_nodes(en, this->mesh)
100 if(en->
ref > 1 || en->
bnd || this->mesh->peek_vertex_node(en->p1, en->
p2) != NULL)
102 int ndofs = this->get_edge_order_internal(en) + 1;
103 this->ndata[en->
id].n = ndofs;
105 if(this->essential_bcs != NULL)
106 if(this->essential_bcs->get_boundary_condition(this->mesh->boundary_markers_conversion.get_user_marker(en->
marker).marker) != NULL)
107 this->ndata[en->
id].dof = this->H2D_CONSTRAINED_DOF;
110 this->ndata[en->
id].dof = this->next_dof;
111 this->next_dof += ndofs * this->stride;
112 this->edge_functions_count += ndofs;
116 this->ndata[en->
id].dof = this->next_dof;
117 this->next_dof += ndofs * this->stride;
118 this->edge_functions_count += ndofs;
122 this->ndata[en->
id].dof = this->next_dof;
123 this->next_dof += ndofs * this->stride;
124 this->edge_functions_count += ndofs;
128 this->ndata[en->
id].n = -1;
132 template<
typename Scalar>
133 void HcurlSpace<Scalar>::assign_bubble_dofs()
136 this->bubble_functions_count = 0;
137 for_all_active_elements(e, this->mesh)
139 typename Space<Scalar>::ElementData* ed = &this->edata[e->id];
140 ed->bdof = this->next_dof;
141 ed->n = this->shapeset->get_num_bubbles(ed->order, e->get_mode());
142 this->next_dof += ed->n * this->stride;
143 this->bubble_functions_count += ed->n;
147 template<
typename Scalar>
148 void HcurlSpace<Scalar>::get_boundary_assembly_list_internal(Element* e,
int surf_num, AsmList<Scalar>* al)
const
150 Node* en = e->en[surf_num];
151 typename Space<Scalar>::NodeData* nd = &this->ndata[en->id];
157 int ori = (e->vn[surf_num]->id < e->vn[e->next_vert(surf_num)]->id) ? 0 : 1;
158 for (
int j = 0, dof = nd->dof; j < nd->n; j++, dof += this->stride)
159 al->add_triplet(this->shapeset->get_edge_index(surf_num, ori, j, e->get_mode()), dof, 1.0);
163 for (
int j = 0; j < nd->n; j++)
164 al->add_triplet(this->shapeset->get_edge_index(surf_num, 0, j, e->get_mode()), -1, nd->edge_bc_proj[j]);
170 int ori = part < 0 ? 1 : 0;
171 if(part < 0) part ^= ~0;
173 nd = &this->ndata[nd->base->id];
174 for (
int j = 0, dof = nd->dof; j < nd->n; j++, dof += this->stride)
175 al->add_triplet(this->shapeset->get_constrained_edge_index(surf_num, j, ori, part, e->get_mode()), dof, 1.0);
181 template<
typename Scalar>
182 Scalar* HcurlSpace<Scalar>::get_bc_projection(SurfPos* surf_pos,
int order, EssentialBoundaryCondition<Scalar> *bc)
185 Scalar* proj =
new Scalar[order + 1];
189 int mo = quad1d.get_max_order();
190 double2* pt = quad1d.get_points(mo);
192 Node* vn1 = this->mesh->get_node(surf_pos->v1);
193 Node* vn2 = this->mesh->get_node(surf_pos->v2);
194 double el = sqrt(sqr(vn1->x - vn2->x) + sqr(vn1->y - vn2->y));
195 el *= 0.5 * (surf_pos->hi - surf_pos->lo);
198 for (
int i = 0; i <= order; i++)
201 int ii = this->shapeset->get_edge_index(0, 0, i, surf_pos->base->get_mode());
202 for (
int j = 0; j < quad1d.get_num_points(mo); j++)
204 double t = (pt[j][0] + 1) * 0.5, s = 1.0 - t;
205 surf_pos->t = surf_pos->lo * s + surf_pos->hi * t;
207 if(bc->get_value_type() == EssentialBoundaryCondition<Scalar>::BC_CONST)
209 rhs[i] += pt[j][1] * this->shapeset->get_fn_value(ii, pt[j][0], -1.0, 0, surf_pos->base->get_mode())
210 * bc->value_const * el;
213 else if(bc->get_value_type() == EssentialBoundaryCondition<Scalar>::BC_FUNCTION)
216 double x, y, n_x, n_y, t_x, t_y;
217 Nurbs* nurbs = surf_pos->base->is_curved() ? surf_pos->base->cm->nurbs[surf_pos->surf_num] : NULL;
218 CurvMap::nurbs_edge(surf_pos->base, nurbs, surf_pos->surf_num, 2.0*surf_pos->t - 1.0, x, y, n_x, n_y, t_x, t_y);
220 rhs[i] += pt[j][1] * this->shapeset->get_fn_value(ii, pt[j][0], -1.0, 0, surf_pos->base->get_mode())
221 * bc->value(x, y, n_x, n_y, t_x, t_y) * el;
227 cholsl(this->proj_mat, order + 1, this->chol_p, rhs, rhs);
231 template<
typename Scalar>
232 void HcurlSpace<Scalar>::update_constrained_nodes(Element* e, EdgeInfo* ei0, EdgeInfo* ei1, EdgeInfo* ei2, EdgeInfo* ei3)
235 EdgeInfo* ei[4] = { ei0, ei1, ei2, ei3 };
236 typename Space<Scalar>::NodeData* nd;
241 for (
unsigned int i = 0; i < e->get_nvert(); i++)
245 nd = &this->ndata[e->en[i]->id];
246 nd->base = ei[i]->node;
247 nd->part = ei[i]->part;
248 if(ei[i]->ori) nd->part ^= ~0;
257 for (
unsigned int i = 0; i < e->get_nvert(); i++)
262 Node* mid_vn = this->get_mid_edge_vertex_node(e, i, j);
263 if(mid_vn != NULL && mid_vn->is_constrained_vertex())
265 Node* mid_en = this->mesh->peek_edge_node(e->vn[i]->id, e->vn[j]->id);
269 ei[i]->node = mid_en;
273 ei[i]->ori = (e->vn[i]->id < e->vn[j]->id) ? 0 : 1;
280 EdgeInfo half_ei_data[4][2];
281 EdgeInfo* half_ei[4][2];
282 for (
unsigned int i = 0; i < e->get_nvert(); i++)
286 half_ei[i][0] = half_ei[i][1] = NULL;
290 EdgeInfo* h0 = half_ei[i][0] = half_ei_data[i];
291 EdgeInfo* h1 = half_ei[i][1] = half_ei_data[i] + 1;
293 h0->node = h1->node = ei[i]->node;
294 h0->part = (ei[i]->part + 1) * 2;
295 h1->part = h0->part + 1;
296 h0->hi = h1->lo = (ei[i]->lo + ei[i]->hi) / 2;
299 h1->ori = h0->ori = ei[i]->ori;
306 update_constrained_nodes(e->sons[0], half_ei[0][0], NULL, half_ei[2][1], NULL);
307 update_constrained_nodes(e->sons[1], half_ei[0][1], half_ei[1][0], NULL, NULL);
308 update_constrained_nodes(e->sons[2], NULL, half_ei[1][1], half_ei[2][0], NULL);
309 update_constrained_nodes(e->sons[3], NULL, NULL, NULL, NULL);
311 else if(e->sons[2] == NULL)
313 update_constrained_nodes(e->sons[0], ei[0], half_ei[1][0], NULL, half_ei[3][1]);
314 update_constrained_nodes(e->sons[1], NULL, half_ei[1][1], ei[2], half_ei[3][0]);
316 else if(e->sons[0] == NULL)
318 update_constrained_nodes(e->sons[2], half_ei[0][0], NULL, half_ei[2][1], ei[3]);
319 update_constrained_nodes(e->sons[3], half_ei[0][1], ei[1], half_ei[2][0], NULL);
323 update_constrained_nodes(e->sons[0], half_ei[0][0], NULL, NULL, half_ei[3][1]);
324 update_constrained_nodes(e->sons[1], half_ei[0][1], half_ei[1][0], NULL, NULL);
325 update_constrained_nodes(e->sons[2], NULL, half_ei[1][1], half_ei[2][0], NULL);
326 update_constrained_nodes(e->sons[3], NULL, NULL, half_ei[2][1], half_ei[3][0]);
331 template<
typename Scalar>
335 for_all_base_elements(e, this->mesh)
336 update_constrained_nodes(e, NULL, NULL, NULL, NULL);