Reference documentation for deal.II version GIT relicensing-1236-g9ea688d5df 2024-07-23 01:40:02+00:00
Searching...
No Matches
polynomials_abf.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
4// Copyright (C) 2005 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
18
19#include <iomanip>
20#include <iostream>
21#include <memory>
22
23
25
26
27
28namespace
29{
30 template <int dim>
31 std::vector<std::vector<Polynomials::Polynomial<double>>>
32 get_abf_polynomials(const unsigned int k)
33 {
34 std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
36
37 if (k == 0)
38 for (unsigned int d = 1; d < dim; ++d)
40 else
41 for (unsigned int d = 1; d < dim; ++d)
43
44 return pols;
45 }
46} // namespace
47
48template <int dim>
50 : TensorPolynomialsBase<dim>(k, n_polynomials(k))
51 , polynomial_space(get_abf_polynomials<dim>(k))
52{
53 // check that the dimensions match. we only store one of the 'dim'
54 // anisotropic polynomials that make up the vector-valued space, so
55 // multiply by 'dim'
57}
58
59
60
61template <int dim>
62void
64 const Point<dim> &unit_point,
65 std::vector<Tensor<1, dim>> &values,
68 std::vector<Tensor<4, dim>> &third_derivatives,
69 std::vector<Tensor<5, dim>> &fourth_derivatives) const
70{
71 Assert(values.size() == this->n() || values.empty(),
72 ExcDimensionMismatch(values.size(), this->n()));
77 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
78 ExcDimensionMismatch(third_derivatives.size(), this->n()));
79 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
80 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
81
82 const unsigned int n_sub = polynomial_space.n();
84 // arrays in the following block
85 // using a mutex to make sure they
86 // are not used by multiple threads
87 // at once
88 std::lock_guard<std::mutex> lock(mutex);
89
90 p_values.resize((values.empty()) ? 0 : n_sub);
93 p_third_derivatives.resize((third_derivatives.empty()) ? 0 : n_sub);
94 p_fourth_derivatives.resize((fourth_derivatives.empty()) ? 0 : n_sub);
95
96 for (unsigned int d = 0; d < dim; ++d)
97 {
98 // First we copy the point. The
99 // polynomial space for
100 // component d consists of
101 // polynomials of degree k+1 in
102 // x_d and degree k in the
103 // other variables. in order to
104 // simplify this, we use the
105 // same AnisotropicPolynomial
106 // space and simply rotate the
107 // coordinates through all
108 // directions.
109 Point<dim> p;
110 for (unsigned int c = 0; c < dim; ++c)
111 p[c] = unit_point[(c + d) % dim];
112
113 polynomial_space.evaluate(p,
114 p_values,
117 p_third_derivatives,
118 p_fourth_derivatives);
119
120 for (unsigned int i = 0; i < p_values.size(); ++i)
121 values[i + d * n_sub][d] = p_values[i];
122
123 for (unsigned int i = 0; i < p_grads.size(); ++i)
124 for (unsigned int d1 = 0; d1 < dim; ++d1)
125 grads[i + d * n_sub][d][(d1 + d) % dim] = p_grads[i][d1];
126
128 for (unsigned int d1 = 0; d1 < dim; ++d1)
129 for (unsigned int d2 = 0; d2 < dim; ++d2)
130 grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
132
133 for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
134 for (unsigned int d1 = 0; d1 < dim; ++d1)
135 for (unsigned int d2 = 0; d2 < dim; ++d2)
136 for (unsigned int d3 = 0; d3 < dim; ++d3)
137 third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
138 [(d2 + d) % dim][(d3 + d) % dim] =
139 p_third_derivatives[i][d1][d2][d3];
140
141 for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
142 for (unsigned int d1 = 0; d1 < dim; ++d1)
143 for (unsigned int d2 = 0; d2 < dim; ++d2)
144 for (unsigned int d3 = 0; d3 < dim; ++d3)
145 for (unsigned int d4 = 0; d4 < dim; ++d4)
146 fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
147 [(d2 + d) % dim][(d3 + d) % dim]
148 [(d4 + d) % dim] =
149 p_fourth_derivatives[i][d1][d2][d3][d4];
150 }
151}
152
153
154template <int dim>
155unsigned int
157{
158 switch (dim)
159 {
160 case 1:
161 // in 1d, we simply have Q_{k+2}, which has dimension k+3
162 return k + 3;
163
164 case 2:
165 // the polynomial space is Q_{k+2,k} \times Q_{k,k+2}, which has
166 // 2(k+3)(k+1) DoFs
167 return 2 * (k + 3) * (k + 1);
168
169 case 3:
170 // the polynomial space is Q_{k+2,k,k} \times Q_{k,k+2,k} \times
171 // Q_{k,k,k+2}, which has 3(k+3)(k+1)(k+1) DoFs
172 return 3 * (k + 3) * (k + 1) * (k + 1);
173
174 default:
176 }
177
178 return 0;
179}
180
181
182template <int dim>
183std::unique_ptr<TensorPolynomialsBase<dim>>
185{
186 return std::make_unique<PolynomialsABF<dim>>(*this);
187}
188
189
190template class PolynomialsABF<1>;
191template class PolynomialsABF<2>;
192template class PolynomialsABF<3>;
193
194
Definition point.h:111
const AnisotropicPolynomials< dim > polynomial_space
PolynomialsABF(const unsigned int k)
static unsigned int n_polynomials(const unsigned int degree)
void evaluate(const Point< dim > &unit_point, std::vector< Tensor< 1, dim > > &values, std::vector< Tensor< 2, dim > > &grads, std::vector< Tensor< 3, dim > > &grad_grads, std::vector< Tensor< 4, dim > > &third_derivatives, std::vector< Tensor< 5, dim > > &fourth_derivatives) const override
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NOT_IMPLEMENTED()
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)