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