Reference documentation for deal.II version 9.0.0
polynomials_rt_bubbles.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 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
9 // Public License as published by the Free Software Foundation; either
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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include "deal.II/base/polynomials_rt_bubbles.h"
18 #include <deal.II/base/quadrature_lib.h>
19 #include <deal.II/base/thread_management.h>
20 #include <iostream>
21 #include <iomanip>
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 
26 template <int dim>
27 PolynomialsRT_Bubbles<dim>::PolynomialsRT_Bubbles (const unsigned int k)
28  :
29  my_degree(k),
30  raviart_thomas_space (k-1),
31  monomials(k+2),
32  n_pols(compute_n_pols(k))
33 {
34  Assert (dim >= 2, ExcImpossibleInDim(dim));
35 
36  for (unsigned int i=0; i<monomials.size(); ++i)
37  monomials[i] = Polynomials::Monomial<double> (i);
38 }
39 
40 
41 
42 template <int dim>
43 void
44 PolynomialsRT_Bubbles<dim>::compute (const Point<dim> &unit_point,
45  std::vector<Tensor<1,dim> > &values,
46  std::vector<Tensor<2,dim> > &grads,
47  std::vector<Tensor<3,dim> > &grad_grads,
48  std::vector<Tensor<4,dim> > &third_derivatives,
49  std::vector<Tensor<5,dim> > &fourth_derivatives) const
50 {
51  Assert(values.size()==n_pols || values.size()==0,
52  ExcDimensionMismatch(values.size(), n_pols));
53  Assert(grads.size()==n_pols || grads.size()==0,
54  ExcDimensionMismatch(grads.size(), n_pols));
55  Assert(grad_grads.size()==n_pols || grad_grads.size()==0,
56  ExcDimensionMismatch(grad_grads.size(), n_pols));
57  Assert(third_derivatives.size()==n_pols || third_derivatives.size()==0,
58  ExcDimensionMismatch(third_derivatives.size(), n_pols));
59  Assert(fourth_derivatives.size()==n_pols || fourth_derivatives.size()==0,
60  ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
61 
62  // Third and fourth derivatives are not implemented
63  (void)third_derivatives;
64  Assert(third_derivatives.size()==0,
65  ExcNotImplemented());
66  (void)fourth_derivatives;
67  Assert(fourth_derivatives.size()==0,
68  ExcNotImplemented());
69 
70  const unsigned int n_sub = raviart_thomas_space.n();
71 
72  // Guard access to the scratch arrays in the following block
73  // using a mutex to make sure they are not used by multiple threads
74  // at once
75  {
76  static Threads::Mutex mutex;
77  Threads::Mutex::ScopedLock lock(mutex);
78 
79  static std::vector<Tensor<1,dim> > p_values;
80  static std::vector<Tensor<2,dim> > p_grads;
81  static std::vector<Tensor<3,dim> > p_grad_grads;
82  static std::vector<Tensor<4,dim> > p_third_derivatives;
83  static std::vector<Tensor<5,dim> > p_fourth_derivatives;
84 
85  p_values.resize((values.size() == 0) ? 0 : n_sub);
86  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
87  p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
88 
89  // This is the Raviart-Thomas part of the space
90  raviart_thomas_space.compute (unit_point, p_values, p_grads, p_grad_grads, p_third_derivatives, p_fourth_derivatives);
91  for (unsigned int i=0; i<p_values.size(); ++i)
92  values[i] = p_values[i];
93  for (unsigned int i=0; i<p_grads.size(); ++i)
94  grads[i] = p_grads[i];
95  for (unsigned int i=0; i<p_grad_grads.size(); ++i)
96  grad_grads[i] = p_grad_grads[i];
97  }
98 
99  // Next we compute the polynomials and derivatives
100  // of the curl part of the space
101  const unsigned int n_derivatives = 3;
102  double monoval_plus[dim][n_derivatives+1];
103  double monoval[dim][n_derivatives+1];
104 
105  double monoval_i[dim][n_derivatives+1];
106  double monoval_j[dim][n_derivatives+1];
107  double monoval_jplus[dim][n_derivatives+1];
108 
109  unsigned int start = n_sub;
110 
111  if (dim == 2)
112  {
113  // In 2d the curl part of the space is spanned by the vectors
114  // of two types. The first one is
115  // [ x^i * [y^(k+1)]' ]
116  // [ -[x^i]' * y^(k+1) ]
117  // The second one can be obtained from the first by a cyclic
118  // rotation of the coordinates.
119  // monoval_i = x^i,
120  // monoval_plus = x^(k+1)
121  for (unsigned int d=0; d<dim; ++d)
122  monomials[my_degree+1].value(unit_point(d), n_derivatives, monoval_plus[d]);
123 
124  for (unsigned int i=0; i<=my_degree; ++i, ++start)
125  {
126  for (unsigned int d=0; d<dim; ++d)
127  monomials[i].value(unit_point(d), n_derivatives, monoval_i[d]);
128 
129  if (values.size() != 0)
130  {
131  values[start][0] = monoval_i[0][0] * monoval_plus[1][1];
132  values[start][1] = -monoval_i[0][1] * monoval_plus[1][0];
133 
134  values[start+my_degree+1][0] = -monoval_plus[0][0] * monoval_i[1][1];
135  values[start+my_degree+1][1] = monoval_plus[0][1] * monoval_i[1][0];
136  }
137 
138  if (grads.size() != 0)
139  {
140  grads[start][0][0] = monoval_i[0][1] * monoval_plus[1][1];
141  grads[start][0][1] = monoval_i[0][0] * monoval_plus[1][2];
142  grads[start][1][0] = -monoval_i[0][2] * monoval_plus[1][0];
143  grads[start][1][1] = -monoval_i[0][1] * monoval_plus[1][1];
144 
145  grads[start+my_degree+1][0][0] = -monoval_plus[0][1] * monoval_i[1][1];
146  grads[start+my_degree+1][0][1] = -monoval_plus[0][0] * monoval_i[1][2];
147  grads[start+my_degree+1][1][0] = monoval_plus[0][2] * monoval_i[1][0];
148  grads[start+my_degree+1][1][1] = monoval_plus[0][1] * monoval_i[1][1];
149  }
150 
151  if (grad_grads.size() != 0)
152  {
153  grad_grads[start][0][0][0] = monoval_i[0][2] * monoval_plus[1][1];
154  grad_grads[start][0][0][1] = monoval_i[0][1] * monoval_plus[1][2];
155  grad_grads[start][0][1][0] = monoval_i[0][1] * monoval_plus[1][2];
156  grad_grads[start][0][1][1] = monoval_i[0][0] * monoval_plus[1][3];
157  grad_grads[start][1][0][0] = -monoval_i[0][3] * monoval_plus[1][0];
158  grad_grads[start][1][0][1] = -monoval_i[0][2] * monoval_plus[1][1];
159  grad_grads[start][1][1][0] = -monoval_i[0][2] * monoval_plus[1][1];
160  grad_grads[start][1][1][1] = -monoval_i[0][1] * monoval_plus[1][2];
161 
162  grad_grads[start+my_degree+1][0][0][0] = -monoval_plus[0][2] * monoval_i[1][1];
163  grad_grads[start+my_degree+1][0][0][1] = -monoval_plus[0][1] * monoval_i[1][2];
164  grad_grads[start+my_degree+1][0][1][0] = -monoval_plus[0][1] * monoval_i[1][2];
165  grad_grads[start+my_degree+1][0][1][1] = -monoval_plus[0][0] * monoval_i[1][3];
166  grad_grads[start+my_degree+1][1][0][0] = monoval_plus[0][3] * monoval_i[1][0];
167  grad_grads[start+my_degree+1][1][0][1] = monoval_plus[0][2] * monoval_i[1][1];
168  grad_grads[start+my_degree+1][1][1][0] = monoval_plus[0][2] * monoval_i[1][1];
169  grad_grads[start+my_degree+1][1][1][1] = monoval_plus[0][1] * monoval_i[1][2];
170  }
171  }
172  Assert(start == n_pols - my_degree - 1, ExcInternalError());
173  }
174  else if (dim == 3)
175  {
176  // In 3d the first type of basis vector is
177  // [ x^i * y^j * z^k * (j+k+2) ]
178  // [ -[x^i]' * y^(j+1) * z^k ]
179  // [ -[x^i]' * y^j * z^(k+1) ],
180  // For the second type of basis vector y and z
181  // are swapped. Then for each of these,
182  // two more are obtained by the cyclic rotation
183  // of the coordinates.
184  // monoval = x^k, monoval_plus = x^(k+1)
185  // monoval_* = x^*, monoval_jplus = x^(j+1)
186  for (unsigned int d=0; d<dim; ++d)
187  {
188  monomials[my_degree+1].value(unit_point(d), n_derivatives, monoval_plus[d]);
189  monomials[my_degree].value(unit_point(d), n_derivatives, monoval[d]);
190  }
191 
192  const unsigned int n_curls = (my_degree+1) * (2*my_degree + 1);
193  // Span of @f$\tilde{B}@f$
194  for (unsigned int i=0; i<=my_degree; ++i)
195  {
196  for (unsigned int d=0; d<dim; ++d)
197  monomials[i].value(unit_point(d), n_derivatives, monoval_i[d]);
198 
199  for (unsigned int j = 0; j<=my_degree; ++j)
200  {
201  for (unsigned int d=0; d<dim; ++d)
202  {
203  monomials[j].value(unit_point(d), n_derivatives, monoval_j[d]);
204  monomials[j+1].value(unit_point(d), n_derivatives, monoval_jplus[d]);
205  }
206 
207  if (values.size() != 0)
208  {
209  values[start][0] = monoval_i[0][0] * monoval_j[1][0] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
210  values[start][1] = -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][0];
211  values[start][2] = -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][0];
212 
213  values[start+n_curls][0] = -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][0];
214  values[start+n_curls][1] = monoval_j[0][0] * monoval_i[1][0] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
215  values[start+n_curls][2] = -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][0];
216 
217  values[start+2*n_curls][0] = -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][1];
218  values[start+2*n_curls][1] = -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][1];
219  values[start+2*n_curls][2] = monoval_j[0][0] * monoval[1][0] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
220 
221  // Only unique triples of powers (i j k)
222  // and (i k j) are allowed, 0 <= i,j <= k
223  if (j != my_degree)
224  {
225  values[start+1][0] = monoval_i[0][0] * monoval[1][0] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
226  values[start+1][1] = -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][0];
227  values[start+1][2] = -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][0];
228 
229  values[start+n_curls+1][0] = -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][0];
230  values[start+n_curls+1][1] = monoval[0][0] * monoval_i[1][0] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
231  values[start+n_curls+1][2] = -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][0];
232 
233  values[start+2*n_curls+1][0] = -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][1];
234  values[start+2*n_curls+1][1] = -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][1];
235  values[start+2*n_curls+1][2] = monoval[0][0] * monoval_j[1][0] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
236  }
237  }
238 
239  if (grads.size() != 0)
240  {
241  grads[start][0][0] = monoval_i[0][1] * monoval_j[1][0] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
242  grads[start][0][1] = monoval_i[0][0] * monoval_j[1][1] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
243  grads[start][0][2] = monoval_i[0][0] * monoval_j[1][0] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
244  grads[start][1][0] = -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][0];
245  grads[start][1][1] = -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][0];
246  grads[start][1][2] = -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][1];
247  grads[start][2][0] = -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][0];
248  grads[start][2][1] = -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][0];
249  grads[start][2][2] = -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][1];
250 
251  grads[start+n_curls][0][0] = -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][0];
252  grads[start+n_curls][0][1] = -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][0];
253  grads[start+n_curls][0][2] = -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][1];
254  grads[start+n_curls][1][0] = monoval_j[0][1] * monoval_i[1][0] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
255  grads[start+n_curls][1][1] = monoval_j[0][0] * monoval_i[1][1] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
256  grads[start+n_curls][1][2] = monoval_j[0][0] * monoval_i[1][0] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
257  grads[start+n_curls][2][0] = -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][0];
258  grads[start+n_curls][2][1] = -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][0];
259  grads[start+n_curls][2][2] = -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][1];
260 
261  grads[start+2*n_curls][0][0] = -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][1];
262  grads[start+2*n_curls][0][1] = -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][1];
263  grads[start+2*n_curls][0][2] = -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][2];
264  grads[start+2*n_curls][1][0] = -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][1];
265  grads[start+2*n_curls][1][1] = -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][1];
266  grads[start+2*n_curls][1][2] = -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][2];
267  grads[start+2*n_curls][2][0] = monoval_j[0][1] * monoval[1][0] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
268  grads[start+2*n_curls][2][1] = monoval_j[0][0] * monoval[1][1] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
269  grads[start+2*n_curls][2][2] = monoval_j[0][0] * monoval[1][0] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
270 
271  if (j != my_degree)
272  {
273  grads[start+1][0][0] = monoval_i[0][1] * monoval[1][0] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
274  grads[start+1][0][1] = monoval_i[0][0] * monoval[1][1] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
275  grads[start+1][0][2] = monoval_i[0][0] * monoval[1][0] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
276  grads[start+1][1][0] = -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][0];
277  grads[start+1][1][1] = -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][0];
278  grads[start+1][1][2] = -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][1];
279  grads[start+1][2][0] = -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][0];
280  grads[start+1][2][1] = -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][0];
281  grads[start+1][2][2] = -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][1];
282 
283  grads[start+n_curls+1][0][0] = -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][0];
284  grads[start+n_curls+1][0][1] = -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][0];
285  grads[start+n_curls+1][0][2] = -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][1];
286  grads[start+n_curls+1][1][0] = monoval[0][1] * monoval_i[1][0] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
287  grads[start+n_curls+1][1][1] = monoval[0][0] * monoval_i[1][1] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
288  grads[start+n_curls+1][1][2] = monoval[0][0] * monoval_i[1][0] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
289  grads[start+n_curls+1][2][0] = -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][0];
290  grads[start+n_curls+1][2][1] = -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][0];
291  grads[start+n_curls+1][2][2] = -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][1];
292 
293  grads[start+2*n_curls+1][0][0] = -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][1];
294  grads[start+2*n_curls+1][0][1] = -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][1];
295  grads[start+2*n_curls+1][0][2] = -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][2];
296  grads[start+2*n_curls+1][1][0] = -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][1];
297  grads[start+2*n_curls+1][1][1] = -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][1];
298  grads[start+2*n_curls+1][1][2] = -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][2];
299  grads[start+2*n_curls+1][2][0] = monoval[0][1] * monoval_j[1][0] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
300  grads[start+2*n_curls+1][2][1] = monoval[0][0] * monoval_j[1][1] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
301  grads[start+2*n_curls+1][2][2] = monoval[0][0] * monoval_j[1][0] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
302  }
303  }
304 
305  if (grad_grads.size() != 0)
306  {
307  grad_grads[start][0][0][0] = monoval_i[0][2] * monoval_j[1][0] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
308  grad_grads[start][0][0][1] = monoval_i[0][1] * monoval_j[1][1] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
309  grad_grads[start][0][0][2] = monoval_i[0][1] * monoval_j[1][0] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
310  grad_grads[start][0][1][0] = monoval_i[0][1] * monoval_j[1][1] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
311  grad_grads[start][0][1][1] = monoval_i[0][0] * monoval_j[1][2] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
312  grad_grads[start][0][1][2] = monoval_i[0][0] * monoval_j[1][1] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
313  grad_grads[start][0][2][0] = monoval_i[0][1] * monoval_j[1][0] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
314  grad_grads[start][0][2][1] = monoval_i[0][0] * monoval_j[1][1] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
315  grad_grads[start][0][2][2] = monoval_i[0][0] * monoval_j[1][0] * monoval[2][2] * static_cast<double>(j + my_degree + 2);
316  grad_grads[start][1][0][0] = -monoval_i[0][3] * monoval_jplus[1][0] * monoval[2][0];
317  grad_grads[start][1][0][1] = -monoval_i[0][2] * monoval_jplus[1][1] * monoval[2][0];
318  grad_grads[start][1][0][2] = -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][1];
319  grad_grads[start][1][1][0] = -monoval_i[0][2] * monoval_jplus[1][1] * monoval[2][0];
320  grad_grads[start][1][1][1] = -monoval_i[0][1] * monoval_jplus[1][2] * monoval[2][0];
321  grad_grads[start][1][1][2] = -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][1];
322  grad_grads[start][1][2][0] = -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][1];
323  grad_grads[start][1][2][1] = -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][1];
324  grad_grads[start][1][2][2] = -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][2];
325  grad_grads[start][2][0][0] = -monoval_i[0][3] * monoval_j[1][0] * monoval_plus[2][0];
326  grad_grads[start][2][0][1] = -monoval_i[0][2] * monoval_j[1][1] * monoval_plus[2][0];
327  grad_grads[start][2][0][2] = -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][1];
328  grad_grads[start][2][1][0] = -monoval_i[0][2] * monoval_j[1][1] * monoval_plus[2][0];
329  grad_grads[start][2][1][1] = -monoval_i[0][1] * monoval_j[1][2] * monoval_plus[2][0];
330  grad_grads[start][2][1][2] = -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][1];
331  grad_grads[start][2][2][0] = -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][1];
332  grad_grads[start][2][2][1] = -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][1];
333  grad_grads[start][2][2][2] = -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][2];
334 
335  grad_grads[start+n_curls][0][0][0] = -monoval_jplus[0][2] * monoval_i[1][1] * monoval[2][0];
336  grad_grads[start+n_curls][0][0][1] = -monoval_jplus[0][1] * monoval_i[1][2] * monoval[2][0];
337  grad_grads[start+n_curls][0][0][2] = -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][1];
338  grad_grads[start+n_curls][0][1][0] = -monoval_jplus[0][1] * monoval_i[1][2] * monoval[2][0];
339  grad_grads[start+n_curls][0][1][1] = -monoval_jplus[0][0] * monoval_i[1][3] * monoval[2][0];
340  grad_grads[start+n_curls][0][1][2] = -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][1];
341  grad_grads[start+n_curls][0][2][0] = -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][1];
342  grad_grads[start+n_curls][0][2][1] = -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][1];
343  grad_grads[start+n_curls][0][2][2] = -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][2];
344  grad_grads[start+n_curls][1][0][0] = monoval_j[0][2] * monoval_i[1][0] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
345  grad_grads[start+n_curls][1][0][1] = monoval_j[0][1] * monoval_i[1][1] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
346  grad_grads[start+n_curls][1][0][2] = monoval_j[0][1] * monoval_i[1][0] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
347  grad_grads[start+n_curls][1][1][0] = monoval_j[0][1] * monoval_i[1][1] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
348  grad_grads[start+n_curls][1][1][1] = monoval_j[0][0] * monoval_i[1][2] * monoval[2][0] * static_cast<double>(j + my_degree + 2);
349  grad_grads[start+n_curls][1][1][2] = monoval_j[0][0] * monoval_i[1][1] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
350  grad_grads[start+n_curls][1][2][0] = monoval_j[0][1] * monoval_i[1][0] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
351  grad_grads[start+n_curls][1][2][1] = monoval_j[0][0] * monoval_i[1][1] * monoval[2][1] * static_cast<double>(j + my_degree + 2);
352  grad_grads[start+n_curls][1][2][2] = monoval_j[0][0] * monoval_i[1][0] * monoval[2][2] * static_cast<double>(j + my_degree + 2);
353  grad_grads[start+n_curls][2][0][0] = -monoval_j[0][2] * monoval_i[1][1] * monoval_plus[2][0];
354  grad_grads[start+n_curls][2][0][1] = -monoval_j[0][1] * monoval_i[1][2] * monoval_plus[2][0];
355  grad_grads[start+n_curls][2][0][2] = -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][1];
356  grad_grads[start+n_curls][2][1][0] = -monoval_j[0][1] * monoval_i[1][2] * monoval_plus[2][0];
357  grad_grads[start+n_curls][2][1][1] = -monoval_j[0][0] * monoval_i[1][3] * monoval_plus[2][0];
358  grad_grads[start+n_curls][2][1][2] = -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][1];
359  grad_grads[start+n_curls][2][2][0] = -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][1];
360  grad_grads[start+n_curls][2][2][1] = -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][1];
361  grad_grads[start+n_curls][2][2][2] = -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][2];
362 
363  grad_grads[start+2*n_curls][0][0][0] = -monoval_jplus[0][2] * monoval[1][0] * monoval_i[2][1];
364  grad_grads[start+2*n_curls][0][0][1] = -monoval_jplus[0][1] * monoval[1][1] * monoval_i[2][1];
365  grad_grads[start+2*n_curls][0][0][2] = -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][2];
366  grad_grads[start+2*n_curls][0][1][0] = -monoval_jplus[0][1] * monoval[1][1] * monoval_i[2][1];
367  grad_grads[start+2*n_curls][0][1][1] = -monoval_jplus[0][0] * monoval[1][2] * monoval_i[2][1];
368  grad_grads[start+2*n_curls][0][1][2] = -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][2];
369  grad_grads[start+2*n_curls][0][2][0] = -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][2];
370  grad_grads[start+2*n_curls][0][2][1] = -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][2];
371  grad_grads[start+2*n_curls][0][2][2] = -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][3];
372  grad_grads[start+2*n_curls][1][0][0] = -monoval_j[0][2] * monoval_plus[1][0] * monoval_i[2][1];
373  grad_grads[start+2*n_curls][1][0][1] = -monoval_j[0][1] * monoval_plus[1][1] * monoval_i[2][1];
374  grad_grads[start+2*n_curls][1][0][2] = -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][2];
375  grad_grads[start+2*n_curls][1][1][0] = -monoval_j[0][1] * monoval_plus[1][1] * monoval_i[2][1];
376  grad_grads[start+2*n_curls][1][1][1] = -monoval_j[0][0] * monoval_plus[1][2] * monoval_i[2][1];
377  grad_grads[start+2*n_curls][1][1][2] = -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][2];
378  grad_grads[start+2*n_curls][1][2][0] = -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][2];
379  grad_grads[start+2*n_curls][1][2][1] = -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][2];
380  grad_grads[start+2*n_curls][1][2][2] = -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][3];
381  grad_grads[start+2*n_curls][2][0][0] = monoval_j[0][2] * monoval[1][0] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
382  grad_grads[start+2*n_curls][2][0][1] = monoval_j[0][1] * monoval[1][1] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
383  grad_grads[start+2*n_curls][2][0][2] = monoval_j[0][1] * monoval[1][0] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
384  grad_grads[start+2*n_curls][2][1][0] = monoval_j[0][1] * monoval[1][1] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
385  grad_grads[start+2*n_curls][2][1][1] = monoval_j[0][0] * monoval[1][2] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
386  grad_grads[start+2*n_curls][2][1][2] = monoval_j[0][0] * monoval[1][1] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
387  grad_grads[start+2*n_curls][2][2][0] = monoval_j[0][1] * monoval[1][0] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
388  grad_grads[start+2*n_curls][2][2][1] = monoval_j[0][0] * monoval[1][1] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
389  grad_grads[start+2*n_curls][2][2][2] = monoval_j[0][0] * monoval[1][0] * monoval_i[2][2] * static_cast<double>(j + my_degree + 2);
390 
391  if (j != my_degree)
392  {
393  grad_grads[start+1][0][0][0] = monoval_i[0][2] * monoval[1][0] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
394  grad_grads[start+1][0][0][1] = monoval_i[0][1] * monoval[1][1] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
395  grad_grads[start+1][0][0][2] = monoval_i[0][1] * monoval[1][0] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
396  grad_grads[start+1][0][1][0] = monoval_i[0][1] * monoval[1][1] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
397  grad_grads[start+1][0][1][1] = monoval_i[0][0] * monoval[1][2] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
398  grad_grads[start+1][0][1][2] = monoval_i[0][0] * monoval[1][1] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
399  grad_grads[start+1][0][2][0] = monoval_i[0][1] * monoval[1][0] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
400  grad_grads[start+1][0][2][1] = monoval_i[0][0] * monoval[1][1] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
401  grad_grads[start+1][0][2][2] = monoval_i[0][0] * monoval[1][0] * monoval_j[2][2] * static_cast<double>(j + my_degree + 2);
402  grad_grads[start+1][1][0][0] = -monoval_i[0][3] * monoval_plus[1][0] * monoval_j[2][0];
403  grad_grads[start+1][1][0][1] = -monoval_i[0][2] * monoval_plus[1][1] * monoval_j[2][0];
404  grad_grads[start+1][1][0][2] = -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][1];
405  grad_grads[start+1][1][1][0] = -monoval_i[0][2] * monoval_plus[1][1] * monoval_j[2][0];
406  grad_grads[start+1][1][1][1] = -monoval_i[0][1] * monoval_plus[1][2] * monoval_j[2][0];
407  grad_grads[start+1][1][1][2] = -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][1];
408  grad_grads[start+1][1][2][0] = -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][1];
409  grad_grads[start+1][1][2][1] = -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][1];
410  grad_grads[start+1][1][2][2] = -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][2];
411  grad_grads[start+1][2][0][0] = -monoval_i[0][3] * monoval[1][0] * monoval_jplus[2][0];
412  grad_grads[start+1][2][0][1] = -monoval_i[0][2] * monoval[1][1] * monoval_jplus[2][0];
413  grad_grads[start+1][2][0][2] = -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][1];
414  grad_grads[start+1][2][1][0] = -monoval_i[0][2] * monoval[1][1] * monoval_jplus[2][0];
415  grad_grads[start+1][2][1][1] = -monoval_i[0][1] * monoval[1][2] * monoval_jplus[2][0];
416  grad_grads[start+1][2][1][2] = -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][1];
417  grad_grads[start+1][2][2][0] = -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][1];
418  grad_grads[start+1][2][2][1] = -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][1];
419  grad_grads[start+1][2][2][2] = -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][2];
420 
421  grad_grads[start+n_curls+1][0][0][0] = -monoval_plus[0][2] * monoval_i[1][1] * monoval_j[2][0];
422  grad_grads[start+n_curls+1][0][0][1] = -monoval_plus[0][1] * monoval_i[1][2] * monoval_j[2][0];
423  grad_grads[start+n_curls+1][0][0][2] = -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][1];
424  grad_grads[start+n_curls+1][0][1][0] = -monoval_plus[0][1] * monoval_i[1][2] * monoval_j[2][0];
425  grad_grads[start+n_curls+1][0][1][1] = -monoval_plus[0][0] * monoval_i[1][3] * monoval_j[2][0];
426  grad_grads[start+n_curls+1][0][1][2] = -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][1];
427  grad_grads[start+n_curls+1][0][2][0] = -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][1];
428  grad_grads[start+n_curls+1][0][2][1] = -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][1];
429  grad_grads[start+n_curls+1][0][2][2] = -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][2];
430  grad_grads[start+n_curls+1][1][0][0] = monoval[0][2] * monoval_i[1][0] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
431  grad_grads[start+n_curls+1][1][0][1] = monoval[0][1] * monoval_i[1][1] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
432  grad_grads[start+n_curls+1][1][0][2] = monoval[0][1] * monoval_i[1][0] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
433  grad_grads[start+n_curls+1][1][1][0] = monoval[0][1] * monoval_i[1][1] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
434  grad_grads[start+n_curls+1][1][1][1] = monoval[0][0] * monoval_i[1][2] * monoval_j[2][0] * static_cast<double>(j + my_degree + 2);
435  grad_grads[start+n_curls+1][1][1][2] = monoval[0][0] * monoval_i[1][1] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
436  grad_grads[start+n_curls+1][1][2][0] = monoval[0][1] * monoval_i[1][0] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
437  grad_grads[start+n_curls+1][1][2][1] = monoval[0][0] * monoval_i[1][1] * monoval_j[2][1] * static_cast<double>(j + my_degree + 2);
438  grad_grads[start+n_curls+1][1][2][2] = monoval[0][0] * monoval_i[1][0] * monoval_j[2][2] * static_cast<double>(j + my_degree + 2);
439  grad_grads[start+n_curls+1][2][0][0] = -monoval[0][2] * monoval_i[1][1] * monoval_jplus[2][0];
440  grad_grads[start+n_curls+1][2][0][1] = -monoval[0][1] * monoval_i[1][2] * monoval_jplus[2][0];
441  grad_grads[start+n_curls+1][2][0][2] = -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][1];
442  grad_grads[start+n_curls+1][2][1][0] = -monoval[0][1] * monoval_i[1][2] * monoval_jplus[2][0];
443  grad_grads[start+n_curls+1][2][1][1] = -monoval[0][0] * monoval_i[1][3] * monoval_jplus[2][0];
444  grad_grads[start+n_curls+1][2][1][2] = -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][1];
445  grad_grads[start+n_curls+1][2][2][0] = -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][1];
446  grad_grads[start+n_curls+1][2][2][1] = -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][1];
447  grad_grads[start+n_curls+1][2][2][2] = -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][2];
448 
449  grad_grads[start+2*n_curls+1][0][0][0] = -monoval_plus[0][2] * monoval_j[1][0] * monoval_i[2][1];
450  grad_grads[start+2*n_curls+1][0][0][1] = -monoval_plus[0][1] * monoval_j[1][1] * monoval_i[2][1];
451  grad_grads[start+2*n_curls+1][0][0][2] = -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][2];
452  grad_grads[start+2*n_curls+1][0][1][0] = -monoval_plus[0][1] * monoval_j[1][1] * monoval_i[2][1];
453  grad_grads[start+2*n_curls+1][0][1][1] = -monoval_plus[0][0] * monoval_j[1][2] * monoval_i[2][1];
454  grad_grads[start+2*n_curls+1][0][1][2] = -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][2];
455  grad_grads[start+2*n_curls+1][0][2][0] = -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][2];
456  grad_grads[start+2*n_curls+1][0][2][1] = -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][2];
457  grad_grads[start+2*n_curls+1][0][2][2] = -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][3];
458  grad_grads[start+2*n_curls+1][1][0][0] = -monoval[0][2] * monoval_jplus[1][0] * monoval_i[2][1];
459  grad_grads[start+2*n_curls+1][1][0][1] = -monoval[0][1] * monoval_jplus[1][1] * monoval_i[2][1];
460  grad_grads[start+2*n_curls+1][1][0][2] = -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][2];
461  grad_grads[start+2*n_curls+1][1][1][0] = -monoval[0][1] * monoval_jplus[1][1] * monoval_i[2][1];
462  grad_grads[start+2*n_curls+1][1][1][1] = -monoval[0][0] * monoval_jplus[1][2] * monoval_i[2][1];
463  grad_grads[start+2*n_curls+1][1][1][2] = -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][2];
464  grad_grads[start+2*n_curls+1][1][2][0] = -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][2];
465  grad_grads[start+2*n_curls+1][1][2][1] = -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][2];
466  grad_grads[start+2*n_curls+1][1][2][2] = -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][3];
467  grad_grads[start+2*n_curls+1][2][0][0] = monoval[0][2] * monoval_j[1][0] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
468  grad_grads[start+2*n_curls+1][2][0][1] = monoval[0][1] * monoval_j[1][1] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
469  grad_grads[start+2*n_curls+1][2][0][2] = monoval[0][1] * monoval_j[1][0] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
470  grad_grads[start+2*n_curls+1][2][1][0] = monoval[0][1] * monoval_j[1][1] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
471  grad_grads[start+2*n_curls+1][2][1][1] = monoval[0][0] * monoval_j[1][2] * monoval_i[2][0] * static_cast<double>(j + my_degree + 2);
472  grad_grads[start+2*n_curls+1][2][1][2] = monoval[0][0] * monoval_j[1][1] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
473  grad_grads[start+2*n_curls+1][2][2][0] = monoval[0][1] * monoval_j[1][0] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
474  grad_grads[start+2*n_curls+1][2][2][1] = monoval[0][0] * monoval_j[1][1] * monoval_i[2][1] * static_cast<double>(j + my_degree + 2);
475  grad_grads[start+2*n_curls+1][2][2][2] = monoval[0][0] * monoval_j[1][0] * monoval_i[2][2] * static_cast<double>(j + my_degree + 2);
476  }
477  }
478 
479  if (j == my_degree)
480  start += 1;
481  else
482  start += 2;
483  }
484  }
485  Assert(start == n_pols - 2*n_curls, ExcInternalError());
486  }
487 }
488 
489 
490 
491 template <int dim>
492 unsigned int
493 PolynomialsRT_Bubbles<dim>::compute_n_pols(const unsigned int k)
494 {
495  if (dim == 1 || dim == 2 || dim == 3)
496  return dim * Utilities::fixed_power<dim>(k+1);
497 
498  Assert(false, ExcNotImplemented());
499  return 0;
500 }
501 
502 
503 template class PolynomialsRT_Bubbles<1>;
504 template class PolynomialsRT_Bubbles<2>;
505 template class PolynomialsRT_Bubbles<3>;
506 
507 
508 DEAL_II_NAMESPACE_CLOSE
PolynomialsRT_Bubbles::compute_n_pols
static unsigned int compute_n_pols(const unsigned int degree)
Definition: polynomials_rt_bubbles.cc:493
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
PolynomialsRT_Bubbles::PolynomialsRT_Bubbles
PolynomialsRT_Bubbles(const unsigned int k)
Definition: polynomials_rt_bubbles.cc:27
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1142
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
PolynomialsRT_Bubbles::compute
void compute(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
Definition: polynomials_rt_bubbles.cc:44
PolynomialsRT_Bubbles::monomials
std::vector< Polynomials::Polynomial< double > > monomials
Definition: polynomials_rt_bubbles.h:149
PolynomialsRT_Bubbles
Definition: polynomials_rt_bubbles.h:83
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.14