Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
polynomials_nedelec.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2013 - 2019 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 #include <deal.II/base/geometry_info.h>
17 #include <deal.II/base/polynomial.h>
18 #include <deal.II/base/polynomials_nedelec.h>
19 #include <deal.II/base/quadrature_lib.h>
20 #include <deal.II/base/std_cxx14/memory.h>
21 
22 #include <iomanip>
23 #include <iostream>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 template <int dim>
29 PolynomialsNedelec<dim>::PolynomialsNedelec(const unsigned int k)
30  : TensorPolynomialsBase<dim>(k, n_polynomials(k))
31  , polynomial_space(create_polynomials(k))
32 {}
33 
34 template <int dim>
35 std::vector<std::vector<Polynomials::Polynomial<double>>>
36 PolynomialsNedelec<dim>::create_polynomials(const unsigned int k)
37 {
38  std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
39 
40  pols[0] = Polynomials::Legendre::generate_complete_basis(k);
41 
42  for (unsigned int i = 1; i < dim; ++i)
43  pols[i] = Polynomials::Lobatto::generate_complete_basis(k + 1);
44 
45  return pols;
46 }
47 
48 
49 // Compute the values, gradients
50 // and double gradients of the
51 // polynomial at the given point.
52 template <int dim>
53 void
54 PolynomialsNedelec<dim>::evaluate(
55  const Point<dim> & unit_point,
56  std::vector<Tensor<1, dim>> &values,
57  std::vector<Tensor<2, dim>> &grads,
58  std::vector<Tensor<3, dim>> &grad_grads,
59  std::vector<Tensor<4, dim>> &third_derivatives,
60  std::vector<Tensor<5, dim>> &fourth_derivatives) const
61 {
62  Assert(values.size() == this->n() || values.size() == 0,
63  ExcDimensionMismatch(values.size(), this->n()));
64  Assert(grads.size() == this->n() || grads.size() == 0,
65  ExcDimensionMismatch(grads.size(), this->n()));
66  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
67  ExcDimensionMismatch(grad_grads.size(), this->n()));
68  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
69  ExcDimensionMismatch(third_derivatives.size(), this->n()));
70  Assert(fourth_derivatives.size() == this->n() ||
71  fourth_derivatives.size() == 0,
72  ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
73 
74  // third and fourth derivatives not implemented
75  (void)third_derivatives;
76  Assert(third_derivatives.size() == 0, ExcNotImplemented());
77  (void)fourth_derivatives;
78  Assert(fourth_derivatives.size() == 0, ExcNotImplemented());
79 
80  // Declare the values, derivatives
81  // and second derivatives vectors of
82  // <tt>polynomial_space</tt> at
83  // <tt>unit_point</tt>
84  const unsigned int n_basis = polynomial_space.n();
85  const unsigned int my_degree = this->degree();
86  std::vector<double> unit_point_values((values.size() == 0) ? 0 : n_basis);
87  std::vector<Tensor<1, dim>> unit_point_grads((grads.size() == 0) ? 0 :
88  n_basis);
89  std::vector<Tensor<2, dim>> unit_point_grad_grads(
90  (grad_grads.size() == 0) ? 0 : n_basis);
91  std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
92  std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
93 
94  switch (dim)
95  {
96  case 1:
97  {
98  polynomial_space.evaluate(unit_point,
99  unit_point_values,
100  unit_point_grads,
101  unit_point_grad_grads,
102  empty_vector_of_3rd_order_tensors,
103  empty_vector_of_4th_order_tensors);
104 
105  // Assign the correct values to the
106  // corresponding shape functions.
107  if (values.size() > 0)
108  for (unsigned int i = 0; i < unit_point_values.size(); ++i)
109  values[i][0] = unit_point_values[i];
110 
111  if (grads.size() > 0)
112  for (unsigned int i = 0; i < unit_point_grads.size(); ++i)
113  grads[i][0][0] = unit_point_grads[i][0];
114 
115  if (grad_grads.size() > 0)
116  for (unsigned int i = 0; i < unit_point_grad_grads.size(); ++i)
117  grad_grads[i][0][0][0] = unit_point_grad_grads[i][0][0];
118 
119  break;
120  }
121 
122  case 2:
123  {
124  polynomial_space.evaluate(unit_point,
125  unit_point_values,
126  unit_point_grads,
127  unit_point_grad_grads,
128  empty_vector_of_3rd_order_tensors,
129  empty_vector_of_4th_order_tensors);
130 
131  // Declare the values, derivatives and
132  // second derivatives vectors of
133  // <tt>polynomial_space</tt> at
134  // <tt>unit_point</tt> with coordinates
135  // shifted one step in positive direction
136  Point<dim> p;
137 
138  p(0) = unit_point(1);
139  p(1) = unit_point(0);
140 
141  std::vector<double> p_values((values.size() == 0) ? 0 : n_basis);
142  std::vector<Tensor<1, dim>> p_grads((grads.size() == 0) ? 0 :
143  n_basis);
144  std::vector<Tensor<2, dim>> p_grad_grads(
145  (grad_grads.size() == 0) ? 0 : n_basis);
146 
147  polynomial_space.evaluate(p,
148  p_values,
149  p_grads,
150  p_grad_grads,
151  empty_vector_of_3rd_order_tensors,
152  empty_vector_of_4th_order_tensors);
153 
154  // Assign the correct values to the
155  // corresponding shape functions.
156  if (values.size() > 0)
157  {
158  for (unsigned int i = 0; i <= my_degree; ++i)
159  for (unsigned int j = 0; j < 2; ++j)
160  {
161  values[i + j * (my_degree + 1)][0] = 0.0;
162  values[i + j * (my_degree + 1)][1] =
163  p_values[i + j * (my_degree + 1)];
164  values[i + (j + 2) * (my_degree + 1)][0] =
165  unit_point_values[i + j * (my_degree + 1)];
166  values[i + (j + 2) * (my_degree + 1)][1] = 0.0;
167  }
168 
169  if (my_degree > 0)
170  for (unsigned int i = 0; i <= my_degree; ++i)
171  for (unsigned int j = 0; j < my_degree; ++j)
172  {
173  values[(i + GeometryInfo<dim>::lines_per_cell) *
174  my_degree +
175  j + GeometryInfo<dim>::lines_per_cell][0] =
176  unit_point_values[i + (j + 2) * (my_degree + 1)];
177  values[(i + GeometryInfo<dim>::lines_per_cell) *
178  my_degree +
179  j + GeometryInfo<dim>::lines_per_cell][1] = 0.0;
180  values[i + (j + my_degree +
181  GeometryInfo<dim>::lines_per_cell) *
182  (my_degree + 1)][0] = 0.0;
183  values[i + (j + my_degree +
184  GeometryInfo<dim>::lines_per_cell) *
185  (my_degree + 1)][1] =
186  p_values[i + (j + 2) * (my_degree + 1)];
187  }
188  }
189 
190  if (grads.size() > 0)
191  {
192  for (unsigned int i = 0; i <= my_degree; ++i)
193  for (unsigned int j = 0; j < 2; ++j)
194  {
195  for (unsigned int k = 0; k < dim; ++k)
196  {
197  grads[i + j * (my_degree + 1)][0][k] = 0.0;
198  grads[i + (j + 2) * (my_degree + 1)][0][k] =
199  unit_point_grads[i + j * (my_degree + 1)][k];
200  grads[i + (j + 2) * (my_degree + 1)][1][k] = 0.0;
201  }
202 
203  grads[i + j * (my_degree + 1)][1][0] =
204  p_grads[i + j * (my_degree + 1)][1];
205  grads[i + j * (my_degree + 1)][1][1] =
206  p_grads[i + j * (my_degree + 1)][0];
207  }
208 
209  if (my_degree > 0)
210  for (unsigned int i = 0; i <= my_degree; ++i)
211  for (unsigned int j = 0; j < my_degree; ++j)
212  {
213  for (unsigned int k = 0; k < dim; ++k)
214  {
215  grads[(i + GeometryInfo<dim>::lines_per_cell) *
216  my_degree +
217  j + GeometryInfo<dim>::lines_per_cell][0][k] =
218  unit_point_grads[i + (j + 2) * (my_degree + 1)][k];
219  grads[(i + GeometryInfo<dim>::lines_per_cell) *
220  my_degree +
221  j + GeometryInfo<dim>::lines_per_cell][1][k] =
222  0.0;
223  grads[i + (j + my_degree +
224  GeometryInfo<dim>::lines_per_cell) *
225  (my_degree + 1)][0][k] = 0.0;
226  }
227 
228  grads[i + (j + my_degree +
229  GeometryInfo<dim>::lines_per_cell) *
230  (my_degree + 1)][1][0] =
231  p_grads[i + (j + 2) * (my_degree + 1)][1];
232  grads[i + (j + my_degree +
233  GeometryInfo<dim>::lines_per_cell) *
234  (my_degree + 1)][1][1] =
235  p_grads[i + (j + 2) * (my_degree + 1)][0];
236  }
237  }
238 
239  if (grad_grads.size() > 0)
240  {
241  for (unsigned int i = 0; i <= my_degree; ++i)
242  for (unsigned int j = 0; j < 2; ++j)
243  {
244  for (unsigned int k = 0; k < dim; ++k)
245  for (unsigned int l = 0; l < dim; ++l)
246  {
247  grad_grads[i + j * (my_degree + 1)][0][k][l] = 0.0;
248  grad_grads[i + (j + 2) * (my_degree + 1)][0][k][l] =
249  unit_point_grad_grads[i + j * (my_degree + 1)][k]
250  [l];
251  grad_grads[i + (j + 2) * (my_degree + 1)][1][k][l] =
252  0.0;
253  }
254 
255  grad_grads[i + j * (my_degree + 1)][1][0][0] =
256  p_grad_grads[i + j * (my_degree + 1)][1][1];
257  grad_grads[i + j * (my_degree + 1)][1][0][1] =
258  p_grad_grads[i + j * (my_degree + 1)][1][0];
259  grad_grads[i + j * (my_degree + 1)][1][1][0] =
260  p_grad_grads[i + j * (my_degree + 1)][0][1];
261  grad_grads[i + j * (my_degree + 1)][1][1][1] =
262  p_grad_grads[i + j * (my_degree + 1)][0][0];
263  }
264 
265  if (my_degree > 0)
266  for (unsigned int i = 0; i <= my_degree; ++i)
267  for (unsigned int j = 0; j < my_degree; ++j)
268  {
269  for (unsigned int k = 0; k < dim; ++k)
270  for (unsigned int l = 0; l < dim; ++l)
271  {
272  grad_grads[(i + GeometryInfo<dim>::lines_per_cell) *
273  my_degree +
274  j + GeometryInfo<dim>::lines_per_cell][0]
275  [k][l] = unit_point_grad_grads
276  [i + (j + 2) * (my_degree + 1)][k][l];
277  grad_grads[(i + GeometryInfo<dim>::lines_per_cell) *
278  my_degree +
279  j + GeometryInfo<dim>::lines_per_cell][1]
280  [k][l] = 0.0;
281  grad_grads[i + (j + my_degree +
282  GeometryInfo<dim>::lines_per_cell) *
283  (my_degree + 1)][0][k][l] = 0.0;
284  }
285 
286  grad_grads[i + (j + my_degree +
287  GeometryInfo<dim>::lines_per_cell) *
288  (my_degree + 1)][1][0][0] =
289  p_grad_grads[i + (j + 2) * (my_degree + 1)][1][1];
290  grad_grads[i + (j + my_degree +
291  GeometryInfo<dim>::lines_per_cell) *
292  (my_degree + 1)][1][0][1] =
293  p_grad_grads[i + (j + 2) * (my_degree + 1)][1][0];
294  grad_grads[i + (j + my_degree +
295  GeometryInfo<dim>::lines_per_cell) *
296  (my_degree + 1)][1][1][0] =
297  p_grad_grads[i + (j + 2) * (my_degree + 1)][0][1];
298  grad_grads[i + (j + my_degree +
299  GeometryInfo<dim>::lines_per_cell) *
300  (my_degree + 1)][1][1][1] =
301  p_grad_grads[i + (j + 2) * (my_degree + 1)][0][0];
302  }
303  }
304 
305  break;
306  }
307 
308  case 3:
309  {
310  polynomial_space.evaluate(unit_point,
311  unit_point_values,
312  unit_point_grads,
313  unit_point_grad_grads,
314  empty_vector_of_3rd_order_tensors,
315  empty_vector_of_4th_order_tensors);
316 
317  // Declare the values, derivatives
318  // and second derivatives vectors of
319  // <tt>polynomial_space</tt> at
320  // <tt>unit_point</tt> with coordinates
321  // shifted two steps in positive
322  // direction
323  Point<dim> p1, p2;
324  std::vector<double> p1_values((values.size() == 0) ? 0 : n_basis);
325  std::vector<Tensor<1, dim>> p1_grads((grads.size() == 0) ? 0 :
326  n_basis);
327  std::vector<Tensor<2, dim>> p1_grad_grads(
328  (grad_grads.size() == 0) ? 0 : n_basis);
329  std::vector<double> p2_values((values.size() == 0) ? 0 : n_basis);
330  std::vector<Tensor<1, dim>> p2_grads((grads.size() == 0) ? 0 :
331  n_basis);
332  std::vector<Tensor<2, dim>> p2_grad_grads(
333  (grad_grads.size() == 0) ? 0 : n_basis);
334 
335  p1(0) = unit_point(1);
336  p1(1) = unit_point(2);
337  p1(2) = unit_point(0);
338  polynomial_space.evaluate(p1,
339  p1_values,
340  p1_grads,
341  p1_grad_grads,
342  empty_vector_of_3rd_order_tensors,
343  empty_vector_of_4th_order_tensors);
344  p2(0) = unit_point(2);
345  p2(1) = unit_point(0);
346  p2(2) = unit_point(1);
347  polynomial_space.evaluate(p2,
348  p2_values,
349  p2_grads,
350  p2_grad_grads,
351  empty_vector_of_3rd_order_tensors,
352  empty_vector_of_4th_order_tensors);
353 
354  // Assign the correct values to the
355  // corresponding shape functions.
356  if (values.size() > 0)
357  {
358  for (unsigned int i = 0; i <= my_degree; ++i)
359  {
360  for (unsigned int j = 0; j < 2; ++j)
361  {
362  for (unsigned int k = 0; k < 2; ++k)
363  {
364  for (unsigned int l = 0; l < 2; ++l)
365  {
366  values[i + (j + 4 * k) * (my_degree + 1)][2 * l] =
367  0.0;
368  values[i + (j + 4 * k + 2) * (my_degree + 1)]
369  [l + 1] = 0.0;
370  values[i + (j + 2 * (k + 4)) * (my_degree + 1)]
371  [l] = 0.0;
372  }
373 
374  values[i + (j + 4 * k + 2) * (my_degree + 1)][0] =
375  unit_point_values[i + (j + k * (my_degree + 2)) *
376  (my_degree + 1)];
377  values[i + (j + 2 * (k + 4)) * (my_degree + 1)][2] =
378  p2_values[i + (j + k * (my_degree + 2)) *
379  (my_degree + 1)];
380  }
381 
382  values[i + j * (my_degree + 1)][1] =
383  p1_values[i + j * (my_degree + 1) * (my_degree + 2)];
384  }
385 
386  values[i + 4 * (my_degree + 1)][1] =
387  p1_values[i + my_degree + 1];
388  values[i + 5 * (my_degree + 1)][1] =
389  p1_values[i + (my_degree + 1) * (my_degree + 3)];
390  }
391 
392  if (my_degree > 0)
393  for (unsigned int i = 0; i <= my_degree; ++i)
394  for (unsigned int j = 0; j < my_degree; ++j)
395  {
396  for (unsigned int k = 0; k < my_degree; ++k)
397  {
398  for (unsigned int l = 0; l < 2; ++l)
399  {
400  values[((i +
401  2 * GeometryInfo<dim>::faces_per_cell) *
402  my_degree +
403  j + GeometryInfo<dim>::lines_per_cell +
404  2 * GeometryInfo<dim>::faces_per_cell) *
405  my_degree +
406  k + GeometryInfo<dim>::lines_per_cell]
407  [l + 1] = 0.0;
408  values[(i +
409  (j +
410  2 * GeometryInfo<dim>::faces_per_cell +
411  my_degree) *
412  (my_degree + 1) +
413  GeometryInfo<dim>::lines_per_cell) *
414  my_degree +
415  k + GeometryInfo<dim>::lines_per_cell]
416  [2 * l] = 0.0;
417  values[i +
418  (j +
419  (k +
420  2 * (GeometryInfo<dim>::faces_per_cell +
421  my_degree)) *
422  my_degree +
423  GeometryInfo<dim>::lines_per_cell) *
424  (my_degree + 1)][l] = 0.0;
425  }
426 
427  values[((i + 2 * GeometryInfo<dim>::faces_per_cell) *
428  my_degree +
429  j + GeometryInfo<dim>::lines_per_cell +
430  2 * GeometryInfo<dim>::faces_per_cell) *
431  my_degree +
432  k + GeometryInfo<dim>::lines_per_cell][0] =
433  unit_point_values[i +
434  (j + (k + 2) * (my_degree + 2) +
435  2) *
436  (my_degree + 1)];
437  values[(i +
438  (j + 2 * GeometryInfo<dim>::faces_per_cell +
439  my_degree) *
440  (my_degree + 1) +
441  GeometryInfo<dim>::lines_per_cell) *
442  my_degree +
443  k + GeometryInfo<dim>::lines_per_cell][1] =
444  p1_values[i + ((j + 2) * (my_degree + 2) + k + 2) *
445  (my_degree + 1)];
446  values[i +
447  (j +
448  (k + 2 * (GeometryInfo<dim>::faces_per_cell +
449  my_degree)) *
450  my_degree +
451  GeometryInfo<dim>::lines_per_cell) *
452  (my_degree + 1)][2] =
453  p2_values[i + (j + (k + 2) * (my_degree + 2) + 2) *
454  (my_degree + 1)];
455  }
456 
457  for (unsigned int k = 0; k < 2; ++k)
458  {
459  for (unsigned int l = 0; l < 2; ++l)
460  {
461  for (unsigned int m = 0; m < 2; ++m)
462  {
463  values[i +
464  (j +
465  (2 * (k + 2 * l) + 1) * my_degree +
466  GeometryInfo<dim>::lines_per_cell) *
467  (my_degree + 1)][m + l] = 0.0;
468  values[(i +
469  2 * (k + 2 * (l + 1)) *
470  (my_degree + 1) +
471  GeometryInfo<dim>::lines_per_cell) *
472  my_degree +
473  j + GeometryInfo<dim>::lines_per_cell]
474  [m + l] = 0.0;
475  }
476 
477  values[(i + 2 * k * (my_degree + 1) +
478  GeometryInfo<dim>::lines_per_cell) *
479  my_degree +
480  j + GeometryInfo<dim>::lines_per_cell]
481  [2 * l] = 0.0;
482  values[i + (j + (2 * k + 9) * my_degree +
483  GeometryInfo<dim>::lines_per_cell) *
484  (my_degree + 1)][2 * l] = 0.0;
485  }
486 
487  values[(i + 2 * k * (my_degree + 1) +
488  GeometryInfo<dim>::lines_per_cell) *
489  my_degree +
490  j + GeometryInfo<dim>::lines_per_cell][1] =
491  p1_values[i + (j + k * (my_degree + 2) + 2) *
492  (my_degree + 1)];
493  values[i + (j + (2 * k + 1) * my_degree +
494  GeometryInfo<dim>::lines_per_cell) *
495  (my_degree + 1)][2] =
496  p2_values[i + ((j + 2) * (my_degree + 2) + k) *
497  (my_degree + 1)];
498  values[(i + 2 * (k + 2) * (my_degree + 1) +
499  GeometryInfo<dim>::lines_per_cell) *
500  my_degree +
501  j + GeometryInfo<dim>::lines_per_cell][2] =
502  p2_values[i + (j + k * (my_degree + 2) + 2) *
503  (my_degree + 1)];
504  values[i + (j + (2 * k + 5) * my_degree +
505  GeometryInfo<dim>::lines_per_cell) *
506  (my_degree + 1)][0] =
507  unit_point_values[i +
508  ((j + 2) * (my_degree + 2) + k) *
509  (my_degree + 1)];
510  values[(i + 2 * (k + 4) * (my_degree + 1) +
511  GeometryInfo<dim>::lines_per_cell) *
512  my_degree +
513  j + GeometryInfo<dim>::lines_per_cell][0] =
514  unit_point_values[i +
515  (j + k * (my_degree + 2) + 2) *
516  (my_degree + 1)];
517  values[i + (j + (2 * k + 9) * my_degree +
518  GeometryInfo<dim>::lines_per_cell) *
519  (my_degree + 1)][1] =
520  p1_values[i + ((j + 2) * (my_degree + 2) + k) *
521  (my_degree + 1)];
522  }
523  }
524  }
525 
526  if (grads.size() > 0)
527  {
528  for (unsigned int i = 0; i <= my_degree; ++i)
529  {
530  for (unsigned int j = 0; j < 2; ++j)
531  {
532  for (unsigned int k = 0; k < 2; ++k)
533  {
534  for (unsigned int l = 0; l < 2; ++l)
535  for (unsigned int m = 0; m < dim; ++m)
536  {
537  grads[i + (j + 4 * k) * (my_degree + 1)][2 * l]
538  [m] = 0.0;
539  grads[i + (j + 4 * k + 2) * (my_degree + 1)]
540  [l + 1][m] = 0.0;
541  grads[i + (j + 2 * (k + 4)) * (my_degree + 1)]
542  [l][m] = 0.0;
543  }
544 
545  for (unsigned int l = 0; l < dim; ++l)
546  grads[i + (j + 4 * k + 2) * (my_degree + 1)][0][l] =
547  unit_point_grads[i + (j + k * (my_degree + 2)) *
548  (my_degree + 1)][l];
549 
550  grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][0] =
551  p2_grads[i + (j + k * (my_degree + 2)) *
552  (my_degree + 1)][1];
553  grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][1] =
554  p2_grads[i + (j + k * (my_degree + 2)) *
555  (my_degree + 1)][2];
556  grads[i + (j + 2 * (k + 4)) * (my_degree + 1)][2][2] =
557  p2_grads[i + (j + k * (my_degree + 2)) *
558  (my_degree + 1)][0];
559  }
560 
561  grads[i + j * (my_degree + 1)][1][0] =
562  p1_grads[i + j * (my_degree + 1) * (my_degree + 2)][2];
563  grads[i + j * (my_degree + 1)][1][1] =
564  p1_grads[i + j * (my_degree + 1) * (my_degree + 2)][0];
565  grads[i + j * (my_degree + 1)][1][2] =
566  p1_grads[i + j * (my_degree + 1) * (my_degree + 2)][1];
567  }
568 
569  grads[i + 4 * (my_degree + 1)][1][0] =
570  p1_grads[i + my_degree + 1][2];
571  grads[i + 4 * (my_degree + 1)][1][1] =
572  p1_grads[i + my_degree + 1][0];
573  grads[i + 4 * (my_degree + 1)][1][2] =
574  p1_grads[i + my_degree + 1][1];
575  grads[i + 5 * (my_degree + 1)][1][0] =
576  p1_grads[i + (my_degree + 1) * (my_degree + 3)][2];
577  grads[i + 5 * (my_degree + 1)][1][1] =
578  p1_grads[i + (my_degree + 1) * (my_degree + 3)][0];
579  grads[i + 5 * (my_degree + 1)][1][2] =
580  p1_grads[i + (my_degree + 1) * (my_degree + 3)][1];
581  }
582 
583  if (my_degree > 0)
584  for (unsigned int i = 0; i <= my_degree; ++i)
585  for (unsigned int j = 0; j < my_degree; ++j)
586  {
587  for (unsigned int k = 0; k < my_degree; ++k)
588  {
589  for (unsigned int l = 0; l < dim; ++l)
590  {
591  for (unsigned int m = 0; m < 2; ++m)
592  {
593  grads
594  [((i +
595  2 * GeometryInfo<dim>::faces_per_cell) *
596  my_degree +
597  j + GeometryInfo<dim>::lines_per_cell +
598  2 * GeometryInfo<dim>::faces_per_cell) *
599  my_degree +
600  k + GeometryInfo<dim>::lines_per_cell]
601  [m + 1][l] = 0.0;
602  grads[(i +
603  (j +
604  2 *
605  GeometryInfo<dim>::faces_per_cell +
606  my_degree) *
607  (my_degree + 1) +
608  GeometryInfo<dim>::lines_per_cell) *
609  my_degree +
610  k + GeometryInfo<dim>::lines_per_cell]
611  [2 * m][l] = 0.0;
612  grads[i +
613  (j +
614  (k +
615  2 *
616  (GeometryInfo<dim>::faces_per_cell +
617  my_degree)) *
618  my_degree +
619  GeometryInfo<dim>::lines_per_cell) *
620  (my_degree + 1)][m][l] = 0.0;
621  }
622 
623  grads[((i +
624  2 * GeometryInfo<dim>::faces_per_cell) *
625  my_degree +
626  j + GeometryInfo<dim>::lines_per_cell +
627  2 * GeometryInfo<dim>::faces_per_cell) *
628  my_degree +
629  k + GeometryInfo<dim>::lines_per_cell][0]
630  [l] = unit_point_grads
631  [i + (j + (k + 2) * (my_degree + 2) + 2) *
632  (my_degree + 1)][l];
633  }
634 
635  grads[(i +
636  (j + 2 * GeometryInfo<dim>::faces_per_cell +
637  my_degree) *
638  (my_degree + 1) +
639  GeometryInfo<dim>::lines_per_cell) *
640  my_degree +
641  k + GeometryInfo<dim>::lines_per_cell][1][0] =
642  p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) *
643  (my_degree + 1)][2];
644  grads[(i +
645  (j + 2 * GeometryInfo<dim>::faces_per_cell +
646  my_degree) *
647  (my_degree + 1) +
648  GeometryInfo<dim>::lines_per_cell) *
649  my_degree +
650  k + GeometryInfo<dim>::lines_per_cell][1][1] =
651  p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) *
652  (my_degree + 1)][0];
653  grads[(i +
654  (j + 2 * GeometryInfo<dim>::faces_per_cell +
655  my_degree) *
656  (my_degree + 1) +
657  GeometryInfo<dim>::lines_per_cell) *
658  my_degree +
659  k + GeometryInfo<dim>::lines_per_cell][1][2] =
660  p1_grads[i + ((j + 2) * (my_degree + 2) + k + 2) *
661  (my_degree + 1)][1];
662  grads[i +
663  (j +
664  (k + 2 * (GeometryInfo<dim>::faces_per_cell +
665  my_degree)) *
666  my_degree +
667  GeometryInfo<dim>::lines_per_cell) *
668  (my_degree + 1)][2][0] =
669  p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) *
670  (my_degree + 1)][1];
671  grads[i +
672  (j +
673  (k + 2 * (GeometryInfo<dim>::faces_per_cell +
674  my_degree)) *
675  my_degree +
676  GeometryInfo<dim>::lines_per_cell) *
677  (my_degree + 1)][2][1] =
678  p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) *
679  (my_degree + 1)][2];
680  grads[i +
681  (j +
682  (k + 2 * (GeometryInfo<dim>::faces_per_cell +
683  my_degree)) *
684  my_degree +
685  GeometryInfo<dim>::lines_per_cell) *
686  (my_degree + 1)][2][2] =
687  p2_grads[i + (j + (k + 2) * (my_degree + 2) + 2) *
688  (my_degree + 1)][0];
689  }
690 
691  for (unsigned int k = 0; k < 2; ++k)
692  {
693  for (unsigned int l = 0; l < 2; ++l)
694  for (unsigned int m = 0; m < dim; ++m)
695  {
696  for (unsigned int n = 0; n < 2; ++n)
697  {
698  grads[i +
699  (j +
700  (2 * (k + 2 * l) + 1) * my_degree +
701  GeometryInfo<dim>::lines_per_cell) *
702  (my_degree + 1)][n + l][m] = 0.0;
703  grads[(i +
704  2 * (k + 2 * (l + 1)) *
705  (my_degree + 1) +
706  GeometryInfo<dim>::lines_per_cell) *
707  my_degree +
708  j + GeometryInfo<dim>::lines_per_cell]
709  [n + l][m] = 0.0;
710  }
711 
712  grads[(i + 2 * k * (my_degree + 1) +
713  GeometryInfo<dim>::lines_per_cell) *
714  my_degree +
715  j + GeometryInfo<dim>::lines_per_cell]
716  [2 * l][m] = 0.0;
717  grads[i + (j + (2 * k + 9) * my_degree +
718  GeometryInfo<dim>::lines_per_cell) *
719  (my_degree + 1)][2 * l][m] = 0.0;
720  }
721 
722  for (unsigned int l = 0; l < dim; ++l)
723  {
724  grads[i + (j + (2 * k + 5) * my_degree +
725  GeometryInfo<dim>::lines_per_cell) *
726  (my_degree + 1)][0][l] =
727  unit_point_grads[i +
728  ((j + 2) * (my_degree + 2) +
729  k) *
730  (my_degree + 1)][l];
731  grads[(i + 2 * (k + 4) * (my_degree + 1) +
732  GeometryInfo<dim>::lines_per_cell) *
733  my_degree +
734  j +
735  GeometryInfo<dim>::lines_per_cell][0][l] =
736  unit_point_grads[i +
737  (j + k * (my_degree + 2) + 2) *
738  (my_degree + 1)][l];
739  }
740 
741  grads[(i + 2 * k * (my_degree + 1) +
742  GeometryInfo<dim>::lines_per_cell) *
743  my_degree +
744  j + GeometryInfo<dim>::lines_per_cell][1][0] =
745  p1_grads[i + (j + k * (my_degree + 2) + 2) *
746  (my_degree + 1)][2];
747  grads[(i + 2 * k * (my_degree + 1) +
748  GeometryInfo<dim>::lines_per_cell) *
749  my_degree +
750  j + GeometryInfo<dim>::lines_per_cell][1][1] =
751  p1_grads[i + (j + k * (my_degree + 2) + 2) *
752  (my_degree + 1)][0];
753  grads[(i + 2 * k * (my_degree + 1) +
754  GeometryInfo<dim>::lines_per_cell) *
755  my_degree +
756  j + GeometryInfo<dim>::lines_per_cell][1][2] =
757  p1_grads[i + (j + k * (my_degree + 2) + 2) *
758  (my_degree + 1)][1];
759  grads[i + (j + (2 * k + 1) * my_degree +
760  GeometryInfo<dim>::lines_per_cell) *
761  (my_degree + 1)][2][0] =
762  p2_grads[i + ((j + 2) * (my_degree + 2) + k) *
763  (my_degree + 1)][1];
764  grads[i + (j + (2 * k + 1) * my_degree +
765  GeometryInfo<dim>::lines_per_cell) *
766  (my_degree + 1)][2][1] =
767  p2_grads[i + ((j + 2) * (my_degree + 2) + k) *
768  (my_degree + 1)][2];
769  grads[i + (j + (2 * k + 1) * my_degree +
770  GeometryInfo<dim>::lines_per_cell) *
771  (my_degree + 1)][2][2] =
772  p2_grads[i + ((j + 2) * (my_degree + 2) + k) *
773  (my_degree + 1)][0];
774  grads[(i + 2 * (k + 2) * (my_degree + 1) +
775  GeometryInfo<dim>::lines_per_cell) *
776  my_degree +
777  j + GeometryInfo<dim>::lines_per_cell][2][0] =
778  p2_grads[i + (j + k * (my_degree + 2) + 2) *
779  (my_degree + 1)][1];
780  grads[(i + 2 * (k + 2) * (my_degree + 1) +
781  GeometryInfo<dim>::lines_per_cell) *
782  my_degree +
783  j + GeometryInfo<dim>::lines_per_cell][2][1] =
784  p2_grads[i + (j + k * (my_degree + 2) + 2) *
785  (my_degree + 1)][2];
786  grads[(i + 2 * (k + 2) * (my_degree + 1) +
787  GeometryInfo<dim>::lines_per_cell) *
788  my_degree +
789  j + GeometryInfo<dim>::lines_per_cell][2][2] =
790  p2_grads[i + (j + k * (my_degree + 2) + 2) *
791  (my_degree + 1)][0];
792  grads[i + (j + (2 * k + 9) * my_degree +
793  GeometryInfo<dim>::lines_per_cell) *
794  (my_degree + 1)][1][0] =
795  p1_grads[i + ((j + 2) * (my_degree + 2) + k) *
796  (my_degree + 1)][2];
797  grads[i + (j + (2 * k + 9) * my_degree +
798  GeometryInfo<dim>::lines_per_cell) *
799  (my_degree + 1)][1][1] =
800  p1_grads[i + ((j + 2) * (my_degree + 2) + k) *
801  (my_degree + 1)][0];
802  grads[i + (j + (2 * k + 9) * my_degree +
803  GeometryInfo<dim>::lines_per_cell) *
804  (my_degree + 1)][1][2] =
805  p1_grads[i + ((j + 2) * (my_degree + 2) + k) *
806  (my_degree + 1)][1];
807  }
808  }
809  }
810 
811  if (grad_grads.size() > 0)
812  {
813  for (unsigned int i = 0; i <= my_degree; ++i)
814  {
815  for (unsigned int j = 0; j < 2; ++j)
816  {
817  for (unsigned int k = 0; k < 2; ++k)
818  {
819  for (unsigned int l = 0; l < dim; ++l)
820  for (unsigned int m = 0; m < dim; ++m)
821  {
822  for (unsigned int n = 0; n < 2; ++n)
823  {
824  grad_grads[i +
825  (j + 4 * k) * (my_degree + 1)]
826  [2 * n][l][m] = 0.0;
827  grad_grads[i + (j + 4 * k + 2) *
828  (my_degree + 1)][n + 1][l]
829  [m] = 0.0;
830  grad_grads[i + (j + 2 * (k + 4)) *
831  (my_degree + 1)][n][l][m] =
832  0.0;
833  }
834 
835  grad_grads[i + (j + 4 * k + 2) *
836  (my_degree + 1)][0][l][m] =
837  unit_point_grad_grads
838  [i + (j + k * (my_degree + 2)) *
839  (my_degree + 1)][l][m];
840  }
841 
842  grad_grads[i + (j + 2 * (k + 4)) *
843  (my_degree + 1)][2][0][0] =
844  p2_grad_grads[i + (j + k * (my_degree + 2)) *
845  (my_degree + 1)][1][1];
846  grad_grads[i + (j + 2 * (k + 4)) *
847  (my_degree + 1)][2][0][1] =
848  p2_grad_grads[i + (j + k * (my_degree + 2)) *
849  (my_degree + 1)][1][2];
850  grad_grads[i + (j + 2 * (k + 4)) *
851  (my_degree + 1)][2][0][2] =
852  p2_grad_grads[i + (j + k * (my_degree + 2)) *
853  (my_degree + 1)][1][0];
854  grad_grads[i + (j + 2 * (k + 4)) *
855  (my_degree + 1)][2][1][0] =
856  p2_grad_grads[i + (j + k * (my_degree + 2)) *
857  (my_degree + 1)][2][1];
858  grad_grads[i + (j + 2 * (k + 4)) *
859  (my_degree + 1)][2][1][1] =
860  p2_grad_grads[i + (j + k * (my_degree + 2)) *
861  (my_degree + 1)][2][2];
862  grad_grads[i + (j + 2 * (k + 4)) *
863  (my_degree + 1)][2][1][2] =
864  p2_grad_grads[i + (j + k * (my_degree + 2)) *
865  (my_degree + 1)][2][0];
866  grad_grads[i + (j + 2 * (k + 4)) *
867  (my_degree + 1)][2][2][0] =
868  p2_grad_grads[i + (j + k * (my_degree + 2)) *
869  (my_degree + 1)][0][1];
870  grad_grads[i + (j + 2 * (k + 4)) *
871  (my_degree + 1)][2][2][1] =
872  p2_grad_grads[i + (j + k * (my_degree + 2)) *
873  (my_degree + 1)][0][2];
874  grad_grads[i + (j + 2 * (k + 4)) *
875  (my_degree + 1)][2][2][2] =
876  p2_grad_grads[i + (j + k * (my_degree + 2)) *
877  (my_degree + 1)][0][0];
878  }
879 
880  grad_grads[i + j * (my_degree + 1)][1][0][0] =
881  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
882  [2][2];
883  grad_grads[i + j * (my_degree + 1)][1][0][1] =
884  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
885  [2][0];
886  grad_grads[i + j * (my_degree + 1)][1][0][2] =
887  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
888  [2][1];
889  grad_grads[i + j * (my_degree + 1)][1][1][0] =
890  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
891  [0][2];
892  grad_grads[i + j * (my_degree + 1)][1][1][1] =
893  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
894  [0][0];
895  grad_grads[i + j * (my_degree + 1)][1][1][2] =
896  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
897  [0][1];
898  grad_grads[i + j * (my_degree + 1)][1][2][0] =
899  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
900  [1][2];
901  grad_grads[i + j * (my_degree + 1)][1][2][1] =
902  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
903  [1][0];
904  grad_grads[i + j * (my_degree + 1)][1][2][2] =
905  p1_grad_grads[i + j * (my_degree + 1) * (my_degree + 2)]
906  [1][1];
907  }
908 
909  grad_grads[i + 4 * (my_degree + 1)][1][0][0] =
910  p1_grad_grads[i + my_degree + 1][2][2];
911  grad_grads[i + 4 * (my_degree + 1)][1][0][1] =
912  p1_grad_grads[i + my_degree + 1][2][0];
913  grad_grads[i + 4 * (my_degree + 1)][1][0][2] =
914  p1_grad_grads[i + my_degree + 1][2][1];
915  grad_grads[i + 4 * (my_degree + 1)][1][1][0] =
916  p1_grad_grads[i + my_degree + 1][0][2];
917  grad_grads[i + 4 * (my_degree + 1)][1][1][1] =
918  p1_grad_grads[i + my_degree + 1][0][0];
919  grad_grads[i + 4 * (my_degree + 1)][1][1][2] =
920  p1_grad_grads[i + my_degree + 1][0][1];
921  grad_grads[i + 4 * (my_degree + 1)][1][2][0] =
922  p1_grad_grads[i + my_degree + 1][1][2];
923  grad_grads[i + 4 * (my_degree + 1)][1][2][1] =
924  p1_grad_grads[i + my_degree + 1][1][0];
925  grad_grads[i + 4 * (my_degree + 1)][1][2][2] =
926  p1_grad_grads[i + my_degree + 1][1][1];
927  grad_grads[i + 5 * (my_degree + 1)][1][0][0] =
928  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2][2];
929  grad_grads[i + 5 * (my_degree + 1)][1][0][1] =
930  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2][0];
931  grad_grads[i + 5 * (my_degree + 1)][1][0][2] =
932  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][2][1];
933  grad_grads[i + 5 * (my_degree + 1)][1][1][0] =
934  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0][2];
935  grad_grads[i + 5 * (my_degree + 1)][1][1][1] =
936  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0][0];
937  grad_grads[i + 5 * (my_degree + 1)][1][1][2] =
938  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][0][1];
939  grad_grads[i + 5 * (my_degree + 1)][1][2][0] =
940  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1][2];
941  grad_grads[i + 5 * (my_degree + 1)][1][2][1] =
942  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1][0];
943  grad_grads[i + 5 * (my_degree + 1)][1][2][2] =
944  p1_grad_grads[i + (my_degree + 1) * (my_degree + 3)][1][1];
945  }
946 
947  if (my_degree > 0)
948  for (unsigned int i = 0; i <= my_degree; ++i)
949  for (unsigned int j = 0; j < my_degree; ++j)
950  {
951  for (unsigned int k = 0; k < my_degree; ++k)
952  {
953  for (unsigned int l = 0; l < dim; ++l)
954  for (unsigned int m = 0; m < dim; ++m)
955  {
956  for (unsigned int n = 0; n < 2; ++n)
957  {
958  grad_grads
959  [((i +
960  2 *
961  GeometryInfo<dim>::faces_per_cell) *
962  my_degree +
963  j + GeometryInfo<dim>::lines_per_cell +
964  2 * GeometryInfo<dim>::faces_per_cell) *
965  my_degree +
966  k + GeometryInfo<dim>::lines_per_cell]
967  [n + 1][l][m] = 0.0;
968  grad_grads
969  [(i +
970  (j +
971  2 * GeometryInfo<dim>::faces_per_cell +
972  my_degree) *
973  (my_degree + 1) +
974  GeometryInfo<dim>::lines_per_cell) *
975  my_degree +
976  k + GeometryInfo<dim>::lines_per_cell]
977  [2 * n][l][m] = 0.0;
978  grad_grads
979  [i + (j +
980  (k + 2 * (GeometryInfo<
981  dim>::faces_per_cell +
982  my_degree)) *
983  my_degree +
984  GeometryInfo<dim>::lines_per_cell) *
985  (my_degree + 1)][n][l][m] = 0.0;
986  }
987 
988  grad_grads
989  [((i +
990  2 * GeometryInfo<dim>::faces_per_cell) *
991  my_degree +
992  j + GeometryInfo<dim>::lines_per_cell +
993  2 * GeometryInfo<dim>::faces_per_cell) *
994  my_degree +
995  k + GeometryInfo<dim>::lines_per_cell][0][l]
996  [m] = unit_point_grad_grads
997  [i + (j + (k + 2) * (my_degree + 2) + 2) *
998  (my_degree + 1)][l][m];
999  }
1000 
1001  grad_grads
1002  [(i +
1003  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1004  my_degree) *
1005  (my_degree + 1) +
1006  GeometryInfo<dim>::lines_per_cell) *
1007  my_degree +
1008  k + GeometryInfo<dim>::lines_per_cell][1][0][0] =
1009  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1010  2) *
1011  (my_degree + 1)][2][2];
1012  grad_grads
1013  [(i +
1014  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1015  my_degree) *
1016  (my_degree + 1) +
1017  GeometryInfo<dim>::lines_per_cell) *
1018  my_degree +
1019  k + GeometryInfo<dim>::lines_per_cell][1][0][1] =
1020  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1021  2) *
1022  (my_degree + 1)][2][0];
1023  grad_grads
1024  [(i +
1025  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1026  my_degree) *
1027  (my_degree + 1) +
1028  GeometryInfo<dim>::lines_per_cell) *
1029  my_degree +
1030  k + GeometryInfo<dim>::lines_per_cell][1][0][2] =
1031  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1032  2) *
1033  (my_degree + 1)][2][1];
1034  grad_grads
1035  [(i +
1036  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1037  my_degree) *
1038  (my_degree + 1) +
1039  GeometryInfo<dim>::lines_per_cell) *
1040  my_degree +
1041  k + GeometryInfo<dim>::lines_per_cell][1][1][0] =
1042  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1043  2) *
1044  (my_degree + 1)][0][2];
1045  grad_grads
1046  [(i +
1047  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1048  my_degree) *
1049  (my_degree + 1) +
1050  GeometryInfo<dim>::lines_per_cell) *
1051  my_degree +
1052  k + GeometryInfo<dim>::lines_per_cell][1][1][1] =
1053  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1054  2) *
1055  (my_degree + 1)][0][0];
1056  grad_grads
1057  [(i +
1058  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1059  my_degree) *
1060  (my_degree + 1) +
1061  GeometryInfo<dim>::lines_per_cell) *
1062  my_degree +
1063  k + GeometryInfo<dim>::lines_per_cell][1][1][2] =
1064  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1065  2) *
1066  (my_degree + 1)][0][1];
1067  grad_grads
1068  [(i +
1069  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1070  my_degree) *
1071  (my_degree + 1) +
1072  GeometryInfo<dim>::lines_per_cell) *
1073  my_degree +
1074  k + GeometryInfo<dim>::lines_per_cell][1][2][0] =
1075  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1076  2) *
1077  (my_degree + 1)][1][2];
1078  grad_grads
1079  [(i +
1080  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1081  my_degree) *
1082  (my_degree + 1) +
1083  GeometryInfo<dim>::lines_per_cell) *
1084  my_degree +
1085  k + GeometryInfo<dim>::lines_per_cell][1][2][1] =
1086  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1087  2) *
1088  (my_degree + 1)][1][0];
1089  grad_grads
1090  [(i +
1091  (j + 2 * GeometryInfo<dim>::faces_per_cell +
1092  my_degree) *
1093  (my_degree + 1) +
1094  GeometryInfo<dim>::lines_per_cell) *
1095  my_degree +
1096  k + GeometryInfo<dim>::lines_per_cell][1][2][2] =
1097  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k +
1098  2) *
1099  (my_degree + 1)][1][1];
1100  grad_grads[i +
1101  (j +
1102  (k +
1103  2 * (GeometryInfo<dim>::faces_per_cell +
1104  my_degree)) *
1105  my_degree +
1106  GeometryInfo<dim>::lines_per_cell) *
1107  (my_degree + 1)][2][0][0] =
1108  p2_grad_grads[i +
1109  (j + (k + 2) * (my_degree + 2) + 2) *
1110  (my_degree + 1)][1][1];
1111  grad_grads[i +
1112  (j +
1113  (k +
1114  2 * (GeometryInfo<dim>::faces_per_cell +
1115  my_degree)) *
1116  my_degree +
1117  GeometryInfo<dim>::lines_per_cell) *
1118  (my_degree + 1)][2][0][1] =
1119  p2_grad_grads[i +
1120  (j + (k + 2) * (my_degree + 2) + 2) *
1121  (my_degree + 1)][1][2];
1122  grad_grads[i +
1123  (j +
1124  (k +
1125  2 * (GeometryInfo<dim>::faces_per_cell +
1126  my_degree)) *
1127  my_degree +
1128  GeometryInfo<dim>::lines_per_cell) *
1129  (my_degree + 1)][2][0][2] =
1130  p2_grad_grads[i +
1131  (j + (k + 2) * (my_degree + 2) + 2) *
1132  (my_degree + 1)][1][0];
1133  grad_grads[i +
1134  (j +
1135  (k +
1136  2 * (GeometryInfo<dim>::faces_per_cell +
1137  my_degree)) *
1138  my_degree +
1139  GeometryInfo<dim>::lines_per_cell) *
1140  (my_degree + 1)][2][1][0] =
1141  p2_grad_grads[i +
1142  (j + (k + 2) * (my_degree + 2) + 2) *
1143  (my_degree + 1)][2][1];
1144  grad_grads[i +
1145  (j +
1146  (k +
1147  2 * (GeometryInfo<dim>::faces_per_cell +
1148  my_degree)) *
1149  my_degree +
1150  GeometryInfo<dim>::lines_per_cell) *
1151  (my_degree + 1)][2][1][1] =
1152  p2_grad_grads[i +
1153  (j + (k + 2) * (my_degree + 2) + 2) *
1154  (my_degree + 1)][2][2];
1155  grad_grads[i +
1156  (j +
1157  (k +
1158  2 * (GeometryInfo<dim>::faces_per_cell +
1159  my_degree)) *
1160  my_degree +
1161  GeometryInfo<dim>::lines_per_cell) *
1162  (my_degree + 1)][2][1][2] =
1163  p2_grad_grads[i +
1164  (j + (k + 2) * (my_degree + 2) + 2) *
1165  (my_degree + 1)][2][0];
1166  grad_grads[i +
1167  (j +
1168  (k +
1169  2 * (GeometryInfo<dim>::faces_per_cell +
1170  my_degree)) *
1171  my_degree +
1172  GeometryInfo<dim>::lines_per_cell) *
1173  (my_degree + 1)][2][2][0] =
1174  p2_grad_grads[i +
1175  (j + (k + 2) * (my_degree + 2) + 2) *
1176  (my_degree + 1)][0][1];
1177  grad_grads[i +
1178  (j +
1179  (k +
1180  2 * (GeometryInfo<dim>::faces_per_cell +
1181  my_degree)) *
1182  my_degree +
1183  GeometryInfo<dim>::lines_per_cell) *
1184  (my_degree + 1)][2][2][1] =
1185  p2_grad_grads[i +
1186  (j + (k + 2) * (my_degree + 2) + 2) *
1187  (my_degree + 1)][0][2];
1188  grad_grads[i +
1189  (j +
1190  (k +
1191  2 * (GeometryInfo<dim>::faces_per_cell +
1192  my_degree)) *
1193  my_degree +
1194  GeometryInfo<dim>::lines_per_cell) *
1195  (my_degree + 1)][2][2][2] =
1196  p2_grad_grads[i +
1197  (j + (k + 2) * (my_degree + 2) + 2) *
1198  (my_degree + 1)][0][0];
1199  }
1200 
1201  for (unsigned int k = 0; k < 2; ++k)
1202  {
1203  for (unsigned int l = 0; l < dim; ++l)
1204  for (unsigned int m = 0; m < dim; ++m)
1205  {
1206  for (unsigned int n = 0; n < 2; ++n)
1207  {
1208  for (unsigned int o = 0; o < 2; ++o)
1209  {
1210  grad_grads
1211  [i +
1212  (j +
1213  (2 * (k + 2 * n) + 1) * my_degree +
1214  GeometryInfo<dim>::lines_per_cell) *
1215  (my_degree + 1)][o + n][l][m] =
1216  0.0;
1217  grad_grads
1218  [(i +
1219  2 * (k + 2 * (n + 1)) *
1220  (my_degree + 1) +
1221  GeometryInfo<dim>::lines_per_cell) *
1222  my_degree +
1223  j +
1224  GeometryInfo<dim>::lines_per_cell]
1225  [o + k][l][m] = 0.0;
1226  }
1227 
1228  grad_grads
1229  [(i + 2 * k * (my_degree + 1) +
1230  GeometryInfo<dim>::lines_per_cell) *
1231  my_degree +
1232  j + GeometryInfo<dim>::lines_per_cell]
1233  [2 * n][l][m] = 0.0;
1234  grad_grads
1235  [i + (j + (2 * k + 9) * my_degree +
1236  GeometryInfo<dim>::lines_per_cell) *
1237  (my_degree + 1)][2 * n][l][m] =
1238  0.0;
1239  }
1240 
1241  grad_grads[i +
1242  (j + (2 * k + 5) * my_degree +
1243  GeometryInfo<dim>::lines_per_cell) *
1244  (my_degree + 1)][0][l][m] =
1245  unit_point_grad_grads
1246  [i + ((j + 2) * (my_degree + 2) + k) *
1247  (my_degree + 1)][l][m];
1248  grad_grads[(i + 2 * (k + 4) * (my_degree + 1) +
1249  GeometryInfo<dim>::lines_per_cell) *
1250  my_degree +
1251  j +
1252  GeometryInfo<dim>::lines_per_cell][0]
1253  [l][m] = unit_point_grad_grads
1254  [i + (j + k * (my_degree + 2) + 2) *
1255  (my_degree + 1)][l][m];
1256  }
1257 
1258  grad_grads
1259  [(i + 2 * k * (my_degree + 1) +
1260  GeometryInfo<dim>::lines_per_cell) *
1261  my_degree +
1262  j + GeometryInfo<dim>::lines_per_cell][1][0][0] =
1263  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1264  (my_degree + 1)][2][2];
1265  grad_grads
1266  [(i + 2 * k * (my_degree + 1) +
1267  GeometryInfo<dim>::lines_per_cell) *
1268  my_degree +
1269  j + GeometryInfo<dim>::lines_per_cell][1][0][1] =
1270  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1271  (my_degree + 1)][2][0];
1272  grad_grads
1273  [(i + 2 * k * (my_degree + 1) +
1274  GeometryInfo<dim>::lines_per_cell) *
1275  my_degree +
1276  j + GeometryInfo<dim>::lines_per_cell][1][0][2] =
1277  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1278  (my_degree + 1)][2][1];
1279  grad_grads
1280  [(i + 2 * k * (my_degree + 1) +
1281  GeometryInfo<dim>::lines_per_cell) *
1282  my_degree +
1283  j + GeometryInfo<dim>::lines_per_cell][1][1][0] =
1284  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1285  (my_degree + 1)][0][2];
1286  grad_grads
1287  [(i + 2 * k * (my_degree + 1) +
1288  GeometryInfo<dim>::lines_per_cell) *
1289  my_degree +
1290  j + GeometryInfo<dim>::lines_per_cell][1][1][1] =
1291  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1292  (my_degree + 1)][0][0];
1293  grad_grads
1294  [(i + 2 * k * (my_degree + 1) +
1295  GeometryInfo<dim>::lines_per_cell) *
1296  my_degree +
1297  j + GeometryInfo<dim>::lines_per_cell][1][1][2] =
1298  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1299  (my_degree + 1)][0][1];
1300  grad_grads
1301  [(i + 2 * k * (my_degree + 1) +
1302  GeometryInfo<dim>::lines_per_cell) *
1303  my_degree +
1304  j + GeometryInfo<dim>::lines_per_cell][1][2][0] =
1305  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1306  (my_degree + 1)][1][2];
1307  grad_grads
1308  [(i + 2 * k * (my_degree + 1) +
1309  GeometryInfo<dim>::lines_per_cell) *
1310  my_degree +
1311  j + GeometryInfo<dim>::lines_per_cell][1][2][1] =
1312  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1313  (my_degree + 1)][1][0];
1314  grad_grads
1315  [(i + 2 * k * (my_degree + 1) +
1316  GeometryInfo<dim>::lines_per_cell) *
1317  my_degree +
1318  j + GeometryInfo<dim>::lines_per_cell][1][2][2] =
1319  p1_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1320  (my_degree + 1)][1][1];
1321  grad_grads[i + (j + (2 * k + 1) * my_degree +
1322  GeometryInfo<dim>::lines_per_cell) *
1323  (my_degree + 1)][2][0][0] =
1324  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1325  (my_degree + 1)][1][1];
1326  grad_grads[i + (j + (2 * k + 1) * my_degree +
1327  GeometryInfo<dim>::lines_per_cell) *
1328  (my_degree + 1)][2][0][1] =
1329  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1330  (my_degree + 1)][1][2];
1331  grad_grads[i + (j + (2 * k + 1) * my_degree +
1332  GeometryInfo<dim>::lines_per_cell) *
1333  (my_degree + 1)][2][0][2] =
1334  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1335  (my_degree + 1)][1][0];
1336  grad_grads[i + (j + (2 * k + 1) * my_degree +
1337  GeometryInfo<dim>::lines_per_cell) *
1338  (my_degree + 1)][2][1][0] =
1339  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1340  (my_degree + 1)][2][1];
1341  grad_grads[i + (j + (2 * k + 1) * my_degree +
1342  GeometryInfo<dim>::lines_per_cell) *
1343  (my_degree + 1)][2][1][1] =
1344  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1345  (my_degree + 1)][2][2];
1346  grad_grads[i + (j + (2 * k + 1) * my_degree +
1347  GeometryInfo<dim>::lines_per_cell) *
1348  (my_degree + 1)][2][1][2] =
1349  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1350  (my_degree + 1)][2][0];
1351  grad_grads[i + (j + (2 * k + 1) * my_degree +
1352  GeometryInfo<dim>::lines_per_cell) *
1353  (my_degree + 1)][2][2][0] =
1354  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1355  (my_degree + 1)][0][1];
1356  grad_grads[i + (j + (2 * k + 1) * my_degree +
1357  GeometryInfo<dim>::lines_per_cell) *
1358  (my_degree + 1)][2][2][1] =
1359  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1360  (my_degree + 1)][0][2];
1361  grad_grads[i + (j + (2 * k + 1) * my_degree +
1362  GeometryInfo<dim>::lines_per_cell) *
1363  (my_degree + 1)][2][2][2] =
1364  p2_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1365  (my_degree + 1)][0][0];
1366  grad_grads
1367  [(i + 2 * (k + 2) * (my_degree + 1) +
1368  GeometryInfo<dim>::lines_per_cell) *
1369  my_degree +
1370  j + GeometryInfo<dim>::lines_per_cell][2][0][0] =
1371  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1372  (my_degree + 1)][1][1];
1373  grad_grads
1374  [(i + 2 * (k + 2) * (my_degree + 1) +
1375  GeometryInfo<dim>::lines_per_cell) *
1376  my_degree +
1377  j + GeometryInfo<dim>::lines_per_cell][2][0][1] =
1378  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1379  (my_degree + 1)][1][2];
1380  grad_grads
1381  [(i + 2 * (k + 2) * (my_degree + 1) +
1382  GeometryInfo<dim>::lines_per_cell) *
1383  my_degree +
1384  j + GeometryInfo<dim>::lines_per_cell][2][0][2] =
1385  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1386  (my_degree + 1)][1][0];
1387  grad_grads
1388  [(i + 2 * (k + 2) * (my_degree + 1) +
1389  GeometryInfo<dim>::lines_per_cell) *
1390  my_degree +
1391  j + GeometryInfo<dim>::lines_per_cell][2][1][0] =
1392  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1393  (my_degree + 1)][2][1];
1394  grad_grads
1395  [(i + 2 * (k + 2) * (my_degree + 1) +
1396  GeometryInfo<dim>::lines_per_cell) *
1397  my_degree +
1398  j + GeometryInfo<dim>::lines_per_cell][2][1][1] =
1399  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1400  (my_degree + 1)][2][2];
1401  grad_grads
1402  [(i + 2 * (k + 2) * (my_degree + 1) +
1403  GeometryInfo<dim>::lines_per_cell) *
1404  my_degree +
1405  j + GeometryInfo<dim>::lines_per_cell][2][1][2] =
1406  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1407  (my_degree + 1)][2][0];
1408  grad_grads
1409  [(i + 2 * (k + 2) * (my_degree + 1) +
1410  GeometryInfo<dim>::lines_per_cell) *
1411  my_degree +
1412  j + GeometryInfo<dim>::lines_per_cell][2][2][0] =
1413  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1414  (my_degree + 1)][0][1];
1415  grad_grads
1416  [(i + 2 * (k + 2) * (my_degree + 1) +
1417  GeometryInfo<dim>::lines_per_cell) *
1418  my_degree +
1419  j + GeometryInfo<dim>::lines_per_cell][2][2][1] =
1420  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1421  (my_degree + 1)][0][2];
1422  grad_grads
1423  [(i + 2 * (k + 2) * (my_degree + 1) +
1424  GeometryInfo<dim>::lines_per_cell) *
1425  my_degree +
1426  j + GeometryInfo<dim>::lines_per_cell][2][2][2] =
1427  p2_grad_grads[i + (j + k * (my_degree + 2) + 2) *
1428  (my_degree + 1)][0][0];
1429  grad_grads[i + (j + (2 * k + 9) * my_degree +
1430  GeometryInfo<dim>::lines_per_cell) *
1431  (my_degree + 1)][1][0][0] =
1432  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1433  (my_degree + 1)][2][2];
1434  grad_grads[i + (j + (2 * k + 9) * my_degree +
1435  GeometryInfo<dim>::lines_per_cell) *
1436  (my_degree + 1)][1][0][1] =
1437  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1438  (my_degree + 1)][2][0];
1439  grad_grads[i + (j + (2 * k + 9) * my_degree +
1440  GeometryInfo<dim>::lines_per_cell) *
1441  (my_degree + 1)][1][0][2] =
1442  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1443  (my_degree + 1)][2][1];
1444  grad_grads[i + (j + (2 * k + 9) * my_degree +
1445  GeometryInfo<dim>::lines_per_cell) *
1446  (my_degree + 1)][1][1][0] =
1447  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1448  (my_degree + 1)][0][2];
1449  grad_grads[i + (j + (2 * k + 9) * my_degree +
1450  GeometryInfo<dim>::lines_per_cell) *
1451  (my_degree + 1)][1][1][1] =
1452  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1453  (my_degree + 1)][0][0];
1454  grad_grads[i + (j + (2 * k + 9) * my_degree +
1455  GeometryInfo<dim>::lines_per_cell) *
1456  (my_degree + 1)][1][1][2] =
1457  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1458  (my_degree + 1)][0][1];
1459  grad_grads[i + (j + (2 * k + 9) * my_degree +
1460  GeometryInfo<dim>::lines_per_cell) *
1461  (my_degree + 1)][1][2][0] =
1462  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1463  (my_degree + 1)][1][2];
1464  grad_grads[i + (j + (2 * k + 9) * my_degree +
1465  GeometryInfo<dim>::lines_per_cell) *
1466  (my_degree + 1)][1][2][1] =
1467  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1468  (my_degree + 1)][1][0];
1469  grad_grads[i + (j + (2 * k + 9) * my_degree +
1470  GeometryInfo<dim>::lines_per_cell) *
1471  (my_degree + 1)][1][2][2] =
1472  p1_grad_grads[i + ((j + 2) * (my_degree + 2) + k) *
1473  (my_degree + 1)][1][1];
1474  }
1475  }
1476  }
1477 
1478  break;
1479  }
1480 
1481  default:
1482  Assert(false, ExcNotImplemented());
1483  }
1484 }
1485 
1486 
1487 template <int dim>
1488 unsigned int
1489 PolynomialsNedelec<dim>::n_polynomials(const unsigned int k)
1490 {
1491  switch (dim)
1492  {
1493  case 1:
1494  return k + 1;
1495 
1496  case 2:
1497  return 2 * (k + 1) * (k + 2);
1498 
1499  case 3:
1500  return 3 * (k + 1) * (k + 2) * (k + 2);
1501 
1502  default:
1503  {
1504  Assert(false, ExcNotImplemented());
1505  return 0;
1506  }
1507  }
1508 }
1509 
1510 
1511 template <int dim>
1512 std::unique_ptr<TensorPolynomialsBase<dim>>
1513 PolynomialsNedelec<dim>::clone() const
1514 {
1515  return std_cxx14::make_unique<PolynomialsNedelec<dim>>(*this);
1516 }
1517 
1518 
1519 template class PolynomialsNedelec<1>;
1520 template class PolynomialsNedelec<2>;
1521 template class PolynomialsNedelec<3>;
1522 
1523 
1524 DEAL_II_NAMESPACE_CLOSE
PolynomialsNedelec::evaluate
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
Definition: polynomials_nedelec.cc:54
PolynomialsNedelec::create_polynomials
static std::vector< std::vector< Polynomials::Polynomial< double > > > create_polynomials(const unsigned int k)
Definition: polynomials_nedelec.cc:36
Polynomials::Lobatto::generate_complete_basis
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:982
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
polynomial.h
GeometryInfo
Definition: geometry_info.h:1224
Polynomials::Legendre::generate_complete_basis
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:877
quadrature_lib.h
geometry_info.h
PolynomialsNedelec::clone
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
Definition: polynomials_nedelec.cc:1513
Tensor< 1, dim >
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TensorPolynomialsBase
Definition: tensor_polynomials_base.h:62
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
PolynomialsNedelec::PolynomialsNedelec
PolynomialsNedelec(const unsigned int k)
Definition: polynomials_nedelec.cc:29
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
memory.h
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
PolynomialsNedelec
Definition: polynomials_nedelec.h:53
polynomials_nedelec.h
PolynomialsNedelec::n_polynomials
static unsigned int n_polynomials(const unsigned int degree)
Definition: polynomials_nedelec.cc:1489
Generated by   doxygen 1.8.17