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