Reference documentation for deal.II version 9.5.0
Searching...
No Matches
polynomials_abf.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2020 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
19
20#include <iomanip>
21#include <iostream>
22#include <memory>
23
24
26
27
28
29namespace
30{
31 template <int dim>
32 std::vector<std::vector<Polynomials::Polynomial<double>>>
33 get_abf_polynomials(const unsigned int k)
34 {
35 std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
37
38 if (k == 0)
39 for (unsigned int d = 1; d < dim; ++d)
41 else
42 for (unsigned int d = 1; d < dim; ++d)
44
45 return pols;
46 }
47} // namespace
48
49template <int dim>
51 : TensorPolynomialsBase<dim>(k, n_polynomials(k))
52 , polynomial_space(get_abf_polynomials<dim>(k))
53{
54 // check that the dimensions match. we only store one of the 'dim'
55 // anisotropic polynomials that make up the vector-valued space, so
56 // multiply by 'dim'
58}
59
60
61
62template <int dim>
63void
65 const Point<dim> & unit_point,
66 std::vector<Tensor<1, dim>> &values,
69 std::vector<Tensor<4, dim>> &third_derivatives,
70 std::vector<Tensor<5, dim>> &fourth_derivatives) const
71{
72 Assert(values.size() == this->n() || values.size() == 0,
73 ExcDimensionMismatch(values.size(), this->n()));
78 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
79 ExcDimensionMismatch(third_derivatives.size(), this->n()));
80 Assert(fourth_derivatives.size() == this->n() ||
81 fourth_derivatives.size() == 0,
82 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
83
84 const unsigned int n_sub = polynomial_space.n();
86 // arrays in the following block
87 // using a mutex to make sure they
88 // are not used by multiple threads
89 // at once
90 std::lock_guard<std::mutex> lock(mutex);
91
92 p_values.resize((values.size() == 0) ? 0 : n_sub);
95 p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
96 p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
97
98 for (unsigned int d = 0; d < dim; ++d)
99 {
100 // First we copy the point. The
101 // polynomial space for
102 // component d consists of
103 // polynomials of degree k+1 in
104 // x_d and degree k in the
105 // other variables. in order to
106 // simplify this, we use the
107 // same AnisotropicPolynomial
108 // space and simply rotate the
109 // coordinates through all
110 // directions.
111 Point<dim> p;
112 for (unsigned int c = 0; c < dim; ++c)
113 p(c) = unit_point((c + d) % dim);
114
115 polynomial_space.evaluate(p,
116 p_values,
119 p_third_derivatives,
120 p_fourth_derivatives);
121
122 for (unsigned int i = 0; i < p_values.size(); ++i)
123 values[i + d * n_sub][d] = p_values[i];
124
125 for (unsigned int i = 0; i < p_grads.size(); ++i)
126 for (unsigned int d1 = 0; d1 < dim; ++d1)
127 grads[i + d * n_sub][d][(d1 + d) % dim] = p_grads[i][d1];
128
130 for (unsigned int d1 = 0; d1 < dim; ++d1)
131 for (unsigned int d2 = 0; d2 < dim; ++d2)
132 grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
134
135 for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
136 for (unsigned int d1 = 0; d1 < dim; ++d1)
137 for (unsigned int d2 = 0; d2 < dim; ++d2)
138 for (unsigned int d3 = 0; d3 < dim; ++d3)
139 third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
140 [(d2 + d) % dim][(d3 + d) % dim] =
141 p_third_derivatives[i][d1][d2][d3];
142
143 for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
144 for (unsigned int d1 = 0; d1 < dim; ++d1)
145 for (unsigned int d2 = 0; d2 < dim; ++d2)
146 for (unsigned int d3 = 0; d3 < dim; ++d3)
147 for (unsigned int d4 = 0; d4 < dim; ++d4)
148 fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
149 [(d2 + d) % dim][(d3 + d) % dim]
150 [(d4 + d) % dim] =
151 p_fourth_derivatives[i][d1][d2][d3][d4];
152 }
153}
154
155
156template <int dim>
157unsigned int
159{
160 switch (dim)
161 {
162 case 1:
163 // in 1d, we simply have Q_{k+2}, which has dimension k+3
164 return k + 3;
165
166 case 2:
167 // the polynomial space is Q_{k+2,k} \times Q_{k,k+2}, which has
168 // 2(k+3)(k+1) DoFs
169 return 2 * (k + 3) * (k + 1);
170
171 case 3:
172 // the polynomial space is Q_{k+2,k,k} \times Q_{k,k+2,k} \times
173 // Q_{k,k,k+2}, which has 3(k+3)(k+1)(k+1) DoFs
174 return 3 * (k + 3) * (k + 1) * (k + 1);
175
176 default:
177 Assert(false, ExcNotImplemented());
178 }
179
180 return 0;
181}
182
183
184template <int dim>
185std::unique_ptr<TensorPolynomialsBase<dim>>
187{
188 return std::make_unique<PolynomialsABF<dim>>(*this);
189}
190
191
192template class PolynomialsABF<1>;
193template class PolynomialsABF<2>;
194template class PolynomialsABF<3>;
195
196
Definition point.h:112
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)