HyperHDG
diffusion_ip.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 DiffusionParametersDefault
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 Neumann boundary.
28  ************************************************************************************************/
29  static constexpr std::array<unsigned int, 0U> neumann_nodes{};
30  /*!***********************************************************************************************
31  * \brief Diffusion coefficient in PDE as analytic function.
32  ************************************************************************************************/
33  static param_float_t diffusion_coeff(const Point<space_dimT, param_float_t>&,
34  const param_float_t = 0.)
35  {
36  return 1.;
37  }
38  /*!***********************************************************************************************
39  * \brief Right-hand side in PDE as analytic function.
40  ************************************************************************************************/
41  static param_float_t right_hand_side(const Point<space_dimT, param_float_t>&,
42  const param_float_t = 0.)
43  {
44  return 0.;
45  }
46  /*!***********************************************************************************************
47  * \brief Dirichlet values of solution as analytic function.
48  ************************************************************************************************/
49  static param_float_t dirichlet_value(const Point<space_dimT, param_float_t>&,
50  const param_float_t = 0.)
51  {
52  return 0.;
53  }
54  /*!***********************************************************************************************
55  * \brief Neumann values of solution as analytic function.
56  ************************************************************************************************/
57  static param_float_t neumann_value(const Point<space_dimT, param_float_t>&,
58  const param_float_t = 0.)
59  {
60  return 0.;
61  }
62  /*!***********************************************************************************************
63  * \brief Analytic result of PDE (for convergence tests).
64  ************************************************************************************************/
65  static param_float_t analytic_result(const Point<space_dimT, param_float_t>&,
66  const param_float_t = 0.)
67  {
68  return 0.;
69  }
70 }; // end of struct DiffusionParametersDefault
71 
72 /*!*************************************************************************************************
73  * \brief Local solver diffusion equation on hypergraph.
74  *
75  * This class contains the local solver for an isotropic diffusion equation, i.e.,
76  * \f[
77  * - \nabla \cdot ( d \nabla u ) = f \quad \text{ in } \Omega, \qquad
78  * u = u_\text D \quad \text{ on } \partial \Omega_\text D, \qquad
79  * - d \nabla u \cdot \nu = g_\text N \quad \text{ on } \partial \Omega_\text N
80  * \f]
81  * in a spatial domain \f$\Omega \subset \mathbb R^d\f$. Here, \f$d\f$ is the spatial dimension
82  * \c space_dim, \f$\Omega\f$ is a regular graph (\c hyEdge_dimT = 1) or hypergraph whose
83  * hyperedges are surfaces (\c hyEdge_dimT = 2) or volumes (\c hyEdge_dimT = 3) or hypervolumes (in
84  * case of \c hyEdge_dimT > 3). \f$f\f$ and \f$d\f$ are scalars defined in the whole domain, the
85  * Dirichlet and Neumann boundary data needs to be defined on their respective hypernodes.
86  *
87  * Moreover, this class provides the functions that approximate the condensed mass matrix. By this,
88  * it is possible to approximate parabolic problems (in a very bad way) and to find good starting
89  * values for the nonlinear eigenvalue problem.
90  *
91  * \tparam hyEdge_dimT Dimension of a hyperedge, i.e., 1 is for PDEs defined on graphs, 2 is for
92  * PDEs defined on surfaces, and 3 is for PDEs defined on volumes.
93  * \tparam poly_deg The polynomial degree of test and trial functions.
94  * \tparam quad_deg The order of the quadrature rule.
95  * \tparam parametersT Struct depending on templates \c space_dimTP and \c lSol_float_TP that
96  * contains static parameter functions.
97  * Defaults to above functions included in \c DiffusionParametersDefault.
98  * \tparam lSol_float_t The floating point type calculations are executed in. Defaults to double.
99  *
100  * \authors Guido Kanschat, Heidelberg University, 2019--2020.
101  * \authors Andreas Rupp, Heidelberg University, 2019--2020.
102  **************************************************************************************************/
103 template <unsigned int hyEdge_dimT,
104  unsigned int poly_deg,
105  unsigned int quad_deg,
106  template <unsigned int, typename> typename parametersT = DiffusionParametersDefault,
107  typename lSol_float_t = double>
108 class DiffusionIP
109 {
110  public:
111  /*!***********************************************************************************************
112  * \brief Define type of (hyperedge related) data that is stored in HyDataContainer.
113  ************************************************************************************************/
114  struct data_type
115  {
116  };
117  /*!***********************************************************************************************
118  * \brief Define type of node elements, especially with respect to nodal shape functions.
119  ************************************************************************************************/
120  struct node_element
121  {
122  typedef std::tuple<TPP::ShapeFunction<
123  TPP::ShapeType::Tensorial<TPP::ShapeType::Legendre<poly_deg>, hyEdge_dimT - 1> > >
124  functions;
125  };
126  /*!***********************************************************************************************
127  * \brief Define how errors are evaluated.
128  ************************************************************************************************/
129  struct error_def
130  {
131  /*!*********************************************************************************************
132  * \brief Define the typename returned by function errors.
133  **********************************************************************************************/
134  typedef std::array<lSol_float_t, 1U> error_t;
135  /*!*********************************************************************************************
136  * \brief Define how initial error is generated.
137  **********************************************************************************************/
138  static error_t initial_error()
139  {
140  std::array<lSol_float_t, 1U> summed_error;
141  summed_error.fill(0.);
142  return summed_error;
143  }
144  /*!*********************************************************************************************
145  * \brief Define how local errors should be accumulated.
146  **********************************************************************************************/
147  static error_t sum_error(error_t& summed_error, const error_t& new_error)
148  {
149  for (unsigned int k = 0; k < summed_error.size(); ++k)
150  summed_error[k] += new_error[k];
151  return summed_error;
152  }
153  /*!*********************************************************************************************
154  * \brief Define how global errors should be postprocessed.
155  **********************************************************************************************/
156  static error_t postprocess_error(error_t& summed_error)
157  {
158  for (unsigned int k = 0; k < summed_error.size(); ++k)
159  summed_error[k] = std::sqrt(summed_error[k]);
160  return summed_error;
161  }
162  };
163 
164  // -----------------------------------------------------------------------------------------------
165  // Public, static constexpr functions
166  // -----------------------------------------------------------------------------------------------
167 
168  /*!***********************************************************************************************
169  * \brief Dimension of hyper edge type that this object solves on.
170  ************************************************************************************************/
171  static constexpr unsigned int hyEdge_dim() { return hyEdge_dimT; }
172  /*!***********************************************************************************************
173  * \brief Evaluate amount of global degrees of freedom per hypernode.
174  *
175  * This number must be equal to HyperNodeFactory::n_dofs_per_node() of the HyperNodeFactory
176  * cooperating with this object.
177  *
178  * \retval n_dofs Number of global degrees of freedom per hypernode.
179  ************************************************************************************************/
180  static constexpr unsigned int n_glob_dofs_per_node()
181  {
182  return Hypercube<hyEdge_dimT - 1>::pow(poly_deg + 1);
183  }
184  /*!***********************************************************************************************
185  * \brief Dimension of of the solution evaluated with respect to a hyperedge.
186  ************************************************************************************************/
187  static constexpr unsigned int system_dimension() { return 1; }
188  /*!***********************************************************************************************
189  * \brief Dimension of of the solution evaluated with respect to a hypernode.
190  ************************************************************************************************/
191  static constexpr unsigned int node_system_dimension() { return 1; }
192 
193  private:
194  // -----------------------------------------------------------------------------------------------
195  // Private, static constexpr functions
196  // -----------------------------------------------------------------------------------------------
197 
198  /*!***********************************************************************************************
199  * \brief Number of local shape functions (with respect to all spatial dimensions).
200  ************************************************************************************************/
201  static constexpr unsigned int n_shape_fct_ = n_glob_dofs_per_node() * (poly_deg + 1);
202  /*!***********************************************************************************************
203  * \brief Number oflocal shape functions (with respect to a face / hypernode).
204  ************************************************************************************************/
205  static constexpr unsigned int n_shape_bdr_ = n_glob_dofs_per_node();
206  /*!***********************************************************************************************
207  * \brief Number of (local) degrees of freedom per hyperedge.
208  ************************************************************************************************/
209  static constexpr unsigned int n_loc_dofs_ = n_shape_fct_;
210  /*!***********************************************************************************************
211  * \brief Dimension of of the solution evaluated with respect to a hypernode.
212  *
213  * This allows to the use of this quantity as template parameter in member functions.
214  ************************************************************************************************/
215  static constexpr unsigned int system_dim = system_dimension();
216  /*!***********************************************************************************************
217  * \brief Dimension of of the solution evaluated with respect to a hypernode.
218  *
219  * This allows to the use of this quantity as template parameter in member functions.
220  ************************************************************************************************/
221  static constexpr unsigned int node_system_dim = node_system_dimension();
222  /*!***********************************************************************************************
223  * \brief Find out whether a node is of Dirichlet type.
224  ************************************************************************************************/
225  template <typename parameters>
226  static constexpr bool is_dirichlet(const unsigned int node_type)
227  {
228  return std::find(parameters::dirichlet_nodes.begin(), parameters::dirichlet_nodes.end(),
229  node_type) != parameters::dirichlet_nodes.end();
230  }
231 
232  // -----------------------------------------------------------------------------------------------
233  // Private, const members: Parameters and auxiliaries that help assembling matrices, etc.
234  // -----------------------------------------------------------------------------------------------
235 
236  /*!***********************************************************************************************
237  * \brief (Globally constant) penalty parameter for HDG scheme.
238  ************************************************************************************************/
239  const lSol_float_t tau_;
240  /*!***********************************************************************************************
241  * \brief An integrator helps to easily evaluate integrals (e.g. via quadrature).
242  ************************************************************************************************/
243  typedef TPP::Quadrature::Tensorial<
244  TPP::Quadrature::GaussLegendre<quad_deg>,
245  TPP::ShapeFunction<TPP::ShapeType::Tensorial<TPP::ShapeType::Legendre<poly_deg>, hyEdge_dimT> >,
246  lSol_float_t>
247  integrator;
248 
249  // -----------------------------------------------------------------------------------------------
250  // Private, internal functions for the local solver
251  // -----------------------------------------------------------------------------------------------
252 
253  /*!***********************************************************************************************
254  * \brief Assemble local matrix for the local solver.
255  *
256  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
257  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
258  * \param time Time, when the local matrix is evaluated.
259  * \retval loc_mat Matrix of the local solver.
260  ************************************************************************************************/
261  template <typename hyEdgeT>
262  inline SmallSquareMat<n_loc_dofs_, lSol_float_t> assemble_loc_matrix(
263  hyEdgeT& hyper_edge,
264  const lSol_float_t time) const;
265  /*!***********************************************************************************************
266  * \brief Assemble local right-hand for the local solver (from skeletal).
267  *
268  * The right hand side needs the values of the global degrees of freedom. Note that we basically
269  * follow the lines of
270  *
271  * B. Cockburn, J. Gopalakrishnan, and R. Lazarov.
272  * Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for
273  * second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009
274  *
275  * and discriminate between local solvers with respect to the skeletal variable and with respect
276  * to the global right-hand side. This assembles the local right-hand side with respect to the
277  * skeletal variable.
278  *
279  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
280  * \tparam SmallMatT The data type of the \c lambda values.
281  * \param lambda_values Global degrees of freedom associated to the hyperedge.
282  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
283  * \retval loc_rhs Local right hand side of the locasl solver.
284  ************************************************************************************************/
285  template <typename hyEdgeT, typename SmallMatT>
286  inline SmallVec<n_loc_dofs_, lSol_float_t> assemble_rhs_from_lambda(
287  const SmallMatT& lambda_values,
288  hyEdgeT& hyper_edge) const;
289  /*!***********************************************************************************************
290  * \brief Assemble local right-hand for the local solver (from global right-hand side).
291  *
292  * Note that we basically follow the lines of
293  *
294  * B. Cockburn, J. Gopalakrishnan, and R. Lazarov.
295  * Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for
296  * second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009
297  *
298  * and discriminate between local solvers with respect to the skeletal variable and with respect
299  * to the global right-hand side. This assembles the local right-hand side with respect to the
300  * global right-hand side. This function implicitly uses the parameters.
301  *
302  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
303  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
304  * \param time Point of time, rhs is evaluated at
305  * \retval loc_rhs Local right hand side of the locasl solver.
306  ************************************************************************************************/
307  template <typename hyEdgeT>
308  inline SmallVec<n_loc_dofs_, lSol_float_t> assemble_rhs_from_global_rhs(
309  hyEdgeT& hyper_edge,
310  const lSol_float_t time) const;
311  /*!***********************************************************************************************
312  * \brief Assemble local right-hand for the local solver (from volume function coefficients).
313  *
314  * Note that we basically follow the lines of
315  *
316  * B. Cockburn, J. Gopalakrishnan, and R. Lazarov.
317  * Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for
318  * second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009
319  *
320  * and discriminate between local solvers with respect to the skeletal variable and with respect
321  * to the global right-hand side. This assembles the local right-hand side with respect to the
322  * global right-hand side. This function implicitly uses the parameters.
323  *
324  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
325  * \tparam SmallVecT The data type of the coefficients.
326  * \param coeffs Local coefficients of bulk variable u.
327  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
328  * \retval loc_rhs Local right hand side of the locasl solver.
329  ************************************************************************************************/
330  template <typename hyEdgeT, typename SmallVecT>
331  inline SmallVec<n_loc_dofs_, lSol_float_t> assemble_rhs_from_coeffs(const SmallVecT& coeffs,
332  hyEdgeT& hyper_edge) const;
333  /*!***********************************************************************************************
334  * \brief Solve local problem (with right-hand side from skeletal).
335  *
336  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
337  * \tparam SmallMatT The data type of \c lambda_values.
338  * \param lambda_values Global degrees of freedom associated to the hyperedge.
339  * \param solution_type Type of local problem that is to be solved.
340  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
341  * \param time Point of time the problem is solved.
342  * \retval loc_sol Solution of the local problem.
343  ************************************************************************************************/
344  template <typename hyEdgeT, typename SmallMatT>
345  inline SmallVec<n_loc_dofs_, lSol_float_t> solve_local_problem(const SmallMatT& lambda_values,
346  const unsigned int solution_type,
347  hyEdgeT& hyper_edge,
348  const lSol_float_t time) const
349  {
350  try
351  {
352  SmallVec<n_loc_dofs_, lSol_float_t> rhs;
353  if (solution_type == 0)
354  rhs = assemble_rhs_from_lambda(lambda_values, hyper_edge);
355  else if (solution_type == 1)
356  rhs = assemble_rhs_from_lambda(lambda_values, hyper_edge) +
357  assemble_rhs_from_global_rhs(hyper_edge, time);
358  else
359  hy_assert(0 == 1, "This has not been implemented!");
360  return rhs / assemble_loc_matrix(hyper_edge, time);
361  }
362  catch (Wrapper::LAPACKexception& exc)
363  {
364  hy_assert(0 == 1, exc.what() << std::endl
365  << "This can happen if quadrature is too inaccurate!");
366  throw exc;
367  }
368  }
369  /*!***********************************************************************************************
370  * \brief Evaluate primal variable at boundary.
371  *
372  * Function to evaluate primal variable of the solution. This function is needed to calculate
373  * the local numerical fluxes.
374  *
375  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
376  * \param coeffs Coefficients of the local solution.
377  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
378  * \retval bdr_coeffs Coefficients of respective (dim-1) dimensional function at boundaries.
379  ************************************************************************************************/
380  template <typename hyEdgeT>
381  inline SmallMat<2 * hyEdge_dimT, n_shape_bdr_, lSol_float_t> primal_at_boundary(
382  const SmallVec<n_loc_dofs_, lSol_float_t>& coeffs,
383  hyEdgeT& hyper_edge) const;
384  /*!***********************************************************************************************
385  * \brief Evaluate dual variable at boundary.
386  *
387  * Function to evaluate dual variable of the solution. This function is needed to calculate the
388  * local numerical fluxes.
389  *
390  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
391  * \param coeffs Coefficients of the local solution.
392  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
393  * \retval bdr_coeffs Coefficients of respective (dim-1) dimensional function at boundaries.
394  ************************************************************************************************/
395  template <typename hyEdgeT>
396  inline SmallMat<2 * hyEdge_dimT, n_shape_bdr_, lSol_float_t> dual_at_boundary(
397  const SmallVec<n_loc_dofs_, lSol_float_t>& coeffs,
398  hyEdgeT& hyper_edge) const;
399 
400  public:
401  // -----------------------------------------------------------------------------------------------
402  // Public functions (and one typedef) to be utilized by external functions.
403  // -----------------------------------------------------------------------------------------------
404 
405  /*!***********************************************************************************************
406  * \brief Class is constructed using a single double indicating the penalty parameter.
407  ************************************************************************************************/
408  typedef lSol_float_t constructor_value_type;
409  /*!***********************************************************************************************
410  * \brief Constructor for local solver.
411  *
412  * \param tau Penalty parameter of HDG scheme.
413  ************************************************************************************************/
414  DiffusionIP(const constructor_value_type& tau = 1.) : tau_(tau) {}
415  /*!***********************************************************************************************
416  * \brief Evaluate local contribution to matrix--vector multiplication.
417  *
418  * Execute matrix--vector multiplication y = A * x, where x represents the vector containing the
419  * skeletal variable (adjacent to one hyperedge), and A is the condensed matrix arising from the
420  * HDG discretization. This function does this multiplication (locally) for one hyperedge. The
421  * hyperedge is no parameter, since all hyperedges are assumed to have the same properties.
422  *
423  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
424  * \tparam SmallMatInT Data type of \c lambda_values_in.
425  * \tparam SmallMatOutT Data type of \c lambda_values_out.
426  * \param lambda_values_in Local part of vector x.
427  * \param lambda_values_out Local part that will be added to A * x.
428  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
429  * \param time Time.
430  * \retval vecAx Local part of vector A * x.
431  ************************************************************************************************/
432  template <typename hyEdgeT, typename SmallMatInT, typename SmallMatOutT>
433  SmallMatOutT& trace_to_flux(const SmallMatInT& lambda_values_in,
434  SmallMatOutT& lambda_values_out,
435  hyEdgeT& hyper_edge,
436  const lSol_float_t time = 0.) const
437  {
438  hy_assert(lambda_values_in.size() == lambda_values_out.size() &&
439  lambda_values_in.size() == 2 * hyEdge_dimT,
440  "Both matrices must be of same size which corresponds to the number of faces!");
441  for (unsigned int i = 0; i < lambda_values_in.size(); ++i)
442  hy_assert(
443  lambda_values_in[i].size() == lambda_values_out[i].size() &&
444  lambda_values_in[i].size() == n_glob_dofs_per_node(),
445  "Both matrices must be of same size which corresponds to the number of dofs per face!");
446 
447  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
448  SmallVec<n_loc_dofs_, lSol_float_t> coeffs =
449  solve_local_problem(lambda_values_in, 0U, hyper_edge, time);
450 
451  SmallMat<2 * hyEdge_dimT, n_shape_bdr_, lSol_float_t> primals(
452  primal_at_boundary(coeffs, hyper_edge)),
453  duals(dual_at_boundary(coeffs, hyper_edge));
454 
455  for (unsigned int i = 0; i < lambda_values_out.size(); ++i)
456  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
457  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
458  lambda_values_out[i][j] = 0.;
459  else
460  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
461  lambda_values_out[i][j] =
462  duals(i, j) + tau_ * primals(i, j) -
463  tau_ * lambda_values_in[i][j] * hyper_edge.geometry.face_area(i);
464 
465  return lambda_values_out;
466  }
467  /*!***********************************************************************************************
468  * \brief Fill an array with 1 if the node is Dirichlet and 0 otherwise.
469  ************************************************************************************************/
470  template <typename hyEdgeT>
471  std::array<unsigned int, 2 * hyEdge_dimT> node_types(hyEdgeT& hyper_edge) const
472  {
473  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
474 
475  std::array<unsigned int, 2 * hyEdge_dimT> result;
476  result.fill(0);
477 
478  for (unsigned int i = 0; i < 2 * hyEdge_dimT; ++i)
479  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
480  result[i] = 1;
481 
482  return result;
483  }
484  /*!***********************************************************************************************
485  * \brief Evaluate local contribution to residual \f$ A x - b \f$.
486  *
487  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
488  * \tparam SmallMatInT Data type of \c lambda_values_in.
489  * \tparam SmallMatOutT Data type of \c lambda_values_out.
490  * \param lambda_values_in Local part of vector x.
491  * \param lambda_values_out Local part that will be added to A * x.
492  * \param hyper_edge The geometry of the considered hyperedge (of typename hyEdgeT).
493  * \param time Time at which analytic functions are evaluated.
494  * \retval vecAx Local part of vector A * x - b.
495  ************************************************************************************************/
496  template <typename hyEdgeT, typename SmallMatInT, typename SmallMatOutT>
497  SmallMatOutT& residual_flux(const SmallMatInT& lambda_values_in,
498  SmallMatOutT& lambda_values_out,
499  hyEdgeT& hyper_edge,
500  const lSol_float_t time = 0.) const
501  {
502  hy_assert(lambda_values_in.size() == lambda_values_out.size() &&
503  lambda_values_in.size() == 2 * hyEdge_dimT,
504  "Both matrices must be of same size which corresponds to the number of faces!");
505  for (unsigned int i = 0; i < lambda_values_in.size(); ++i)
506  hy_assert(
507  lambda_values_in[i].size() == lambda_values_out[i].size() &&
508  lambda_values_in[i].size() == n_glob_dofs_per_node(),
509  "Both matrices must be of same size which corresponds to the number of dofs per face!");
510 
511  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
512  SmallVec<n_loc_dofs_, lSol_float_t> coeffs =
513  solve_local_problem(lambda_values_in, 1U, hyper_edge, time);
514 
515  SmallMat<2 * hyEdge_dimT, n_shape_bdr_, lSol_float_t> primals(
516  primal_at_boundary(coeffs, hyper_edge)),
517  duals(dual_at_boundary(coeffs, hyper_edge));
518 
519  for (unsigned int i = 0; i < lambda_values_out.size(); ++i)
520  {
521  if (is_dirichlet<parameters>(hyper_edge.node_descriptor[i]))
522  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
523  lambda_values_out[i][j] = 0.;
524  else
525  for (unsigned int j = 0; j < lambda_values_out[i].size(); ++j)
526  lambda_values_out[i][j] =
527  duals(i, j) + tau_ * primals(i, j) -
528  tau_ * lambda_values_in[i][j] * hyper_edge.geometry.face_area(i);
529  }
530 
531  return lambda_values_out;
532  }
533  /*!***********************************************************************************************
534  * \brief Calculate squared local contribution of L2 error.
535  *
536  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
537  * \tparam SmallMatT Data type of \c lambda_values.
538  * \param lambda_values The values of the skeletal variable's coefficients.
539  * \param hy_edge The geometry of the considered hyperedge (of typename GeomT).
540  * \param time Time at which anaytic functions are evaluated.
541  * \retval vec_b Local part of vector b.
542  ************************************************************************************************/
543  template <typename hyEdgeT, typename SmallMatT>
544  std::array<lSol_float_t, 1U> errors(const SmallMatT& lambda_values,
545  hyEdgeT& hy_edge,
546  const lSol_float_t time = 0.) const
547  {
548  hy_assert(lambda_values.size() == 2 * hyEdge_dimT, "Matrix must have appropriate size!");
549  for (unsigned int i = 0; i < lambda_values.size(); ++i)
550  hy_assert(lambda_values[i].size() == n_glob_dofs_per_node(),
551  "Matrix must have appropriate size!");
552 
553  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
554  SmallVec<n_loc_dofs_, lSol_float_t> coefficients =
555  solve_local_problem(lambda_values, 1U, hy_edge, time);
556  return std::array<lSol_float_t, 1U>({integrator::template integrate_vol_diffsquare_discana<
557  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>, decltype(hyEdgeT::geometry),
558  parameters::analytic_result, Point<hyEdge_dimT, lSol_float_t> >(coefficients.data(),
559  hy_edge.geometry, time)});
560  }
561  /*!***********************************************************************************************
562  * \brief Evaluate local reconstruction at tensorial products of abscissas.
563  *
564  * \tparam absc_float_t Floating type for the abscissa values.
565  * \tparam abscissas_sizeT Size of the array of array of abscissas.
566  * \tparam input_array_t Type of input array.
567  * \tparam hyEdgeT The geometry type / typename of the considered hyEdge's geometry.
568  * \param abscissas Abscissas of the supporting points.
569  * \param lambda_values The values of the skeletal variable's coefficients.
570  * \param hyper_edge The geometry of the considered hyperedge (of typename GeomT).
571  * \param time Time at which analytic functions are evaluated.
572  * \retval function_values Function values at quadrature points.
573  ************************************************************************************************/
574  template <typename abscissa_float_t,
575  std::size_t abscissas_sizeT,
576  class input_array_t,
577  class hyEdgeT>
578  std::array<std::array<lSol_float_t, Hypercube<hyEdge_dimT>::pow(abscissas_sizeT)>, system_dim>
579  bulk_values(const std::array<abscissa_float_t, abscissas_sizeT>& abscissas,
580  const input_array_t& lambda_values,
581  hyEdgeT& hyper_edge,
582  const lSol_float_t time = 0.) const;
583 }; // end of class DiffusionIP
584 
585 // -------------------------------------------------------------------------------------------------
587 // -------------------------------------------------------------------------------------------------
588 
589 // -------------------------------------------------------------------------------------------------
590 // -------------------------------------------------------------------------------------------------
591 //
592 // IMPLEMENTATION OF MEMBER FUNCTIONS OF Diffusion
593 //
594 // -------------------------------------------------------------------------------------------------
595 // -------------------------------------------------------------------------------------------------
596 
597 // -------------------------------------------------------------------------------------------------
598 // assemble_loc_matrix
599 // -------------------------------------------------------------------------------------------------
600 
601 template <unsigned int hyEdge_dimT,
602  unsigned int poly_deg,
603  unsigned int quad_deg,
604  template <unsigned int, typename>
605  typename parametersT,
606  typename lSol_float_t>
607 template <typename hyEdgeT>
608 inline SmallSquareMat<
609  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
610  lSol_float_t>
611 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::assemble_loc_matrix(
612  hyEdgeT& hyper_edge,
613  const lSol_float_t time) const
614 {
615  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
616  SmallSquareMat<n_loc_dofs_, lSol_float_t> local_mat;
617 
618  for (unsigned int i = 0; i < n_shape_fct_; ++i)
619  {
620  for (unsigned int j = 0; j < n_shape_fct_; ++j)
621  {
622  // Integral_element phi_i phi_j dx in diagonal blocks
623  local_mat(i, j) += integrator::template integrate_vol_nablaphinablaphifunc<
624  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>, decltype(hyEdgeT::geometry),
625  parameters::diffusion_coeff, Point<hyEdge_dimT, lSol_float_t> >(i, j, hyper_edge.geometry,
626  time);
627 
628  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
629  {
630  local_mat(i, j) +=
631  tau_ * integrator::template integrate_bdr_phiphi<decltype(hyEdgeT::geometry)>(
632  i, j, face, hyper_edge.geometry);
633  local_mat(i, j) -= integrator::template integrate_bdr_nablaphiphinufunc<
634  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>,
635  decltype(hyEdgeT::geometry), parameters::diffusion_coeff,
636  Point<hyEdge_dimT, lSol_float_t> >(j, i, face, hyper_edge.geometry);
637  }
638  }
639  }
640 
641  return local_mat;
642 } // end of Diffusion::assemble_loc_matrix
643 
644 // -------------------------------------------------------------------------------------------------
645 // assemble_rhs_from_lambda
646 // -------------------------------------------------------------------------------------------------
647 
648 template <unsigned int hyEdge_dimT,
649  unsigned int poly_deg,
650  unsigned int quad_deg,
651  template <unsigned int, typename>
652  typename parametersT,
653  typename lSol_float_t>
654 template <typename hyEdgeT, typename SmallMatT>
655 inline SmallVec<
656  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
657  lSol_float_t>
658 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::assemble_rhs_from_lambda(
659  const SmallMatT& lambda_values,
660  hyEdgeT& hyper_edge) const
661 {
662  static_assert(std::is_same<typename SmallMatT::value_type::value_type, lSol_float_t>::value,
663  "Lambda values should have same floating point arithmetics as local solver!");
664  hy_assert(lambda_values.size() == 2 * hyEdge_dimT,
665  "The size of the lambda values should be twice the dimension of a hyperedge.");
666  for (unsigned int i = 0; i < 2 * hyEdge_dimT; ++i)
667  hy_assert(lambda_values[i].size() == n_shape_bdr_,
668  "The size of lambda should be the amount of ansatz functions at boundary.");
669 
670  SmallVec<n_loc_dofs_, lSol_float_t> right_hand_side;
671 
672  for (unsigned int i = 0; i < n_shape_fct_; ++i)
673  for (unsigned int j = 0; j < n_shape_bdr_; ++j)
674  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
675  right_hand_side[i] +=
676  tau_ * lambda_values[face][j] *
677  integrator::template integrate_bdr_phipsi<decltype(hyEdgeT::geometry)>(
678  i, j, face, hyper_edge.geometry);
679 
680  return right_hand_side;
681 } // end of Diffusion::assemble_rhs_from_lambda
682 
683 // -------------------------------------------------------------------------------------------------
684 // assemble_rhs_from_global_rhs
685 // -------------------------------------------------------------------------------------------------
686 
687 template <unsigned int hyEdge_dimT,
688  unsigned int poly_deg,
689  unsigned int quad_deg,
690  template <unsigned int, typename>
691  typename parametersT,
692  typename lSol_float_t>
693 template <typename hyEdgeT>
694 inline SmallVec<
695  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
696  lSol_float_t>
697 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::
698  assemble_rhs_from_global_rhs(hyEdgeT& hyper_edge, const lSol_float_t time) const
699 {
700  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
701  SmallVec<n_loc_dofs_, lSol_float_t> right_hand_side;
702 
703  for (unsigned int i = 0; i < n_shape_fct_; ++i)
704  {
705  right_hand_side[i] = integrator::template integrate_vol_phifunc<
706  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>, decltype(hyEdgeT::geometry),
707  parameters::right_hand_side, Point<hyEdge_dimT, lSol_float_t> >(i, hyper_edge.geometry, time);
708  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
709  {
710  if (!is_dirichlet<parameters>(hyper_edge.node_descriptor[face]))
711  continue;
712  right_hand_side[i] +=
713  tau_ * integrator::template integrate_bdr_phifunc<
714  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>,
715  decltype(hyEdgeT::geometry), parameters::dirichlet_value,
716  Point<hyEdge_dimT, lSol_float_t> >(i, face, hyper_edge.geometry, time);
717  }
718  }
719 
720  return right_hand_side;
721 } // end of Diffusion::assemble_rhs_from_global_rhs
722 
723 // -------------------------------------------------------------------------------------------------
724 // assemble_rhs_from_coeffs
725 // -------------------------------------------------------------------------------------------------
726 
727 template <unsigned int hyEdge_dimT,
728  unsigned int poly_deg,
729  unsigned int quad_deg,
730  template <unsigned int, typename>
731  typename parametersT,
732  typename lSol_float_t>
733 template <typename hyEdgeT, typename SmallVecT>
734 inline SmallVec<
735  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_loc_dofs_,
736  lSol_float_t>
737 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::assemble_rhs_from_coeffs(
738  const SmallVecT& coeffs,
739  hyEdgeT& hyper_edge) const
740 {
741  SmallVec<n_loc_dofs_, lSol_float_t> right_hand_side;
742  for (unsigned int i = 0; i < n_shape_fct_; ++i)
743  for (unsigned int j = 0; j < n_shape_fct_; ++j)
744  right_hand_side[hyEdge_dimT * n_shape_fct_ + i] +=
745  coeffs[j] * integrator::template integrate_vol_phiphi(i, j, hyper_edge.geometry);
746 
747  return right_hand_side;
748 } // end of Diffusion::assemble_rhs_from_coeffs
749 
750 // -------------------------------------------------------------------------------------------------
751 // primal_at_boundary
752 // -------------------------------------------------------------------------------------------------
753 
754 template <unsigned int hyEdge_dimT,
755  unsigned int poly_deg,
756  unsigned int quad_deg,
757  template <unsigned int, typename>
758  typename parametersT,
759  typename lSol_float_t>
760 template <typename hyEdgeT>
761 inline SmallMat<
762  2 * hyEdge_dimT,
763  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_shape_bdr_,
764  lSol_float_t>
765 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::primal_at_boundary(
766  const SmallVec<n_loc_dofs_, lSol_float_t>& coeffs,
767  hyEdgeT& hyper_edge) const
768 {
769  SmallMat<2 * hyEdge_dimT, n_shape_bdr_, lSol_float_t> bdr_values;
770 
771  for (unsigned int i = 0; i < n_shape_fct_; ++i)
772  for (unsigned int j = 0; j < n_shape_bdr_; ++j)
773  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
774  bdr_values(face, j) +=
775  coeffs[i] * integrator::template integrate_bdr_phipsi<decltype(hyEdgeT::geometry)>(
776  i, j, face, hyper_edge.geometry);
777 
778  return bdr_values;
779 } // end of Diffusion::primal_at_boundary
780 
781 // -------------------------------------------------------------------------------------------------
782 // dual_at_boundary
783 // -------------------------------------------------------------------------------------------------
784 
785 template <unsigned int hyEdge_dimT,
786  unsigned int poly_deg,
787  unsigned int quad_deg,
788  template <unsigned int, typename>
789  typename parametersT,
790  typename lSol_float_t>
791 template <typename hyEdgeT>
792 inline SmallMat<
793  2 * hyEdge_dimT,
794  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::n_shape_bdr_,
795  lSol_float_t>
796 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::dual_at_boundary(
797  const SmallVec<n_loc_dofs_, lSol_float_t>& coeffs,
798  hyEdgeT& hyper_edge) const
799 {
800  using parameters = parametersT<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>;
801  SmallMat<2 * hyEdge_dimT, n_shape_bdr_, lSol_float_t> bdr_values;
802 
803  for (unsigned int i = 0; i < n_shape_fct_; ++i)
804  for (unsigned int j = 0; j < n_shape_bdr_; ++j)
805  for (unsigned int face = 0; face < 2 * hyEdge_dimT; ++face)
806  bdr_values(face, j) -=
807  integrator::template integrate_bdr_nablaphipsinufunc<
808  Point<decltype(hyEdgeT::geometry)::space_dim(), lSol_float_t>,
809  decltype(hyEdgeT::geometry), parameters::diffusion_coeff,
810  Point<hyEdge_dimT, lSol_float_t> >(i, j, face, hyper_edge.geometry) *
811  coeffs[i];
812 
813  return bdr_values;
814 } // end of Diffusion::dual_at_boundary
815 
816 // -------------------------------------------------------------------------------------------------
817 // bulk_values
818 // -------------------------------------------------------------------------------------------------
819 
820 template <unsigned int hyEdge_dimT,
821  unsigned int poly_deg,
822  unsigned int quad_deg,
823  template <unsigned int, typename>
824  typename parametersT,
825  typename lSol_float_t>
826 template <typename abscissa_float_t,
827  std::size_t abscissas_sizeT,
828  class input_array_t,
829  typename hyEdgeT>
830 std::array<std::array<lSol_float_t, Hypercube<hyEdge_dimT>::pow(abscissas_sizeT)>,
831  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::system_dim>
832 DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::bulk_values(
833  const std::array<abscissa_float_t, abscissas_sizeT>& abscissas,
834  const input_array_t& lambda_values,
835  hyEdgeT& hyper_edge,
836  const lSol_float_t time) const
837 {
838  SmallVec<n_loc_dofs_, lSol_float_t> coefficients =
839  solve_local_problem(lambda_values, 1U, hyper_edge, time);
840  SmallVec<n_shape_fct_, lSol_float_t> coeffs;
841  SmallVec<static_cast<unsigned int>(abscissas_sizeT), abscissa_float_t> helper(abscissas);
842 
843  std::array<std::array<lSol_float_t, Hypercube<hyEdge_dimT>::pow(abscissas_sizeT)>,
844  DiffusionIP<hyEdge_dimT, poly_deg, quad_deg, parametersT, lSol_float_t>::system_dim>
845  point_vals;
846 
847  for (unsigned int d = 0; d < system_dim; ++d)
848  {
849  for (unsigned int i = 0; i < coeffs.size(); ++i)
850  coeffs[i] = coefficients[d * n_shape_fct_ + i];
851  for (unsigned int pt = 0; pt < Hypercube<hyEdge_dimT>::pow(abscissas_sizeT); ++pt)
852  point_vals[d][pt] = integrator::shape_fun_t::template lin_comb_fct_val<float>(
853  coeffs, Hypercube<hyEdge_dimT>::template tensorial_pt<Point<hyEdge_dimT> >(pt, helper));
854  }
855 
856  return point_vals;
857 }
858 // end of Diffusion::bulk_values
859 
860 // -------------------------------------------------------------------------------------------------
862 // -------------------------------------------------------------------------------------------------
863 
864 } // namespace LocalSolver
LocalSolver::DiffusionIP::n_glob_dofs_per_node
static constexpr unsigned int n_glob_dofs_per_node()
Evaluate amount of global degrees of freedom per hypernode.
Definition: diffusion_ip.hxx:180
shape_function.hxx
LocalSolver::DiffusionIP::dual_at_boundary
SmallMat< 2 *hyEdge_dimT, n_shape_bdr_, lSol_float_t > dual_at_boundary(const SmallVec< n_loc_dofs_, lSol_float_t > &coeffs, hyEdgeT &hyper_edge) const
Evaluate dual variable at boundary.
LocalSolver::DiffusionIP::constructor_value_type
lSol_float_t constructor_value_type
Class is constructed using a single double indicating the penalty parameter.
Definition: diffusion_ip.hxx:408
Wrapper::LAPACKexception::what
const char * what() const
Definition: lapack.hxx:27
LocalSolver::DiffusionIP::error_def::postprocess_error
static error_t postprocess_error(error_t &summed_error)
Define how global errors should be postprocessed.
Definition: diffusion_ip.hxx:156
LocalSolver
A namespace for local solvers.
Definition: advection_parab_ldgh.hxx:11
tensorial.hxx
LocalSolver::DiffusionIP::data_type
Define type of (hyperedge related) data that is stored in HyDataContainer.
Definition: diffusion_ip.hxx:114
diffusion_uniform.geometry
geometry
Definition: diffusion_uniform.py:19
LocalSolver::DiffusionIP::assemble_rhs_from_coeffs
SmallVec< n_loc_dofs_, lSol_float_t > assemble_rhs_from_coeffs(const SmallVecT &coeffs, hyEdgeT &hyper_edge) const
Assemble local right-hand for the local solver (from volume function coefficients).
LocalSolver::DiffusionIP::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::DiffusionIP::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: diffusion_ip.hxx:471
LocalSolver::DiffusionIP::n_shape_fct_
static constexpr unsigned int n_shape_fct_
Number of local shape functions (with respect to all spatial dimensions).
Definition: diffusion_ip.hxx:201
LocalSolver::DiffusionParametersDefault::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: diffusion_ip.hxx:49
LocalSolver::DiffusionParametersDefault::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: diffusion_ip.hxx:57
LocalSolver::DiffusionIP::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: diffusion_ip.hxx:544
LocalSolver::DiffusionIP
Local solver diffusion equation on hypergraph.
Definition: diffusion_ip.hxx:108
LocalSolver::DiffusionParametersDefault::diffusion_coeff
static param_float_t diffusion_coeff(const Point< space_dimT, param_float_t > &, const param_float_t=0.)
Diffusion coefficient in PDE as analytic function.
Definition: diffusion_ip.hxx:33
TPP::Quadrature::GaussLegendre
Gauss–Legendre quadrature points and weights in one spatial dimension.
Definition: one_dimensional.hxx:19
SmallMat::size
static constexpr unsigned int size()
Return size a SmallMat.
Definition: dense_la.hxx:54
TPP::ShapeType::Tensorial
Struct that handles the evaluation of tensorial shape functions.
Definition: tensorial.hxx:22
LocalSolver::DiffusionIP::n_loc_dofs_
static constexpr unsigned int n_loc_dofs_
Number of (local) degrees of freedom per hyperedge.
Definition: diffusion_ip.hxx:209
LocalSolver::DiffusionIP::n_shape_bdr_
static constexpr unsigned int n_shape_bdr_
Number oflocal shape functions (with respect to a face / hypernode).
Definition: diffusion_ip.hxx:205
LocalSolver::DiffusionIP::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: diffusion_ip.hxx:247
LocalSolver::DiffusionIP::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::DiffusionIP::solve_local_problem
SmallVec< n_loc_dofs_, lSol_float_t > 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: diffusion_ip.hxx:345
LocalSolver::DiffusionIP::assemble_loc_matrix
SmallSquareMat< n_loc_dofs_, lSol_float_t > assemble_loc_matrix(hyEdgeT &hyper_edge, const lSol_float_t time) const
Assemble local matrix for the local solver.
TPP::ShapeFunction
Struct that handles different types of evaluation of shape functions.
Definition: shape_function.hxx:15
LocalSolver::DiffusionIP::error_def::error_t
std::array< lSol_float_t, 1U > error_t
Define the typename returned by function errors.
Definition: diffusion_ip.hxx:134
LocalSolver::DiffusionParametersDefault::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: diffusion_ip.hxx:65
LocalSolver::DiffusionIP::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: diffusion_ip.hxx:433
TPP::Quadrature::Tensorial
General integrator class on tensorial hypergraphs.
Definition: tensorial.hxx:24
LocalSolver::DiffusionIP::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: diffusion_ip.hxx:147
LocalSolver::DiffusionIP::error_def
Define how errors are evaluated.
Definition: diffusion_ip.hxx:129
LocalSolver::DiffusionParametersDefault::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: diffusion_ip.hxx:41
LocalSolver::DiffusionParametersDefault::neumann_nodes
static constexpr std::array< unsigned int, 0U > neumann_nodes
Array containing hypernode types corresponding to Neumann boundary.
Definition: advection_parab_ldgh.hxx:29
SmallMat::data
std::array< mat_entry_t, size()> & data()
Return data array that allows to manipulate the SmallMat.
Definition: dense_la.hxx:190
LocalSolver::DiffusionIP::node_element
Define type of node elements, especially with respect to nodal shape functions.
Definition: diffusion_ip.hxx:120
LocalSolver::DiffusionIP::node_system_dimension
static constexpr unsigned int node_system_dimension()
Dimension of of the solution evaluated with respect to a hypernode.
Definition: diffusion_ip.hxx:191
hypercube.hxx
LocalSolver::DiffusionIP::primal_at_boundary
SmallMat< 2 *hyEdge_dimT, n_shape_bdr_, lSol_float_t > primal_at_boundary(const SmallVec< n_loc_dofs_, lSol_float_t > &coeffs, hyEdgeT &hyper_edge) const
Evaluate primal variable at boundary.
Wrapper::LAPACKexception
Exception to be thrown if LAPACK's solve fails.
Definition: lapack.hxx:25
LocalSolver::DiffusionIP::system_dim
static constexpr unsigned int system_dim
Dimension of of the solution evaluated with respect to a hypernode.
Definition: diffusion_ip.hxx:215
LocalSolver::DiffusionParametersDefault::dirichlet_nodes
static constexpr std::array< unsigned int, 0U > dirichlet_nodes
Array containing hypernode types corresponding to Dirichlet boundary.
Definition: advection_parab_ldgh.hxx:25
LocalSolver::DiffusionIP::system_dimension
static constexpr unsigned int system_dimension()
Dimension of of the solution evaluated with respect to a hyperedge.
Definition: diffusion_ip.hxx:187
Hypercube::pow
static constexpr unsigned int pow(unsigned int base)
Return n to the power dim.
Definition: hypercube.hxx:25
LocalSolver::DiffusionIP::is_dirichlet
static constexpr bool is_dirichlet(const unsigned int node_type)
Find out whether a node is of Dirichlet type.
Definition: diffusion_ip.hxx:226
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
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::DiffusionIP::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: diffusion_ip.hxx:497
LocalSolver::DiffusionIP::node_system_dim
static constexpr unsigned int node_system_dim
Dimension of of the solution evaluated with respect to a hypernode.
Definition: diffusion_ip.hxx:221
LocalSolver::DiffusionIP::bulk_values
std::array< std::array< lSol_float_t, Hypercube< hyEdge_dimT >::pow(abscissas_sizeT)>, system_dim > bulk_values(const std::array< abscissa_float_t, abscissas_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::DiffusionIP::error_def::initial_error
static error_t initial_error()
Define how initial error is generated.
Definition: diffusion_ip.hxx:138
LocalSolver::DiffusionIP::DiffusionIP
DiffusionIP(const constructor_value_type &tau=1.)
Constructor for local solver.
Definition: diffusion_ip.hxx:414
LocalSolver::DiffusionIP::hyEdge_dim
static constexpr unsigned int hyEdge_dim()
Dimension of hyper edge type that this object solves on.
Definition: diffusion_ip.hxx:171
LocalSolver::DiffusionIP::tau_
const lSol_float_t tau_
(Globally constant) penalty parameter for HDG scheme.
Definition: diffusion_ip.hxx:239
LocalSolver::DiffusionIP::node_element::functions
std::tuple< TPP::ShapeFunction< TPP::ShapeType::Tensorial< TPP::ShapeType::Legendre< poly_deg >, hyEdge_dimT - 1 > > > functions
Definition: diffusion_ip.hxx:124
SmallMat
This class implements a small/dense matrix.
Definition: dense_la.hxx:34
Hypercube
Helper class containing numbers and functions related to hypercubes.
Definition: hypercube.hxx:12
Generated by   doxygen 1.8.17