17 #include "space_hdiv.h"
20 #include "shapeset/shapeset_hd_all.h"
21 #include "boundary_conditions/essential_boundary_conditions.h"
26 template<
typename Scalar>
27 HdivSpace<Scalar>::HdivSpace() : Space<Scalar>()
31 template<
typename Scalar>
39 if(this->shapeset->
get_num_components() < 2)
throw Hermes::Exceptions::Exception(
"HdivSpace requires a vector shapeset.");
41 this->precalculate_projection_matrix(0, this->proj_mat, this->chol_p);
44 if(p_init < 0)
throw Hermes::Exceptions::Exception(
"P_INIT must be >= 0 in an Hdiv space.");
51 template<
typename Scalar>
53 :
Space<Scalar>(mesh, shapeset, essential_bcs)
55 init(shapeset, p_init);
58 template<
typename Scalar>
59 HdivSpace<Scalar>::HdivSpace(
const Mesh* mesh,
int p_init, Shapeset* shapeset)
60 : Space<Scalar>(mesh, shapeset, NULL)
62 init(shapeset, p_init);
65 template<
typename Scalar>
66 HdivSpace<Scalar>::~HdivSpace()
68 if(this->own_shapeset)
69 delete this->shapeset;
72 template<
typename Scalar>
77 this->precalculate_projection_matrix(0, this->proj_mat, this->chol_p);
80 template<
typename Scalar>
85 this->shapeset = shapeset;
86 this->own_shapeset =
false;
89 throw Hermes::Exceptions::Exception(
"Wrong shapeset type in HdivSpace<Scalar>::set_shapeset()");
92 template<
typename Scalar>
96 this->edge_functions_count = 0;
97 for_all_edge_nodes(en, this->mesh)
99 if(en->
ref > 1 || en->
bnd || this->mesh->peek_vertex_node(en->p1, en->
p2) != NULL)
101 int ndofs = this->get_edge_order_internal(en) + 1;
102 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;
129 this->ndata[en->
id].n = -1;
134 template<
typename Scalar>
135 void HdivSpace<Scalar>::assign_bubble_dofs()
138 this->bubble_functions_count = 0;
139 for_all_active_elements(e, this->mesh)
141 typename Space<Scalar>::ElementData* ed = &this->edata[e->id];
142 ed->bdof = this->next_dof;
143 ed->n = this->shapeset->get_num_bubbles(ed->order, e->get_mode());
144 this->next_dof += ed->n * this->stride;
145 this->bubble_functions_count += ed->n;
149 template<
typename Scalar>
150 void HdivSpace<Scalar>::get_boundary_assembly_list_internal(Element* e,
int surf_num, AsmList<Scalar>* al)
const
152 Node* en = e->en[surf_num];
153 typename Space<Scalar>::NodeData* nd = &this->ndata[en->id];
159 int ori = (e->vn[surf_num]->id < e->vn[e->next_vert(surf_num)]->id) ? 0 : 1;
160 for (
int j = 0, dof = nd->dof; j < nd->n; j++, dof += this->stride)
161 al->add_triplet(this->shapeset->get_edge_index(surf_num, ori, j, e->get_mode()), dof, 1.0);
165 for (
int j = 0; j < nd->n; j++)
166 al->add_triplet(this->shapeset->get_edge_index(surf_num, 0, j, e->get_mode()), -1, nd->edge_bc_proj[j]);
172 int ori = part < 0 ? 1 : 0;
173 if(part < 0) part ^= ~0;
175 nd = &this->ndata[nd->base->id];
176 for (
int j = 0, dof = nd->dof; j < nd->n; j++, dof += this->stride)
177 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 void HdivSpace<Scalar>::get_bubble_assembly_list(Element* e, AsmList<Scalar>* al)
const
184 typename Space<Scalar>::ElementData* ed = &this->edata[e->id];
187 int* indices = this->shapeset->get_bubble_indices(ed->order, e->get_mode());
188 for (
int i = 0, dof = ed->bdof; i < ed->n; i++, dof += this->stride)
189 al->add_triplet(*indices++, dof, 1.0);
194 template<
typename Scalar>
195 Scalar* HdivSpace<Scalar>::get_bc_projection(SurfPos* surf_pos,
int order, EssentialBoundaryCondition<Scalar> *bc)
198 Scalar* proj =
new Scalar[order + 1];
202 int mo = quad1d.get_max_order();
203 double2* pt = quad1d.get_points(mo);
205 Node* vn1 = this->mesh->get_node(surf_pos->v1);
206 Node* vn2 = this->mesh->get_node(surf_pos->v2);
207 double el = sqrt(sqr(vn1->x - vn2->x) + sqr(vn1->y - vn2->y));
208 el *= 0.5 * (surf_pos->hi - surf_pos->lo);
211 for (
int i = 0; i <= order; i++)
214 int ii = this->shapeset->get_edge_index(0, 0, i, surf_pos->base->get_mode());
215 for (
int j = 0; j < quad1d.get_num_points(mo); j++)
217 double t = (pt[j][0] + 1) * 0.5, s = 1.0 - t;
218 surf_pos->t = surf_pos->lo * s + surf_pos->hi * t;
221 EssentialBoundaryCondition<Scalar> *bc = this->essential_bcs->get_boundary_condition(this->mesh->boundary_markers_conversion.get_user_marker(surf_pos->marker).marker);
223 if(bc->get_value_type() == EssentialBoundaryCondition<Scalar>::BC_CONST)
225 rhs[i] += pt[j][1] * this->shapeset->get_fn_value(ii, pt[j][0], -1.0, 1, surf_pos->base->get_mode())
226 * bc->value_const * el;
229 else if(bc->get_value_type() == EssentialBoundaryCondition<Scalar>::BC_FUNCTION)
232 double x, y, n_x, n_y, t_x, t_y;
233 Nurbs* nurbs = surf_pos->base->is_curved() ? surf_pos->base->cm->nurbs[surf_pos->surf_num] : NULL;
234 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);
236 rhs[i] += pt[j][1] * this->shapeset->get_fn_value(ii, pt[j][0], -1.0, 1, surf_pos->base->get_mode())
237 * bc->value(x, y, n_x, n_y, t_x, t_y) * el;
243 cholsl(this->proj_mat, order + 1, this->chol_p, rhs, rhs);
245 for (
int i = 0; i < order + 1; i++)
251 template<
typename Scalar>
252 void HdivSpace<Scalar>::update_constrained_nodes(Element* e, EdgeInfo* ei0, EdgeInfo* ei1, EdgeInfo* ei2, EdgeInfo* ei3)
255 EdgeInfo* ei[4] = { ei0, ei1, ei2, ei3 };
256 typename Space<Scalar>::NodeData* nd;
261 for (
unsigned int i = 0; i < e->get_nvert(); i++)
265 nd = &this->ndata[e->en[i]->id];
266 nd->base = ei[i]->node;
267 nd->part = ei[i]->part;
268 if(ei[i]->ori) nd->part ^= ~0;
277 for (
unsigned int i = 0; i < e->get_nvert(); i++)
282 Node* mid_vn = this->get_mid_edge_vertex_node(e, i, j);
283 if(mid_vn != NULL && mid_vn->is_constrained_vertex())
285 Node* mid_en = this->mesh->peek_edge_node(e->vn[i]->id, e->vn[j]->id);
289 ei[i]->node = mid_en;
293 ei[i]->ori = (e->vn[i]->id < e->vn[j]->id) ? 0 : 1;
300 EdgeInfo half_ei_data[4][2];
301 EdgeInfo* half_ei[4][2];
302 for (
unsigned int i = 0; i < e->get_nvert(); i++)
306 half_ei[i][0] = half_ei[i][1] = NULL;
310 EdgeInfo* h0 = half_ei[i][0] = half_ei_data[i];
311 EdgeInfo* h1 = half_ei[i][1] = half_ei_data[i] + 1;
313 h0->node = h1->node = ei[i]->node;
314 h0->part = (ei[i]->part + 1) * 2;
315 h1->part = h0->part + 1;
316 h0->hi = h1->lo = (ei[i]->lo + ei[i]->hi) / 2;
319 h1->ori = h0->ori = ei[i]->ori;
326 update_constrained_nodes(e->sons[0], half_ei[0][0], NULL, half_ei[2][1], NULL);
327 update_constrained_nodes(e->sons[1], half_ei[0][1], half_ei[1][0], NULL, NULL);
328 update_constrained_nodes(e->sons[2], NULL, half_ei[1][1], half_ei[2][0], NULL);
329 update_constrained_nodes(e->sons[3], NULL, NULL, NULL, NULL);
331 else if(e->sons[2] == NULL)
333 update_constrained_nodes(e->sons[0], ei[0], half_ei[1][0], NULL, half_ei[3][1]);
334 update_constrained_nodes(e->sons[1], NULL, half_ei[1][1], ei[2], half_ei[3][0]);
336 else if(e->sons[0] == NULL)
338 update_constrained_nodes(e->sons[2], half_ei[0][0], NULL, half_ei[2][1], ei[3]);
339 update_constrained_nodes(e->sons[3], half_ei[0][1], ei[1], half_ei[2][0], NULL);
343 update_constrained_nodes(e->sons[0], half_ei[0][0], NULL, NULL, half_ei[3][1]);
344 update_constrained_nodes(e->sons[1], half_ei[0][1], half_ei[1][0], NULL, NULL);
345 update_constrained_nodes(e->sons[2], NULL, half_ei[1][1], half_ei[2][0], NULL);
346 update_constrained_nodes(e->sons[3], NULL, NULL, half_ei[2][1], half_ei[3][0]);
351 template<
typename Scalar>
355 for_all_base_elements(e, this->mesh)
356 update_constrained_nodes(e, NULL, NULL, NULL, NULL);