Doxygen API reference documentation for ideal.II
slab_dof_tools.hh
1 /*
2  * slab_dof_tools.hh
3  *
4  * Created on: Nov 22, 2022
5  * Author: thiele
6  */
7 
8 #ifndef INCLUDE_IDEAL_II_DOFS_SLAB_DOF_TOOLS_HH_
9 #define INCLUDE_IDEAL_II_DOFS_SLAB_DOF_TOOLS_HH_
10 
11 #include <ideal.II/dofs/slab_dof_handler.hh>
12 
13 #include <ideal.II/fe/fe_dg.hh>
14 
15 #include <deal.II/base/types.h>
16 
17 #include <deal.II/dofs/dof_tools.h>
18 
19 #include <deal.II/lac/dynamic_sparsity_pattern.h>
20 
21 namespace idealii::slab::DoFTools
22 {
23  namespace internal
24  {
25  template <int dim>
26  void
27  upwind_temporal_pattern(DoFHandler<dim> &dof,
28  dealii::DynamicSparsityPattern &time_dsp)
29  {
30  dealii::DoFTools::make_sparsity_pattern(*dof.temporal(), time_dsp);
31  if (dof.fe_support_type() ==
32  spacetime::DG_FiniteElement<dim>::support_type::Lobatto)
33  {
34  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
35  ii < dof.n_dofs_time();
36  ii += dof.dofs_per_cell_time())
37  {
38  time_dsp.add(ii, ii - 1);
39  }
40  }
41  else if (dof.fe_support_type() ==
42  spacetime::DG_FiniteElement<dim>::support_type::RadauLeft)
43  {
44  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
45  ii < dof.n_dofs_time();
46  ii += dof.dofs_per_cell_time())
47  {
48  for (dealii::types::global_dof_index l = 0;
49  l < dof.dofs_per_cell_time();
50  l++)
51  {
52  time_dsp.add(ii, ii - l - 1);
53  }
54  }
55  }
56  else if (dof.fe_support_type() ==
57  spacetime::DG_FiniteElement<dim>::support_type::RadauRight)
58  {
59  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
60  ii < dof.n_dofs_time();
61  ii += dof.dofs_per_cell_time())
62  {
63  for (dealii::types::global_dof_index k = 0;
64  k < dof.dofs_per_cell_time();
65  k++)
66  {
67  time_dsp.add(ii + k, ii - 1);
68  }
69  }
70  }
71  else
72  {
73  // go over first DoF of each cell
74  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
75  ii < dof.n_dofs_time();
76  ii += dof.dofs_per_cell_time())
77  {
78  // row offset
79  for (dealii::types::global_dof_index k = 0;
80  k < dof.dofs_per_cell_time();
81  k++)
82  {
83  for (dealii::types::global_dof_index l = 0;
84  l < dof.dofs_per_cell_time();
85  l++)
86  {
87  time_dsp.add(ii + k, ii - l - 1);
88  }
89  }
90  }
91  }
92  }
93  template <int dim>
94  void
95  downwind_temporal_pattern(DoFHandler<dim> &dof,
96  dealii::DynamicSparsityPattern &time_dsp)
97  {
98  dealii::DoFTools::make_sparsity_pattern(*dof.temporal(), time_dsp);
99  if (dof.fe_support_type() ==
100  spacetime::DG_FiniteElement<dim>::support_type::Lobatto)
101  {
102  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
103  ii < dof.n_dofs_time();
104  ii += dof.dofs_per_cell_time())
105  {
106  time_dsp.add(ii - 1, ii);
107  }
108  }
109  else if (dof.fe_support_type() ==
110  spacetime::DG_FiniteElement<dim>::support_type::RadauRight)
111  {
112  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
113  ii < dof.n_dofs_time();
114  ii += dof.dofs_per_cell_time())
115  {
116  for (dealii::types::global_dof_index k = 0;
117  k < dof.dofs_per_cell_time();
118  k++)
119  {
120  time_dsp.add(ii - 1, ii + k);
121  }
122  }
123  }
124  else if (dof.fe_support_type() ==
125  spacetime::DG_FiniteElement<dim>::support_type::RadauLeft)
126  {
127  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
128  ii < dof.n_dofs_time();
129  ii += dof.dofs_per_cell_time())
130  {
131  for (dealii::types::global_dof_index l = 0;
132  l < dof.dofs_per_cell_time();
133  l++)
134  {
135  time_dsp.add(ii - l - 1, ii);
136  }
137  }
138  }
139  else
140  {
141  // go over first DoF of each cell
142  for (dealii::types::global_dof_index ii = dof.dofs_per_cell_time();
143  ii < dof.n_dofs_time();
144  ii += dof.dofs_per_cell_time())
145  {
146  // row offset
147  for (dealii::types::global_dof_index k = 0;
148  k < dof.dofs_per_cell_time();
149  k++)
150  {
151  for (dealii::types::global_dof_index l = 0;
152  l < dof.dofs_per_cell_time();
153  l++)
154  {
155  time_dsp.add(ii - l - 1, ii + k);
156  }
157  }
158  }
159  }
160  }
161  } // namespace internal
162 
195  template <int dim>
196  void
197  make_upwind_sparsity_pattern(
198  DoFHandler<dim> &dof,
199  dealii::DynamicSparsityPattern &st_dsp,
200  std::shared_ptr<dealii::AffineConstraints<double>> space_constraints =
201  std::make_shared<dealii::AffineConstraints<double>>(),
202  const bool keep_constrained_dofs = true)
203  {
204  dealii::DynamicSparsityPattern space_dsp(dof.n_dofs_space());
205  dealii::DoFTools::make_sparsity_pattern(*dof.spatial(),
206  space_dsp,
207  *space_constraints,
208  keep_constrained_dofs);
209 
210  dealii::DynamicSparsityPattern time_dsp(dof.n_dofs_time());
211  internal::upwind_temporal_pattern(dof, time_dsp);
212 
213  for (auto &space_entry : space_dsp)
214  {
215  for (auto &time_entry : time_dsp)
216  {
217  st_dsp.add(time_entry.row() * dof.n_dofs_space() +
218  space_entry.row(), // test function
219  time_entry.column() * dof.n_dofs_space() +
220  space_entry.column() // trial function
221  );
222  }
223  }
224  }
238  template <int dim>
239  void
240  make_upwind_sparsity_pattern(
241  DoFHandler<dim> &dof,
242  const dealii::Table<2, dealii::DoFTools::Coupling> &space_couplings,
243  dealii::DynamicSparsityPattern &st_dsp,
244  std::shared_ptr<dealii::AffineConstraints<double>> space_constraints =
245  std::make_shared<dealii::AffineConstraints<double>>(),
246  const bool keep_constrained_dofs = true)
247  {
248  dealii::DynamicSparsityPattern space_dsp(dof.n_dofs_space());
249  dealii::DoFTools::make_sparsity_pattern(*dof.spatial(),
250  space_couplings,
251  space_dsp,
252  *space_constraints,
253  keep_constrained_dofs);
254 
255  dealii::DynamicSparsityPattern time_dsp(dof.n_dofs_time());
256  internal::upwind_temporal_pattern(dof, time_dsp);
257 
258  for (auto &space_entry : space_dsp)
259  {
260  for (auto &time_entry : time_dsp)
261  {
262  st_dsp.add(time_entry.row() * dof.n_dofs_space() +
263  space_entry.row(), // test function
264  time_entry.column() * dof.n_dofs_space() +
265  space_entry.column() // trial function
266  );
267  }
268  }
269  }
270 
303  template <int dim>
304  void
305  make_downwind_sparsity_pattern(
306  DoFHandler<dim> &dof,
307  dealii::DynamicSparsityPattern &st_dsp,
308  std::shared_ptr<dealii::AffineConstraints<double>> space_constraints =
309  std::make_shared<dealii::AffineConstraints<double>>(),
310  const bool keep_constrained_dofs = true)
311  {
312  dealii::DynamicSparsityPattern space_dsp(dof.n_dofs_space());
313  dealii::DoFTools::make_sparsity_pattern(*dof.spatial(),
314  space_dsp,
315  *space_constraints,
316  keep_constrained_dofs);
317 
318  dealii::DynamicSparsityPattern time_dsp(dof.n_dofs_time());
319  internal::downwind_temporal_pattern(dof, time_dsp);
320 
321  for (auto &space_entry : space_dsp)
322  {
323  for (auto &time_entry : time_dsp)
324  {
325  st_dsp.add(time_entry.row() * dof.n_dofs_space() +
326  space_entry.row(), // test function
327  time_entry.column() * dof.n_dofs_space() +
328  space_entry.column() // trial function
329  );
330  }
331  }
332  }
346  template <int dim>
347  void
348  make_downwind_sparsity_pattern(
349  DoFHandler<dim> &dof,
350  const dealii::Table<2, dealii::DoFTools::Coupling> &space_couplings,
351  dealii::DynamicSparsityPattern &st_dsp,
352  std::shared_ptr<dealii::AffineConstraints<double>> space_constraints =
353  std::make_shared<dealii::AffineConstraints<double>>(),
354  const bool keep_constrained_dofs = true)
355  {
356  dealii::DynamicSparsityPattern space_dsp(dof.n_dofs_space());
357  dealii::DoFTools::make_sparsity_pattern(*dof.spatial(),
358  space_couplings,
359  space_dsp,
360  *space_constraints,
361  keep_constrained_dofs);
362 
363  dealii::DynamicSparsityPattern time_dsp(dof.n_dofs_time());
364  internal::downwind_temporal_pattern(dof, time_dsp);
365 
366  for (auto &space_entry : space_dsp)
367  {
368  for (auto &time_entry : time_dsp)
369  {
370  st_dsp.add(time_entry.row() * dof.n_dofs_space() +
371  space_entry.row(), // test function
372  time_entry.column() * dof.n_dofs_space() +
373  space_entry.column() // trial function
374  );
375  }
376  }
377  }
378 
391  template <int dim, typename Number>
392  void
393  make_hanging_node_constraints(
394  idealii::slab::DoFHandler<dim> &dof_handler,
395  std::shared_ptr<dealii::AffineConstraints<Number>> spacetime_constraints)
396  {
397  auto space_constraints =
398  std::make_shared<dealii::AffineConstraints<Number>>();
399  dealii::DoFTools::make_hanging_node_constraints(*dof_handler.spatial(),
400  *space_constraints);
401 
402  unsigned int n_space_dofs = dof_handler.n_dofs_space();
403  for (unsigned int time_dof = 0; time_dof < dof_handler.n_dofs_time();
404  time_dof++)
405  {
406  for (unsigned int i = 0; i < n_space_dofs; i++)
407  {
408  if (space_constraints->is_constrained(i))
409  {
410  const std::vector<
411  std::pair<dealii::types::global_dof_index, double>> *entries =
412  space_constraints->get_constraint_entries(i);
413  spacetime_constraints->add_line(i + time_dof * n_space_dofs);
414  // non Dirichlet constraint
415  if (entries->size() > 0)
416  {
417  for (auto entry : *entries)
418  {
419  spacetime_constraints->add_entry(
420  i + time_dof * n_space_dofs,
421  entry.first + time_dof * n_space_dofs,
422  entry.second);
423  }
424  }
425  else
426  {
427  spacetime_constraints->set_inhomogeneity(
428  i + time_dof * n_space_dofs,
429  space_constraints->get_inhomogeneity(i));
430  }
431  }
432  }
433  }
434  }
435 
448  template <int dim, typename Number>
449  void
450  make_hanging_node_constraints(
451  dealii::IndexSet &space_relevant_dofs,
452  idealii::slab::DoFHandler<dim> &dof_handler,
453  std::shared_ptr<dealii::AffineConstraints<Number>> spacetime_constraints)
454  {
455  auto space_constraints =
456  std::make_shared<dealii::AffineConstraints<Number>>();
457  space_constraints->reinit(space_relevant_dofs);
458  dealii::DoFTools::make_hanging_node_constraints(*dof_handler.spatial(),
459  *space_constraints);
460 
461  unsigned int n_space_dofs = dof_handler.n_dofs_space();
462  for (unsigned int time_dof = 0; time_dof < dof_handler.n_dofs_time();
463  time_dof++)
464  {
465  for (auto id = space_relevant_dofs.begin();
466  id != space_relevant_dofs.end();
467  id++)
468  {
469  if (space_constraints->is_constrained(*id))
470  {
471  const std::vector<
472  std::pair<dealii::types::global_dof_index, double>> *entries =
473  space_constraints->get_constraint_entries(*id);
474 
475  spacetime_constraints->add_line(*id + time_dof * n_space_dofs);
476  // non Dirichlet constraint
477  if (entries->size() > 0)
478  {
479  for (auto entry : *entries)
480  {
481  spacetime_constraints->add_entry(
482  *id + time_dof * n_space_dofs,
483  entry.first + time_dof * n_space_dofs,
484  entry.second);
485  }
486  }
487  else
488  {
489  spacetime_constraints->set_inhomogeneity(
490  *id + time_dof * n_space_dofs,
491  space_constraints->get_inhomogeneity(*id));
492  }
493  }
494  }
495  }
496  }
497 
505  template <int dim>
506  dealii::IndexSet
507  extract_locally_relevant_dofs(slab::DoFHandler<dim> &dof_handler)
508  {
509  dealii::IndexSet space_relevant_dofs;
510  dealii::DoFTools::extract_locally_relevant_dofs(*dof_handler.spatial(),
511  space_relevant_dofs);
512 
513  dealii::IndexSet spacetime_relevant_dofs(space_relevant_dofs.size() *
514  dof_handler.n_dofs_time());
515 
516  for (dealii::types::global_dof_index time_dof_index = 0;
517  time_dof_index < dof_handler.n_dofs_time();
518  time_dof_index++)
519  {
520  spacetime_relevant_dofs.add_indices(
521  space_relevant_dofs,
522  time_dof_index * dof_handler.n_dofs_space() // offset
523  );
524  }
525 
526  return spacetime_relevant_dofs;
527  }
528 
529 } // namespace idealii::slab::DoFTools
530 
531 #endif /* INCLUDE_IDEAL_II_DOFS_SLAB_DOF_TOOLS_HH_ */
idealii::slab::DoFHandler
Actual DoFHandler for a specific slab.
Definition: slab_dof_handler.hh:43
idealii::slab::DoFHandler::n_dofs_time
unsigned int n_dofs_time()
Number of temporal degrees of fredom on this slab.
idealii::slab::DoFHandler::spatial
std::shared_ptr< dealii::DoFHandler< dim > > spatial()
The underlying spatial dof handler.
idealii::slab::DoFHandler::n_dofs_space
unsigned int n_dofs_space()
Number of spatial degrees of fredom on this slab.
idealii::slab::DoFTools
Collection of functions working on degrees of freedom.
Definition: idealii.hh:58
idealii::slab::DoFTools::extract_locally_relevant_dofs
dealii::IndexSet extract_locally_relevant_dofs(slab::DoFHandler< dim > &dof_handler)
Extract global space-time dof indices active on the current DoFHandler.
Definition: slab_dof_tools.hh:507
idealii::slab::DoFTools::make_downwind_sparsity_pattern
void make_downwind_sparsity_pattern(DoFHandler< dim > &dof, dealii::DynamicSparsityPattern &st_dsp, std::shared_ptr< dealii::AffineConstraints< double >> space_constraints=std::make_shared< dealii::AffineConstraints< double >>(), const bool keep_constrained_dofs=true)
Construction of a sparsity pattern with upper off diagonal jump terms.
Definition: slab_dof_tools.hh:305
idealii::slab::DoFTools::make_upwind_sparsity_pattern
void make_upwind_sparsity_pattern(DoFHandler< dim > &dof, dealii::DynamicSparsityPattern &st_dsp, std::shared_ptr< dealii::AffineConstraints< double >> space_constraints=std::make_shared< dealii::AffineConstraints< double >>(), const bool keep_constrained_dofs=true)
Construction of a sparsity pattern with lower off diagonal jump terms.
Definition: slab_dof_tools.hh:197
idealii::slab::DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(idealii::slab::DoFHandler< dim > &dof_handler, std::shared_ptr< dealii::AffineConstraints< Number >> spacetime_constraints)
Construction of space-time hanging node constraints.
Definition: slab_dof_tools.hh:393
Generated by doxygen 1.9.1