Reference documentation for deal.II version 9.1.0-pre
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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include "deal.II/base/polynomials_rt_bubbles.h"
18 
19 #include <deal.II/base/quadrature_lib.h>
20 #include <deal.II/base/thread_management.h>
21 
22 #include <iomanip>
23 #include <iostream>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 template <int dim>
30  : my_degree(k)
31  , raviart_thomas_space(k - 1)
32  , monomials(k + 2)
33  , n_pols(compute_n_pols(k))
34 {
35  Assert(dim >= 2, ExcImpossibleInDim(dim));
36 
37  for (unsigned int i = 0; i < monomials.size(); ++i)
39 }
40 
41 
42 
43 template <int dim>
44 void
46  const Point<dim> & unit_point,
47  std::vector<Tensor<1, dim>> &values,
48  std::vector<Tensor<2, dim>> &grads,
49  std::vector<Tensor<3, dim>> &grad_grads,
50  std::vector<Tensor<4, dim>> &third_derivatives,
51  std::vector<Tensor<5, dim>> &fourth_derivatives) const
52 {
53  Assert(values.size() == n_pols || values.size() == 0,
54  ExcDimensionMismatch(values.size(), n_pols));
55  Assert(grads.size() == n_pols || grads.size() == 0,
56  ExcDimensionMismatch(grads.size(), n_pols));
57  Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
58  ExcDimensionMismatch(grad_grads.size(), n_pols));
59  Assert(third_derivatives.size() == n_pols || third_derivatives.size() == 0,
60  ExcDimensionMismatch(third_derivatives.size(), n_pols));
61  Assert(fourth_derivatives.size() == n_pols || fourth_derivatives.size() == 0,
62  ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
63 
64  // Third and fourth derivatives are not implemented
65  (void)third_derivatives;
66  Assert(third_derivatives.size() == 0, ExcNotImplemented());
67  (void)fourth_derivatives;
68  Assert(fourth_derivatives.size() == 0, 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,
91  p_values,
92  p_grads,
93  p_grad_grads,
94  p_third_derivatives,
95  p_fourth_derivatives);
96  for (unsigned int i = 0; i < p_values.size(); ++i)
97  values[i] = p_values[i];
98  for (unsigned int i = 0; i < p_grads.size(); ++i)
99  grads[i] = p_grads[i];
100  for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
101  grad_grads[i] = p_grad_grads[i];
102  }
103 
104  // Next we compute the polynomials and derivatives
105  // of the curl part of the space
106  const unsigned int n_derivatives = 3;
107  double monoval_plus[dim][n_derivatives + 1];
108  double monoval[dim][n_derivatives + 1];
109 
110  double monoval_i[dim][n_derivatives + 1];
111  double monoval_j[dim][n_derivatives + 1];
112  double monoval_jplus[dim][n_derivatives + 1];
113 
114  unsigned int start = n_sub;
115 
116  if (dim == 2)
117  {
118  // In 2d the curl part of the space is spanned by the vectors
119  // of two types. The first one is
120  // [ x^i * [y^(k+1)]' ]
121  // [ -[x^i]' * y^(k+1) ]
122  // The second one can be obtained from the first by a cyclic
123  // rotation of the coordinates.
124  // monoval_i = x^i,
125  // monoval_plus = x^(k+1)
126  for (unsigned int d = 0; d < dim; ++d)
127  monomials[my_degree + 1].value(unit_point(d),
128  n_derivatives,
129  monoval_plus[d]);
130 
131  for (unsigned int i = 0; i <= my_degree; ++i, ++start)
132  {
133  for (unsigned int d = 0; d < dim; ++d)
134  monomials[i].value(unit_point(d), n_derivatives, monoval_i[d]);
135 
136  if (values.size() != 0)
137  {
138  values[start][0] = monoval_i[0][0] * monoval_plus[1][1];
139  values[start][1] = -monoval_i[0][1] * monoval_plus[1][0];
140 
141  values[start + my_degree + 1][0] =
142  -monoval_plus[0][0] * monoval_i[1][1];
143  values[start + my_degree + 1][1] =
144  monoval_plus[0][1] * monoval_i[1][0];
145  }
146 
147  if (grads.size() != 0)
148  {
149  grads[start][0][0] = monoval_i[0][1] * monoval_plus[1][1];
150  grads[start][0][1] = monoval_i[0][0] * monoval_plus[1][2];
151  grads[start][1][0] = -monoval_i[0][2] * monoval_plus[1][0];
152  grads[start][1][1] = -monoval_i[0][1] * monoval_plus[1][1];
153 
154  grads[start + my_degree + 1][0][0] =
155  -monoval_plus[0][1] * monoval_i[1][1];
156  grads[start + my_degree + 1][0][1] =
157  -monoval_plus[0][0] * monoval_i[1][2];
158  grads[start + my_degree + 1][1][0] =
159  monoval_plus[0][2] * monoval_i[1][0];
160  grads[start + my_degree + 1][1][1] =
161  monoval_plus[0][1] * monoval_i[1][1];
162  }
163 
164  if (grad_grads.size() != 0)
165  {
166  grad_grads[start][0][0][0] = monoval_i[0][2] * monoval_plus[1][1];
167  grad_grads[start][0][0][1] = monoval_i[0][1] * monoval_plus[1][2];
168  grad_grads[start][0][1][0] = monoval_i[0][1] * monoval_plus[1][2];
169  grad_grads[start][0][1][1] = monoval_i[0][0] * monoval_plus[1][3];
170  grad_grads[start][1][0][0] =
171  -monoval_i[0][3] * monoval_plus[1][0];
172  grad_grads[start][1][0][1] =
173  -monoval_i[0][2] * monoval_plus[1][1];
174  grad_grads[start][1][1][0] =
175  -monoval_i[0][2] * monoval_plus[1][1];
176  grad_grads[start][1][1][1] =
177  -monoval_i[0][1] * monoval_plus[1][2];
178 
179  grad_grads[start + my_degree + 1][0][0][0] =
180  -monoval_plus[0][2] * monoval_i[1][1];
181  grad_grads[start + my_degree + 1][0][0][1] =
182  -monoval_plus[0][1] * monoval_i[1][2];
183  grad_grads[start + my_degree + 1][0][1][0] =
184  -monoval_plus[0][1] * monoval_i[1][2];
185  grad_grads[start + my_degree + 1][0][1][1] =
186  -monoval_plus[0][0] * monoval_i[1][3];
187  grad_grads[start + my_degree + 1][1][0][0] =
188  monoval_plus[0][3] * monoval_i[1][0];
189  grad_grads[start + my_degree + 1][1][0][1] =
190  monoval_plus[0][2] * monoval_i[1][1];
191  grad_grads[start + my_degree + 1][1][1][0] =
192  monoval_plus[0][2] * monoval_i[1][1];
193  grad_grads[start + my_degree + 1][1][1][1] =
194  monoval_plus[0][1] * monoval_i[1][2];
195  }
196  }
197  Assert(start == n_pols - my_degree - 1, ExcInternalError());
198  }
199  else if (dim == 3)
200  {
201  // In 3d the first type of basis vector is
202  // [ x^i * y^j * z^k * (j+k+2) ]
203  // [ -[x^i]' * y^(j+1) * z^k ]
204  // [ -[x^i]' * y^j * z^(k+1) ],
205  // For the second type of basis vector y and z
206  // are swapped. Then for each of these,
207  // two more are obtained by the cyclic rotation
208  // of the coordinates.
209  // monoval = x^k, monoval_plus = x^(k+1)
210  // monoval_* = x^*, monoval_jplus = x^(j+1)
211  for (unsigned int d = 0; d < dim; ++d)
212  {
213  monomials[my_degree + 1].value(unit_point(d),
214  n_derivatives,
215  monoval_plus[d]);
216  monomials[my_degree].value(unit_point(d), n_derivatives, monoval[d]);
217  }
218 
219  const unsigned int n_curls = (my_degree + 1) * (2 * my_degree + 1);
220  // Span of @f$\tilde{B}@f$
221  for (unsigned int i = 0; i <= my_degree; ++i)
222  {
223  for (unsigned int d = 0; d < dim; ++d)
224  monomials[i].value(unit_point(d), n_derivatives, monoval_i[d]);
225 
226  for (unsigned int j = 0; j <= my_degree; ++j)
227  {
228  for (unsigned int d = 0; d < dim; ++d)
229  {
230  monomials[j].value(unit_point(d),
231  n_derivatives,
232  monoval_j[d]);
233  monomials[j + 1].value(unit_point(d),
234  n_derivatives,
235  monoval_jplus[d]);
236  }
237 
238  if (values.size() != 0)
239  {
240  values[start][0] = monoval_i[0][0] * monoval_j[1][0] *
241  monoval[2][0] *
242  static_cast<double>(j + my_degree + 2);
243  values[start][1] =
244  -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][0];
245  values[start][2] =
246  -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][0];
247 
248  values[start + n_curls][0] =
249  -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][0];
250  values[start + n_curls][1] =
251  monoval_j[0][0] * monoval_i[1][0] * monoval[2][0] *
252  static_cast<double>(j + my_degree + 2);
253  values[start + n_curls][2] =
254  -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][0];
255 
256  values[start + 2 * n_curls][0] =
257  -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][1];
258  values[start + 2 * n_curls][1] =
259  -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][1];
260  values[start + 2 * n_curls][2] =
261  monoval_j[0][0] * monoval[1][0] * monoval_i[2][0] *
262  static_cast<double>(j + my_degree + 2);
263 
264  // Only unique triples of powers (i j k)
265  // and (i k j) are allowed, 0 <= i,j <= k
266  if (j != my_degree)
267  {
268  values[start + 1][0] =
269  monoval_i[0][0] * monoval[1][0] * monoval_j[2][0] *
270  static_cast<double>(j + my_degree + 2);
271  values[start + 1][1] =
272  -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][0];
273  values[start + 1][2] =
274  -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][0];
275 
276  values[start + n_curls + 1][0] =
277  -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][0];
278  values[start + n_curls + 1][1] =
279  monoval[0][0] * monoval_i[1][0] * monoval_j[2][0] *
280  static_cast<double>(j + my_degree + 2);
281  values[start + n_curls + 1][2] =
282  -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][0];
283 
284  values[start + 2 * n_curls + 1][0] =
285  -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][1];
286  values[start + 2 * n_curls + 1][1] =
287  -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][1];
288  values[start + 2 * n_curls + 1][2] =
289  monoval[0][0] * monoval_j[1][0] * monoval_i[2][0] *
290  static_cast<double>(j + my_degree + 2);
291  }
292  }
293 
294  if (grads.size() != 0)
295  {
296  grads[start][0][0] = monoval_i[0][1] * monoval_j[1][0] *
297  monoval[2][0] *
298  static_cast<double>(j + my_degree + 2);
299  grads[start][0][1] = monoval_i[0][0] * monoval_j[1][1] *
300  monoval[2][0] *
301  static_cast<double>(j + my_degree + 2);
302  grads[start][0][2] = monoval_i[0][0] * monoval_j[1][0] *
303  monoval[2][1] *
304  static_cast<double>(j + my_degree + 2);
305  grads[start][1][0] =
306  -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][0];
307  grads[start][1][1] =
308  -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][0];
309  grads[start][1][2] =
310  -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][1];
311  grads[start][2][0] =
312  -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][0];
313  grads[start][2][1] =
314  -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][0];
315  grads[start][2][2] =
316  -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][1];
317 
318  grads[start + n_curls][0][0] =
319  -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][0];
320  grads[start + n_curls][0][1] =
321  -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][0];
322  grads[start + n_curls][0][2] =
323  -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][1];
324  grads[start + n_curls][1][0] =
325  monoval_j[0][1] * monoval_i[1][0] * monoval[2][0] *
326  static_cast<double>(j + my_degree + 2);
327  grads[start + n_curls][1][1] =
328  monoval_j[0][0] * monoval_i[1][1] * monoval[2][0] *
329  static_cast<double>(j + my_degree + 2);
330  grads[start + n_curls][1][2] =
331  monoval_j[0][0] * monoval_i[1][0] * monoval[2][1] *
332  static_cast<double>(j + my_degree + 2);
333  grads[start + n_curls][2][0] =
334  -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][0];
335  grads[start + n_curls][2][1] =
336  -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][0];
337  grads[start + n_curls][2][2] =
338  -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][1];
339 
340  grads[start + 2 * n_curls][0][0] =
341  -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][1];
342  grads[start + 2 * n_curls][0][1] =
343  -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][1];
344  grads[start + 2 * n_curls][0][2] =
345  -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][2];
346  grads[start + 2 * n_curls][1][0] =
347  -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][1];
348  grads[start + 2 * n_curls][1][1] =
349  -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][1];
350  grads[start + 2 * n_curls][1][2] =
351  -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][2];
352  grads[start + 2 * n_curls][2][0] =
353  monoval_j[0][1] * monoval[1][0] * monoval_i[2][0] *
354  static_cast<double>(j + my_degree + 2);
355  grads[start + 2 * n_curls][2][1] =
356  monoval_j[0][0] * monoval[1][1] * monoval_i[2][0] *
357  static_cast<double>(j + my_degree + 2);
358  grads[start + 2 * n_curls][2][2] =
359  monoval_j[0][0] * monoval[1][0] * monoval_i[2][1] *
360  static_cast<double>(j + my_degree + 2);
361 
362  if (j != my_degree)
363  {
364  grads[start + 1][0][0] =
365  monoval_i[0][1] * monoval[1][0] * monoval_j[2][0] *
366  static_cast<double>(j + my_degree + 2);
367  grads[start + 1][0][1] =
368  monoval_i[0][0] * monoval[1][1] * monoval_j[2][0] *
369  static_cast<double>(j + my_degree + 2);
370  grads[start + 1][0][2] =
371  monoval_i[0][0] * monoval[1][0] * monoval_j[2][1] *
372  static_cast<double>(j + my_degree + 2);
373  grads[start + 1][1][0] =
374  -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][0];
375  grads[start + 1][1][1] =
376  -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][0];
377  grads[start + 1][1][2] =
378  -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][1];
379  grads[start + 1][2][0] =
380  -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][0];
381  grads[start + 1][2][1] =
382  -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][0];
383  grads[start + 1][2][2] =
384  -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][1];
385 
386  grads[start + n_curls + 1][0][0] =
387  -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][0];
388  grads[start + n_curls + 1][0][1] =
389  -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][0];
390  grads[start + n_curls + 1][0][2] =
391  -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][1];
392  grads[start + n_curls + 1][1][0] =
393  monoval[0][1] * monoval_i[1][0] * monoval_j[2][0] *
394  static_cast<double>(j + my_degree + 2);
395  grads[start + n_curls + 1][1][1] =
396  monoval[0][0] * monoval_i[1][1] * monoval_j[2][0] *
397  static_cast<double>(j + my_degree + 2);
398  grads[start + n_curls + 1][1][2] =
399  monoval[0][0] * monoval_i[1][0] * monoval_j[2][1] *
400  static_cast<double>(j + my_degree + 2);
401  grads[start + n_curls + 1][2][0] =
402  -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][0];
403  grads[start + n_curls + 1][2][1] =
404  -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][0];
405  grads[start + n_curls + 1][2][2] =
406  -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][1];
407 
408  grads[start + 2 * n_curls + 1][0][0] =
409  -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][1];
410  grads[start + 2 * n_curls + 1][0][1] =
411  -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][1];
412  grads[start + 2 * n_curls + 1][0][2] =
413  -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][2];
414  grads[start + 2 * n_curls + 1][1][0] =
415  -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][1];
416  grads[start + 2 * n_curls + 1][1][1] =
417  -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][1];
418  grads[start + 2 * n_curls + 1][1][2] =
419  -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][2];
420  grads[start + 2 * n_curls + 1][2][0] =
421  monoval[0][1] * monoval_j[1][0] * monoval_i[2][0] *
422  static_cast<double>(j + my_degree + 2);
423  grads[start + 2 * n_curls + 1][2][1] =
424  monoval[0][0] * monoval_j[1][1] * monoval_i[2][0] *
425  static_cast<double>(j + my_degree + 2);
426  grads[start + 2 * n_curls + 1][2][2] =
427  monoval[0][0] * monoval_j[1][0] * monoval_i[2][1] *
428  static_cast<double>(j + my_degree + 2);
429  }
430  }
431 
432  if (grad_grads.size() != 0)
433  {
434  grad_grads[start][0][0][0] =
435  monoval_i[0][2] * monoval_j[1][0] * monoval[2][0] *
436  static_cast<double>(j + my_degree + 2);
437  grad_grads[start][0][0][1] =
438  monoval_i[0][1] * monoval_j[1][1] * monoval[2][0] *
439  static_cast<double>(j + my_degree + 2);
440  grad_grads[start][0][0][2] =
441  monoval_i[0][1] * monoval_j[1][0] * monoval[2][1] *
442  static_cast<double>(j + my_degree + 2);
443  grad_grads[start][0][1][0] =
444  monoval_i[0][1] * monoval_j[1][1] * monoval[2][0] *
445  static_cast<double>(j + my_degree + 2);
446  grad_grads[start][0][1][1] =
447  monoval_i[0][0] * monoval_j[1][2] * monoval[2][0] *
448  static_cast<double>(j + my_degree + 2);
449  grad_grads[start][0][1][2] =
450  monoval_i[0][0] * monoval_j[1][1] * monoval[2][1] *
451  static_cast<double>(j + my_degree + 2);
452  grad_grads[start][0][2][0] =
453  monoval_i[0][1] * monoval_j[1][0] * monoval[2][1] *
454  static_cast<double>(j + my_degree + 2);
455  grad_grads[start][0][2][1] =
456  monoval_i[0][0] * monoval_j[1][1] * monoval[2][1] *
457  static_cast<double>(j + my_degree + 2);
458  grad_grads[start][0][2][2] =
459  monoval_i[0][0] * monoval_j[1][0] * monoval[2][2] *
460  static_cast<double>(j + my_degree + 2);
461  grad_grads[start][1][0][0] =
462  -monoval_i[0][3] * monoval_jplus[1][0] * monoval[2][0];
463  grad_grads[start][1][0][1] =
464  -monoval_i[0][2] * monoval_jplus[1][1] * monoval[2][0];
465  grad_grads[start][1][0][2] =
466  -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][1];
467  grad_grads[start][1][1][0] =
468  -monoval_i[0][2] * monoval_jplus[1][1] * monoval[2][0];
469  grad_grads[start][1][1][1] =
470  -monoval_i[0][1] * monoval_jplus[1][2] * monoval[2][0];
471  grad_grads[start][1][1][2] =
472  -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][1];
473  grad_grads[start][1][2][0] =
474  -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][1];
475  grad_grads[start][1][2][1] =
476  -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][1];
477  grad_grads[start][1][2][2] =
478  -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][2];
479  grad_grads[start][2][0][0] =
480  -monoval_i[0][3] * monoval_j[1][0] * monoval_plus[2][0];
481  grad_grads[start][2][0][1] =
482  -monoval_i[0][2] * monoval_j[1][1] * monoval_plus[2][0];
483  grad_grads[start][2][0][2] =
484  -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][1];
485  grad_grads[start][2][1][0] =
486  -monoval_i[0][2] * monoval_j[1][1] * monoval_plus[2][0];
487  grad_grads[start][2][1][1] =
488  -monoval_i[0][1] * monoval_j[1][2] * monoval_plus[2][0];
489  grad_grads[start][2][1][2] =
490  -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][1];
491  grad_grads[start][2][2][0] =
492  -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][1];
493  grad_grads[start][2][2][1] =
494  -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][1];
495  grad_grads[start][2][2][2] =
496  -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][2];
497 
498  grad_grads[start + n_curls][0][0][0] =
499  -monoval_jplus[0][2] * monoval_i[1][1] * monoval[2][0];
500  grad_grads[start + n_curls][0][0][1] =
501  -monoval_jplus[0][1] * monoval_i[1][2] * monoval[2][0];
502  grad_grads[start + n_curls][0][0][2] =
503  -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][1];
504  grad_grads[start + n_curls][0][1][0] =
505  -monoval_jplus[0][1] * monoval_i[1][2] * monoval[2][0];
506  grad_grads[start + n_curls][0][1][1] =
507  -monoval_jplus[0][0] * monoval_i[1][3] * monoval[2][0];
508  grad_grads[start + n_curls][0][1][2] =
509  -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][1];
510  grad_grads[start + n_curls][0][2][0] =
511  -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][1];
512  grad_grads[start + n_curls][0][2][1] =
513  -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][1];
514  grad_grads[start + n_curls][0][2][2] =
515  -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][2];
516  grad_grads[start + n_curls][1][0][0] =
517  monoval_j[0][2] * monoval_i[1][0] * monoval[2][0] *
518  static_cast<double>(j + my_degree + 2);
519  grad_grads[start + n_curls][1][0][1] =
520  monoval_j[0][1] * monoval_i[1][1] * monoval[2][0] *
521  static_cast<double>(j + my_degree + 2);
522  grad_grads[start + n_curls][1][0][2] =
523  monoval_j[0][1] * monoval_i[1][0] * monoval[2][1] *
524  static_cast<double>(j + my_degree + 2);
525  grad_grads[start + n_curls][1][1][0] =
526  monoval_j[0][1] * monoval_i[1][1] * monoval[2][0] *
527  static_cast<double>(j + my_degree + 2);
528  grad_grads[start + n_curls][1][1][1] =
529  monoval_j[0][0] * monoval_i[1][2] * monoval[2][0] *
530  static_cast<double>(j + my_degree + 2);
531  grad_grads[start + n_curls][1][1][2] =
532  monoval_j[0][0] * monoval_i[1][1] * monoval[2][1] *
533  static_cast<double>(j + my_degree + 2);
534  grad_grads[start + n_curls][1][2][0] =
535  monoval_j[0][1] * monoval_i[1][0] * monoval[2][1] *
536  static_cast<double>(j + my_degree + 2);
537  grad_grads[start + n_curls][1][2][1] =
538  monoval_j[0][0] * monoval_i[1][1] * monoval[2][1] *
539  static_cast<double>(j + my_degree + 2);
540  grad_grads[start + n_curls][1][2][2] =
541  monoval_j[0][0] * monoval_i[1][0] * monoval[2][2] *
542  static_cast<double>(j + my_degree + 2);
543  grad_grads[start + n_curls][2][0][0] =
544  -monoval_j[0][2] * monoval_i[1][1] * monoval_plus[2][0];
545  grad_grads[start + n_curls][2][0][1] =
546  -monoval_j[0][1] * monoval_i[1][2] * monoval_plus[2][0];
547  grad_grads[start + n_curls][2][0][2] =
548  -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][1];
549  grad_grads[start + n_curls][2][1][0] =
550  -monoval_j[0][1] * monoval_i[1][2] * monoval_plus[2][0];
551  grad_grads[start + n_curls][2][1][1] =
552  -monoval_j[0][0] * monoval_i[1][3] * monoval_plus[2][0];
553  grad_grads[start + n_curls][2][1][2] =
554  -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][1];
555  grad_grads[start + n_curls][2][2][0] =
556  -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][1];
557  grad_grads[start + n_curls][2][2][1] =
558  -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][1];
559  grad_grads[start + n_curls][2][2][2] =
560  -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][2];
561 
562  grad_grads[start + 2 * n_curls][0][0][0] =
563  -monoval_jplus[0][2] * monoval[1][0] * monoval_i[2][1];
564  grad_grads[start + 2 * n_curls][0][0][1] =
565  -monoval_jplus[0][1] * monoval[1][1] * monoval_i[2][1];
566  grad_grads[start + 2 * n_curls][0][0][2] =
567  -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][2];
568  grad_grads[start + 2 * n_curls][0][1][0] =
569  -monoval_jplus[0][1] * monoval[1][1] * monoval_i[2][1];
570  grad_grads[start + 2 * n_curls][0][1][1] =
571  -monoval_jplus[0][0] * monoval[1][2] * monoval_i[2][1];
572  grad_grads[start + 2 * n_curls][0][1][2] =
573  -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][2];
574  grad_grads[start + 2 * n_curls][0][2][0] =
575  -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][2];
576  grad_grads[start + 2 * n_curls][0][2][1] =
577  -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][2];
578  grad_grads[start + 2 * n_curls][0][2][2] =
579  -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][3];
580  grad_grads[start + 2 * n_curls][1][0][0] =
581  -monoval_j[0][2] * monoval_plus[1][0] * monoval_i[2][1];
582  grad_grads[start + 2 * n_curls][1][0][1] =
583  -monoval_j[0][1] * monoval_plus[1][1] * monoval_i[2][1];
584  grad_grads[start + 2 * n_curls][1][0][2] =
585  -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][2];
586  grad_grads[start + 2 * n_curls][1][1][0] =
587  -monoval_j[0][1] * monoval_plus[1][1] * monoval_i[2][1];
588  grad_grads[start + 2 * n_curls][1][1][1] =
589  -monoval_j[0][0] * monoval_plus[1][2] * monoval_i[2][1];
590  grad_grads[start + 2 * n_curls][1][1][2] =
591  -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][2];
592  grad_grads[start + 2 * n_curls][1][2][0] =
593  -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][2];
594  grad_grads[start + 2 * n_curls][1][2][1] =
595  -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][2];
596  grad_grads[start + 2 * n_curls][1][2][2] =
597  -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][3];
598  grad_grads[start + 2 * n_curls][2][0][0] =
599  monoval_j[0][2] * monoval[1][0] * monoval_i[2][0] *
600  static_cast<double>(j + my_degree + 2);
601  grad_grads[start + 2 * n_curls][2][0][1] =
602  monoval_j[0][1] * monoval[1][1] * monoval_i[2][0] *
603  static_cast<double>(j + my_degree + 2);
604  grad_grads[start + 2 * n_curls][2][0][2] =
605  monoval_j[0][1] * monoval[1][0] * monoval_i[2][1] *
606  static_cast<double>(j + my_degree + 2);
607  grad_grads[start + 2 * n_curls][2][1][0] =
608  monoval_j[0][1] * monoval[1][1] * monoval_i[2][0] *
609  static_cast<double>(j + my_degree + 2);
610  grad_grads[start + 2 * n_curls][2][1][1] =
611  monoval_j[0][0] * monoval[1][2] * monoval_i[2][0] *
612  static_cast<double>(j + my_degree + 2);
613  grad_grads[start + 2 * n_curls][2][1][2] =
614  monoval_j[0][0] * monoval[1][1] * monoval_i[2][1] *
615  static_cast<double>(j + my_degree + 2);
616  grad_grads[start + 2 * n_curls][2][2][0] =
617  monoval_j[0][1] * monoval[1][0] * monoval_i[2][1] *
618  static_cast<double>(j + my_degree + 2);
619  grad_grads[start + 2 * n_curls][2][2][1] =
620  monoval_j[0][0] * monoval[1][1] * monoval_i[2][1] *
621  static_cast<double>(j + my_degree + 2);
622  grad_grads[start + 2 * n_curls][2][2][2] =
623  monoval_j[0][0] * monoval[1][0] * monoval_i[2][2] *
624  static_cast<double>(j + my_degree + 2);
625 
626  if (j != my_degree)
627  {
628  grad_grads[start + 1][0][0][0] =
629  monoval_i[0][2] * monoval[1][0] * monoval_j[2][0] *
630  static_cast<double>(j + my_degree + 2);
631  grad_grads[start + 1][0][0][1] =
632  monoval_i[0][1] * monoval[1][1] * monoval_j[2][0] *
633  static_cast<double>(j + my_degree + 2);
634  grad_grads[start + 1][0][0][2] =
635  monoval_i[0][1] * monoval[1][0] * monoval_j[2][1] *
636  static_cast<double>(j + my_degree + 2);
637  grad_grads[start + 1][0][1][0] =
638  monoval_i[0][1] * monoval[1][1] * monoval_j[2][0] *
639  static_cast<double>(j + my_degree + 2);
640  grad_grads[start + 1][0][1][1] =
641  monoval_i[0][0] * monoval[1][2] * monoval_j[2][0] *
642  static_cast<double>(j + my_degree + 2);
643  grad_grads[start + 1][0][1][2] =
644  monoval_i[0][0] * monoval[1][1] * monoval_j[2][1] *
645  static_cast<double>(j + my_degree + 2);
646  grad_grads[start + 1][0][2][0] =
647  monoval_i[0][1] * monoval[1][0] * monoval_j[2][1] *
648  static_cast<double>(j + my_degree + 2);
649  grad_grads[start + 1][0][2][1] =
650  monoval_i[0][0] * monoval[1][1] * monoval_j[2][1] *
651  static_cast<double>(j + my_degree + 2);
652  grad_grads[start + 1][0][2][2] =
653  monoval_i[0][0] * monoval[1][0] * monoval_j[2][2] *
654  static_cast<double>(j + my_degree + 2);
655  grad_grads[start + 1][1][0][0] =
656  -monoval_i[0][3] * monoval_plus[1][0] * monoval_j[2][0];
657  grad_grads[start + 1][1][0][1] =
658  -monoval_i[0][2] * monoval_plus[1][1] * monoval_j[2][0];
659  grad_grads[start + 1][1][0][2] =
660  -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][1];
661  grad_grads[start + 1][1][1][0] =
662  -monoval_i[0][2] * monoval_plus[1][1] * monoval_j[2][0];
663  grad_grads[start + 1][1][1][1] =
664  -monoval_i[0][1] * monoval_plus[1][2] * monoval_j[2][0];
665  grad_grads[start + 1][1][1][2] =
666  -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][1];
667  grad_grads[start + 1][1][2][0] =
668  -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][1];
669  grad_grads[start + 1][1][2][1] =
670  -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][1];
671  grad_grads[start + 1][1][2][2] =
672  -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][2];
673  grad_grads[start + 1][2][0][0] =
674  -monoval_i[0][3] * monoval[1][0] * monoval_jplus[2][0];
675  grad_grads[start + 1][2][0][1] =
676  -monoval_i[0][2] * monoval[1][1] * monoval_jplus[2][0];
677  grad_grads[start + 1][2][0][2] =
678  -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][1];
679  grad_grads[start + 1][2][1][0] =
680  -monoval_i[0][2] * monoval[1][1] * monoval_jplus[2][0];
681  grad_grads[start + 1][2][1][1] =
682  -monoval_i[0][1] * monoval[1][2] * monoval_jplus[2][0];
683  grad_grads[start + 1][2][1][2] =
684  -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][1];
685  grad_grads[start + 1][2][2][0] =
686  -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][1];
687  grad_grads[start + 1][2][2][1] =
688  -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][1];
689  grad_grads[start + 1][2][2][2] =
690  -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][2];
691 
692  grad_grads[start + n_curls + 1][0][0][0] =
693  -monoval_plus[0][2] * monoval_i[1][1] * monoval_j[2][0];
694  grad_grads[start + n_curls + 1][0][0][1] =
695  -monoval_plus[0][1] * monoval_i[1][2] * monoval_j[2][0];
696  grad_grads[start + n_curls + 1][0][0][2] =
697  -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][1];
698  grad_grads[start + n_curls + 1][0][1][0] =
699  -monoval_plus[0][1] * monoval_i[1][2] * monoval_j[2][0];
700  grad_grads[start + n_curls + 1][0][1][1] =
701  -monoval_plus[0][0] * monoval_i[1][3] * monoval_j[2][0];
702  grad_grads[start + n_curls + 1][0][1][2] =
703  -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][1];
704  grad_grads[start + n_curls + 1][0][2][0] =
705  -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][1];
706  grad_grads[start + n_curls + 1][0][2][1] =
707  -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][1];
708  grad_grads[start + n_curls + 1][0][2][2] =
709  -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][2];
710  grad_grads[start + n_curls + 1][1][0][0] =
711  monoval[0][2] * monoval_i[1][0] * monoval_j[2][0] *
712  static_cast<double>(j + my_degree + 2);
713  grad_grads[start + n_curls + 1][1][0][1] =
714  monoval[0][1] * monoval_i[1][1] * monoval_j[2][0] *
715  static_cast<double>(j + my_degree + 2);
716  grad_grads[start + n_curls + 1][1][0][2] =
717  monoval[0][1] * monoval_i[1][0] * monoval_j[2][1] *
718  static_cast<double>(j + my_degree + 2);
719  grad_grads[start + n_curls + 1][1][1][0] =
720  monoval[0][1] * monoval_i[1][1] * monoval_j[2][0] *
721  static_cast<double>(j + my_degree + 2);
722  grad_grads[start + n_curls + 1][1][1][1] =
723  monoval[0][0] * monoval_i[1][2] * monoval_j[2][0] *
724  static_cast<double>(j + my_degree + 2);
725  grad_grads[start + n_curls + 1][1][1][2] =
726  monoval[0][0] * monoval_i[1][1] * monoval_j[2][1] *
727  static_cast<double>(j + my_degree + 2);
728  grad_grads[start + n_curls + 1][1][2][0] =
729  monoval[0][1] * monoval_i[1][0] * monoval_j[2][1] *
730  static_cast<double>(j + my_degree + 2);
731  grad_grads[start + n_curls + 1][1][2][1] =
732  monoval[0][0] * monoval_i[1][1] * monoval_j[2][1] *
733  static_cast<double>(j + my_degree + 2);
734  grad_grads[start + n_curls + 1][1][2][2] =
735  monoval[0][0] * monoval_i[1][0] * monoval_j[2][2] *
736  static_cast<double>(j + my_degree + 2);
737  grad_grads[start + n_curls + 1][2][0][0] =
738  -monoval[0][2] * monoval_i[1][1] * monoval_jplus[2][0];
739  grad_grads[start + n_curls + 1][2][0][1] =
740  -monoval[0][1] * monoval_i[1][2] * monoval_jplus[2][0];
741  grad_grads[start + n_curls + 1][2][0][2] =
742  -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][1];
743  grad_grads[start + n_curls + 1][2][1][0] =
744  -monoval[0][1] * monoval_i[1][2] * monoval_jplus[2][0];
745  grad_grads[start + n_curls + 1][2][1][1] =
746  -monoval[0][0] * monoval_i[1][3] * monoval_jplus[2][0];
747  grad_grads[start + n_curls + 1][2][1][2] =
748  -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][1];
749  grad_grads[start + n_curls + 1][2][2][0] =
750  -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][1];
751  grad_grads[start + n_curls + 1][2][2][1] =
752  -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][1];
753  grad_grads[start + n_curls + 1][2][2][2] =
754  -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][2];
755 
756  grad_grads[start + 2 * n_curls + 1][0][0][0] =
757  -monoval_plus[0][2] * monoval_j[1][0] * monoval_i[2][1];
758  grad_grads[start + 2 * n_curls + 1][0][0][1] =
759  -monoval_plus[0][1] * monoval_j[1][1] * monoval_i[2][1];
760  grad_grads[start + 2 * n_curls + 1][0][0][2] =
761  -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][2];
762  grad_grads[start + 2 * n_curls + 1][0][1][0] =
763  -monoval_plus[0][1] * monoval_j[1][1] * monoval_i[2][1];
764  grad_grads[start + 2 * n_curls + 1][0][1][1] =
765  -monoval_plus[0][0] * monoval_j[1][2] * monoval_i[2][1];
766  grad_grads[start + 2 * n_curls + 1][0][1][2] =
767  -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][2];
768  grad_grads[start + 2 * n_curls + 1][0][2][0] =
769  -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][2];
770  grad_grads[start + 2 * n_curls + 1][0][2][1] =
771  -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][2];
772  grad_grads[start + 2 * n_curls + 1][0][2][2] =
773  -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][3];
774  grad_grads[start + 2 * n_curls + 1][1][0][0] =
775  -monoval[0][2] * monoval_jplus[1][0] * monoval_i[2][1];
776  grad_grads[start + 2 * n_curls + 1][1][0][1] =
777  -monoval[0][1] * monoval_jplus[1][1] * monoval_i[2][1];
778  grad_grads[start + 2 * n_curls + 1][1][0][2] =
779  -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][2];
780  grad_grads[start + 2 * n_curls + 1][1][1][0] =
781  -monoval[0][1] * monoval_jplus[1][1] * monoval_i[2][1];
782  grad_grads[start + 2 * n_curls + 1][1][1][1] =
783  -monoval[0][0] * monoval_jplus[1][2] * monoval_i[2][1];
784  grad_grads[start + 2 * n_curls + 1][1][1][2] =
785  -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][2];
786  grad_grads[start + 2 * n_curls + 1][1][2][0] =
787  -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][2];
788  grad_grads[start + 2 * n_curls + 1][1][2][1] =
789  -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][2];
790  grad_grads[start + 2 * n_curls + 1][1][2][2] =
791  -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][3];
792  grad_grads[start + 2 * n_curls + 1][2][0][0] =
793  monoval[0][2] * monoval_j[1][0] * monoval_i[2][0] *
794  static_cast<double>(j + my_degree + 2);
795  grad_grads[start + 2 * n_curls + 1][2][0][1] =
796  monoval[0][1] * monoval_j[1][1] * monoval_i[2][0] *
797  static_cast<double>(j + my_degree + 2);
798  grad_grads[start + 2 * n_curls + 1][2][0][2] =
799  monoval[0][1] * monoval_j[1][0] * monoval_i[2][1] *
800  static_cast<double>(j + my_degree + 2);
801  grad_grads[start + 2 * n_curls + 1][2][1][0] =
802  monoval[0][1] * monoval_j[1][1] * monoval_i[2][0] *
803  static_cast<double>(j + my_degree + 2);
804  grad_grads[start + 2 * n_curls + 1][2][1][1] =
805  monoval[0][0] * monoval_j[1][2] * monoval_i[2][0] *
806  static_cast<double>(j + my_degree + 2);
807  grad_grads[start + 2 * n_curls + 1][2][1][2] =
808  monoval[0][0] * monoval_j[1][1] * monoval_i[2][1] *
809  static_cast<double>(j + my_degree + 2);
810  grad_grads[start + 2 * n_curls + 1][2][2][0] =
811  monoval[0][1] * monoval_j[1][0] * monoval_i[2][1] *
812  static_cast<double>(j + my_degree + 2);
813  grad_grads[start + 2 * n_curls + 1][2][2][1] =
814  monoval[0][0] * monoval_j[1][1] * monoval_i[2][1] *
815  static_cast<double>(j + my_degree + 2);
816  grad_grads[start + 2 * n_curls + 1][2][2][2] =
817  monoval[0][0] * monoval_j[1][0] * monoval_i[2][2] *
818  static_cast<double>(j + my_degree + 2);
819  }
820  }
821 
822  if (j == my_degree)
823  start += 1;
824  else
825  start += 2;
826  }
827  }
828  Assert(start == n_pols - 2 * n_curls, ExcInternalError());
829  }
830 }
831 
832 
833 
834 template <int dim>
835 unsigned int
837 {
838  if (dim == 1 || dim == 2 || dim == 3)
839  return dim * Utilities::fixed_power<dim>(k + 1);
840 
841  Assert(false, ExcNotImplemented());
842  return 0;
843 }
844 
845 
846 template class PolynomialsRT_Bubbles<1>;
847 template class PolynomialsRT_Bubbles<2>;
848 template class PolynomialsRT_Bubbles<3>;
849 
850 
851 DEAL_II_NAMESPACE_CLOSE
const unsigned int my_degree
static unsigned int compute_n_pols(const unsigned int degree)
static::ExceptionBase & ExcImpossibleInDim(int arg1)
PolynomialsRT_Bubbles(const unsigned int k)
#define Assert(cond, exc)
Definition: exceptions.h:1227
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
const PolynomialsRaviartThomas< dim > raviart_thomas_space
static::ExceptionBase & ExcNotImplemented()
std::vector< Polynomials::Polynomial< double > > monomials
static::ExceptionBase & ExcInternalError()