HyperHDG
one_dimensional.hxx
Go to the documentation of this file.
1 #pragma once // Ensure that file is included only once in a single compilation.
2 
4 
5 #include <array>
6 
7 namespace TPP
8 {
9 namespace Quadrature
10 {
11 /*!*************************************************************************************************
12  * \brief Gauss--Legendre quadrature points and weights in one spatial dimension.
13  *
14  * \tparam quad_deg Desired degree of accuracy.
15  *
16  * \authors Andreas Rupp, Heidelberg University, 2021.
17  **************************************************************************************************/
18 template <unsigned int quad_deg>
20 {
21  public:
22  /*!***********************************************************************************************
23  * \brief Calculate the amount of quadrature points at compile time.
24  *
25  * Naive implementation to calculate the amount of needed quadrature points if a rule of accuracy
26  * \c quad_deg is desired in one dimension.
27  *
28  * \retval n_quad_points Amount of needed quadrature points.
29  ************************************************************************************************/
30  static constexpr unsigned int n_points()
31  {
32  unsigned int amount = 1;
33  for (; 2 * amount - 1 < quad_deg; ++amount)
34  ;
35  return amount;
36  }
37 
38  private:
39  /*!***********************************************************************************************
40  * \brief Transform array of quadrature points in interval \f$[-1,1]\f$ to \f$[0,1]\f$.
41  ************************************************************************************************/
42  template <typename return_t>
43  static constexpr std::array<return_t, n_points()> transform_points(
44  std::array<return_t, n_points()> points)
45  {
46  for (unsigned int index = 0; index < n_points(); ++index)
47  points[index] = 0.5 * (points[index] + 1.);
48  return points;
49  }
50  /*!***********************************************************************************************
51  * \brief Transform array of quadrature weights for interval \f$[-1,1]\f$ to \f$[0,1]\f$.
52  ************************************************************************************************/
53  template <typename return_t>
54  static constexpr std::array<return_t, n_points()> transform_weights(
55  std::array<return_t, n_points()> weights)
56  {
57  for (unsigned int index = 0; index < n_points(); ++index)
58  weights[index] *= 0.5;
59  return weights;
60  }
61 
62  public:
63  /*!***********************************************************************************************
64  * \brief Gauss--Legendre quadrature points on one-dimensional unit interval.
65  *
66  * Returns the quadrature points of the quadrature rule with accuracy order \c quad_deg on the
67  * one-dimensional unit interval \f$[0,1]\f$.
68  *
69  * \tparam return_t Floating type specification. Default is double.
70  * \retval quad_points \c std::array containing the quadrature points.
71  ************************************************************************************************/
72  template <typename return_t = double>
73  static constexpr std::array<return_t, n_points()> points()
74  {
75  static_assert(n_points() < 10, "Amount of points needs to be smaller than 10!");
76 
77  if constexpr (n_points() == 1)
78  return transform_points(std::array<return_t, n_points()>({static_cast<return_t>(0.)}));
79  if constexpr (n_points() == 2)
80  return transform_points(std::array<return_t, n_points()>(
81  {-heron_root<return_t>(1. / 3.), heron_root<return_t>(1. / 3.)}));
82  if constexpr (n_points() == 3)
83  return transform_points(
84  std::array<return_t, n_points()>({-heron_root<return_t>(3. / 5.), static_cast<return_t>(0.),
85  heron_root<return_t>(3. / 5.)}));
86  if constexpr (n_points() == 4)
87  return transform_points(std::array<return_t, n_points()>(
88  {-heron_root<return_t>(3. / 7. + 2. / 7. * heron_root<return_t>(6. / 5.)),
89  -heron_root<return_t>(3. / 7. - 2. / 7. * heron_root<return_t>(6. / 5.)),
90  heron_root<return_t>(3. / 7. - 2. / 7. * heron_root<return_t>(6. / 5.)),
91  heron_root<return_t>(3. / 7. + 2. / 7. * heron_root<return_t>(6. / 5.))}));
92  if constexpr (n_points() == 5)
93  return transform_points(std::array<return_t, n_points()>(
94  {static_cast<return_t>(-heron_root<return_t>(5. + 2. * heron_root<return_t>(10. / 7.)) /
95  3.),
96  static_cast<return_t>(-heron_root<return_t>(5. - 2. * heron_root<return_t>(10. / 7.)) /
97  3.),
98  static_cast<return_t>(0.),
99  static_cast<return_t>(heron_root<return_t>(5. - 2. * heron_root<return_t>(10. / 7.)) / 3.),
100  static_cast<return_t>(heron_root<return_t>(5. + 2. * heron_root<return_t>(10. / 7.)) /
101  3.)}));
102  if constexpr (n_points() == 6)
103  return transform_points(std::array<return_t, n_points()>(
104  {static_cast<return_t>(0.6612093864662645), static_cast<return_t>(-0.6612093864662645),
105  static_cast<return_t>(-0.2386191860831969), static_cast<return_t>(0.2386191860831969),
106  static_cast<return_t>(-0.9324695142031521), static_cast<return_t>(0.9324695142031521)}));
107  if constexpr (n_points() == 7)
108  return transform_points(std::array<return_t, n_points()>(
109  {static_cast<return_t>(0.0000000000000000), static_cast<return_t>(0.4058451513773972),
110  static_cast<return_t>(-0.4058451513773972), static_cast<return_t>(-0.7415311855993945),
111  static_cast<return_t>(0.7415311855993945), static_cast<return_t>(-0.9491079123427585),
112  static_cast<return_t>(0.9491079123427585)}));
113  if constexpr (n_points() == 8)
114  return transform_points(std::array<return_t, n_points()>(
115  {static_cast<return_t>(-0.1834346424956498), static_cast<return_t>(0.1834346424956498),
116  static_cast<return_t>(-0.5255324099163290), static_cast<return_t>(0.5255324099163290),
117  static_cast<return_t>(-0.7966664774136267), static_cast<return_t>(0.7966664774136267),
118  static_cast<return_t>(-0.9602898564975363), static_cast<return_t>(0.9602898564975363)}));
119  if constexpr (n_points() == 9)
120  return transform_points(std::array<return_t, n_points()>(
121  {static_cast<return_t>(0.0000000000000000), static_cast<return_t>(-0.8360311073266358),
122  static_cast<return_t>(0.8360311073266358), static_cast<return_t>(-0.9681602395076261),
123  static_cast<return_t>(0.9681602395076261), static_cast<return_t>(-0.3242534234038089),
124  static_cast<return_t>(0.3123470770400029), static_cast<return_t>(0.2606106964029354),
125  static_cast<return_t>(0.2606106964029354)}));
126  }
127  /*!***********************************************************************************************
128  * \brief Gauss--Legendre quadrature weights on one-dimensional unit interval.
129  *
130  * Returns the quadrature weights of the quadrature rule with accuracy order \c quad_deg on the
131  * one-dimensional unit interval \f$[0,1]\f$.
132  *
133  * \tparam return_t Floating type specification. Default is double.
134  * \retval quad_weights \c std::array containing the quadrature weights.
135  ************************************************************************************************/
136  template <typename return_t = double>
137  static constexpr std::array<return_t, n_points()> weights()
138  {
139  static_assert(n_points() < 10, "Amount of points needs to be smaller than 10!");
140 
141  if constexpr (n_points() == 1)
142  return transform_weights(std::array<return_t, n_points()>({static_cast<return_t>(2.)}));
143  if constexpr (n_points() == 2)
144  return transform_weights(
145  std::array<return_t, n_points()>({static_cast<return_t>(1.), static_cast<return_t>(1.)}));
146  if constexpr (n_points() == 3)
147  return transform_weights(std::array<return_t, n_points()>({static_cast<return_t>(5. / 9.),
148  static_cast<return_t>(8. / 9.),
149  static_cast<return_t>(5. / 9.)}));
150  if constexpr (n_points() == 4)
151  return transform_weights(std::array<return_t, n_points()>(
152  {static_cast<return_t>(1. / 36. * (18. - heron_root<return_t>(30.))),
153  static_cast<return_t>(1. / 36. * (18. + heron_root<return_t>(30.))),
154  static_cast<return_t>(1. / 36. * (18. + heron_root<return_t>(30.))),
155  static_cast<return_t>(1. / 36. * (18. - heron_root<return_t>(30.)))}));
156  if constexpr (n_points() == 5)
157  return transform_weights(std::array<return_t, n_points()>(
158  {static_cast<return_t>(1. / 900. * (322. - 13. * heron_root<return_t>(70.))),
159  static_cast<return_t>(1. / 900. * (322. + 13. * heron_root<return_t>(70.))),
160  static_cast<return_t>(1. / 900. * (322. + 190.)),
161  static_cast<return_t>(1. / 900. * (322. + 13. * heron_root<return_t>(70.))),
162  static_cast<return_t>(1. / 900. * (322. - 13. * heron_root<return_t>(70.)))}));
163  if constexpr (n_points() == 6)
164  return transform_weights(std::array<return_t, n_points()>(
165  {static_cast<return_t>(0.3607615730481386), static_cast<return_t>(0.3607615730481386),
166  static_cast<return_t>(0.4679139345726910), static_cast<return_t>(0.4679139345726910),
167  static_cast<return_t>(0.1713244923791704), static_cast<return_t>(0.1713244923791700)}));
168  if constexpr (n_points() == 7)
169  return transform_weights(std::array<return_t, n_points()>(
170  {static_cast<return_t>(0.4179591836734694), static_cast<return_t>(0.3818300505051189),
171  static_cast<return_t>(0.3818300505051189), static_cast<return_t>(0.2797053914892766),
172  static_cast<return_t>(0.2797053914892766), static_cast<return_t>(0.1294849661688697),
173  static_cast<return_t>(0.1294849661688697)}));
174  if constexpr (n_points() == 8)
175  return transform_weights(std::array<return_t, n_points()>(
176  {static_cast<return_t>(0.3626837833783620), static_cast<return_t>(0.3626837833783620),
177  static_cast<return_t>(0.3137066458778873), static_cast<return_t>(0.3137066458778873),
178  static_cast<return_t>(0.2223810344533745), static_cast<return_t>(0.2223810344533745),
179  static_cast<return_t>(0.1012285362903763), static_cast<return_t>(0.1012285362903763)}));
180  if constexpr (n_points() == 9)
181  return transform_weights(std::array<return_t, n_points()>(
182  {static_cast<return_t>(0.3302393550012598), static_cast<return_t>(0.1806481606948574),
183  static_cast<return_t>(0.1806481606948574), static_cast<return_t>(0.0812743883615744),
184  static_cast<return_t>(0.0812743883615744), static_cast<return_t>(0.3123470770400029),
185  static_cast<return_t>(0.3123470770400029), static_cast<return_t>(0.2606106964029354),
186  static_cast<return_t>(0.2606106964029354)}));
187  }
188 }; // end of struct GaussLegendre
189 
190 } // end of namespace Quadrature
191 
192 } // end of namespace TPP
dim3
constexpr std::array< unsigned short, 3 > dim3
Definition: lexicographic_03.cc:13
TPCC::Element
Tensor coordinates for a facet of dimension k in the complex of dimension n.
Definition: element.h:38
dim4
constexpr std::array< unsigned short, 4 > dim4
Definition: lexicographic_02.cc:13
test
void test(const TPCC::Slab< n, k, B, S, T > &slab)
Definition: slab_01.cc:16
test
void test(const TPCC::Slab< n, k, B, S, T > &slab)
Definition: slab_02.cc:13
dim4
constexpr std::array< unsigned short, 4 > dim4
Definition: lexicographic_03.cc:14
test
void test(const TPCC::Slab< n, k, B, S, T > &slab)
Definition: slab_03.cc:13
main
int main()
Definition: element_01.cc:49
TPCC::Combination::eliminate
constexpr std::enable_if<(kk > 0), Combination< n, k - 1, T > >::type eliminate(unsigned int i) const
The combination obtained by eliminating the ith element.
Definition: combinations.h:98
check_push_test.index
index
Definition: check_push_test.py:10
dim2
constexpr std::array< unsigned short, 2 > dim2
Definition: slab_03.cc:9
TPCC::Combination::print_debug
void print_debug(std::ostream &os) const
Print the content of this object for debugging.
Definition: combinations.h:176
combinations.h
dim2
constexpr std::array< unsigned short, 2 > dim2
Definition: slab_02.cc:9
block_sizes21
const unsigned int block_sizes21[]
Definition: lexicographic_01.cc:14
dim2
constexpr std::array< unsigned short, 2 > dim2
Definition: slab_01.cc:12
output_2
const char * output_2[]
Definition: combination_02.cc:19
output_2
const char * output_2[]
Definition: combination_01.cc:15
lexicographic.h
pascal
void pascal(unsigned int padding=n)
Definition: binomial.cc:92
main
int main()
Definition: combination_01.cc:58
main
int main()
Definition: combination_02.cc:61
test
void test(const TPCC::Element< n, k, T1, T2 > &e)
Definition: element_01.cc:7
block_sizes30
const unsigned int block_sizes30[]
Definition: lexicographic_01.cc:17
combination_data_1
const char * combination_data_1[]
Definition: binomial.cc:23
block_sizes31
const unsigned int block_sizes31[]
Definition: lexicographic_01.cc:18
dim2
constexpr std::array< unsigned short, 2 > dim2
Definition: lexicographic_02.cc:11
dim2
constexpr std::array< unsigned short, 2 > dim2
Definition: lexicographic_01.cc:23
TPP::Quadrature::GaussLegendre::transform_weights
static constexpr std::array< return_t, n_points()> transform_weights(std::array< return_t, n_points()> weights)
Transform array of quadrature weights for interval to .
Definition: one_dimensional.hxx:54
dim2
constexpr std::array< unsigned short, 2 > dim2
Definition: lexicographic_03.cc:12
TPCC::Combinations::index
static constexpr unsigned int index(const Combination< n, k, T > &combi)
The index of a combination within the lexicographic enumeration.
Definition: combinations.h:272
combination_data
const char ** combination_data[]
Definition: binomial.cc:35
TPCC::Lexicographic::block_size
constexpr Bint block_size(Tint block) const
The number of elements in one direction.
Definition: lexicographic.h:128
main
int main()
Definition: binomial.cc:118
output_data
const char ** output_data[]
Definition: combination_02.cc:33
output_data
const char ** output_data[]
Definition: combination_01.cc:32
TPCC::Combination::add
constexpr std::enable_if<(kk< n), Combination< n, k+1, T > >::type add(unsigned int i) const
The combination obtained by adding the element i.
Definition: combinations.h:130
block_sizes32
const unsigned int block_sizes32[]
Definition: lexicographic_01.cc:19
pascal_data
unsigned int pascal_data[11][11]
Definition: binomial.cc:9
TPP::Quadrature::GaussLegendre::points
static constexpr std::array< return_t, n_points()> points()
Gauss–Legendre quadrature points on one-dimensional unit interval.
Definition: one_dimensional.hxx:73
TPCC::Lexicographic::size
constexpr Bint size() const
The number of elements in this set.
Definition: lexicographic.h:117
size2
const unsigned int size2[]
Definition: lexicographic_01.cc:11
TPP::Quadrature::GaussLegendre
Gauss–Legendre quadrature points and weights in one spatial dimension.
Definition: one_dimensional.hxx:19
print_combination::doit
static void doit()
Definition: binomial.cc:57
main
int main()
Definition: lexicographic_03.cc:68
combination_data_3
const char * combination_data_3[]
Definition: binomial.cc:28
main
int main()
Definition: lexicographic_02.cc:61
output_4
const char * output_4[]
Definition: combination_01.cc:26
main
int main()
Definition: lexicographic_01.cc:58
block_sizes2
const unsigned int * block_sizes2[]
Definition: lexicographic_01.cc:16
output_4
const char * output_4[]
Definition: combination_02.cc:29
block_sizes20
const unsigned int block_sizes20[]
Definition: lexicographic_01.cc:13
TPP
Contains the functionalities to evaluate and integrate tensor product type polynomials.
Definition: compile_time_tricks.hxx:7
operator<<
std::ostream & operator<<(std::ostream &os, std::array< T, N > a)
Definition: binomial.cc:39
size3
const unsigned int size3[]
Definition: lexicographic_01.cc:12
TPP::Quadrature::GaussLegendre::transform_points
static constexpr std::array< return_t, n_points()> transform_points(std::array< return_t, n_points()> points)
Transform array of quadrature points in interval to .
Definition: one_dimensional.hxx:43
print_mesh
void print_mesh(const MESH &mesh)
Definition: lexicographic_03.cc:17
TPCC::Combinations::dual
static std::array< unsigned int, n - k > dual(unsigned int index)
The array of numbers (of length n-k) of numbers not in the combination of the given index.
Definition: combinations.h:296
dim1
constexpr std::array< unsigned short, 1 > dim1
Definition: lexicographic_02.cc:10
output_0
const char * output_0[]
Definition: combination_02.cc:13
dim1
constexpr std::array< unsigned short, 1 > dim1
Definition: lexicographic_03.cc:11
print_mesh
void print_mesh(const MESH &mesh)
Definition: lexicographic_02.cc:16
element.h
TPCC::Lexicographic
Lexicographic enumeration of the k-dimensional faces in a tensor product chain complex of dimension n...
Definition: lexicographic.h:44
output_5
const char * output_5[]
Definition: combination_01.cc:30
TPCC::Combinations::size
static constexpr unsigned int size()
The number of such combinations.
Definition: combinations.h:240
test_3
void test_3()
Definition: lexicographic_02.cc:46
test_3
void test_3()
Definition: lexicographic_01.cc:43
test_3
void test_3()
Definition: lexicographic_03.cc:53
TPP::Quadrature::GaussLegendre::n_points
static constexpr unsigned int n_points()
Calculate the amount of quadrature points at compile time.
Definition: one_dimensional.hxx:30
test
void test()
Definition: combination_01.cc:35
test
void test()
Definition: combination_02.cc:36
test_2
void test_2()
Definition: lexicographic_03.cc:45
combination_data_5
const char * combination_data_5[]
Definition: binomial.cc:33
test_2
void test_2()
Definition: lexicographic_02.cc:38
test_2
void test_2()
Definition: lexicographic_01.cc:27
combination_data_4
const char * combination_data_4[]
Definition: binomial.cc:31
dim1
constexpr std::array< unsigned short, 1 > dim1
Definition: slab_03.cc:8
dim1
constexpr std::array< unsigned short, 1 > dim1
Definition: slab_02.cc:8
dim1
constexpr std::array< unsigned short, 1 > dim1
Definition: slab_01.cc:11
combination_data_2
const char * combination_data_2[]
Definition: binomial.cc:25
output_3
const char * output_3[]
Definition: combination_01.cc:20
TPCC::Combination< n, k >
TPCC::Combination::out
T out(unsigned int i) const
The ith element which is not part of the combination in descending order.
Definition: combinations.h:84
output_3
const char * output_3[]
Definition: combination_02.cc:25
slab.h
block_sizes3
const unsigned int * block_sizes3[]
Definition: lexicographic_01.cc:21
TPCC::Combinations::value
static std::array< unsigned int, k > value(unsigned int index)
The array of numbers (of length k) in the combination with given index.
Definition: combinations.h:284
combination_data_0
const char * combination_data_0[]
Definition: binomial.cc:21
test_4
void test_4()
Definition: lexicographic_02.cc:54
TPCC::Element::print_debug
void print_debug(std::ostream &os) const
Function for printing the data stored in the element.
Definition: element.h:112
test_4
void test_4()
Definition: lexicographic_03.cc:61
TPCC::Slab::block_size
constexpr Bint block_size(Tint block) const
The number of elements in one direction.
Definition: slab.h:73
test_1
void test_1()
Definition: lexicographic_02.cc:30
fibonacci
constexpr unsigned int fibonacci(unsigned int i)
Definition: element_01.cc:20
test_1
void test_1()
Definition: lexicographic_03.cc:37
TPCC::Combinations
Definition: combinations.h:195
block_sizes22
const unsigned int block_sizes22[]
Definition: lexicographic_01.cc:15
diffusion_uniform.A
A
Definition: diffusion_uniform.py:39
TPCC::Element::facet
constexpr Element< n, k - 1, Sint, Tint > facet(Tint index) const
Enumeration of the boundary facets of the element.
Definition: element.h:136
TPP::Quadrature::GaussLegendre::weights
static constexpr std::array< return_t, n_points()> weights()
Gauss–Legendre quadrature weights on one-dimensional unit interval.
Definition: one_dimensional.hxx:137
output_1
const char * output_1[]
Definition: combination_02.cc:15
TPCC::binomial
constexpr T binomial(T n, T k)
Compute the binomial coefficient n over k.
Definition: combinations.h:17
output_1
const char * output_1[]
Definition: combination_01.cc:13
print_combination
Definition: binomial.cc:55
compile_time_tricks.hxx
dim3
constexpr std::array< unsigned short, 3 > dim3
Definition: slab_03.cc:10
dim3
constexpr std::array< unsigned short, 3 > dim3
Definition: slab_02.cc:10
dim3
constexpr std::array< unsigned short, 3 > dim3
Definition: slab_01.cc:13
main
int main(int, char *[])
Definition: slab_01.cc:69
main
int main(int, char *[])
Definition: slab_02.cc:63
TPCC::Slab::size
constexpr Bint size() const
Definition: slab.h:68
dim3
constexpr std::array< unsigned short, 3 > dim3
Definition: lexicographic_01.cc:24
block_sizes33
const unsigned int block_sizes33[]
Definition: lexicographic_01.cc:20
main
int main(int, char *[])
Definition: slab_03.cc:64
TPCC::Slab
A slab of thickness one cell cut out of a tensor product chain complex.
Definition: element.h:9
dim3
constexpr std::array< unsigned short, 3 > dim3
Definition: lexicographic_02.cc:12