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