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