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