HyperHDG
bilaplacian_ldgh.hxx
Go to the documentation of this file.
1 #pragma once // Ensure that file is included only once in a single compilation.
2 
3 #include <HyperHDG/dense_la.hxx>
4 #include <HyperHDG/hypercube.hxx>
5 #include <tpp/quadrature/tensorial.hxx>
6 #include <tpp/shape_function/shape_function.hxx>
7 
8 #include <algorithm>
9 #include <tuple>
10 
11 namespace LocalSolver
12 {
13 /*!*************************************************************************************************
14  * \brief Default parameters for the diffusion equation, cf. below.
15  *
16  * \authors Guido Kanschat, Heidelberg University, 2019--2020.
17  * \authors Andreas Rupp, Heidelberg University, 2019--2020.
18  **************************************************************************************************/
19 template <unsigned int space_dimT, typename param_float_t = double>
20 struct Bilaplacian_parameters_default
21 {
22  /*!***********************************************************************************************
23  * \brief Array containing hypernode types corresponding to Dirichlet boundary.
24  ************************************************************************************************/
25  static constexpr std::array<unsigned int, 0U> dirichlet_nodes{};
26  /*!***********************************************************************************************
27  * \brief Array containing hypernode types corresponding to Dirichlet boundary of Laplacian.
28  ************************************************************************************************/
29  static constexpr std::array<unsigned int, 0U> dirichlet_laplacian_nodes{};
30  /*!***********************************************************************************************
31  * \brief Array containing hypernode types corresponding to Neumann boundary.
32  ************************************************************************************************/
33  static constexpr std::array<unsigned int, 0U> neumann_nodes{};
34  /*!***********************************************************************************************
35  * \brief Array containing hypernode types corresponding to Neumann boundary of Laplacian.
36  ************************************************************************************************/
37  static constexpr std::array<unsigned int, 0U> neumann_laplacian_nodes{};
38  /*!***********************************************************************************************
39  * \brief Inverse bilaplacian coefficient in PDE as analytic function.
40  ************************************************************************************************/
41  static param_float_t inverse_bilaplacian_coefficient(const Point<space_dimT, param_float_t>&,
42  const param_float_t = 0.)
43  {
44  return 1.;
45  }
46  /*!***********************************************************************************************
47  * \brief Right-hand side in PDE as analytic function.
48  ************************************************************************************************/
49  static param_float_t right_hand_side(const Point<space_dimT, param_float_t>&,
50  const param_float_t = 0.)
51  {
52  return 0.;
53  }
54  /*!***********************************************************************************************
55  * \brief Dirichlet values of solution as analytic function.
56  ************************************************************************************************/
57  static param_float_t dirichlet_value(const Point<space_dimT, param_float_t>&,
58  const param_float_t = 0.)
59  {
60  return 0.;
61  }
62  /*!***********************************************************************************************
63  * \brief Dirichlet values of solution's Laplacian as analytic function.
64  ************************************************************************************************/
65  static param_float_t dirichlet_laplace_value(const Point<space_dimT, param_float_t>&,
66  const param_float_t = 0.)
67  {
68  return 0.;
69  }
70  /*!***********************************************************************************************
71  * \brief Neumann values of solution as analytic function.
72  ************************************************************************************************/
73  static param_float_t neumann_value(const Point<space_dimT, param_float_t>&,
74  const param_float_t = 0.)
75  {
76  return 0.;
77  }
78  /*!***********************************************************************************************
79  * \brief Neumann values of solution's Laplacian as analytic function.
80  ************************************************************************************************/
81  static param_float_t neumann_laplace_value(const Point<space_dimT, param_float_t>&,
82  const param_float_t = 0.)
83  {
84  return 0.;
85  }
86  /*!***********************************************************************************************
87  * \brief Analytic result of PDE (for convergence tests).
88  ************************************************************************************************/
89  static param_float_t analytic_result(const Point<space_dimT, param_float_t>&,
90  const param_float_t = 0.)
91  {
92  return 0.;
93  }
94 }; // end of struct Bilaplacian_parameters_default
95 
96 /*!*************************************************************************************************
97  * \brief Local solver for bilaplacian equation on uniform hypergraph.
98  *
99  * Moreover, this class provides the functions that approximate the condensed mass matrix. By this,
100  * it is possible to approximate parabolic problems (in a very bad way) and to find good starting
101  * values for the nonlinear eigenvalue problem.
102  *
103  * \tparam hyEdge_dimT Dimension of a hyperedge, i.e., 1 is for PDEs defined on graphs, 2 is for
104  * PDEs defined on surfaces, and 3 is for PDEs defined on volumes.
105  * \tparam poly_deg The polynomial degree of test and trial functions.
106  * \tparam quad_deg The order of the quadrature rule.
107  * \tparam parametersT Struct depending on templates \c space_dimTP and \c lSol_float_TP that
108  * contains static parameter functions.
109  * Defaults to above functions included in \c BilaplacianParametersDefault.
110  * \tparam lSol_float_t The floating point type calculations are executed in. Defaults to double.
111  *
112  * \authors Guido Kanschat, Heidelberg University, 2019--2020.
113  * \authors Andreas Rupp, Heidelberg University, 2019--2020.
114  **************************************************************************************************/
115 template <unsigned int hyEdge_dimT,
116  unsigned int poly_deg,
117  unsigned int quad_deg,
118  template <unsigned int, typename> typename parametersT = Bilaplacian_parameters_default,
119  typename lSol_float_t = double>
120 class Bilaplacian
121 {
122  public:
123  /*!***********************************************************************************************
124  * \brief Define type of (hyperedge related) data that is stored in HyDataContainer.
125  ************************************************************************************************/
126  struct data_type
127  {
128  };
129  /*!***********************************************************************************************
130  * \brief Define type of node elements, especially with respect to nodal shape functions.
131  ************************************************************************************************/
132  struct node_element
133  {
134  typedef std::tuple<TPP::ShapeFunction<
135  TPP::ShapeType::Tensorial<TPP::ShapeType::Legendre<poly_deg>, hyEdge_dimT - 1>>>
136  functions;
137  };
138  /*!***********************************************************************************************
139  * \brief Define how errors are evaluated.
140  ************************************************************************************************/
141  struct error_def
142  {
143  /*!*********************************************************************************************
144  * \brief Define the typename returned by function errors.
145  **********************************************************************************************/
146  typedef std::array<lSol_float_t, 1U> error_t;
147  /*!*********************************************************************************************
148  * \brief Define how initial error is generated.
149  **********************************************************************************************/
150  static error_t initial_error()
151  {
152  std::array<lSol_float_t, 1U> summed_error;
153  summed_error.fill(0.);
154  return summed_error;
155  }
156  /*!*********************************************************************************************
157  * \brief Define how local errors should be accumulated.
158  **********************************************************************************************/
159  static error_t sum_error(error_t& summed_error, const error_t& new_error)
160  {
161  for (unsigned int k = 0; k < summed_error.size(); ++k)
162  summed_error[k] += new_error[k];
163  return summed_error;
164  }
165  /*!*********************************************************************************************
166  * \brief Define how global errors should be postprocessed.
167  **********************************************************************************************/
168  static error_t postprocess_error(error_t& summed_error)
169  {
170  for (unsigned int k = 0; k < summed_error.size(); ++k)
171  summed_error[k] = std::sqrt(summed_error[k]);
172  return summed_error;
173  }
174  };
175 
176  // -----------------------------------------------------------------------------------------------
177  // Public, static constexpr functions
178  // -----------------------------------------------------------------------------------------------
179 
180  /*!***********************************************************************************************
181  * \brief Dimension of hyper edge type that this object solves on.
182  ************************************************************************************************/
183  static constexpr unsigned int hyEdge_dim() { return hyEdge_dimT; }
184  /*!***********************************************************************************************
185  * \brief Evaluate amount of global degrees of freedom per hypernode.
186  *
187  * This number must be equal to HyperNodeFactory::n_dofs_per_node() of the HyperNodeFactory
188  * cooperating with this object.
189  *
190  * \retval n_dofs Number of global degrees of freedom per hypernode.
191  ************************************************************************************************/
192  static constexpr unsigned int n_glob_dofs_per_node()
193  {
194  return 2 * Hypercube<hyEdge_dimT - 1>::pow(poly_deg + 1);
195  }
196  /*!***********************************************************************************************
197  * \brief Dimension of of the solution evaluated with respect to a hyperedge.
198  ************************************************************************************************/
199  static constexpr unsigned int system_dimension() { return hyEdge_dimT + 1; }
200  /*!***********************************************************************************************
201  * \brief Dimension of of the solution evaluated with respect to a hypernode.
202  ************************************************************************************************/
203  static constexpr unsigned int node_system_dimension() { return 1; }
204 
205  private:
206  // -----------------------------------------------------------------------------------------------
207  // Private, static constexpr functions
208  // -----------------------------------------------------------------------------------------------
209 
210  /*!***********************************************************************************************
211  * \brief Number of local shape functions (with respect to all spatial dimensions).
212  ************************************************************************************************/
213  static constexpr unsigned int n_shape_fct_ = Hypercube<hyEdge_dimT>::pow(poly_deg + 1);
214  /*!***********************************************************************************************
215  * \brief Number oflocal shape functions (with respect to a face / hypernode).
216  ************************************************************************************************/
217  static constexpr unsigned int n_shape_bdr_ = Hypercube<hyEdge_dimT - 1>::pow(poly_deg + 1);
218  /*!***********************************************************************************************
219  * \brief Number of (local) degrees of freedom per hyperedge.
220  ************************************************************************************************/
221  static constexpr unsigned int n_loc_dofs_ = 2 * (hyEdge_dimT + 1) * n_shape_fct_;
222  /*!***********************************************************************************************
223  * \brief Dimension of of the solution evaluated with respect to a hypernode.
224  *
225  * This allows to the use of this quantity as template parameter in member functions.
226  ************************************************************************************************/
227  static constexpr unsigned int system_dim = system_dimension();
228  /*!***********************************************************************************************
229  * \brief Dimension of of the solution evaluated with respect to a hypernode.
230  *
231  * This allows to the use of this quantity as template parameter in member functions.
232  ************************************************************************************************/
233  static constexpr unsigned int node_system_dim = node_system_dimension();
234  /*!***********************************************************************************************
235  * \brief Find out whether a node is of Dirichlet type.
236  ************************************************************************************************/
237  template <typename parameters>
238  static constexpr bool is_dirichlet(const unsigned int node_type)
239  {
240  return std::find(parameters::dirichlet_nodes.begin(), parameters::dirichlet_nodes.end(),
241  node_type) != parameters::dirichlet_nodes.end();
242  }
243  template <typename parameters>
244  static constexpr bool is_dirichlet_laplacian(const unsigned int node_type)
245  {
246  return std::find(parameters::dirichlet_laplacian_nodes.begin(),
247  parameters::dirichlet_laplacian_nodes.end(),
248  node_type) != parameters::dirichlet_laplacian_nodes.end();
249  }
250 
251  // -----------------------------------------------------------------------------------------------
252  // Private, const members: Parameters and auxiliaries that help assembling matrices, etc.
253  // -----------------------------------------------------------------------------------------------
254 
255  /*!***********************************************************************************************
256  * \brief (Globally constant) penalty parameter for HDG scheme.
257  ************************************************************************************************/
258  const lSol_float_t tau_;
259  /*!***********************************************************************************************
260  * \brief An integrator helps to easily evaluate integrals (e.g. via quadrature).
261  ************************************************************************************************/
262  typedef TPP::Quadrature::Tensorial<
263  TPP::Quadrature::GaussLegendre<quad_deg>,
264  TPP::ShapeFunction<TPP::ShapeType::Tensorial<TPP::ShapeType::Legendre<poly_deg>, hyEdge_dimT>>,
265  lSol_float_t>
266  integrator;
267 
268  // -----------------------------------------------------------------------------------------------
269  // Private, internal functions for the local solver
270  // -----------------------------------------------------------------------------------------------
271 
272  /*!***********************************************************************************************
273  * \brief Assemble local matrix for the local solver.
274  *
275  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
276  * \param tau Penalty parameter.
277  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
278  * \param time Time, when the local matrix is evaluated.
279  * \retval loc_mat Matrix of the local solver.
280  ************************************************************************************************/
281  template <typename hyEdgeT>
282  inline SmallSquareMat<n_loc_dofs_, lSol_float_t>
283  assemble_loc_matrix(const lSol_float_t tau, hyEdgeT& hyper_edge, const lSol_float_t time) const;
284  /*!***********************************************************************************************
285  * \brief Assemble local right-hand for the local solver (from skeletal).
286  *
287  * The right hand side needs the values of the global degrees of freedom. Note that we basically
288  * follow the lines of
289  *
290  * B. Cockburn, J. Gopalakrishnan, and R. Lazarov.
291  * Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for
292  * second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009
293  *
294  * and discriminate between local solvers with respect to the skeletal variable and with respect
295  * to the global right-hand side. This assembles the local right-hand side with respect to the
296  * skeletal variable.
297  *
298  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
299  * \tparam SmallMatT The data type of the \c lambda values.
300  * \param lambda_values Global degrees of freedom associated to the hyperedge.
301  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
302  * \retval loc_rhs Local right hand side of the locasl solver.
303  ************************************************************************************************/
304  template <typename hyEdgeT, typename SmallMatT>
305  inline SmallVec<n_loc_dofs_, lSol_float_t> assemble_rhs_from_lambda(
306  const SmallMatT& lambda_values,
307  hyEdgeT& hyper_edge) const;
308  /*!***********************************************************************************************
309  * \brief Assemble local right-hand for the local solver (from global right-hand side).
310  *
311  * Note that we basically follow the lines of
312  *
313  * B. Cockburn, J. Gopalakrishnan, and R. Lazarov.
314  * Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for
315  * second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009
316  *
317  * and discriminate between local solvers with respect to the skeletal variable and with respect
318  * to the global right-hand side. This assembles the local right-hand side with respect to the
319  * global right-hand side. This function implicitly uses the parameters.
320  *
321  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
322  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
323  * \param time Point of time, rhs is evaluated at
324  * \retval loc_rhs Local right hand side of the locasl solver.
325  ************************************************************************************************/
326  template <typename hyEdgeT>
327  inline SmallVec<n_loc_dofs_, lSol_float_t> assemble_rhs_from_global_rhs(
328  hyEdgeT& hyper_edge,
329  const lSol_float_t time) const;
330  /*!***********************************************************************************************
331  * \brief Assemble local right-hand for the local solver (from volume function coefficients).
332  *
333  * Note that we basically follow the lines of
334  *
335  * B. Cockburn, J. Gopalakrishnan, and R. Lazarov.
336  * Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for
337  * second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009
338  *
339  * and discriminate between local solvers with respect to the skeletal variable and with respect
340  * to the global right-hand side. This assembles the local right-hand side with respect to the
341  * global right-hand side. This function implicitly uses the parameters.
342  *
343  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
344  * \param coeffs Local coefficients of bulk variable u.
345  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
346  * \retval loc_rhs Local right hand side of the locasl solver.
347  ************************************************************************************************/
348  template <typename hyEdgeT>
349  inline SmallVec<n_loc_dofs_, lSol_float_t> assemble_rhs_from_coeffs(
350  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
351  hyEdgeT& hyper_edge) const;
352  /*!***********************************************************************************************
353  * \brief Solve local problem (with right-hand side from skeletal).
354  *
355  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
356  * \tparam SmallMatT The data type of \c lambda_values.
357  * \param lambda_values Global degrees of freedom associated to the hyperedge.
358  * \param solution_type Type of local problem that is to be solved.
359  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
360  * \param time Point of time the problem is solved.
361  * \retval loc_sol Solution of the local problem.
362  ************************************************************************************************/
363  template <typename hyEdgeT, typename SmallMatT>
364  inline std::array<lSol_float_t, n_loc_dofs_> solve_local_problem(const SmallMatT& lambda_values,
365  const unsigned int solution_type,
366  hyEdgeT& hyper_edge,
367  const lSol_float_t time) const
368  {
369  try
370  {
371  SmallVec<n_loc_dofs_, lSol_float_t> rhs;
372  if (solution_type == 0)
373  rhs = assemble_rhs_from_lambda(lambda_values, hyper_edge);
374  else if (solution_type == 1)
375  rhs = assemble_rhs_from_lambda(lambda_values, hyper_edge) +
376  assemble_rhs_from_global_rhs(hyper_edge, time);
377  else
378  hy_assert(0 == 1, "This has not been implemented!");
379  return (rhs / assemble_loc_matrix(tau_, hyper_edge, time)).data();
380  }
381  catch (Wrapper::LAPACKexception& exc)
382  {
383  hy_assert(0 == 1, exc.what() << std::endl
384  << "This can happen if quadrature is too inaccurate!");
385  throw exc;
386  }
387  }
388  /*!***********************************************************************************************
389  * \brief Solve local problem for mass matrix approximation.
390  *
391  * \note This function is not use for elliptic problems.
392  *
393  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
394  * \param coeffs Global degrees of freedom associated to the hyperedge.
395  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
396  * \param time Point of time the local problem is solved at.
397  * \retval loc_sol Solution of the local problem.
398  ************************************************************************************************/
399  template <typename hyEdgeT>
400  inline std::array<lSol_float_t, n_loc_dofs_> solve_mass_problem(
401  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
402  hyEdgeT& hyper_edge,
403  const lSol_float_t time) const
404  {
405  try
406  {
407  SmallVec<n_loc_dofs_, lSol_float_t> rhs = assemble_rhs_from_coeffs(coeffs, hyper_edge);
408  return (rhs / assemble_loc_matrix(tau_, hyper_edge, time)).data();
409  }
410  catch (Wrapper::LAPACKexception& exc)
411  {
412  hy_assert(0 == 1, exc.what() << std::endl
413  << "This can happen if quadrature is too inaccurate!");
414  throw exc;
415  }
416  }
417  /*!***********************************************************************************************
418  * \brief Solve local problem for parabolic approximations.
419  *
420  * \note This function is not use for elliptic problems.
421  ************************************************************************************************/
422  template <typename hyEdgeT>
423  inline std::array<lSol_float_t, n_loc_dofs_> solve_loc_prob_cor(
424  const std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT>& lambda_values,
425  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
426  hyEdgeT& hyper_edge,
427  const lSol_float_t delta_time,
428  const lSol_float_t time) const
429  {
430  try
431  {
432  SmallVec<n_loc_dofs_, lSol_float_t> rhs;
433  rhs = delta_time * assemble_rhs_from_lambda(lambda_values, hyper_edge) +
434  delta_time * assemble_rhs_from_global_rhs(hyper_edge, time) +
435  assemble_rhs_from_coeffs(coeffs, hyper_edge);
436  return (rhs / (delta_time * assemble_loc_matrix(tau_, hyper_edge, time))).data();
437  }
438  catch (Wrapper::LAPACKexception& exc)
439  {
440  hy_assert(0 == 1, exc.what() << std::endl
441  << "This can happen if quadrature is too inaccurate!");
442  throw exc;
443  }
444  }
445  /*!***********************************************************************************************
446  * \brief Evaluate primal variable at boundary.
447  *
448  * Function to evaluate primal variable of the solution. This function is needed to calculate
449  * the local numerical fluxes.
450  *
451  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
452  * \param coeffs Coefficients of the local solution.
453  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
454  * \retval bdr_coeffs Coefficients of respective (dim-1) dimensional function at boundaries.
455  ************************************************************************************************/
456  template <typename hyEdgeT>
457  inline std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> primal_at_boundary(
458  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
459  hyEdgeT& hyper_edge) const;
460  /*!***********************************************************************************************
461  * \brief Evaluate dual variable at boundary.
462  *
463  * Function to evaluate dual variable of the solution. This function is needed to calculate the
464  * local numerical fluxes.
465  *
466  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
467  * \param coeffs Coefficients of the local solution.
468  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
469  * \retval bdr_coeffs Coefficients of respective (dim-1) dimensional function at boundaries.
470  ************************************************************************************************/
471  template <typename hyEdgeT>
472  inline std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> dual_at_boundary(
473  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
474  hyEdgeT& hyper_edge) const;
475 
476  public:
477  // -----------------------------------------------------------------------------------------------
478  // Public functions (and one typedef) to be utilized by external functions.
479  // -----------------------------------------------------------------------------------------------
480 
481  /*!***********************************************************************************************
482  * \brief Class is constructed using a single double indicating the penalty parameter.
483  ************************************************************************************************/
484  typedef lSol_float_t constructor_value_type;
485  /*!***********************************************************************************************
486  * \brief Constructor for local solver.
487  *
488  * \param tau Penalty parameter of HDG scheme.
489  ************************************************************************************************/
490  Bilaplacian(const constructor_value_type& tau = 1.) : tau_(tau) {}
491  /*!***********************************************************************************************
492  * \brief Evaluate local contribution to matrix--vector multiplication.
493  *
494  * Execute matrix--vector multiplication y = A * x, where x represents the vector containing the
495  * skeletal variable (adjacent to one hyperedge), and A is the condensed matrix arising from the
496  * HDG discretization. This function does this multiplication (locally) for one hyperedge. The
497  * hyperedge is no parameter, since all hyperedges are assumed to have the same properties.
498  *
499  **\tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
500  * \tparam SmallMatInT Data type of \c lambda_values_in.
501  * \tparam SmallMatOutT Data type of \c lambda_values_out.
502  * \param lambda_values_in Local part of vector x.
503  * \param lambda_values_out Local part that will be added to A * x.
504  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
505  * \param time Time.
506  * \retval vecAx Local part of vector A * x.
507  ************************************************************************************************/
508  template <class hyEdgeT, typename SmallMatInT, typename SmallMatOutT>
509  SmallMatOutT& trace_to_flux(const SmallMatInT& lambda_values_in,
510  SmallMatOutT& lambda_values_out,
511  hyEdgeT& hyper_edge,
512  const lSol_float_t time = 0.) const
513  {
514  hy_assert(lambda_values_in.size() == lambda_values_out.size() &&
515  lambda_values_in.size() == 2 * hyEdge_dimT,
516  "Both matrices must be of same size which corresponds to the number of faces!");
517  for (unsigned int i = 0; i < lambda_values_in.size(); ++i)
518  hy_assert(
519  lambda_values_in[i].size() == lambda_values_out[i].size() &&
520  lambda_values_in[i].size() == n_glob_dofs_per_node(),
521  "Both matrices must be of same size which corresponds to the number of dofs per face!");
522 
523  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
524  std::array<lSol_float_t, n_loc_dofs_> coeffs =
525  solve_local_problem(lambda_values_in, 0U, hyper_edge, time);
526 
527  std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> primals(
528  primal_at_boundary(coeffs, hyper_edge)),
529  duals(dual_at_boundary(coeffs, hyper_edge));
530 
531  for (unsigned int i = 0; i < lambda_values_out.size(); ++i)
532  {
533  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
534  lambda_values_out[i][j] =
535  duals[i][j] + tau_ * primals[i][j] -
536  tau_ * lambda_values_in[i][j < n_shape_bdr_ ? j + n_shape_bdr_ : j - n_shape_bdr_] *
537  hyper_edge.geometry.face_area(i);
538  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
539  for (unsigned int j = 0; j < lambda_values_out[i].size() / 2; ++j)
540  lambda_values_out[i][j] = 0.;
541  if (is_dirichlet_laplacian<parameters>(hyper_edge.node_descriptor[i]))
542  for (unsigned int j = lambda_values_out[i].size() / 2; j < lambda_values_out[i].size(); ++j)
543  lambda_values_out[i][j] = 0.;
544  }
545 
546  return lambda_values_out;
547  }
548  /*!***********************************************************************************************
549  * \brief Fill an array with 1 if the node is Dirichlet and 0 otherwise.
550  ************************************************************************************************/
551  template <class hyEdgeT>
552  std::array<unsigned int, 2 * hyEdge_dimT> node_types(hyEdgeT& hyper_edge) const
553  {
554  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
555 
556  std::array<unsigned int, 2 * hyEdge_dimT> result;
557  result.fill(0);
558 
559  for (unsigned int i = 0; i < 2 * hyEdge_dimT; ++i)
560  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
561  result[i] = 1;
562 
563  return result;
564  }
565  /*!***********************************************************************************************
566  * \brief Evaluate local contribution to residual \f$ A x - b \f$.
567  *
568  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
569  * \tparam SmallMatInT Data type of \c lambda_values_in.
570  * \tparam SmallMatOutT Data type of \c lambda_values_out.
571  * \param lambda_values_in Local part of vector x.
572  * \param lambda_values_out Local part that will be added to A * x.
573  * \param hyper_edge The geometry of the considered hyperedge (of typename hyEdgeT).
574  * \param time Time at which analytic functions are evaluated.
575  * \retval vecAx Local part of vector A * x - b.
576  ************************************************************************************************/
577  template <class hyEdgeT, typename SmallMatInT, typename SmallMatOutT>
578  SmallMatOutT& residual_flux(const SmallMatInT& lambda_values_in,
579  SmallMatOutT& lambda_values_out,
580  hyEdgeT& hyper_edge,
581  const lSol_float_t time = 0.) const
582  {
583  hy_assert(lambda_values_in.size() == lambda_values_out.size() &&
584  lambda_values_in.size() == 2 * hyEdge_dimT,
585  "Both matrices must be of same size which corresponds to the number of faces!");
586  for (unsigned int i = 0; i < lambda_values_in.size(); ++i)
587  hy_assert(
588  lambda_values_in[i].size() == lambda_values_out[i].size() &&
589  lambda_values_in[i].size() == n_glob_dofs_per_node(),
590  "Both matrices must be of same size which corresponds to the number of dofs per face!");
591 
592  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
593  std::array<lSol_float_t, n_loc_dofs_> coeffs =
594  solve_local_problem(lambda_values_in, 1U, hyper_edge, time);
595 
596  std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> primals(
597  primal_at_boundary(coeffs, hyper_edge)),
598  duals(dual_at_boundary(coeffs, hyper_edge));
599 
600  for (unsigned int i = 0; i < lambda_values_out.size(); ++i)
601  {
602  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
603  lambda_values_out[i][j] = duals[i][j] + tau_ * primals[i][j] -
604  tau_ * lambda_values_in[i][j] * hyper_edge.geometry.face_area(i);
605  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
606  for (unsigned int j = 0; j < lambda_values_out[i].size() / 2; ++j)
607  lambda_values_out[i][j] = 0.;
608  if (is_dirichlet_laplacian<parameters>(hyper_edge.node_descriptor[i]))
609  for (unsigned int j = lambda_values_out[i].size() / 2; j < lambda_values_out[i].size(); ++j)
610  lambda_values_out[i][j] = 0.;
611  }
612 
613  return lambda_values_out;
614  }
615  /*!***********************************************************************************************
616  * \brief Evaluate L2 projected lambda values at initial state.
617  *
618  * \note This function is not used for elliptic problems.
619  *
620  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
621  * \tparam SmallMatT Data type of \c lambda_values.
622  * \param lambda_values Local part of vector x.
623  * \param hyper_edge The geometry of the considered hyperedge (of typename hyEdgeT).
624  * \param time Initial time.
625  * \retval lambda_vakues L2 projected lambda ar initial state.
626  ************************************************************************************************/
627  template <class hyEdgeT, typename SmallMatT>
628  SmallMatT& make_initial(SmallMatT& lambda_values,
629  hyEdgeT& hyper_edge,
630  const lSol_float_t time = 0.) const
631  {
632  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
633 
634  for (unsigned int i = 0; i < lambda_values.size(); ++i)
635  {
636  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
637  for (unsigned int j = 0; j < lambda_values[i].size(); ++j)
638  lambda_values[i][j] = 0.;
639  else
640  {
641  for (unsigned int j = 0; j < lambda_values[i].size() / 2; ++j)
642  lambda_values[i][j] =
643  integrator ::template integrate_bdrUni_psifunc<decltype(hyEdgeT::geometry),
644  parameters::initial>(
645  i, j, hyper_edge.geometry, time);
646  for (unsigned int j = lambda_values[i].size() / 2; j < lambda_values[i].size(); ++j)
647  lambda_values[i][j] =
648  integrator::template integrate_bdrUni_psifunc<decltype(hyEdgeT::geometry),
649  parameters::initial_laplace>(
650  i, j, hyper_edge.geometry, time);
651  }
652  }
653 
654  return lambda_values;
655  }
656  /*!*********************************************************************************************
657  * \brief Evaluate local contribution to mass matrix--vector multiplication.
658  *
659  * Execute matrix--vector multiplication y = A * x, where x represents the vector containing the
660  * skeletal variable (adjacent to one hyperedge), and A is the condensed matrix arising from the
661  * HDG discretization. This function does this multiplication (locally) for one hyperedge. The
662  * hyperedge is no parameter, since all hyperedges are assumed to have the same properties.
663  *
664  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
665  * \tparam SmallMatInT Data type of \c lambda_values_in.
666  * \tparam SmallMatOutT Data type of \c lambda_values_out.
667  * \param lambda_values_in Local part of vector x.
668  * \param lambda_values_out Local part that will be added to A * x.
669  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
670  * \param time Time.
671  * \retval vecAx Local part of vector A * x.
672  **********************************************************************************************/
673  template <typename hyEdgeT, typename SmallMatInT, typename SmallMatOutT>
674  SmallMatOutT& trace_to_mass_flux(const SmallMatInT& lambda_values_in,
675  SmallMatOutT& lambda_values_out,
676  hyEdgeT& hyper_edge,
677  const lSol_float_t time = 0.) const
678  {
679  hy_assert(lambda_values_in.size() == lambda_values_out.size() &&
680  lambda_values_in.size() == 2 * hyEdge_dimT,
681  "Both matrices must be of same size which corresponds to the number of faces!");
682  for (unsigned int i = 0; i < lambda_values_in.size(); ++i)
683  hy_assert(
684  lambda_values_in[i].size() == lambda_values_out[i].size() &&
685  lambda_values_in[i].size() == n_glob_dofs_per_node(),
686  "Both matrices must be of same size which corresponds to the number of dofs per face!");
687 
688  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
689  std::array<lSol_float_t, n_loc_dofs_> coeffs =
690  solve_local_problem(lambda_values_in, 0U, hyper_edge, time);
691  coeffs = solve_mass_problem(coeffs, hyper_edge, time);
692 
693  std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> primals(
694  primal_at_boundary(coeffs, hyper_edge)),
695  duals(dual_at_boundary(coeffs, hyper_edge));
696 
697  for (unsigned int i = 0; i < lambda_values_out.size(); ++i)
698  {
699  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
700  lambda_values_out[i][j] = duals[i][j] + tau_ * primals[i][j];
701  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
702  for (unsigned int j = 0; j < lambda_values_out[i].size() / 2; ++j)
703  lambda_values_out[i][j] = 0.;
704  if (is_dirichlet_laplacian<parameters>(hyper_edge.node_descriptor[i]))
705  for (unsigned int j = lambda_values_out[i].size() / 2; j < lambda_values_out[i].size(); ++j)
706  lambda_values_out[i][j] = 0.;
707  }
708 
709  return lambda_values_out;
710  }
711  /*!*********************************************************************************************
712  * \brief Evaluate local contribution to mass matrix--vector residual.
713  *
714  * Execute matrix--vector multiplication y = A * x, where x represents the vector containing the
715  * skeletal variable (adjacent to one hyperedge), and A is the condensed matrix arising from the
716  * HDG discretization. This function does this multiplication (locally) for one hyperedge. The
717  * hyperedge is no parameter, since all hyperedges are assumed to have the same properties.
718  *
719  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
720  * \tparam SmallMatInT Data type of \c lambda_values_in.
721  * \tparam SmallMatOutT Data type of \c lambda_values_out.
722  * \param lambda_values_in Local part of vector x.
723  * \param lambda_values_out Local part that will be added to A * x.
724  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
725  * \param time Time.
726  * \retval vecAx Local part of vector A * x.
727  **********************************************************************************************/
728  template <class hyEdgeT, typename SmallMatInT, typename SmallMatOutT>
729  SmallMatOutT& total_numerical_flux_mass(const SmallMatInT& lambda_values_in,
730  SmallMatOutT& lambda_values_out,
731  hyEdgeT& hyper_edge,
732  const lSol_float_t time = 0.) const
733  {
734  hy_assert(lambda_values_in.size() == lambda_values_out.size() &&
735  lambda_values_in.size() == 2 * hyEdge_dimT,
736  "Both matrices must be of same size which corresponds to the number of faces!");
737  for (unsigned int i = 0; i < lambda_values_in.size(); ++i)
738  hy_assert(
739  lambda_values_in[i].size() == lambda_values_out[i].size() &&
740  lambda_values_in[i].size() == n_glob_dofs_per_node(),
741  "Both matrices must be of same size which corresponds to the number of dofs per face!");
742 
743  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
744  std::array<lSol_float_t, n_loc_dofs_> coeffs =
745  solve_local_problem(lambda_values_in, 1U, hyper_edge, time);
746  coeffs = solve_mass_problem(coeffs, hyper_edge, time);
747 
748  std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> primals(
749  primal_at_boundary(coeffs, hyper_edge)),
750  duals(dual_at_boundary(coeffs, hyper_edge));
751 
752  for (unsigned int i = 0; i < lambda_values_out.size(); ++i)
753  {
754  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
755  lambda_values_out[i][j] = duals[i][j] + tau_ * primals[i][j];
756  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
757  for (unsigned int j = 0; j < lambda_values_out[i].size() / 2; ++j)
758  lambda_values_out[i][j] = 0.;
759  if (is_dirichlet_laplacian<parameters>(hyper_edge.node_descriptor[i]))
760  for (unsigned int j = lambda_values_out[i].size() / 2; j < lambda_values_out[i].size(); ++j)
761  lambda_values_out[i][j] = 0.;
762  }
763 
764  return lambda_values_out;
765  }
766 
767  /*template < class hyEdgeT >
768  std::array< std::array<lSol_float_t, 2*n_shape_bdr_>, 2*hyEdge_dimT > trace_to_mass_flux
769  (
770  const std::array< std::array<lSol_float_t, 2*n_shape_bdr_>, 2*hyEdge_dimT > & lambda_values,
771  hyEdgeT & hyper_edge,
772  const lSol_float_t time = 0.
773  ) const
774  {
775  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
776  std::array<lSol_float_t, n_loc_dofs_> coeffs
777  = solve_local_problem(lambda_values, 0U, hyper_edge, time);
778  std::array< std::array<lSol_float_t, 2*n_shape_bdr_> , 2 * hyEdge_dimT > bdr_values;
779 
780  SmallVec<n_shape_fct_, lSol_float_t> u_coeffs, test_coeffs;
781  for (unsigned int i = 0; i < n_shape_fct_; ++i)
782  u_coeffs[i] = coeffs[hyEdge_dimT*n_shape_fct_+i];
783 
784  std::array< std::array<lSol_float_t, 2*n_shape_bdr_>, 2*hyEdge_dimT > lambda_values_uni;
785  for (unsigned int i = 0; i < lambda_values.size(); ++i) lambda_values_uni[i].fill(0.);
786  for (unsigned int i = 0; i < lambda_values.size(); ++i)
787  if ( is_dirichlet<parameters>(hyper_edge.node_descriptor[i]) )
788  for (unsigned int j = 0; j < lambda_values[i].size(); ++j) bdr_values[i][j] = 0.;
789  else
790  for (unsigned int j = 0; j < lambda_values[i].size(); ++j)
791  {
792  lambda_values_uni[i][j] = 1.;
793  coeffs = solve_local_problem(lambda_values_uni, 0U, hyper_edge, time);
794  for (unsigned int k = 0; k < n_shape_fct_; ++k)
795  test_coeffs[k] = coeffs[hyEdge_dimT*n_shape_fct_+k];
796  bdr_values[i][j] = integrator::integrate_vol_phiphi
797  (u_coeffs.data(), test_coeffs.data(), hyper_edge.geometry);
798  lambda_values_uni[i][j] = 0.;
799  }
800 
801  return bdr_values;
802  }
803 
804  template < class hyEdgeT >
805  std::array< std::array<lSol_float_t, 2*n_shape_bdr_>, 2*hyEdge_dimT > total_numerical_flux_mass
806  (
807  const std::array< std::array<lSol_float_t, 2*n_shape_bdr_>, 2*hyEdge_dimT > & lambda_values,
808  hyEdgeT & hyper_edge,
809  const lSol_float_t time = 0.
810  ) const
811  {
812  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
813  std::array<lSol_float_t, n_loc_dofs_> coeffs
814  = solve_local_problem(lambda_values, 1U, hyper_edge, time);
815  std::array< std::array<lSol_float_t, 2*n_shape_bdr_> , 2 * hyEdge_dimT > bdr_values;
816 
817  SmallVec<n_shape_fct_, lSol_float_t> u_coeffs, test_coeffs;
818  for (unsigned int i = 0; i < n_shape_fct_; ++i)
819  u_coeffs[i] = coeffs[hyEdge_dimT*n_shape_fct_+i];
820 
821  std::array< std::array<lSol_float_t, 2*n_shape_bdr_>, 2*hyEdge_dimT > lambda_values_uni;
822  for (unsigned int i = 0; i < lambda_values.size(); ++i) lambda_values_uni[i].fill(0.);
823  for (unsigned int i = 0; i < lambda_values.size(); ++i)
824  if ( is_dirichlet<parameters>(hyper_edge.node_descriptor[i]) )
825  for (unsigned int j = 0; j < lambda_values[i].size(); ++j) bdr_values[i][j] = 0.;
826  else
827  for (unsigned int j = 0; j < lambda_values[i].size(); ++j)
828  {
829  lambda_values_uni[i][j] = 1.;
830  coeffs = solve_local_problem(lambda_values_uni, 0U, hyper_edge, time);
831  for (unsigned int k = 0; k < n_shape_fct_; ++k)
832  test_coeffs[k] = coeffs[hyEdge_dimT*n_shape_fct_+k];
833  bdr_values[i][j] = integrator::integrate_vol_phiphi
834  (u_coeffs.data(), test_coeffs.data(), hyper_edge.geometry);
835  lambda_values_uni[i][j] = 0.;
836  }
837 
838  return bdr_values;
839  }*/
840 
841  /*!***********************************************************************************************
842  * \brief Calculate squared local contribution of L2 error.
843  *
844  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
845  * \tparam SmallMatT Data type of \c lambda_values.
846  * \param lambda_values The values of the skeletal variable's coefficients.
847  * \param hy_edge The geometry of the considered hyperedge (of typename GeomT).
848  * \param time Time at which anaytic functions are evaluated.
849  * \retval vec_b Local part of vector b.
850  ************************************************************************************************/
851  template <class hyEdgeT, typename SmallMatT>
852  std::array<lSol_float_t, 1U> errors(const SmallMatT& lambda_values,
853  hyEdgeT& hy_edge,
854  const lSol_float_t time = 0.) const
855  {
856  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
857 
858  for (unsigned int i = 0; i < lambda_values.size(); ++i)
859  for (unsigned int j = 0; j < lambda_values[i].size(); ++j)
860  hy_assert(lambda_values[i][j] == lambda_values[i][j],
861  "Lambda value wit index " << i << "," << j << " is NaN!");
862 
863  std::array<lSol_float_t, n_loc_dofs_> coefficients =
864  solve_local_problem(lambda_values, 1U, hy_edge, time);
865  std::array<lSol_float_t, n_shape_fct_> coeffs;
866  for (unsigned int i = 0; i < coeffs.size(); ++i)
867  coeffs[i] = coefficients[i + hyEdge_dimT * n_shape_fct_];
868 
869  for (unsigned int i = 0; i < coeffs.size(); ++i)
870  hy_assert(coeffs[i] == coeffs[i], "The " << i << "-th coeff is NaN!");
871 
872  lSol_float_t result = integrator::template integrate_vol_diffsquare_discana<
873  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>, decltype(hyEdgeT::geometry),
874  parameters::analytic_result, Point<hyEdge_dimT, lSol_float_t>>(coeffs, hy_edge.geometry,
875  time);
876  hy_assert(result >= 0., "The squared error must be non-negative, but was " << result);
877  return std::array<lSol_float_t, 1U>({result});
878  }
879  /*!***********************************************************************************************
880  * \brief Parabolic approximation version of local squared L2 error.
881  *
882  * \note This function is not used for elliptic problems.
883  *
884  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
885  * \tparam SmallMatT Data type of \c lambda_values.
886  * \param lambda_values_new Abscissas of the supporting points.
887  * \param lambda_values_old The values of the skeletal variable's coefficients.
888  * \param hy_edge The geometry of the considered hyperedge (of typename GeomT).
889  * \param delta_t Time step size.
890  * \param time Time.
891  * \retval vec_b Local part of vector b.
892  ************************************************************************************************/
893  template <class hyEdgeT, typename SmallMatT>
894  lSol_float_t errors_temp(const SmallMatT& lambda_values_new,
895  const SmallMatT& lambda_values_old,
896  hyEdgeT& hy_edge,
897  const lSol_float_t delta_t,
898  const lSol_float_t time) const
899  {
900  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
901 
902  std::array<lSol_float_t, n_loc_dofs_> coeffs_new =
903  solve_local_problem(lambda_values_new, 1U, hy_edge, time);
904  std::array<lSol_float_t, n_loc_dofs_> coeffs_old =
905  solve_local_problem(lambda_values_old, 1U, hy_edge, time - delta_t);
906  for (unsigned int i = 0; i < coeffs_old.size(); ++i)
907  coeffs_old[i] = (coeffs_old[i] - coeffs_new[i]);
908 
909  std::array<lSol_float_t, n_loc_dofs_> coefficients =
910  solve_loc_prob_cor(lambda_values_new, coeffs_old, hy_edge, delta_t, time);
911  std::array<lSol_float_t, n_shape_fct_> coeffs;
912  for (unsigned int i = 0; i < coeffs.size(); ++i)
913  coeffs[i] = coefficients[i + hyEdge_dimT * n_shape_fct_];
914  return integrator::template integrate_vol_diffsquare_discana<decltype(hyEdgeT::geometry),
915  parameters::analytic_result>(
916  coeffs, hy_edge.geometry, time);
917  }
918  /*!***********************************************************************************************
919  * \brief Evaluate local reconstruction at tensorial products of abscissas.
920  *
921  * \tparam absc_float_t Floating type for the abscissa values.
922  * \tparam abscissas_sizeT Size of the array of array of abscissas.
923  * \tparam input_array_t Type of input array.
924  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
925  * \param abscissas Abscissas of the supporting points.
926  * \param lambda_values The values of the skeletal variable's coefficients.
927  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
928  * \param time Time at which analytic functions are evaluated.
929  * \retval function_values Function values at quadrature points.
930  ************************************************************************************************/
931  template <typename abscissa_float_t, std::size_t sizeT, class input_array_t, class hyEdgeT>
932  std::array<std::array<lSol_float_t, Hypercube<hyEdge_dimT>::pow(sizeT)>, system_dim> bulk_values(
933  const std::array<abscissa_float_t, sizeT>& abscissas,
934  const input_array_t& lambda_values,
935  hyEdgeT& hyper_edge,
936  const lSol_float_t time = 0.) const;
937 }; // end of class Bilaplacian
938 
939 // -------------------------------------------------------------------------------------------------
941 // -------------------------------------------------------------------------------------------------
942 
943 // -------------------------------------------------------------------------------------------------
944 // -------------------------------------------------------------------------------------------------
945 //
946 // IMPLEMENTATION OF MEMBER FUNCTIONS OF Bilaplacian
947 //
948 // -------------------------------------------------------------------------------------------------
949 // -------------------------------------------------------------------------------------------------
950 
951 // -------------------------------------------------------------------------------------------------
952 // assemble_loc_matrix
953 // -------------------------------------------------------------------------------------------------
954 
955 template <unsigned int hyEdge_dimT,
956  unsigned int poly_deg,
957  unsigned int quad_deg,
958  template <unsigned int, typename>
959  typename parametersT,
960  typename lSol_float_t>
961 template <typename hyEdgeT>
962 inline SmallSquareMat<
963  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
964  lSol_float_t>
965 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::assemble_loc_matrix(
966  const lSol_float_t tau,
967  hyEdgeT& hyper_edge,
968  const lSol_float_t time) const
969 {
970  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
971  constexpr unsigned int n_dofs_lap = n_loc_dofs_ / 2;
972 
973  SmallSquareMat<n_loc_dofs_, lSol_float_t> local_mat;
974  lSol_float_t vol_integral, vol_func_integral, face_integral, helper;
975  SmallVec<hyEdge_dimT, lSol_float_t> grad_int_vec, normal_int_vec;
976 
977  for (unsigned int i = 0; i < n_shape_fct_; ++i)
978  {
979  for (unsigned int j = 0; j < n_shape_fct_; ++j)
980  {
981  // Integral_element phi_i phi_j dx in diagonal blocks
982  vol_integral = integrator::template integrate_vol_phiphi(i, j, hyper_edge.geometry);
983  vol_func_integral = integrator::template integrate_vol_phiphifunc<
984  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>, decltype(hyEdgeT::geometry),
985  parameters::inverse_bilaplacian_coefficient, Point<hyEdge_dimT, lSol_float_t>>(
986  i, j, hyper_edge.geometry, time);
987  // Integral_element - nabla phi_i \vec phi_j dx
988  // = Integral_element - div \vec phi_i phi_j dx in right upper and left lower blocks
989  grad_int_vec =
990  integrator::template integrate_vol_nablaphiphi<SmallVec<hyEdge_dimT, lSol_float_t>,
991  decltype(hyEdgeT::geometry)>(
992  i, j, hyper_edge.geometry);
993 
994  face_integral = 0.;
995  normal_int_vec = 0.;
996  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
997  {
998  helper = integrator::template integrate_bdr_phiphi<decltype(hyEdgeT::geometry)>(
999  i, j, face, hyper_edge.geometry);
1000  face_integral += helper;
1001  for (unsigned int dim = 0; dim < hyEdge_dimT; ++dim)
1002  normal_int_vec[dim] += hyper_edge.geometry.local_normal(face).operator[](dim) * helper;
1003  }
1004 
1005  local_mat(hyEdge_dimT * n_shape_fct_ + i, n_dofs_lap + hyEdge_dimT * n_shape_fct_ + j) -=
1006  vol_func_integral;
1007 
1008  local_mat(hyEdge_dimT * n_shape_fct_ + i, hyEdge_dimT * n_shape_fct_ + j) +=
1009  tau * face_integral;
1010  local_mat(n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i,
1011  n_dofs_lap + hyEdge_dimT * n_shape_fct_ + j) += tau * face_integral;
1012 
1013  for (unsigned int dim = 0; dim < hyEdge_dimT; ++dim)
1014  {
1015  local_mat(dim * n_shape_fct_ + i, dim * n_shape_fct_ + j) += vol_integral;
1016  local_mat(n_dofs_lap + dim * n_shape_fct_ + i, n_dofs_lap + dim * n_shape_fct_ + j) +=
1017  vol_integral;
1018 
1019  local_mat(hyEdge_dimT * n_shape_fct_ + i, dim * n_shape_fct_ + j) -= grad_int_vec[dim];
1020  local_mat(dim * n_shape_fct_ + i, hyEdge_dimT * n_shape_fct_ + j) -= grad_int_vec[dim];
1021  local_mat(n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i,
1022  n_dofs_lap + dim * n_shape_fct_ + j) -= grad_int_vec[dim];
1023  local_mat(n_dofs_lap + dim * n_shape_fct_ + i,
1024  n_dofs_lap + hyEdge_dimT * n_shape_fct_ + j) -= grad_int_vec[dim];
1025 
1026  local_mat(hyEdge_dimT * n_shape_fct_ + i, dim * n_shape_fct_ + j) += normal_int_vec[dim];
1027  local_mat(n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i,
1028  n_dofs_lap + dim * n_shape_fct_ + j) += normal_int_vec[dim];
1029  }
1030  }
1031  }
1032 
1033  return local_mat;
1034 } // end of Bilaplacian::assemble_loc_matrix
1035 
1036 // -------------------------------------------------------------------------------------------------
1037 // assemble_rhs_from_lambda
1038 // -------------------------------------------------------------------------------------------------
1039 
1040 template <unsigned int hyEdge_dimT,
1041  unsigned int poly_deg,
1042  unsigned int quad_deg,
1043  template <unsigned int, typename>
1044  typename parametersT,
1045  typename lSol_float_t>
1046 template <typename hyEdgeT, typename SmallMatT>
1047 inline SmallVec<
1048  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
1049  lSol_float_t>
1050 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::assemble_rhs_from_lambda(
1051  const SmallMatT& lambda_values,
1052  hyEdgeT& hyper_edge) const
1053 {
1054  constexpr unsigned int n_dofs_lap = n_loc_dofs_ / 2;
1055  hy_assert(lambda_values.size() == 2 * hyEdge_dimT,
1056  "The size of the lambda values should be twice the dimension of a hyperedge.");
1057  for (unsigned int i = 0; i < 2 * hyEdge_dimT; ++i)
1058  hy_assert(lambda_values[i].size() == 2 * n_shape_bdr_,
1059  "The size of lambda should be the amount of ansatz functions at boundary.");
1060 
1061  SmallVec<n_loc_dofs_, lSol_float_t> right_hand_side;
1062  lSol_float_t integral;
1063 
1064  for (unsigned int i = 0; i < n_shape_fct_; ++i)
1065  for (unsigned int j = 0; j < n_shape_bdr_; ++j)
1066  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
1067  {
1068  integral = integrator::template integrate_bdr_phipsi<decltype(hyEdgeT::geometry)>(
1069  i, j, face, hyper_edge.geometry);
1070  right_hand_side[hyEdge_dimT * n_shape_fct_ + i] += tau_ * lambda_values[face][j] * integral;
1071  right_hand_side[n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i] +=
1072  tau_ * lambda_values[face][n_shape_bdr_ + j] * integral;
1073 
1074  for (unsigned int dim = 0; dim < hyEdge_dimT; ++dim)
1075  {
1076  right_hand_side[dim * n_shape_fct_ + i] -=
1077  hyper_edge.geometry.local_normal(face).operator[](dim) * lambda_values[face][j] *
1078  integral;
1079  right_hand_side[n_dofs_lap + dim * n_shape_fct_ + i] -=
1080  hyper_edge.geometry.local_normal(face).operator[](dim) *
1081  lambda_values[face][n_shape_bdr_ + j] * integral;
1082  }
1083  }
1084 
1085  return right_hand_side;
1086 } // end of Bilaplacian::assemble_rhs_from_lambda
1087 
1088 // -------------------------------------------------------------------------------------------------
1089 // assemble_rhs_from_global_rhs
1090 // -------------------------------------------------------------------------------------------------
1091 
1092 template <unsigned int hyEdge_dimT,
1093  unsigned int poly_deg,
1094  unsigned int quad_deg,
1095  template <unsigned int, typename>
1096  typename parametersT,
1097  typename lSol_float_t>
1098 template <typename hyEdgeT>
1099 inline SmallVec<
1100  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
1101  lSol_float_t>
1102 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::
1103  assemble_rhs_from_global_rhs(hyEdgeT& hyper_edge, const lSol_float_t time) const
1104 {
1105  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
1106  constexpr unsigned int n_dofs_lap = n_loc_dofs_ / 2;
1107  SmallVec<n_loc_dofs_, lSol_float_t> right_hand_side;
1108  lSol_float_t integral;
1109  for (unsigned int i = 0; i < n_shape_fct_; ++i)
1110  {
1111  right_hand_side[n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i] =
1112  integrator::template integrate_vol_phifunc<
1113  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>, decltype(hyEdgeT::geometry),
1114  parameters::right_hand_side, Point<hyEdge_dimT, lSol_float_t>>(i, hyper_edge.geometry,
1115  time);
1116  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
1117  {
1118  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[face]))
1119  {
1120  integral = integrator::template integrate_bdr_phifunc<
1121  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>,
1122  decltype(hyEdgeT::geometry), parameters::dirichlet_value,
1123  Point<hyEdge_dimT, lSol_float_t>>(i, face, hyper_edge.geometry, time);
1124  right_hand_side[hyEdge_dimT * n_shape_fct_ + i] += tau_ * integral;
1125  for (unsigned int dim = 0; dim < hyEdge_dimT; ++dim)
1126  right_hand_side[dim * n_shape_fct_ + i] -=
1127  hyper_edge.geometry.local_normal(face).operator[](dim) * integral;
1128  }
1129  if (is_dirichlet_laplacian<parameters>(hyper_edge.node_descriptor[face]))
1130  {
1131  integral = integrator::template integrate_bdr_phifunc<
1132  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>,
1133  decltype(hyEdgeT::geometry), parameters::dirichlet_laplace_value,
1134  Point<hyEdge_dimT, lSol_float_t>>(i, face, hyper_edge.geometry, time);
1135  right_hand_side[n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i] += tau_ * integral;
1136  for (unsigned int dim = 0; dim < hyEdge_dimT; ++dim)
1137  right_hand_side[n_dofs_lap + dim * n_shape_fct_ + i] -=
1138  hyper_edge.geometry.local_normal(face).operator[](dim) * integral;
1139  }
1140  }
1141  }
1142  return right_hand_side;
1143 } // end of Bilaplacian::assemble_rhs_from_global_rhs
1144 
1145 // -------------------------------------------------------------------------------------------------
1146 // assemble_rhs_from_coeffs
1147 // -------------------------------------------------------------------------------------------------
1148 
1149 template <unsigned int hyEdge_dimT,
1150  unsigned int poly_deg,
1151  unsigned int quad_deg,
1152  template <unsigned int, typename>
1153  typename parametersT,
1154  typename lSol_float_t>
1155 template <typename hyEdgeT>
1156 inline SmallVec<
1157  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
1158  lSol_float_t>
1159 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::assemble_rhs_from_coeffs(
1160  const std::array<
1161  lSol_float_t,
1162  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_>& coeffs,
1163  hyEdgeT& hyper_edge) const
1164 {
1165  SmallVec<n_loc_dofs_, lSol_float_t> right_hand_side;
1166  for (unsigned int i = 0; i < n_shape_fct_; ++i)
1167  for (unsigned int j = 0; j < n_shape_fct_; ++j)
1168  {
1169  right_hand_side[n_loc_dofs_ / 2 + hyEdge_dimT * n_shape_fct_ + i] +=
1170  coeffs[hyEdge_dimT * n_shape_fct_ + j] *
1171  integrator::template integrate_vol_phiphi(i, j, hyper_edge.geometry);
1172  }
1173 
1174  return right_hand_side;
1175 } // end of Diffusion::assemble_rhs_from_coeffs
1176 
1177 // -------------------------------------------------------------------------------------------------
1178 // primal_at_boundary
1179 // -------------------------------------------------------------------------------------------------
1180 
1181 template <unsigned int hyEdge_dimT,
1182  unsigned int poly_deg,
1183  unsigned int quad_deg,
1184  template <unsigned int, typename>
1185  typename parametersT,
1186  typename lSol_float_t>
1187 template <typename hyEdgeT>
1188 inline std::array<
1189  std::array<
1190  lSol_float_t,
1191  2 * Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_shape_bdr_>,
1192  2 * hyEdge_dimT>
1193 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::primal_at_boundary(
1194  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
1195  hyEdgeT& hyper_edge) const
1196 {
1197  constexpr unsigned int n_dofs_lap = n_loc_dofs_ / 2;
1198  std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> bdr_values;
1199 
1200  for (unsigned int dim_n = 0; dim_n < 2 * hyEdge_dimT; ++dim_n)
1201  bdr_values[dim_n].fill(0.);
1202 
1203  for (unsigned int i = 0; i < n_shape_fct_; ++i)
1204  for (unsigned int j = 0; j < n_shape_bdr_; ++j)
1205  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
1206  {
1207  bdr_values[face][n_shape_bdr_ + j] +=
1208  coeffs[hyEdge_dimT * n_shape_fct_ + i] *
1209  integrator::template integrate_bdr_phipsi<decltype(hyEdgeT::geometry)>(
1210  i, j, face, hyper_edge.geometry);
1211 
1212  bdr_values[face][j] +=
1213  coeffs[n_dofs_lap + hyEdge_dimT * n_shape_fct_ + i] *
1214  integrator::template integrate_bdr_phipsi<decltype(hyEdgeT::geometry)>(
1215  i, j, face, hyper_edge.geometry);
1216  }
1217 
1218  return bdr_values;
1219 } // end of Bilaplacian::primal_at_boundary
1220 
1221 // -------------------------------------------------------------------------------------------------
1222 // dual_at_boundary
1223 // -------------------------------------------------------------------------------------------------
1224 
1225 template <unsigned int hyEdge_dimT,
1226  unsigned int poly_deg,
1227  unsigned int quad_deg,
1228  template <unsigned int, typename>
1229  typename parametersT,
1230  typename lSol_float_t>
1231 template <typename hyEdgeT>
1232 inline std::array<
1233  std::array<
1234  lSol_float_t,
1235  2 * Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_shape_bdr_>,
1236  2 * hyEdge_dimT>
1237 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::dual_at_boundary(
1238  const std::array<lSol_float_t, n_loc_dofs_>& coeffs,
1239  hyEdgeT& hyper_edge) const
1240 {
1241  constexpr unsigned int n_dofs_lap = n_loc_dofs_ / 2;
1242  std::array<std::array<lSol_float_t, 2 * n_shape_bdr_>, 2 * hyEdge_dimT> bdr_values;
1243  lSol_float_t integral;
1244 
1245  for (unsigned int dim_n = 0; dim_n < 2 * hyEdge_dimT; ++dim_n)
1246  bdr_values[dim_n].fill(0.);
1247 
1248  for (unsigned int i = 0; i < n_shape_fct_; ++i)
1249  for (unsigned int j = 0; j < n_shape_bdr_; ++j)
1250  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
1251  {
1252  integral = integrator::template integrate_bdr_phipsi<decltype(hyEdgeT::geometry)>(
1253  i, j, face, hyper_edge.geometry);
1254  for (unsigned int dim = 0; dim < hyEdge_dimT; ++dim)
1255  {
1256  bdr_values[face][n_shape_bdr_ + j] +=
1257  hyper_edge.geometry.local_normal(face).operator[](dim) * integral *
1258  coeffs[dim * n_shape_fct_ + i];
1259 
1260  bdr_values[face][j] += hyper_edge.geometry.local_normal(face).operator[](dim) * integral *
1261  coeffs[n_dofs_lap + dim * n_shape_fct_ + i];
1262  }
1263  }
1264 
1265  return bdr_values;
1266 } // end of Bilaplacian::dual_at_boundary
1267 
1268 // -------------------------------------------------------------------------------------------------
1269 // bulk_values
1270 // -------------------------------------------------------------------------------------------------
1271 
1272 template <unsigned int hyEdge_dimT,
1273  unsigned int poly_deg,
1274  unsigned int quad_deg,
1275  template <unsigned int, typename>
1276  typename parametersT,
1277  typename lSol_float_t>
1278 template <typename abscissa_float_t, std::size_t sizeT, class input_array_t, typename hyEdgeT>
1279 std::array<std::array<lSol_float_t, Hypercube<hyEdge_dimT>::pow(sizeT)>,
1280  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::system_dim>
1281 Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::bulk_values(
1282  const std::array<abscissa_float_t, sizeT>& abscissas,
1283  const input_array_t& lambda_values,
1284  hyEdgeT& hyper_edge,
1285  const lSol_float_t time) const
1286 {
1287  SmallVec<n_loc_dofs_, lSol_float_t> coefficients =
1288  solve_local_problem(lambda_values, 1U, hyper_edge, time);
1289  SmallVec<n_shape_fct_, lSol_float_t> coeffs;
1290  SmallVec<static_cast<unsigned int>(sizeT), abscissa_float_t> helper(abscissas);
1291 
1292  std::array<std::array<lSol_float_t, Hypercube<hyEdge_dimT>::pow(sizeT)>,
1293  Bilaplacian<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::system_dim>
1294  point_vals;
1295 
1296  for (unsigned int d = 0; d < system_dim; ++d)
1297  {
1298  for (unsigned int i = 0; i < coeffs.size(); ++i)
1299  coeffs[i] = coefficients[d * n_shape_fct_ + i];
1300  for (unsigned int pt = 0; pt < Hypercube<hyEdge_dimT>::pow(sizeT); ++pt)
1301  point_vals[d][pt] = integrator::shape_fun_t::template lin_comb_fct_val<float>(
1302  coeffs, Hypercube<hyEdge_dimT>::template tensorial_pt<Point<hyEdge_dimT>>(pt, helper));
1303  }
1304 
1305  return point_vals;
1306 } // end of Bilaplacian::bulk_values
1307 
1308 // -------------------------------------------------------------------------------------------------
1310 // -------------------------------------------------------------------------------------------------
1311 
1312 } // namespace LocalSolver
LocalSolver::Bilaplacian::data_type
Define type of (hyperedge related) data that is stored in HyDataContainer.
Definition: bilaplacian_ldgh.hxx:126
LocalSolver::Bilaplacian::assemble_rhs_from_global_rhs
SmallVec< n_loc_dofs_, lSol_float_t > assemble_rhs_from_global_rhs(hyEdgeT &hyper_edge, const lSol_float_t time) const
Assemble local right-hand for the local solver (from global right-hand side).
LocalSolver::Bilaplacian_parameters_default::inverse_bilaplacian_coefficient
static param_float_t inverse_bilaplacian_coefficient(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Inverse bilaplacian coefficient in PDE as analytic function.
Definition: bilaplacian_ldgh.hxx:41
LocalSolver::Bilaplacian::make_initial
SmallMatT & make_initial(SmallMatT &lambda_values, hyEdgeT &hyper_edge, const lSol_float_t time=0.) const
Evaluate L2 projected lambda values at initial state.
Definition: bilaplacian_ldgh.hxx:628
shape_function.hxx
LocalSolver::Bilaplacian::hyEdge_dim
static constexpr unsigned int hyEdge_dim()
Dimension of hyper edge type that this object solves on.
Definition: bilaplacian_ldgh.hxx:183
LocalSolver::Bilaplacian::trace_to_flux
SmallMatOutT & trace_to_flux(const SmallMatInT &lambda_values_in, SmallMatOutT &lambda_values_out, hyEdgeT &hyper_edge, const lSol_float_t time=0.) const
Evaluate local contribution to matrix–vector multiplication.
Definition: bilaplacian_ldgh.hxx:509
LocalSolver::Bilaplacian::solve_mass_problem
std::array< lSol_float_t, n_loc_dofs_ > solve_mass_problem(const std::array< lSol_float_t, n_loc_dofs_ > &coeffs, hyEdgeT &hyper_edge, const lSol_float_t time) const
Solve local problem for mass matrix approximation.
Definition: bilaplacian_ldgh.hxx:400
LocalSolver::Bilaplacian_parameters_default::dirichlet_laplacian_nodes
static constexpr std::array< unsigned int, 0U > dirichlet_laplacian_nodes
Array containing hypernode types corresponding to Dirichlet boundary of Laplacian.
Definition: bilaplacian_eigs_ldgh.hxx:29
LocalSolver::Bilaplacian::errors
std::array< lSol_float_t, 1U > errors(const SmallMatT &lambda_values, hyEdgeT &hy_edge, const lSol_float_t time=0.) const
Calculate squared local contribution of L2 error.
Definition: bilaplacian_ldgh.hxx:852
LocalSolver::Bilaplacian::node_system_dim
static constexpr unsigned int node_system_dim
Dimension of of the solution evaluated with respect to a hypernode.
Definition: bilaplacian_ldgh.hxx:233
LocalSolver::Bilaplacian::n_loc_dofs_
static constexpr unsigned int n_loc_dofs_
Number of (local) degrees of freedom per hyperedge.
Definition: bilaplacian_ldgh.hxx:221
Wrapper::LAPACKexception::what
const char * what() const
Definition: lapack.hxx:27
LocalSolver
A namespace for local solvers.
Definition: advection_parab_ldgh.hxx:11
tensorial.hxx
diffusion_uniform.geometry
geometry
Definition: diffusion_uniform.py:19
LocalSolver::Bilaplacian_parameters_default::dirichlet_nodes
static constexpr std::array< unsigned int, 0U > dirichlet_nodes
Array containing hypernode types corresponding to Dirichlet boundary.
Definition: bilaplacian_eigs_ldgh.hxx:25
LocalSolver::Bilaplacian::assemble_rhs_from_coeffs
SmallVec< n_loc_dofs_, lSol_float_t > assemble_rhs_from_coeffs(const std::array< lSol_float_t, n_loc_dofs_ > &coeffs, hyEdgeT &hyper_edge) const
Assemble local right-hand for the local solver (from volume function coefficients).
LocalSolver::Bilaplacian::n_shape_fct_
static constexpr unsigned int n_shape_fct_
Number of local shape functions (with respect to all spatial dimensions).
Definition: bilaplacian_ldgh.hxx:213
LocalSolver::Bilaplacian_parameters_default::neumann_laplacian_nodes
static constexpr std::array< unsigned int, 0U > neumann_laplacian_nodes
Array containing hypernode types corresponding to Neumann boundary of Laplacian.
Definition: bilaplacian_eigs_ldgh.hxx:37
LocalSolver::Bilaplacian::error_def::sum_error
static error_t sum_error(error_t &summed_error, const error_t &new_error)
Define how local errors should be accumulated.
Definition: bilaplacian_ldgh.hxx:159
LocalSolver::Bilaplacian::constructor_value_type
lSol_float_t constructor_value_type
Class is constructed using a single double indicating the penalty parameter.
Definition: bilaplacian_ldgh.hxx:484
LocalSolver::Bilaplacian::trace_to_mass_flux
SmallMatOutT & trace_to_mass_flux(const SmallMatInT &lambda_values_in, SmallMatOutT &lambda_values_out, hyEdgeT &hyper_edge, const lSol_float_t time=0.) const
Evaluate local contribution to mass matrix–vector multiplication.
Definition: bilaplacian_ldgh.hxx:674
TPP::Quadrature::GaussLegendre
Gauss–Legendre quadrature points and weights in one spatial dimension.
Definition: one_dimensional.hxx:19
LocalSolver::Bilaplacian_parameters_default::dirichlet_value
static param_float_t dirichlet_value(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Dirichlet values of solution as analytic function.
Definition: bilaplacian_ldgh.hxx:57
SmallMat::size
static constexpr unsigned int size()
Return size a SmallMat.
Definition: dense_la.hxx:54
LocalSolver::Bilaplacian_parameters_default::neumann_nodes
static constexpr std::array< unsigned int, 0U > neumann_nodes
Array containing hypernode types corresponding to Neumann boundary.
Definition: bilaplacian_eigs_ldgh.hxx:33
LocalSolver::Bilaplacian_parameters_default::dirichlet_laplace_value
static param_float_t dirichlet_laplace_value(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Dirichlet values of solution's Laplacian as analytic function.
Definition: bilaplacian_ldgh.hxx:65
LocalSolver::Bilaplacian::n_shape_bdr_
static constexpr unsigned int n_shape_bdr_
Number oflocal shape functions (with respect to a face / hypernode).
Definition: bilaplacian_ldgh.hxx:217
LocalSolver::Bilaplacian::is_dirichlet_laplacian
static constexpr bool is_dirichlet_laplacian(const unsigned int node_type)
Definition: bilaplacian_ldgh.hxx:244
TPP::ShapeType::Tensorial
Struct that handles the evaluation of tensorial shape functions.
Definition: tensorial.hxx:22
LocalSolver::Bilaplacian::node_types
std::array< unsigned int, 2 *hyEdge_dimT > node_types(hyEdgeT &hyper_edge) const
Fill an array with 1 if the node is Dirichlet and 0 otherwise.
Definition: bilaplacian_ldgh.hxx:552
LocalSolver::Bilaplacian::node_system_dimension
static constexpr unsigned int node_system_dimension()
Dimension of of the solution evaluated with respect to a hypernode.
Definition: bilaplacian_ldgh.hxx:203
LocalSolver::Bilaplacian::system_dim
static constexpr unsigned int system_dim
Dimension of of the solution evaluated with respect to a hypernode.
Definition: bilaplacian_ldgh.hxx:227
TPP::ShapeFunction
Struct that handles different types of evaluation of shape functions.
Definition: shape_function.hxx:15
TPP::Quadrature::Tensorial
General integrator class on tensorial hypergraphs.
Definition: tensorial.hxx:24
LocalSolver::Bilaplacian::dual_at_boundary
std::array< std::array< lSol_float_t, 2 *n_shape_bdr_ >, 2 *hyEdge_dimT > dual_at_boundary(const std::array< lSol_float_t, n_loc_dofs_ > &coeffs, hyEdgeT &hyper_edge) const
Evaluate dual variable at boundary.
LocalSolver::Bilaplacian::integrator
TPP::Quadrature::Tensorial< TPP::Quadrature::GaussLegendre< quad_deg >, TPP::ShapeFunction< TPP::ShapeType::Tensorial< TPP::ShapeType::Legendre< poly_deg >, hyEdge_dimT > >, lSol_float_t > integrator
An integrator helps to easily evaluate integrals (e.g. via quadrature).
Definition: bilaplacian_ldgh.hxx:266
LocalSolver::Bilaplacian::solve_loc_prob_cor
std::array< lSol_float_t, n_loc_dofs_ > solve_loc_prob_cor(const std::array< std::array< lSol_float_t, 2 *n_shape_bdr_ >, 2 *hyEdge_dimT > &lambda_values, const std::array< lSol_float_t, n_loc_dofs_ > &coeffs, hyEdgeT &hyper_edge, const lSol_float_t delta_time, const lSol_float_t time) const
Solve local problem for parabolic approximations.
Definition: bilaplacian_ldgh.hxx:423
LocalSolver::Bilaplacian::error_def::error_t
std::array< lSol_float_t, 1U > error_t
Define the typename returned by function errors.
Definition: bilaplacian_ldgh.hxx:146
LocalSolver::Bilaplacian::is_dirichlet
static constexpr bool is_dirichlet(const unsigned int node_type)
Find out whether a node is of Dirichlet type.
Definition: bilaplacian_ldgh.hxx:238
LocalSolver::Bilaplacian::error_def::initial_error
static error_t initial_error()
Define how initial error is generated.
Definition: bilaplacian_ldgh.hxx:150
LocalSolver::Bilaplacian::node_element
Define type of node elements, especially with respect to nodal shape functions.
Definition: bilaplacian_ldgh.hxx:132
hypercube.hxx
LocalSolver::Bilaplacian::system_dimension
static constexpr unsigned int system_dimension()
Dimension of of the solution evaluated with respect to a hyperedge.
Definition: bilaplacian_ldgh.hxx:199
Wrapper::LAPACKexception
Exception to be thrown if LAPACK's solve fails.
Definition: lapack.hxx:25
Hypercube::pow
static constexpr unsigned int pow(unsigned int base)
Return n to the power dim.
Definition: hypercube.hxx:25
LocalSolver::Bilaplacian::tau_
const lSol_float_t tau_
(Globally constant) penalty parameter for HDG scheme.
Definition: bilaplacian_ldgh.hxx:258
LocalSolver::Bilaplacian::residual_flux
SmallMatOutT & residual_flux(const SmallMatInT &lambda_values_in, SmallMatOutT &lambda_values_out, hyEdgeT &hyper_edge, const lSol_float_t time=0.) const
Evaluate local contribution to residual .
Definition: bilaplacian_ldgh.hxx:578
LocalSolver::Bilaplacian::errors_temp
lSol_float_t errors_temp(const SmallMatT &lambda_values_new, const SmallMatT &lambda_values_old, hyEdgeT &hy_edge, const lSol_float_t delta_t, const lSol_float_t time) const
Parabolic approximation version of local squared L2 error.
Definition: bilaplacian_ldgh.hxx:894
LocalSolver::Bilaplacian::n_glob_dofs_per_node
static constexpr unsigned int n_glob_dofs_per_node()
Evaluate amount of global degrees of freedom per hypernode.
Definition: bilaplacian_ldgh.hxx:192
LocalSolver::Bilaplacian_parameters_default::analytic_result
static param_float_t analytic_result(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Analytic result of PDE (for convergence tests).
Definition: bilaplacian_ldgh.hxx:89
LocalSolver::Bilaplacian::bulk_values
std::array< std::array< lSol_float_t, Hypercube< hyEdge_dimT >::pow(sizeT)>, system_dim > bulk_values(const std::array< abscissa_float_t, sizeT > &abscissas, const input_array_t &lambda_values, hyEdgeT &hyper_edge, const lSol_float_t time=0.) const
Evaluate local reconstruction at tensorial products of abscissas.
LocalSolver::Bilaplacian
Local solver for bilaplacian equation on uniform hypergraph.
Definition: bilaplacian_ldgh.hxx:120
parameters
the intent is to exercise the right to control the distribution of derivative or collective works based on the Library In mere aggregation of another work not based on the Library with the you must alter all the notices that refer to this so that they refer to the ordinary GNU General Public instead of to this it is irreversible for that so the ordinary GNU General Public License applies to all subsequent copies and derivative works made from that copy This option is useful when you wish to copy part of the code of the Library into a program that is not a library You may copy and distribute the which must be distributed under the terms of Sections and above on a medium customarily used for software interchange If distribution of object code is made by offering access to copy from a designated then offering equivalent access to copy the source code from the same place satisfies the requirement to distribute the source even though third parties are not compelled to copy the source along with the object code A program that contains no derivative of any portion of the but is designed to work with the Library by being compiled or linked with is called a work that uses the Library Such a in is not a derivative work of the and therefore falls outside the scope of this License linking a work that uses the Library with the Library creates an executable that is a derivative of the rather than a work that uses the library The executable is therefore covered by this License Section states terms for distribution of such executables When a work that uses the Library uses material from a header file that is part of the the object code for the work may be a derivative work of the Library even though the source code is not Whether this is true is especially significant if the work can be linked without the or if the work is itself a library The threshold for this to be true is not precisely defined by law If such an object file uses only numerical parameters
Definition: License.txt:259
LocalSolver::Bilaplacian_parameters_default::neumann_laplace_value
static param_float_t neumann_laplace_value(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Neumann values of solution's Laplacian as analytic function.
Definition: bilaplacian_ldgh.hxx:81
LocalSolver::Bilaplacian::primal_at_boundary
std::array< std::array< lSol_float_t, 2 *n_shape_bdr_ >, 2 *hyEdge_dimT > primal_at_boundary(const std::array< lSol_float_t, n_loc_dofs_ > &coeffs, hyEdgeT &hyper_edge) const
Evaluate primal variable at boundary.
hy_assert
#define hy_assert(Expr, Msg)
The assertion to be used within HyperHDG — deactivate using -DNDEBUG compile flag.
Definition: hy_assert.hxx:38
dense_la.hxx
LocalSolver::Bilaplacian::assemble_loc_matrix
SmallSquareMat< n_loc_dofs_, lSol_float_t > assemble_loc_matrix(const lSol_float_t tau, hyEdgeT &hyper_edge, const lSol_float_t time) const
Assemble local matrix for the local solver.
LocalSolver::Bilaplacian_parameters_default::neumann_value
static param_float_t neumann_value(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Neumann values of solution as analytic function.
Definition: bilaplacian_ldgh.hxx:73
LocalSolver::Bilaplacian::assemble_rhs_from_lambda
SmallVec< n_loc_dofs_, lSol_float_t > assemble_rhs_from_lambda(const SmallMatT &lambda_values, hyEdgeT &hyper_edge) const
Assemble local right-hand for the local solver (from skeletal).
LocalSolver::Bilaplacian::total_numerical_flux_mass
SmallMatOutT & total_numerical_flux_mass(const SmallMatInT &lambda_values_in, SmallMatOutT &lambda_values_out, hyEdgeT &hyper_edge, const lSol_float_t time=0.) const
Evaluate local contribution to mass matrix–vector residual.
Definition: bilaplacian_ldgh.hxx:729
LocalSolver::Bilaplacian::node_element::functions
std::tuple< TPP::ShapeFunction< TPP::ShapeType::Tensorial< TPP::ShapeType::Legendre< poly_deg >, hyEdge_dimT - 1 > > > functions
Definition: bilaplacian_ldgh.hxx:136
LocalSolver::Bilaplacian_parameters_default::right_hand_side
static param_float_t right_hand_side(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Right-hand side in PDE as analytic function.
Definition: bilaplacian_ldgh.hxx:49
LocalSolver::Bilaplacian::Bilaplacian
Bilaplacian(const constructor_value_type &tau=1.)
Constructor for local solver.
Definition: bilaplacian_ldgh.hxx:490
LocalSolver::Bilaplacian::error_def::postprocess_error
static error_t postprocess_error(error_t &summed_error)
Define how global errors should be postprocessed.
Definition: bilaplacian_ldgh.hxx:168
SmallMat
This class implements a small/dense matrix.
Definition: dense_la.hxx:34
LocalSolver::Bilaplacian::error_def
Define how errors are evaluated.
Definition: bilaplacian_ldgh.hxx:141
LocalSolver::Bilaplacian::solve_local_problem
std::array< lSol_float_t, n_loc_dofs_ > solve_local_problem(const SmallMatT &lambda_values, const unsigned int solution_type, hyEdgeT &hyper_edge, const lSol_float_t time) const
Solve local problem (with right-hand side from skeletal).
Definition: bilaplacian_ldgh.hxx:364
Hypercube
Helper class containing numbers and functions related to hypercubes.
Definition: hypercube.hxx:12
Generated by   doxygen 1.8.17