Hermes2D  3.0
space_hdiv.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 "global.h"
17 #include "space_hdiv.h"
18 #include "matrix.h"
19 #include "quad_all.h"
20 #include "shapeset/shapeset_hd_all.h"
21 #include "boundary_conditions/essential_boundary_conditions.h"
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  template<typename Scalar>
27  HdivSpace<Scalar>::HdivSpace() : Space<Scalar>()
28  {
29  }
30 
31  template<typename Scalar>
32  HdivSpace<Scalar>::HdivSpace(MeshSharedPtr mesh, EssentialBCs<Scalar>* essential_bcs, int p_init, Shapeset* shapeset)
33  : Space<Scalar>(mesh, shapeset, essential_bcs)
34  {
35  init(shapeset, p_init);
36  }
37 
38  template<typename Scalar>
39  HdivSpace<Scalar>::HdivSpace(MeshSharedPtr mesh, int p_init, Shapeset* shapeset)
40  : Space<Scalar>(mesh, shapeset, nullptr)
41  {
42  init(shapeset, p_init);
43  }
44 
45  template<typename Scalar>
46  void HdivSpace<Scalar>::init(Shapeset* shapeset, int p_init, bool assign_dofs_init)
47  {
48  if (shapeset == nullptr)
49  {
50  this->shapeset = new HdivShapeset;
51  this->own_shapeset = true;
52  }
53  if (this->shapeset->get_num_components() < 2) throw Hermes::Exceptions::Exception("HdivSpace requires a vector shapeset.");
54 
55  this->precalculate_projection_matrix(0, this->proj_mat, this->chol_p);
56 
57  // enumerate basis functions
58  if (assign_dofs_init)
59  {
60  // set uniform poly order in elements
61  if (p_init < 0)
62  throw Hermes::Exceptions::Exception("P_INIT must be >= 0 in an Hdiv space.");
63  else
64  this->set_uniform_order_internal(p_init, HERMES_ANY_INT);
65 
66  this->assign_dofs();
67  }
68  }
69 
70  template<typename Scalar>
72  {
73  }
74 
75  template<typename Scalar>
76  void HdivSpace<Scalar>::copy(SpaceSharedPtr<Scalar> space, MeshSharedPtr new_mesh)
77  {
78  this->set_shapeset(space->get_shapeset(), true);
79 
80  this->precalculate_projection_matrix(0, this->proj_mat, this->chol_p);
81 
82  Space<Scalar>::copy(space, new_mesh);
83  }
84 
85  template<typename Scalar>
86  void HdivSpace<Scalar>::set_shapeset(Shapeset *shapeset, bool clone)
87  {
88  if (!(shapeset->get_id() < 30 && shapeset->get_id() > 19))
89  throw Hermes::Exceptions::Exception("Wrong shapeset type in HdivSpace<Scalar>::set_shapeset()");
90 
91  if (clone)
92  {
93  if (this->own_shapeset)
94  delete this->shapeset;
95 
96  this->shapeset = shapeset->clone();
97  }
98  else
99  {
100  this->shapeset = shapeset;
101  this->own_shapeset = false;
102  }
103  }
104 
105  template<typename Scalar>
107  {
108  Node* en;
109  this->edge_functions_count = 0;
110  for_all_edge_nodes(en, this->mesh)
111  {
112  if (en->ref > 1 || en->bnd || this->mesh->peek_vertex_node(en->p1, en->p2) != nullptr)
113  {
114  int ndofs = this->get_edge_order_internal(en) + 1;
115  this->ndata[en->id].n = ndofs;
116 
117  if (en->bnd)
118  if (this->essential_bcs != nullptr)
119  if (this->essential_bcs->get_boundary_condition(this->mesh->boundary_markers_conversion.get_user_marker(en->marker).marker) != nullptr)
120  this->ndata[en->id].dof = this->H2D_CONSTRAINED_DOF;
121  else
122  {
123  this->ndata[en->id].dof = this->next_dof;
124  this->next_dof += ndofs;
125  this->edge_functions_count += ndofs;
126  }
127  else
128  {
129  this->ndata[en->id].dof = this->next_dof;
130  this->next_dof += ndofs;
131  this->edge_functions_count += ndofs;
132  }
133  else
134  {
135  this->ndata[en->id].dof = this->next_dof;
136  this->next_dof += ndofs;
137  this->edge_functions_count += ndofs;
138  }
139  }
140  else
141  {
142  this->ndata[en->id].n = -1;
143  }
144  }
145  }
146 
147  template<typename Scalar>
148  void HdivSpace<Scalar>::assign_bubble_dofs()
149  {
150  Element* e;
151  this->bubble_functions_count = 0;
152  for_all_active_elements(e, this->mesh)
153  {
154  typename Space<Scalar>::ElementData* ed = &this->edata[e->id];
155  ed->bdof = this->next_dof;
156  ed->n = this->shapeset->get_num_bubbles(ed->order, e->get_mode());
157  this->next_dof += ed->n;
158  this->bubble_functions_count += ed->n;
159  }
160  }
161 
162  template<typename Scalar>
163  void HdivSpace<Scalar>::get_boundary_assembly_list_internal(Element* e, int surf_num, AsmList<Scalar>* al) const
164  {
165  Node* en = e->en[surf_num];
166  typename Space<Scalar>::NodeData* nd = &this->ndata[en->id];
167 
168  if (nd->n >= 0) // unconstrained
169  {
170  if (nd->dof >= 0)
171  {
172  int ori = (e->vn[surf_num]->id < e->vn[e->next_vert(surf_num)]->id) ? 0 : 1;
173  for (int j = 0, dof = nd->dof; j < nd->n; j++, dof++)
174  al->add_triplet(this->shapeset->get_edge_index(surf_num, ori, j, e->get_mode()), dof, 1.0);
175  }
176  else
177  {
178  for (int j = 0; j < nd->n; j++)
179  al->add_triplet(this->shapeset->get_edge_index(surf_num, 0, j, e->get_mode()), -1, nd->edge_bc_proj[j]);
180  }
181  }
182  else // constrained
183  {
184  int part = nd->part;
185  int ori = part < 0 ? 1 : 0;
186  if (part < 0) part ^= ~0;
187 
188  // ccc
189  nd = &this->ndata[nd->base->id];
190  for (int j = 0, dof = nd->dof; j < nd->n; j++, dof++)
191  al->add_triplet(this->shapeset->get_constrained_edge_index(surf_num, j, ori, part, e->get_mode()), dof, 1.0);
192  }
193  }
194 
195  template<typename Scalar>
196  void HdivSpace<Scalar>::get_bubble_assembly_list(Element* e, AsmList<Scalar>* al) const
197  {
198  typename Space<Scalar>::ElementData* ed = &this->edata[e->id];
199  if (!ed->n) return;
200 
201  short* indices = this->shapeset->get_bubble_indices(ed->order, e->get_mode());
202  for (unsigned short i = 0, dof = ed->bdof; i < ed->n; i++, dof++)
203  al->add_triplet(*indices++, dof, 1.0);
204  }
205 
207 
208  template<typename Scalar>
209  Scalar* HdivSpace<Scalar>::get_bc_projection(SurfPos* surf_pos, int order, EssentialBoundaryCondition<Scalar> *bc)
210  {
211  assert(order >= 0);
212  Scalar* proj = malloc_with_check<HdivSpace<Scalar>, Scalar>(order + 1, this);
213 
214  Quad1DStd quad1d;
215  Scalar* rhs = proj;
216  int mo = quad1d.get_max_order();
217  double2* pt = quad1d.get_points(mo);
218 
219  Node* vn1 = this->mesh->get_node(surf_pos->v1);
220  Node* vn2 = this->mesh->get_node(surf_pos->v2);
221  double el = sqrt(sqr(vn1->x - vn2->x) + sqr(vn1->y - vn2->y));
222  el *= 0.5 * (surf_pos->hi - surf_pos->lo);
223 
224  // get boundary values at integration points, construct rhs
225  for (int i = 0; i <= order; i++)
226  {
227  rhs[i] = 0.0;
228  int ii = this->shapeset->get_edge_index(0, 0, i, surf_pos->base->get_mode());
229  for (int j = 0; j < quad1d.get_num_points(mo); j++)
230  {
231  double t = (pt[j][0] + 1) * 0.5, s = 1.0 - t;
232  surf_pos->t = surf_pos->lo * s + surf_pos->hi * t;
233 
234  // If the BC on this part of the boundary is constant.
235  EssentialBoundaryCondition<Scalar> *bc = this->essential_bcs->get_boundary_condition(this->mesh->boundary_markers_conversion.get_user_marker(surf_pos->marker).marker);
236 
237  if (bc->get_value_type() == BC_CONST)
238  {
239  rhs[i] += pt[j][1] * this->shapeset->get_fn_value(ii, pt[j][0], -1.0, 1, surf_pos->base->get_mode())
240  * bc->value_const * el;
241  }
242  // If the BC is not constant.
243  else if (bc->get_value_type() == BC_FUNCTION)
244  {
245  // Find out the (x, y) coordinate.
246  double x, y;
247  Curve* curve = surf_pos->base->is_curved() ? surf_pos->base->cm->curves[surf_pos->surf_num] : nullptr;
248  CurvMap::nurbs_edge(surf_pos->base, curve, surf_pos->surf_num, 2.0*surf_pos->t - 1.0, x, y);
249  // Calculate.
250  rhs[i] += pt[j][1] * this->shapeset->get_fn_value(ii, pt[j][0], -1.0, 1, surf_pos->base->get_mode())
251  * bc->value(x, y) * el;
252  }
253  }
254  }
255 
256  // solve the system using a precalculated Cholesky decomposed projection matrix
257  cholsl(this->proj_mat, order + 1, this->chol_p, rhs, rhs);
258 
259  for (int i = 0; i < order + 1; i++)
260  proj[i] = 0.0;
261 
262  return proj;
263  }
264 
265  template<typename Scalar>
266  void HdivSpace<Scalar>::update_constrained_nodes(Element* e, typename Space<Scalar>::EdgeInfo* ei0, typename Space<Scalar>::EdgeInfo* ei1, typename Space<Scalar>::EdgeInfo* ei2, typename Space<Scalar>::EdgeInfo* ei3)
267  {
268  int j;
269  typename Space<Scalar>::EdgeInfo* ei[4] = { ei0, ei1, ei2, ei3 };
270  typename Space<Scalar>::NodeData* nd;
271 
272  // on non-refined elements all we have to do is update edge nodes lying on constrained edges
273  if (e->active)
274  {
275  for (unsigned char i = 0; i < e->get_nvert(); i++)
276  {
277  if (ei[i] != nullptr)
278  {
279  nd = &this->ndata[e->en[i]->id];
280  nd->base = ei[i]->node;
281  nd->part = ei[i]->part;
282  if (ei[i]->ori) nd->part ^= ~0;
283  }
284  }
285  }
286  // the element has sons - update mid-edge constrained vertex nodes
287  else
288  {
289  // create new_ edge infos where we don't have them yet
290  typename Space<Scalar>::EdgeInfo ei_data[4];
291  for (unsigned char i = 0; i < e->get_nvert(); i++)
292  {
293  if (ei[i] == nullptr)
294  {
295  j = e->next_vert(i);
296  Node* mid_vn = this->get_mid_edge_vertex_node(e, i, j);
297  if (mid_vn != nullptr && mid_vn->is_constrained_vertex())
298  {
299  Node* mid_en = this->mesh->peek_edge_node(e->vn[i]->id, e->vn[j]->id);
300  if (mid_en != nullptr)
301  {
302  ei[i] = ei_data + i;
303  ei[i]->node = mid_en;
304  ei[i]->part = -1;
305  ei[i]->lo = -1.0;
306  ei[i]->hi = 1.0;
307  ei[i]->ori = (e->vn[i]->id < e->vn[j]->id) ? 0 : 1;
308  }
309  }
310  }
311  }
312 
313  // create edge infos for half-edges
314  typename Space<Scalar>::EdgeInfo half_ei_data[4][2];
315  typename Space<Scalar>::EdgeInfo* half_ei[4][2];
316  for (unsigned char i = 0; i < e->get_nvert(); i++)
317  {
318  if (ei[i] == nullptr)
319  {
320  half_ei[i][0] = half_ei[i][1] = nullptr;
321  }
322  else
323  {
324  typename Space<Scalar>::EdgeInfo* h0 = half_ei[i][0] = half_ei_data[i];
325  typename Space<Scalar>::EdgeInfo* h1 = half_ei[i][1] = half_ei_data[i] + 1;
326 
327  h0->node = h1->node = ei[i]->node;
328  h0->part = (ei[i]->part + 1) * 2;
329  h1->part = h0->part + 1;
330  h0->hi = h1->lo = (ei[i]->lo + ei[i]->hi) / 2;
331  h0->lo = ei[i]->lo;
332  h1->hi = ei[i]->hi;
333  h1->ori = h0->ori = ei[i]->ori;
334  }
335  }
336 
337  // recur to sons
338  if (e->is_triangle())
339  {
340  update_constrained_nodes(e->sons[0], half_ei[0][0], nullptr, half_ei[2][1], nullptr);
341  update_constrained_nodes(e->sons[1], half_ei[0][1], half_ei[1][0], nullptr, nullptr);
342  update_constrained_nodes(e->sons[2], nullptr, half_ei[1][1], half_ei[2][0], nullptr);
343  update_constrained_nodes(e->sons[3], nullptr, nullptr, nullptr, nullptr);
344  }
345  else if (e->sons[2] == nullptr) // 'horizontally' split quad
346  {
347  update_constrained_nodes(e->sons[0], ei[0], half_ei[1][0], nullptr, half_ei[3][1]);
348  update_constrained_nodes(e->sons[1], nullptr, half_ei[1][1], ei[2], half_ei[3][0]);
349  }
350  else if (e->sons[0] == nullptr) // 'vertically' split quad
351  {
352  update_constrained_nodes(e->sons[2], half_ei[0][0], nullptr, half_ei[2][1], ei[3]);
353  update_constrained_nodes(e->sons[3], half_ei[0][1], ei[1], half_ei[2][0], nullptr);
354  }
355  else // fully split quad
356  {
357  update_constrained_nodes(e->sons[0], half_ei[0][0], nullptr, nullptr, half_ei[3][1]);
358  update_constrained_nodes(e->sons[1], half_ei[0][1], half_ei[1][0], nullptr, nullptr);
359  update_constrained_nodes(e->sons[2], nullptr, half_ei[1][1], half_ei[2][0], nullptr);
360  update_constrained_nodes(e->sons[3], nullptr, nullptr, half_ei[2][1], half_ei[3][0]);
361  }
362  }
363  }
364 
365  template<typename Scalar>
367  {
368  Element* e;
369  for_all_base_elements(e, this->mesh)
370  update_constrained_nodes(e, nullptr, nullptr, nullptr, nullptr);
371  }
372 
373  template HERMES_API class HdivSpace < double > ;
374  template HERMES_API class HdivSpace < std::complex<double> > ;
375  }
376 }
int id
node id number
Definition: element.h:48
Definition: adapt.h:24
Stores one element of a mesh.
Definition: element.h:107
HdivShapesetLegendre HdivShapeset
This is the default Hdiv shapeset typedef.
virtual void copy(SpaceSharedPtr< Scalar > space, MeshSharedPtr new_mesh)
Copy from Space instance 'space'.
Definition: space_hdiv.cpp:76
void init()
Common code for constructors.
Definition: space.cpp:85
Used to pass the instances of Space around.
Definition: space.h:34
Common definitions for Hermes2D.
Stores one node of a mesh.
Definition: element.h:45
virtual void copy(SpaceSharedPtr< Scalar > space, MeshSharedPtr new_mesh)
Definition: space.cpp:243
virtual unsigned char get_id() const =0
Should be exactly the same as is the count of enum ShapesetType.
Definition: shapeset.h:95
int p1
parent id numbers
Definition: element.h:75
unsigned bnd
1 = boundary node; 0 = inner node
Definition: element.h:54
virtual void set_shapeset(Shapeset *shapeset, bool clone=false)
Sets the shapeset.
Definition: space_hdiv.cpp:86
unsigned char get_num_components() const
Returns 2 if this is a vector shapeset, 1 otherwise.
Definition: shapeset.cpp:195
unsigned ref
the number of elements using the node
Definition: element.h:50
int marker
edge marker
Definition: element.h:68
virtual void update_constraints()
Definition: space_hdiv.cpp:366