Reference documentation for deal.II version GIT relicensing-12-g89ba5bf9db 2024-02-24 00:00:03+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
fe_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
20
22
25#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/mapping.h>
28
29#include <deal.II/grid/tria.h>
31
33#include <deal.II/lac/vector.h>
34
35#include <iostream>
36#include <memory>
37#include <sstream>
38
40
41// #define DEBUG_NEDELEC
42
43namespace internal
44{
45 namespace FE_Nedelec
46 {
47 namespace
48 {
49 double
50 get_embedding_computation_tolerance(const unsigned int p)
51 {
52 // This heuristic was computed by monitoring the worst residual
53 // resulting from the least squares computation when computing
54 // the face embedding matrices in the FE_Nedelec constructor.
55 // The residual growth is exponential, but is bounded by this
56 // function up to degree 12.
57 return 1.e-15 * std::exp(std::pow(p, 1.075));
58 }
59 } // namespace
60 } // namespace FE_Nedelec
61} // namespace internal
62
63
64// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
65// adjust_line_dof_index_for_line_orientation_table fields, and write tests
66// similar to bits/face_orientation_and_fe_q_*
67
68template <int dim>
69FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
70 : FE_PolyTensor<dim>(
71 PolynomialsNedelec<dim>(order),
72 FiniteElementData<dim>(get_dpo_vector(order),
73 dim,
74 order + 1,
75 FiniteElementData<dim>::Hcurl),
76 std::vector<bool>(PolynomialsNedelec<dim>::n_polynomials(order), true),
77 std::vector<ComponentMask>(PolynomialsNedelec<dim>::n_polynomials(order),
78 ComponentMask(std::vector<bool>(dim, true))))
79{
80#ifdef DEBUG_NEDELEC
81 deallog << get_name() << std::endl;
82#endif
83
84 Assert(dim >= 2, ExcImpossibleInDim(dim));
85
87 // First, initialize the
88 // generalized support points and
89 // quadrature weights, since they
90 // are required for interpolation.
92
93 // We already use the correct basis, so no basis transformation is required
94 // from the polynomial space we have described above to the one that is dual
95 // to the node functionals. As documented in the base class, this is
96 // expressed by setting the inverse node matrix to the empty matrix.
97 this->inverse_node_matrix.clear();
98
99 // do not initialize embedding and restriction here. these matrices are
100 // initialized on demand in get_restriction_matrix and
101 // get_prolongation_matrix
102
103#ifdef DEBUG_NEDELEC
104 deallog << "Face Embedding" << std::endl;
105#endif
107
108 // TODO: the implementation makes the assumption that all faces have the
109 // same number of dofs
111 const unsigned int face_no = 0;
112
113 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
114 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
115 this->n_dofs_per_face(face_no));
116
117 FETools::compute_face_embedding_matrices<dim, double>(
118 *this,
119 face_embeddings,
120 0,
121 0,
122 internal::FE_Nedelec::get_embedding_computation_tolerance(order));
123
124 switch (dim)
125 {
126 case 1:
127 {
128 this->interface_constraints.reinit(0, 0);
129 break;
130 }
131
132 case 2:
133 {
134 this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no),
135 this->n_dofs_per_face(face_no));
136
137 for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
138 ++i)
139 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
140 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
141 this->interface_constraints(i * this->n_dofs_per_face(face_no) +
142 j,
143 k) = face_embeddings[i](j, k);
144
145 break;
146 }
147
148 case 3:
149 {
150 this->interface_constraints.reinit(
151 4 * (this->n_dofs_per_face(face_no) - this->degree),
152 this->n_dofs_per_face(face_no));
153
154 unsigned int target_row = 0;
155
156 for (unsigned int i = 0; i < 2; ++i)
157 for (unsigned int j = this->degree; j < 2 * this->degree;
158 ++j, ++target_row)
159 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
160 this->interface_constraints(target_row, k) =
161 face_embeddings[2 * i](j, k);
162
163 for (unsigned int i = 0; i < 2; ++i)
164 for (unsigned int j = 3 * this->degree;
165 j < GeometryInfo<3>::lines_per_face * this->degree;
166 ++j, ++target_row)
167 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
168 this->interface_constraints(target_row, k) =
169 face_embeddings[i](j, k);
170
171 for (unsigned int i = 0; i < 2; ++i)
172 for (unsigned int j = 0; j < 2; ++j)
173 for (unsigned int k = i * this->degree;
174 k < (i + 1) * this->degree;
175 ++k, ++target_row)
176 for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
177 ++l)
178 this->interface_constraints(target_row, l) =
179 face_embeddings[i + 2 * j](k, l);
180
181 for (unsigned int i = 0; i < 2; ++i)
182 for (unsigned int j = 0; j < 2; ++j)
183 for (unsigned int k = (i + 2) * this->degree;
184 k < (i + 3) * this->degree;
185 ++k, ++target_row)
186 for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
187 ++l)
188 this->interface_constraints(target_row, l) =
189 face_embeddings[2 * i + j](k, l);
190
191 for (unsigned int i = 0; i < GeometryInfo<3>::max_children_per_face;
192 ++i)
193 for (unsigned int j =
194 GeometryInfo<3>::lines_per_face * this->degree;
195 j < this->n_dofs_per_face(face_no);
196 ++j, ++target_row)
197 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
198 this->interface_constraints(target_row, k) =
199 face_embeddings[i](j, k);
200
201 break;
202 }
203
204 default:
206 }
207
208 // We need to initialize the dof permutation table and the one for the sign
209 // change.
211}
212
213
214template <int dim>
215void
217{
218 // for 1d and 2d, do nothing
219 if (dim < 3)
220 return;
221
222 // TODO: Implement this for this class
223 return;
224}
225
226
227template <int dim>
228std::string
230{
231 // note that the
232 // FETools::get_fe_by_name
233 // function depends on the
234 // particular format of the string
235 // this function returns, so they
236 // have to be kept in synch
237
238 std::ostringstream namebuf;
239 namebuf << "FE_Nedelec<" << dim << ">(" << this->degree - 1 << ")";
240
241 return namebuf.str();
242}
243
244
245template <int dim>
246std::unique_ptr<FiniteElement<dim, dim>>
248{
249 return std::make_unique<FE_Nedelec<dim>>(*this);
250}
251
252//---------------------------------------------------------------------------
253// Auxiliary and internal functions
254//---------------------------------------------------------------------------
255
256
257
258// Set the generalized support
259// points and precompute the
260// parts of the projection-based
261// interpolation, which does
262// not depend on the interpolated
263// function.
264template <>
265void
270
271
272
273template <>
274void
276{
277 const int dim = 2;
278
279 // TODO: the implementation makes the assumption that all faces have the
280 // same number of dofs
281 AssertDimension(this->n_unique_faces(), 1);
282 const unsigned int face_no = 0;
283
284 // Create polynomial basis.
285 const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
287 std::vector<Polynomials::Polynomial<double>> lobatto_polynomials_grad(order +
288 1);
289
290 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
291 lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative();
292
293 // Initialize quadratures to obtain
294 // quadrature points later on.
295 const QGauss<dim - 1> reference_edge_quadrature(order + 1);
296 const unsigned int n_edge_points = reference_edge_quadrature.size();
297 const unsigned int n_boundary_points =
298 GeometryInfo<dim>::lines_per_cell * n_edge_points;
299 const Quadrature<dim> edge_quadrature =
300 QProjector<dim>::project_to_all_faces(this->reference_cell(),
301 reference_edge_quadrature);
302
303 this->generalized_face_support_points[face_no].resize(n_edge_points);
304
305 // Create face support points.
306 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
307 this->generalized_face_support_points[face_no][q_point] =
308 reference_edge_quadrature.point(q_point);
309
310 if (order > 0)
311 {
312 // If the polynomial degree is positive
313 // we have support points on the faces
314 // and in the interior of a cell.
315 const QGauss<dim> quadrature(order + 1);
316 const unsigned int n_interior_points = quadrature.size();
317
318 this->generalized_support_points.resize(n_boundary_points +
319 n_interior_points);
320 boundary_weights.reinit(n_edge_points, order);
321
322 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
323 {
324 for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
325 ++line)
326 this->generalized_support_points[line * n_edge_points + q_point] =
327 edge_quadrature.point(
328 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
329 line,
330 true,
331 false,
332 false,
333 n_edge_points) +
334 q_point);
335
336 for (unsigned int i = 0; i < order; ++i)
337 boundary_weights(q_point, i) =
338 reference_edge_quadrature.weight(q_point) *
339 lobatto_polynomials_grad[i + 1].value(
340 this->generalized_face_support_points[face_no][q_point][0]);
341 }
342
343 for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
344 this->generalized_support_points[q_point + n_boundary_points] =
345 quadrature.point(q_point);
346 }
347
348 else
349 {
350 // In this case we only need support points
351 // on the faces of a cell.
352 this->generalized_support_points.resize(n_boundary_points);
353
354 for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
355 ++line)
356 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
357 this->generalized_support_points[line * n_edge_points + q_point] =
358 edge_quadrature.point(
359 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
360 line,
361 true,
362 false,
363 false,
364 n_edge_points) +
365 q_point);
366 }
367}
368
369
370
371template <>
372void
374{
375 const int dim = 3;
376
377 // TODO: the implementation makes the assumption that all faces have the
378 // same number of dofs
379 AssertDimension(this->n_unique_faces(), 1);
380 const unsigned int face_no = 0;
381
382 // Create polynomial basis.
383 const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
385 std::vector<Polynomials::Polynomial<double>> lobatto_polynomials_grad(order +
386 1);
387
388 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
389 lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative();
390
391 // Initialize quadratures to obtain
392 // quadrature points later on.
393 const QGauss<1> reference_edge_quadrature(order + 1);
394 const unsigned int n_edge_points = reference_edge_quadrature.size();
395 const Quadrature<dim - 1> &edge_quadrature =
397 ReferenceCells::get_hypercube<dim - 1>(), reference_edge_quadrature);
398
399 if (order > 0)
400 {
401 // If the polynomial order is positive
402 // we have support points on the edges,
403 // faces and in the interior of a cell.
404 const QGauss<dim - 1> reference_face_quadrature(order + 1);
405 const unsigned int n_face_points = reference_face_quadrature.size();
406 const unsigned int n_boundary_points =
407 GeometryInfo<dim>::lines_per_cell * n_edge_points +
408 GeometryInfo<dim>::faces_per_cell * n_face_points;
409 const QGauss<dim> quadrature(order + 1);
410 const unsigned int n_interior_points = quadrature.size();
411
412 boundary_weights.reinit(n_edge_points + n_face_points,
413 2 * (order + 1) * order);
414 this->generalized_face_support_points[face_no].resize(4 * n_edge_points +
415 n_face_points);
416 this->generalized_support_points.resize(n_boundary_points +
417 n_interior_points);
418
419 // Create support points on edges.
420 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
421 {
422 for (unsigned int line = 0;
423 line < GeometryInfo<dim - 1>::lines_per_cell;
424 ++line)
425 this
426 ->generalized_face_support_points[face_no][line * n_edge_points +
427 q_point] =
428 edge_quadrature.point(
430 ReferenceCells::get_hypercube<dim - 1>(),
431 line,
432 true,
433 false,
434 false,
435 n_edge_points) +
436 q_point);
437
438 for (unsigned int i = 0; i < 2; ++i)
439 for (unsigned int j = 0; j < 2; ++j)
440 {
441 this->generalized_support_points[q_point +
442 (i + 4 * j) * n_edge_points] =
443 Point<dim>(i, reference_edge_quadrature.point(q_point)[0], j);
444 this->generalized_support_points[q_point + (i + 4 * j + 2) *
445 n_edge_points] =
446 Point<dim>(reference_edge_quadrature.point(q_point)[0], i, j);
447 this->generalized_support_points[q_point + (i + 2 * (j + 4)) *
448 n_edge_points] =
449 Point<dim>(i, j, reference_edge_quadrature.point(q_point)[0]);
450 }
451
452 for (unsigned int i = 0; i < order; ++i)
453 boundary_weights(q_point, i) =
454 reference_edge_quadrature.weight(q_point) *
455 lobatto_polynomials_grad[i + 1].value(
456 this->generalized_face_support_points[face_no][q_point][1]);
457 }
458
459 // Create support points on faces.
460 for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
461 {
462 this->generalized_face_support_points[face_no]
463 [q_point + 4 * n_edge_points] =
464 reference_face_quadrature.point(q_point);
465
466 for (unsigned int i = 0; i <= order; ++i)
467 for (unsigned int j = 0; j < order; ++j)
468 {
469 boundary_weights(q_point + n_edge_points, 2 * (i * order + j)) =
470 reference_face_quadrature.weight(q_point) *
471 lobatto_polynomials_grad[i].value(
472 this->generalized_face_support_points
473 [face_no][q_point + 4 * n_edge_points][0]) *
474 lobatto_polynomials[j + 2].value(
475 this->generalized_face_support_points
476 [face_no][q_point + 4 * n_edge_points][1]);
477 boundary_weights(q_point + n_edge_points,
478 2 * (i * order + j) + 1) =
479 reference_face_quadrature.weight(q_point) *
480 lobatto_polynomials_grad[i].value(
481 this->generalized_face_support_points
482 [face_no][q_point + 4 * n_edge_points][1]) *
483 lobatto_polynomials[j + 2].value(
484 this->generalized_face_support_points
485 [face_no][q_point + 4 * n_edge_points][0]);
486 }
487 }
488
489 const Quadrature<dim> &face_quadrature =
490 QProjector<dim>::project_to_all_faces(this->reference_cell(),
491 reference_face_quadrature);
492
493 for (const unsigned int face : GeometryInfo<dim>::face_indices())
494 for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
495 {
496 this->generalized_support_points[face * n_face_points + q_point +
498 n_edge_points] =
499 face_quadrature.point(
500 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
501 face,
502 true,
503 false,
504 false,
505 n_face_points) +
506 q_point);
507 }
508
509 // Create support points in the interior.
510 for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
511 this->generalized_support_points[q_point + n_boundary_points] =
512 quadrature.point(q_point);
513 }
514
515 else
516 {
517 this->generalized_face_support_points[face_no].resize(4 * n_edge_points);
518 this->generalized_support_points.resize(
519 GeometryInfo<dim>::lines_per_cell * n_edge_points);
520
521 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
522 {
523 for (unsigned int line = 0;
524 line < GeometryInfo<dim - 1>::lines_per_cell;
525 ++line)
526 this
527 ->generalized_face_support_points[face_no][line * n_edge_points +
528 q_point] =
529 edge_quadrature.point(
531 ReferenceCells::get_hypercube<dim - 1>(),
532 line,
533 true,
534 false,
535 false,
536 n_edge_points) +
537 q_point);
538
539 for (unsigned int i = 0; i < 2; ++i)
540 for (unsigned int j = 0; j < 2; ++j)
541 {
542 this->generalized_support_points[q_point +
543 (i + 4 * j) * n_edge_points] =
544 Point<dim>(i, reference_edge_quadrature.point(q_point)[0], j);
545 this->generalized_support_points[q_point + (i + 4 * j + 2) *
546 n_edge_points] =
547 Point<dim>(reference_edge_quadrature.point(q_point)[0], i, j);
548 this->generalized_support_points[q_point + (i + 2 * (j + 4)) *
549 n_edge_points] =
550 Point<dim>(i, j, reference_edge_quadrature.point(q_point)[0]);
551 }
552 }
553 }
554}
555
556
557
558// Set the restriction matrices.
559template <>
560void
562{
563 // there is only one refinement case in 1d,
564 // which is the isotropic one
565 for (unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
566 this->restriction[0][i].reinit(0, 0);
567}
568
569
570
571// Restriction operator
572template <int dim>
573void
575{
576 // This function does the same as the
577 // function interpolate further below.
578 // But since the functions, which we
579 // interpolate here, are discontinuous
580 // we have to use more quadrature
581 // points as in interpolate.
582 const QGauss<1> edge_quadrature(2 * this->degree);
583 const std::vector<Point<1>> &edge_quadrature_points =
584 edge_quadrature.get_points();
585 const unsigned int n_edge_quadrature_points = edge_quadrature.size();
586 const unsigned int index = RefinementCase<dim>::isotropic_refinement - 1;
587
588 switch (dim)
589 {
590 case 2:
591 {
592 // First interpolate the shape
593 // functions of the child cells
594 // to the lowest order shape
595 // functions of the parent cell.
596 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
597 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
598 ++q_point)
599 {
600 const double weight = 2.0 * edge_quadrature.weight(q_point);
601
602 if (edge_quadrature_points[q_point][0] < 0.5)
603 {
604 Point<dim> quadrature_point(
605 0.0, 2.0 * edge_quadrature_points[q_point][0]);
606
607 this->restriction[index][0](0, dof) +=
608 weight *
609 this->shape_value_component(dof, quadrature_point, 1);
610 quadrature_point[0] = 1.0;
611 this->restriction[index][1](this->degree, dof) +=
612 weight *
613 this->shape_value_component(dof, quadrature_point, 1);
614 quadrature_point[0] = quadrature_point[1];
615 quadrature_point[1] = 0.0;
616 this->restriction[index][0](2 * this->degree, dof) +=
617 weight *
618 this->shape_value_component(dof, quadrature_point, 0);
619 quadrature_point[1] = 1.0;
620 this->restriction[index][2](3 * this->degree, dof) +=
621 weight *
622 this->shape_value_component(dof, quadrature_point, 0);
623 }
624
625 else
626 {
627 Point<dim> quadrature_point(
628 0.0, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
629
630 this->restriction[index][2](0, dof) +=
631 weight *
632 this->shape_value_component(dof, quadrature_point, 1);
633 quadrature_point[0] = 1.0;
634 this->restriction[index][3](this->degree, dof) +=
635 weight *
636 this->shape_value_component(dof, quadrature_point, 1);
637 quadrature_point[0] = quadrature_point[1];
638 quadrature_point[1] = 0.0;
639 this->restriction[index][1](2 * this->degree, dof) +=
640 weight *
641 this->shape_value_component(dof, quadrature_point, 0);
642 quadrature_point[1] = 1.0;
643 this->restriction[index][3](3 * this->degree, dof) +=
644 weight *
645 this->shape_value_component(dof, quadrature_point, 0);
646 }
647 }
648
649 // Then project the shape functions
650 // of the child cells to the higher
651 // order shape functions of the
652 // parent cell.
653 if (this->degree > 1)
654 {
655 const unsigned int deg = this->degree - 1;
656 const std::vector<Polynomials::Polynomial<double>>
657 &legendre_polynomials =
659 FullMatrix<double> system_matrix_inv(deg, deg);
660
661 {
662 FullMatrix<double> assembling_matrix(deg,
663 n_edge_quadrature_points);
664
665 for (unsigned int q_point = 0;
666 q_point < n_edge_quadrature_points;
667 ++q_point)
668 {
669 const double weight =
670 std::sqrt(edge_quadrature.weight(q_point));
671
672 for (unsigned int i = 0; i < deg; ++i)
673 assembling_matrix(i, q_point) =
674 weight * legendre_polynomials[i + 1].value(
675 edge_quadrature_points[q_point][0]);
676 }
677
678 FullMatrix<double> system_matrix(deg, deg);
679
680 assembling_matrix.mTmult(system_matrix, assembling_matrix);
681 system_matrix_inv.invert(system_matrix);
682 }
683
684 FullMatrix<double> solution(this->degree - 1, 4);
685 FullMatrix<double> system_rhs(this->degree - 1, 4);
686 Vector<double> tmp(4);
687
688 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
689 for (unsigned int i = 0; i < 2; ++i)
690 {
691 system_rhs = 0.0;
692
693 for (unsigned int q_point = 0;
694 q_point < n_edge_quadrature_points;
695 ++q_point)
696 {
697 const double weight = edge_quadrature.weight(q_point);
698 const Point<dim> quadrature_point_0(
699 i, edge_quadrature_points[q_point][0]);
700 const Point<dim> quadrature_point_1(
701 edge_quadrature_points[q_point][0], i);
702
703 if (edge_quadrature_points[q_point][0] < 0.5)
704 {
705 Point<dim> quadrature_point_2(
706 i, 2.0 * edge_quadrature_points[q_point][0]);
707
708 tmp(0) =
709 weight *
710 (2.0 * this->shape_value_component(
711 dof, quadrature_point_2, 1) -
712 this->restriction[index][i](i * this->degree,
713 dof) *
714 this->shape_value_component(i * this->degree,
715 quadrature_point_0,
716 1));
717 tmp(1) =
718 -1.0 * weight *
719 this->restriction[index][i + 2](i * this->degree,
720 dof) *
721 this->shape_value_component(i * this->degree,
722 quadrature_point_0,
723 1);
724 quadrature_point_2 = Point<dim>(
725 2.0 * edge_quadrature_points[q_point][0], i);
726 tmp(2) =
727 weight *
728 (2.0 * this->shape_value_component(
729 dof, quadrature_point_2, 0) -
730 this->restriction[index][2 * i]((i + 2) *
731 this->degree,
732 dof) *
733 this->shape_value_component((i + 2) *
734 this->degree,
735 quadrature_point_1,
736 0));
737 tmp(3) =
738 -1.0 * weight *
739 this->restriction[index][2 * i + 1](
740 (i + 2) * this->degree, dof) *
741 this->shape_value_component(
742 (i + 2) * this->degree, quadrature_point_1, 0);
743 }
744
745 else
746 {
747 tmp(0) =
748 -1.0 * weight *
749 this->restriction[index][i](i * this->degree,
750 dof) *
751 this->shape_value_component(i * this->degree,
752 quadrature_point_0,
753 1);
754
755 Point<dim> quadrature_point_2(
756 i,
757 2.0 * edge_quadrature_points[q_point][0] - 1.0);
758
759 tmp(1) =
760 weight *
761 (2.0 * this->shape_value_component(
762 dof, quadrature_point_2, 1) -
763 this->restriction[index][i + 2](i * this->degree,
764 dof) *
765 this->shape_value_component(i * this->degree,
766 quadrature_point_0,
767 1));
768 tmp(2) =
769 -1.0 * weight *
770 this->restriction[index][2 * i]((i + 2) *
771 this->degree,
772 dof) *
773 this->shape_value_component(
774 (i + 2) * this->degree, quadrature_point_1, 0);
775 quadrature_point_2 = Point<dim>(
776 2.0 * edge_quadrature_points[q_point][0] - 1.0,
777 i);
778 tmp(3) =
779 weight *
780 (2.0 * this->shape_value_component(
781 dof, quadrature_point_2, 0) -
782 this->restriction[index][2 * i + 1](
783 (i + 2) * this->degree, dof) *
784 this->shape_value_component((i + 2) *
785 this->degree,
786 quadrature_point_1,
787 0));
788 }
789
790 for (unsigned int j = 0; j < this->degree - 1; ++j)
791 {
792 const double L_j =
793 legendre_polynomials[j + 1].value(
794 edge_quadrature_points[q_point][0]);
795
796 for (unsigned int k = 0; k < tmp.size(); ++k)
797 system_rhs(j, k) += tmp(k) * L_j;
798 }
799 }
800
801 system_matrix_inv.mmult(solution, system_rhs);
802
803 for (unsigned int j = 0; j < this->degree - 1; ++j)
804 for (unsigned int k = 0; k < 2; ++k)
805 {
806 if (std::abs(solution(j, k)) > 1e-14)
807 this->restriction[index][i + 2 * k](
808 i * this->degree + j + 1, dof) = solution(j, k);
809
810 if (std::abs(solution(j, k + 2)) > 1e-14)
811 this->restriction[index][2 * i + k](
812 (i + 2) * this->degree + j + 1, dof) =
813 solution(j, k + 2);
814 }
815 }
816
817 const QGauss<dim> quadrature(2 * this->degree);
818 const std::vector<Point<dim>> &quadrature_points =
819 quadrature.get_points();
820 const std::vector<Polynomials::Polynomial<double>>
821 &lobatto_polynomials =
823 const unsigned int n_boundary_dofs =
825 const unsigned int n_quadrature_points = quadrature.size();
826
827 {
828 FullMatrix<double> assembling_matrix((this->degree - 1) *
829 this->degree,
830 n_quadrature_points);
831
832 for (unsigned int q_point = 0; q_point < n_quadrature_points;
833 ++q_point)
834 {
835 const double weight = std::sqrt(quadrature.weight(q_point));
836
837 for (unsigned int i = 0; i < this->degree; ++i)
838 {
839 const double L_i =
840 weight * legendre_polynomials[i].value(
841 quadrature_points[q_point][0]);
842
843 for (unsigned int j = 0; j < this->degree - 1; ++j)
844 assembling_matrix(i * (this->degree - 1) + j,
845 q_point) =
846 L_i * lobatto_polynomials[j + 2].value(
847 quadrature_points[q_point][1]);
848 }
849 }
850
851 FullMatrix<double> system_matrix(assembling_matrix.m(),
852 assembling_matrix.m());
853
854 assembling_matrix.mTmult(system_matrix, assembling_matrix);
855 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
856 system_matrix_inv.invert(system_matrix);
857 }
858
859 solution.reinit(system_matrix_inv.m(), 8);
860 system_rhs.reinit(system_matrix_inv.m(), 8);
861 tmp.reinit(8);
862
863 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
864 {
865 system_rhs = 0.0;
866
867 for (unsigned int q_point = 0; q_point < n_quadrature_points;
868 ++q_point)
869 {
870 tmp = 0.0;
871
872 if (quadrature_points[q_point][0] < 0.5)
873 {
874 if (quadrature_points[q_point][1] < 0.5)
875 {
876 const Point<dim> quadrature_point(
877 2.0 * quadrature_points[q_point][0],
878 2.0 * quadrature_points[q_point][1]);
879
880 tmp(0) += 2.0 * this->shape_value_component(
881 dof, quadrature_point, 0);
882 tmp(1) += 2.0 * this->shape_value_component(
883 dof, quadrature_point, 1);
884 }
885
886 else
887 {
888 const Point<dim> quadrature_point(
889 2.0 * quadrature_points[q_point][0],
890 2.0 * quadrature_points[q_point][1] - 1.0);
891
892 tmp(4) += 2.0 * this->shape_value_component(
893 dof, quadrature_point, 0);
894 tmp(5) += 2.0 * this->shape_value_component(
895 dof, quadrature_point, 1);
896 }
897 }
898
899 else if (quadrature_points[q_point][1] < 0.5)
900 {
901 const Point<dim> quadrature_point(
902 2.0 * quadrature_points[q_point][0] - 1.0,
903 2.0 * quadrature_points[q_point][1]);
904
905 tmp(2) +=
906 2.0 * this->shape_value_component(dof,
907 quadrature_point,
908 0);
909 tmp(3) +=
910 2.0 * this->shape_value_component(dof,
911 quadrature_point,
912 1);
913 }
914
915 else
916 {
917 const Point<dim> quadrature_point(
918 2.0 * quadrature_points[q_point][0] - 1.0,
919 2.0 * quadrature_points[q_point][1] - 1.0);
920
921 tmp(6) +=
922 2.0 * this->shape_value_component(dof,
923 quadrature_point,
924 0);
925 tmp(7) +=
926 2.0 * this->shape_value_component(dof,
927 quadrature_point,
928 1);
929 }
930
931 for (unsigned int i = 0; i < 2; ++i)
932 for (unsigned int j = 0; j < this->degree; ++j)
933 {
934 tmp(2 * i) -=
935 this->restriction[index][i](j + 2 * this->degree,
936 dof) *
937 this->shape_value_component(
938 j + 2 * this->degree,
939 quadrature_points[q_point],
940 0);
941 tmp(2 * i + 1) -=
942 this->restriction[index][i](i * this->degree + j,
943 dof) *
944 this->shape_value_component(
945 i * this->degree + j,
946 quadrature_points[q_point],
947 1);
948 tmp(2 * (i + 2)) -= this->restriction[index][i + 2](
949 j + 3 * this->degree, dof) *
950 this->shape_value_component(
951 j + 3 * this->degree,
952 quadrature_points[q_point],
953 0);
954 tmp(2 * i + 5) -= this->restriction[index][i + 2](
955 i * this->degree + j, dof) *
956 this->shape_value_component(
957 i * this->degree + j,
958 quadrature_points[q_point],
959 1);
960 }
961
962 tmp *= quadrature.weight(q_point);
963
964 for (unsigned int i = 0; i < this->degree; ++i)
965 {
966 const double L_i_0 = legendre_polynomials[i].value(
967 quadrature_points[q_point][0]);
968 const double L_i_1 = legendre_polynomials[i].value(
969 quadrature_points[q_point][1]);
970
971 for (unsigned int j = 0; j < this->degree - 1; ++j)
972 {
973 const double l_j_0 =
974 L_i_0 * lobatto_polynomials[j + 2].value(
975 quadrature_points[q_point][1]);
976 const double l_j_1 =
977 L_i_1 * lobatto_polynomials[j + 2].value(
978 quadrature_points[q_point][0]);
979
980 for (unsigned int k = 0; k < 4; ++k)
981 {
982 system_rhs(i * (this->degree - 1) + j,
983 2 * k) += tmp(2 * k) * l_j_0;
984 system_rhs(i * (this->degree - 1) + j,
985 2 * k + 1) +=
986 tmp(2 * k + 1) * l_j_1;
987 }
988 }
989 }
990 }
991
992 system_matrix_inv.mmult(solution, system_rhs);
993
994 for (unsigned int i = 0; i < this->degree; ++i)
995 for (unsigned int j = 0; j < this->degree - 1; ++j)
996 for (unsigned int k = 0; k < 4; ++k)
997 {
998 if (std::abs(solution(i * (this->degree - 1) + j,
999 2 * k)) > 1e-14)
1000 this->restriction[index][k](i * (this->degree - 1) +
1001 j + n_boundary_dofs,
1002 dof) =
1003 solution(i * (this->degree - 1) + j, 2 * k);
1004
1005 if (std::abs(solution(i * (this->degree - 1) + j,
1006 2 * k + 1)) > 1e-14)
1007 this->restriction[index][k](
1008 i + (this->degree - 1 + j) * this->degree +
1009 n_boundary_dofs,
1010 dof) =
1011 solution(i * (this->degree - 1) + j, 2 * k + 1);
1012 }
1013 }
1014 }
1015
1016 break;
1017 }
1018
1019 case 3:
1020 {
1021 // First interpolate the shape
1022 // functions of the child cells
1023 // to the lowest order shape
1024 // functions of the parent cell.
1025 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1026 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
1027 ++q_point)
1028 {
1029 const double weight = 2.0 * edge_quadrature.weight(q_point);
1030
1031 if (edge_quadrature_points[q_point][0] < 0.5)
1032 for (unsigned int i = 0; i < 2; ++i)
1033 for (unsigned int j = 0; j < 2; ++j)
1034 {
1035 Point<dim> quadrature_point(
1036 i, 2.0 * edge_quadrature_points[q_point][0], j);
1037
1038 this->restriction[index][i + 4 * j]((i + 4 * j) *
1039 this->degree,
1040 dof) +=
1041 weight *
1042 this->shape_value_component(dof, quadrature_point, 1);
1043 quadrature_point =
1044 Point<dim>(2.0 * edge_quadrature_points[q_point][0],
1045 i,
1046 j);
1047 this->restriction[index][2 * (i + 2 * j)](
1048 (i + 4 * j + 2) * this->degree, dof) +=
1049 weight *
1050 this->shape_value_component(dof, quadrature_point, 0);
1051 quadrature_point =
1052 Point<dim>(i,
1053 j,
1054 2.0 * edge_quadrature_points[q_point][0]);
1055 this->restriction[index][i + 2 * j]((i + 2 * (j + 4)) *
1056 this->degree,
1057 dof) +=
1058 weight *
1059 this->shape_value_component(dof, quadrature_point, 2);
1060 }
1061
1062 else
1063 for (unsigned int i = 0; i < 2; ++i)
1064 for (unsigned int j = 0; j < 2; ++j)
1065 {
1066 Point<dim> quadrature_point(
1067 i, 2.0 * edge_quadrature_points[q_point][0] - 1.0, j);
1068
1069 this->restriction[index][i + 4 * j + 2]((i + 4 * j) *
1070 this->degree,
1071 dof) +=
1072 weight *
1073 this->shape_value_component(dof, quadrature_point, 1);
1074 quadrature_point = Point<dim>(
1075 2.0 * edge_quadrature_points[q_point][0] - 1.0, i, j);
1076 this->restriction[index][2 * (i + 2 * j) + 1](
1077 (i + 4 * j + 2) * this->degree, dof) +=
1078 weight *
1079 this->shape_value_component(dof, quadrature_point, 0);
1080 quadrature_point = Point<dim>(
1081 i, j, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1082 this->restriction[index][i + 2 * (j + 2)](
1083 (i + 2 * (j + 4)) * this->degree, dof) +=
1084 weight *
1085 this->shape_value_component(dof, quadrature_point, 2);
1086 }
1087 }
1088
1089 // Then project the shape functions
1090 // of the child cells to the higher
1091 // order shape functions of the
1092 // parent cell.
1093 if (this->degree > 1)
1094 {
1095 const unsigned int deg = this->degree - 1;
1096 const std::vector<Polynomials::Polynomial<double>>
1097 &legendre_polynomials =
1099 FullMatrix<double> system_matrix_inv(deg, deg);
1100
1101 {
1102 FullMatrix<double> assembling_matrix(deg,
1103 n_edge_quadrature_points);
1104
1105 for (unsigned int q_point = 0;
1106 q_point < n_edge_quadrature_points;
1107 ++q_point)
1108 {
1109 const double weight =
1110 std::sqrt(edge_quadrature.weight(q_point));
1111
1112 for (unsigned int i = 0; i < deg; ++i)
1113 assembling_matrix(i, q_point) =
1114 weight * legendre_polynomials[i + 1].value(
1115 edge_quadrature_points[q_point][0]);
1116 }
1117
1118 FullMatrix<double> system_matrix(deg, deg);
1119
1120 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1121 system_matrix_inv.invert(system_matrix);
1122 }
1123
1124 FullMatrix<double> solution(deg, 6);
1125 FullMatrix<double> system_rhs(deg, 6);
1126 Vector<double> tmp(6);
1127
1128 for (unsigned int i = 0; i < 2; ++i)
1129 for (unsigned int j = 0; j < 2; ++j)
1130 for (unsigned int dof = 0; dof < this->n_dofs_per_cell();
1131 ++dof)
1132 {
1133 system_rhs = 0.0;
1134
1135 for (unsigned int q_point = 0;
1136 q_point < n_edge_quadrature_points;
1137 ++q_point)
1138 {
1139 const double weight = edge_quadrature.weight(q_point);
1140 const Point<dim> quadrature_point_0(
1141 i, edge_quadrature_points[q_point][0], j);
1142 const Point<dim> quadrature_point_1(
1143 edge_quadrature_points[q_point][0], i, j);
1144 const Point<dim> quadrature_point_2(
1145 i, j, edge_quadrature_points[q_point][0]);
1146
1147 if (edge_quadrature_points[q_point][0] < 0.5)
1148 {
1149 Point<dim> quadrature_point_3(
1150 i, 2.0 * edge_quadrature_points[q_point][0], j);
1151
1152 tmp(0) =
1153 weight * (2.0 * this->shape_value_component(
1154 dof, quadrature_point_3, 1) -
1155 this->restriction[index][i + 4 * j](
1156 (i + 4 * j) * this->degree, dof) *
1157 this->shape_value_component(
1158 (i + 4 * j) * this->degree,
1159 quadrature_point_0,
1160 1));
1161 tmp(1) =
1162 -1.0 * weight *
1163 this->restriction[index][i + 4 * j + 2](
1164 (i + 4 * j) * this->degree, dof) *
1165 this->shape_value_component((i + 4 * j) *
1166 this->degree,
1167 quadrature_point_0,
1168 1);
1169 quadrature_point_3 = Point<dim>(
1170 2.0 * edge_quadrature_points[q_point][0], i, j);
1171 tmp(2) =
1172 weight *
1173 (2.0 * this->shape_value_component(
1174 dof, quadrature_point_3, 0) -
1175 this->restriction[index][2 * (i + 2 * j)](
1176 (i + 4 * j + 2) * this->degree, dof) *
1177 this->shape_value_component(
1178 (i + 4 * j + 2) * this->degree,
1179 quadrature_point_1,
1180 0));
1181 tmp(3) =
1182 -1.0 * weight *
1183 this->restriction[index][2 * (i + 2 * j) + 1](
1184 (i + 4 * j + 2) * this->degree, dof) *
1185 this->shape_value_component((i + 4 * j + 2) *
1186 this->degree,
1187 quadrature_point_1,
1188 0);
1189 quadrature_point_3 = Point<dim>(
1190 i, j, 2.0 * edge_quadrature_points[q_point][0]);
1191 tmp(4) =
1192 weight *
1193 (2.0 * this->shape_value_component(
1194 dof, quadrature_point_3, 2) -
1195 this->restriction[index][i + 2 * j](
1196 (i + 2 * (j + 4)) * this->degree, dof) *
1197 this->shape_value_component(
1198 (i + 2 * (j + 4)) * this->degree,
1199 quadrature_point_2,
1200 2));
1201 tmp(5) =
1202 -1.0 * weight *
1203 this->restriction[index][i + 2 * (j + 2)](
1204 (i + 2 * (j + 4)) * this->degree, dof) *
1205 this->shape_value_component((i + 2 * (j + 4)) *
1206 this->degree,
1207 quadrature_point_2,
1208 2);
1209 }
1210
1211 else
1212 {
1213 tmp(0) =
1214 -1.0 * weight *
1215 this->restriction[index][i + 4 * j](
1216 (i + 4 * j) * this->degree, dof) *
1217 this->shape_value_component((i + 4 * j) *
1218 this->degree,
1219 quadrature_point_0,
1220 1);
1221
1222 Point<dim> quadrature_point_3(
1223 i,
1224 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1225 j);
1226
1227 tmp(1) = weight *
1228 (2.0 * this->shape_value_component(
1229 dof, quadrature_point_3, 1) -
1230 this->restriction[index][i + 4 * j + 2](
1231 (i + 4 * j) * this->degree, dof) *
1232 this->shape_value_component(
1233 (i + 4 * j) * this->degree,
1234 quadrature_point_0,
1235 1));
1236 tmp(2) =
1237 -1.0 * weight *
1238 this->restriction[index][2 * (i + 2 * j)](
1239 (i + 4 * j + 2) * this->degree, dof) *
1240 this->shape_value_component((i + 4 * j + 2) *
1241 this->degree,
1242 quadrature_point_1,
1243 0);
1244 quadrature_point_3 = Point<dim>(
1245 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1246 i,
1247 j);
1248 tmp(3) =
1249 weight *
1250 (2.0 * this->shape_value_component(
1251 dof, quadrature_point_3, 0) -
1252 this->restriction[index][2 * (i + 2 * j) + 1](
1253 (i + 4 * j + 2) * this->degree, dof) *
1254 this->shape_value_component(
1255 (i + 4 * j + 2) * this->degree,
1256 quadrature_point_1,
1257 0));
1258 tmp(4) =
1259 -1.0 * weight *
1260 this->restriction[index][i + 2 * j](
1261 (i + 2 * (j + 4)) * this->degree, dof) *
1262 this->shape_value_component((i + 2 * (j + 4)) *
1263 this->degree,
1264 quadrature_point_2,
1265 2);
1266 quadrature_point_3 = Point<dim>(
1267 i,
1268 j,
1269 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1270 tmp(5) =
1271 weight *
1272 (2.0 * this->shape_value_component(
1273 dof, quadrature_point_3, 2) -
1274 this->restriction[index][i + 2 * (j + 2)](
1275 (i + 2 * (j + 4)) * this->degree, dof) *
1276 this->shape_value_component(
1277 (i + 2 * (j + 4)) * this->degree,
1278 quadrature_point_2,
1279 2));
1280 }
1281
1282 for (unsigned int k = 0; k < deg; ++k)
1283 {
1284 const double L_k =
1285 legendre_polynomials[k + 1].value(
1286 edge_quadrature_points[q_point][0]);
1287
1288 for (unsigned int l = 0; l < tmp.size(); ++l)
1289 system_rhs(k, l) += tmp(l) * L_k;
1290 }
1291 }
1292
1293 system_matrix_inv.mmult(solution, system_rhs);
1294
1295 for (unsigned int k = 0; k < 2; ++k)
1296 for (unsigned int l = 0; l < deg; ++l)
1297 {
1298 if (std::abs(solution(l, k)) > 1e-14)
1299 this->restriction[index][i + 2 * (2 * j + k)](
1300 (i + 4 * j) * this->degree + l + 1, dof) =
1301 solution(l, k);
1302
1303 if (std::abs(solution(l, k + 2)) > 1e-14)
1304 this->restriction[index][2 * (i + 2 * j) + k](
1305 (i + 4 * j + 2) * this->degree + l + 1, dof) =
1306 solution(l, k + 2);
1307
1308 if (std::abs(solution(l, k + 4)) > 1e-14)
1309 this->restriction[index][i + 2 * (j + 2 * k)](
1310 (i + 2 * (j + 4)) * this->degree + l + 1, dof) =
1311 solution(l, k + 4);
1312 }
1313 }
1314
1315 const QGauss<2> face_quadrature(2 * this->degree);
1316 const std::vector<Point<2>> &face_quadrature_points =
1317 face_quadrature.get_points();
1318 const std::vector<Polynomials::Polynomial<double>>
1319 &lobatto_polynomials =
1321 const unsigned int n_edge_dofs =
1322 GeometryInfo<dim>::lines_per_cell * this->degree;
1323 const unsigned int n_face_quadrature_points =
1324 face_quadrature.size();
1325
1326 {
1327 FullMatrix<double> assembling_matrix(deg * this->degree,
1328 n_face_quadrature_points);
1329
1330 for (unsigned int q_point = 0;
1331 q_point < n_face_quadrature_points;
1332 ++q_point)
1333 {
1334 const double weight =
1335 std::sqrt(face_quadrature.weight(q_point));
1336
1337 for (unsigned int i = 0; i <= deg; ++i)
1338 {
1339 const double L_i =
1340 weight * legendre_polynomials[i].value(
1341 face_quadrature_points[q_point][0]);
1342
1343 for (unsigned int j = 0; j < deg; ++j)
1344 assembling_matrix(i * deg + j, q_point) =
1345 L_i * lobatto_polynomials[j + 2].value(
1346 face_quadrature_points[q_point][1]);
1347 }
1348 }
1349
1350 FullMatrix<double> system_matrix(assembling_matrix.m(),
1351 assembling_matrix.m());
1352
1353 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1354 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
1355 system_matrix_inv.invert(system_matrix);
1356 }
1357
1358 solution.reinit(system_matrix_inv.m(), 24);
1359 system_rhs.reinit(system_matrix_inv.m(), 24);
1360 tmp.reinit(24);
1361
1362 for (unsigned int i = 0; i < 2; ++i)
1363 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1364 {
1365 system_rhs = 0.0;
1366
1367 for (unsigned int q_point = 0;
1368 q_point < n_face_quadrature_points;
1369 ++q_point)
1370 {
1371 tmp = 0.0;
1372
1373 if (face_quadrature_points[q_point][0] < 0.5)
1374 {
1375 if (face_quadrature_points[q_point][1] < 0.5)
1376 {
1377 Point<dim> quadrature_point_0(
1378 i,
1379 2.0 * face_quadrature_points[q_point][0],
1380 2.0 * face_quadrature_points[q_point][1]);
1381
1382 tmp(0) += 2.0 * this->shape_value_component(
1383 dof, quadrature_point_0, 1);
1384 tmp(1) += 2.0 * this->shape_value_component(
1385 dof, quadrature_point_0, 2);
1386 quadrature_point_0 = Point<dim>(
1387 2.0 * face_quadrature_points[q_point][0],
1388 i,
1389 2.0 * face_quadrature_points[q_point][1]);
1390 tmp(8) += 2.0 * this->shape_value_component(
1391 dof, quadrature_point_0, 2);
1392 tmp(9) += 2.0 * this->shape_value_component(
1393 dof, quadrature_point_0, 0);
1394 quadrature_point_0 = Point<dim>(
1395 2.0 * face_quadrature_points[q_point][0],
1396 2.0 * face_quadrature_points[q_point][1],
1397 i);
1398 tmp(16) += 2.0 * this->shape_value_component(
1399 dof, quadrature_point_0, 0);
1400 tmp(17) += 2.0 * this->shape_value_component(
1401 dof, quadrature_point_0, 1);
1402 }
1403
1404 else
1405 {
1406 Point<dim> quadrature_point_0(
1407 i,
1408 2.0 * face_quadrature_points[q_point][0],
1409 2.0 * face_quadrature_points[q_point][1] -
1410 1.0);
1411
1412 tmp(2) += 2.0 * this->shape_value_component(
1413 dof, quadrature_point_0, 1);
1414 tmp(3) += 2.0 * this->shape_value_component(
1415 dof, quadrature_point_0, 2);
1416 quadrature_point_0 = Point<dim>(
1417 2.0 * face_quadrature_points[q_point][0],
1418 i,
1419 2.0 * face_quadrature_points[q_point][1] -
1420 1.0);
1421 tmp(10) += 2.0 * this->shape_value_component(
1422 dof, quadrature_point_0, 2);
1423 tmp(11) += 2.0 * this->shape_value_component(
1424 dof, quadrature_point_0, 0);
1425 quadrature_point_0 = Point<dim>(
1426 2.0 * face_quadrature_points[q_point][0],
1427 2.0 * face_quadrature_points[q_point][1] -
1428 1.0,
1429 i);
1430 tmp(18) += 2.0 * this->shape_value_component(
1431 dof, quadrature_point_0, 0);
1432 tmp(19) += 2.0 * this->shape_value_component(
1433 dof, quadrature_point_0, 1);
1434 }
1435 }
1436
1437 else if (face_quadrature_points[q_point][1] < 0.5)
1438 {
1439 Point<dim> quadrature_point_0(
1440 i,
1441 2.0 * face_quadrature_points[q_point][0] - 1.0,
1442 2.0 * face_quadrature_points[q_point][1]);
1443
1444 tmp(4) += 2.0 * this->shape_value_component(
1445 dof, quadrature_point_0, 1);
1446 tmp(5) += 2.0 * this->shape_value_component(
1447 dof, quadrature_point_0, 2);
1448 quadrature_point_0 = Point<dim>(
1449 2.0 * face_quadrature_points[q_point][0] - 1.0,
1450 i,
1451 2.0 * face_quadrature_points[q_point][1]);
1452 tmp(12) += 2.0 * this->shape_value_component(
1453 dof, quadrature_point_0, 2);
1454 tmp(13) += 2.0 * this->shape_value_component(
1455 dof, quadrature_point_0, 0);
1456 quadrature_point_0 = Point<dim>(
1457 2.0 * face_quadrature_points[q_point][0] - 1.0,
1458 2.0 * face_quadrature_points[q_point][1],
1459 i);
1460 tmp(20) += 2.0 * this->shape_value_component(
1461 dof, quadrature_point_0, 0);
1462 tmp(21) += 2.0 * this->shape_value_component(
1463 dof, quadrature_point_0, 1);
1464 }
1465
1466 else
1467 {
1468 Point<dim> quadrature_point_0(
1469 i,
1470 2.0 * face_quadrature_points[q_point][0] - 1.0,
1471 2.0 * face_quadrature_points[q_point][1] - 1.0);
1472
1473 tmp(6) += 2.0 * this->shape_value_component(
1474 dof, quadrature_point_0, 1);
1475 tmp(7) += 2.0 * this->shape_value_component(
1476 dof, quadrature_point_0, 2);
1477 quadrature_point_0 = Point<dim>(
1478 2.0 * face_quadrature_points[q_point][0] - 1.0,
1479 i,
1480 2.0 * face_quadrature_points[q_point][1] - 1.0);
1481 tmp(14) += 2.0 * this->shape_value_component(
1482 dof, quadrature_point_0, 2);
1483 tmp(15) += 2.0 * this->shape_value_component(
1484 dof, quadrature_point_0, 0);
1485 quadrature_point_0 = Point<dim>(
1486 2.0 * face_quadrature_points[q_point][0] - 1.0,
1487 2.0 * face_quadrature_points[q_point][1] - 1.0,
1488 i);
1489 tmp(22) += 2.0 * this->shape_value_component(
1490 dof, quadrature_point_0, 0);
1491 tmp(23) += 2.0 * this->shape_value_component(
1492 dof, quadrature_point_0, 1);
1493 }
1494
1495 const Point<dim> quadrature_point_0(
1496 i,
1497 face_quadrature_points[q_point][0],
1498 face_quadrature_points[q_point][1]);
1499 const Point<dim> quadrature_point_1(
1500 face_quadrature_points[q_point][0],
1501 i,
1502 face_quadrature_points[q_point][1]);
1503 const Point<dim> quadrature_point_2(
1504 face_quadrature_points[q_point][0],
1505 face_quadrature_points[q_point][1],
1506 i);
1507
1508 for (unsigned int j = 0; j < 2; ++j)
1509 for (unsigned int k = 0; k < 2; ++k)
1510 for (unsigned int l = 0; l <= deg; ++l)
1511 {
1512 tmp(2 * (j + 2 * k)) -=
1513 this->restriction[index][i + 2 * (2 * j + k)](
1514 (i + 4 * j) * this->degree + l, dof) *
1515 this->shape_value_component(
1516 (i + 4 * j) * this->degree + l,
1517 quadrature_point_0,
1518 1);
1519 tmp(2 * (j + 2 * k) + 1) -=
1520 this->restriction[index][i + 2 * (2 * j + k)](
1521 (i + 2 * (k + 4)) * this->degree + l, dof) *
1522 this->shape_value_component(
1523 (i + 2 * (k + 4)) * this->degree + l,
1524 quadrature_point_0,
1525 2);
1526 tmp(2 * (j + 2 * (k + 2))) -=
1527 this->restriction[index][2 * (i + 2 * j) + k](
1528 (2 * (i + 4) + k) * this->degree + l, dof) *
1529 this->shape_value_component(
1530 (2 * (i + 4) + k) * this->degree + l,
1531 quadrature_point_1,
1532 2);
1533 tmp(2 * (j + 2 * k) + 9) -=
1534 this->restriction[index][2 * (i + 2 * j) + k](
1535 (i + 4 * j + 2) * this->degree + l, dof) *
1536 this->shape_value_component(
1537 (i + 4 * j + 2) * this->degree + l,
1538 quadrature_point_1,
1539 0);
1540 tmp(2 * (j + 2 * (k + 4))) -=
1541 this->restriction[index][2 * (2 * i + j) + k](
1542 (4 * i + j + 2) * this->degree + l, dof) *
1543 this->shape_value_component(
1544 (4 * i + j + 2) * this->degree + l,
1545 quadrature_point_2,
1546 0);
1547 tmp(2 * (j + 2 * k) + 17) -=
1548 this->restriction[index][2 * (2 * i + j) + k](
1549 (4 * i + k) * this->degree + l, dof) *
1550 this->shape_value_component(
1551 (4 * i + k) * this->degree + l,
1552 quadrature_point_2,
1553 1);
1554 }
1555
1556 tmp *= face_quadrature.weight(q_point);
1557
1558 for (unsigned int j = 0; j <= deg; ++j)
1559 {
1560 const double L_j_0 = legendre_polynomials[j].value(
1561 face_quadrature_points[q_point][0]);
1562 const double L_j_1 = legendre_polynomials[j].value(
1563 face_quadrature_points[q_point][1]);
1564
1565 for (unsigned int k = 0; k < deg; ++k)
1566 {
1567 const double l_k_0 =
1568 L_j_0 * lobatto_polynomials[k + 2].value(
1569 face_quadrature_points[q_point][1]);
1570 const double l_k_1 =
1571 L_j_1 * lobatto_polynomials[k + 2].value(
1572 face_quadrature_points[q_point][0]);
1573
1574 for (unsigned int l = 0; l < 4; ++l)
1575 {
1576 system_rhs(j * deg + k, 2 * l) +=
1577 tmp(2 * l) * l_k_0;
1578 system_rhs(j * deg + k, 2 * l + 1) +=
1579 tmp(2 * l + 1) * l_k_1;
1580 system_rhs(j * deg + k, 2 * (l + 4)) +=
1581 tmp(2 * (l + 4)) * l_k_1;
1582 system_rhs(j * deg + k, 2 * l + 9) +=
1583 tmp(2 * l + 9) * l_k_0;
1584 system_rhs(j * deg + k, 2 * (l + 8)) +=
1585 tmp(2 * (l + 8)) * l_k_0;
1586 system_rhs(j * deg + k, 2 * l + 17) +=
1587 tmp(2 * l + 17) * l_k_1;
1588 }
1589 }
1590 }
1591 }
1592
1593 system_matrix_inv.mmult(solution, system_rhs);
1594
1595 for (unsigned int j = 0; j < 2; ++j)
1596 for (unsigned int k = 0; k < 2; ++k)
1597 for (unsigned int l = 0; l <= deg; ++l)
1598 for (unsigned int m = 0; m < deg; ++m)
1599 {
1600 if (std::abs(solution(l * deg + m,
1601 2 * (j + 2 * k))) > 1e-14)
1602 this->restriction[index][i + 2 * (2 * j + k)](
1603 (2 * i * this->degree + l) * deg + m +
1604 n_edge_dofs,
1605 dof) = solution(l * deg + m, 2 * (j + 2 * k));
1606
1607 if (std::abs(solution(l * deg + m,
1608 2 * (j + 2 * k) + 1)) >
1609 1e-14)
1610 this->restriction[index][i + 2 * (2 * j + k)](
1611 ((2 * i + 1) * deg + m) * this->degree + l +
1612 n_edge_dofs,
1613 dof) =
1614 solution(l * deg + m, 2 * (j + 2 * k) + 1);
1615
1616 if (std::abs(solution(l * deg + m,
1617 2 * (j + 2 * (k + 2)))) >
1618 1e-14)
1619 this->restriction[index][2 * (i + 2 * j) + k](
1620 (2 * (i + 2) * this->degree + l) * deg + m +
1621 n_edge_dofs,
1622 dof) =
1623 solution(l * deg + m, 2 * (j + 2 * (k + 2)));
1624
1625 if (std::abs(solution(l * deg + m,
1626 2 * (j + 2 * k) + 9)) >
1627 1e-14)
1628 this->restriction[index][2 * (i + 2 * j) + k](
1629 ((2 * i + 5) * deg + m) * this->degree + l +
1630 n_edge_dofs,
1631 dof) =
1632 solution(l * deg + m, 2 * (j + 2 * k) + 9);
1633
1634 if (std::abs(solution(l * deg + m,
1635 2 * (j + 2 * (k + 4)))) >
1636 1e-14)
1637 this->restriction[index][2 * (2 * i + j) + k](
1638 (2 * (i + 4) * this->degree + l) * deg + m +
1639 n_edge_dofs,
1640 dof) =
1641 solution(l * deg + m, 2 * (j + 2 * (k + 4)));
1642
1643 if (std::abs(solution(l * deg + m,
1644 2 * (j + 2 * k) + 17)) >
1645 1e-14)
1646 this->restriction[index][2 * (2 * i + j) + k](
1647 ((2 * i + 9) * deg + m) * this->degree + l +
1648 n_edge_dofs,
1649 dof) =
1650 solution(l * deg + m, 2 * (j + 2 * k) + 17);
1651 }
1652 }
1653
1654 const QGauss<dim> quadrature(2 * this->degree);
1655 const std::vector<Point<dim>> &quadrature_points =
1656 quadrature.get_points();
1657 const unsigned int n_boundary_dofs =
1658 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree +
1659 n_edge_dofs;
1660 const unsigned int n_quadrature_points = quadrature.size();
1661
1662 {
1663 FullMatrix<double> assembling_matrix(deg * deg * this->degree,
1664 n_quadrature_points);
1665
1666 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1667 ++q_point)
1668 {
1669 const double weight = std::sqrt(quadrature.weight(q_point));
1670
1671 for (unsigned int i = 0; i <= deg; ++i)
1672 {
1673 const double L_i =
1674 weight * legendre_polynomials[i].value(
1675 quadrature_points[q_point][0]);
1676
1677 for (unsigned int j = 0; j < deg; ++j)
1678 {
1679 const double l_j =
1680 L_i * lobatto_polynomials[j + 2].value(
1681 quadrature_points[q_point][1]);
1682
1683 for (unsigned int k = 0; k < deg; ++k)
1684 assembling_matrix((i * deg + j) * deg + k,
1685 q_point) =
1686 l_j * lobatto_polynomials[k + 2].value(
1687 quadrature_points[q_point][2]);
1688 }
1689 }
1690 }
1691
1692 FullMatrix<double> system_matrix(assembling_matrix.m(),
1693 assembling_matrix.m());
1694
1695 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1696 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
1697 system_matrix_inv.invert(system_matrix);
1698 }
1699
1700 solution.reinit(system_matrix_inv.m(), 24);
1701 system_rhs.reinit(system_matrix_inv.m(), 24);
1702 tmp.reinit(24);
1703
1704 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1705 {
1706 system_rhs = 0.0;
1707
1708 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1709 ++q_point)
1710 {
1711 tmp = 0.0;
1712
1713 if (quadrature_points[q_point][0] < 0.5)
1714 {
1715 if (quadrature_points[q_point][1] < 0.5)
1716 {
1717 if (quadrature_points[q_point][2] < 0.5)
1718 {
1719 const Point<dim> quadrature_point(
1720 2.0 * quadrature_points[q_point][0],
1721 2.0 * quadrature_points[q_point][1],
1722 2.0 * quadrature_points[q_point][2]);
1723
1724 tmp(0) += 2.0 * this->shape_value_component(
1725 dof, quadrature_point, 0);
1726 tmp(1) += 2.0 * this->shape_value_component(
1727 dof, quadrature_point, 1);
1728 tmp(2) += 2.0 * this->shape_value_component(
1729 dof, quadrature_point, 2);
1730 }
1731
1732 else
1733 {
1734 const Point<dim> quadrature_point(
1735 2.0 * quadrature_points[q_point][0],
1736 2.0 * quadrature_points[q_point][1],
1737 2.0 * quadrature_points[q_point][2] - 1.0);
1738
1739 tmp(3) += 2.0 * this->shape_value_component(
1740 dof, quadrature_point, 0);
1741 tmp(4) += 2.0 * this->shape_value_component(
1742 dof, quadrature_point, 1);
1743 tmp(5) += 2.0 * this->shape_value_component(
1744 dof, quadrature_point, 2);
1745 }
1746 }
1747
1748 else if (quadrature_points[q_point][2] < 0.5)
1749 {
1750 const Point<dim> quadrature_point(
1751 2.0 * quadrature_points[q_point][0],
1752 2.0 * quadrature_points[q_point][1] - 1.0,
1753 2.0 * quadrature_points[q_point][2]);
1754
1755 tmp(6) += 2.0 * this->shape_value_component(
1756 dof, quadrature_point, 0);
1757 tmp(7) += 2.0 * this->shape_value_component(
1758 dof, quadrature_point, 1);
1759 tmp(8) += 2.0 * this->shape_value_component(
1760 dof, quadrature_point, 2);
1761 }
1762
1763 else
1764 {
1765 const Point<dim> quadrature_point(
1766 2.0 * quadrature_points[q_point][0],
1767 2.0 * quadrature_points[q_point][1] - 1.0,
1768 2.0 * quadrature_points[q_point][2] - 1.0);
1769
1770 tmp(9) += 2.0 * this->shape_value_component(
1771 dof, quadrature_point, 0);
1772 tmp(10) += 2.0 * this->shape_value_component(
1773 dof, quadrature_point, 1);
1774 tmp(11) += 2.0 * this->shape_value_component(
1775 dof, quadrature_point, 2);
1776 }
1777 }
1778
1779 else if (quadrature_points[q_point][1] < 0.5)
1780 {
1781 if (quadrature_points[q_point][2] < 0.5)
1782 {
1783 const Point<dim> quadrature_point(
1784 2.0 * quadrature_points[q_point][0] - 1.0,
1785 2.0 * quadrature_points[q_point][1],
1786 2.0 * quadrature_points[q_point][2]);
1787
1788 tmp(12) += 2.0 * this->shape_value_component(
1789 dof, quadrature_point, 0);
1790 tmp(13) += 2.0 * this->shape_value_component(
1791 dof, quadrature_point, 1);
1792 tmp(14) += 2.0 * this->shape_value_component(
1793 dof, quadrature_point, 2);
1794 }
1795
1796 else
1797 {
1798 const Point<dim> quadrature_point(
1799 2.0 * quadrature_points[q_point][0] - 1.0,
1800 2.0 * quadrature_points[q_point][1],
1801 2.0 * quadrature_points[q_point][2] - 1.0);
1802
1803 tmp(15) += 2.0 * this->shape_value_component(
1804 dof, quadrature_point, 0);
1805 tmp(16) += 2.0 * this->shape_value_component(
1806 dof, quadrature_point, 1);
1807 tmp(17) += 2.0 * this->shape_value_component(
1808 dof, quadrature_point, 2);
1809 }
1810 }
1811
1812 else if (quadrature_points[q_point][2] < 0.5)
1813 {
1814 const Point<dim> quadrature_point(
1815 2.0 * quadrature_points[q_point][0] - 1.0,
1816 2.0 * quadrature_points[q_point][1] - 1.0,
1817 2.0 * quadrature_points[q_point][2]);
1818
1819 tmp(18) +=
1820 2.0 * this->shape_value_component(dof,
1821 quadrature_point,
1822 0);
1823 tmp(19) +=
1824 2.0 * this->shape_value_component(dof,
1825 quadrature_point,
1826 1);
1827 tmp(20) +=
1828 2.0 * this->shape_value_component(dof,
1829 quadrature_point,
1830 2);
1831 }
1832
1833 else
1834 {
1835 const Point<dim> quadrature_point(
1836 2.0 * quadrature_points[q_point][0] - 1.0,
1837 2.0 * quadrature_points[q_point][1] - 1.0,
1838 2.0 * quadrature_points[q_point][2] - 1.0);
1839
1840 tmp(21) +=
1841 2.0 * this->shape_value_component(dof,
1842 quadrature_point,
1843 0);
1844 tmp(22) +=
1845 2.0 * this->shape_value_component(dof,
1846 quadrature_point,
1847 1);
1848 tmp(23) +=
1849 2.0 * this->shape_value_component(dof,
1850 quadrature_point,
1851 2);
1852 }
1853
1854 for (unsigned int i = 0; i < 2; ++i)
1855 for (unsigned int j = 0; j < 2; ++j)
1856 for (unsigned int k = 0; k < 2; ++k)
1857 for (unsigned int l = 0; l <= deg; ++l)
1858 {
1859 tmp(3 * (i + 2 * (j + 2 * k))) -=
1860 this->restriction[index][2 * (2 * i + j) + k](
1861 (4 * i + j + 2) * this->degree + l, dof) *
1862 this->shape_value_component(
1863 (4 * i + j + 2) * this->degree + l,
1864 quadrature_points[q_point],
1865 0);
1866 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1867 this->restriction[index][2 * (2 * i + j) + k](
1868 (4 * i + k) * this->degree + l, dof) *
1869 this->shape_value_component(
1870 (4 * i + k) * this->degree + l,
1871 quadrature_points[q_point],
1872 1);
1873 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1874 this->restriction[index][2 * (2 * i + j) + k](
1875 (2 * (j + 4) + k) * this->degree + l, dof) *
1876 this->shape_value_component(
1877 (2 * (j + 4) + k) * this->degree + l,
1878 quadrature_points[q_point],
1879 2);
1880
1881 for (unsigned int m = 0; m < deg; ++m)
1882 {
1883 tmp(3 * (i + 2 * (j + 2 * k))) -=
1884 this->restriction[index][2 * (2 * i + j) +
1885 k](
1886 ((2 * j + 5) * deg + m) * this->degree +
1887 l + n_edge_dofs,
1888 dof) *
1889 this->shape_value_component(
1890 ((2 * j + 5) * deg + m) * this->degree +
1891 l + n_edge_dofs,
1892 quadrature_points[q_point],
1893 0);
1894 tmp(3 * (i + 2 * (j + 2 * k))) -=
1895 this->restriction[index][2 * (2 * i + j) +
1896 k](
1897 (2 * (i + 4) * this->degree + l) * deg +
1898 m + n_edge_dofs,
1899 dof) *
1900 this->shape_value_component(
1901 (2 * (i + 4) * this->degree + l) * deg +
1902 m + n_edge_dofs,
1903 quadrature_points[q_point],
1904 0);
1905 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1906 this->restriction[index][2 * (2 * i + j) +
1907 k](
1908 (2 * k * this->degree + l) * deg + m +
1909 n_edge_dofs,
1910 dof) *
1911 this->shape_value_component(
1912 (2 * k * this->degree + l) * deg + m +
1913 n_edge_dofs,
1914 quadrature_points[q_point],
1915 1);
1916 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1917 this->restriction[index][2 * (2 * i + j) +
1918 k](
1919 ((2 * i + 9) * deg + m) * this->degree +
1920 l + n_edge_dofs,
1921 dof) *
1922 this->shape_value_component(
1923 ((2 * i + 9) * deg + m) * this->degree +
1924 l + n_edge_dofs,
1925 quadrature_points[q_point],
1926 1);
1927 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1928 this->restriction[index][2 * (2 * i + j) +
1929 k](
1930 ((2 * k + 1) * deg + m) * this->degree +
1931 l + n_edge_dofs,
1932 dof) *
1933 this->shape_value_component(
1934 ((2 * k + 1) * deg + m) * this->degree +
1935 l + n_edge_dofs,
1936 quadrature_points[q_point],
1937 2);
1938 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1939 this->restriction[index][2 * (2 * i + j) +
1940 k](
1941 (2 * (j + 2) * this->degree + l) * deg +
1942 m + n_edge_dofs,
1943 dof) *
1944 this->shape_value_component(
1945 (2 * (j + 2) * this->degree + l) * deg +
1946 m + n_edge_dofs,
1947 quadrature_points[q_point],
1948 2);
1949 }
1950 }
1951
1952 tmp *= quadrature.weight(q_point);
1953
1954 for (unsigned int i = 0; i <= deg; ++i)
1955 {
1956 const double L_i_0 = legendre_polynomials[i].value(
1957 quadrature_points[q_point][0]);
1958 const double L_i_1 = legendre_polynomials[i].value(
1959 quadrature_points[q_point][1]);
1960 const double L_i_2 = legendre_polynomials[i].value(
1961 quadrature_points[q_point][2]);
1962
1963 for (unsigned int j = 0; j < deg; ++j)
1964 {
1965 const double l_j_0 =
1966 L_i_0 * lobatto_polynomials[j + 2].value(
1967 quadrature_points[q_point][1]);
1968 const double l_j_1 =
1969 L_i_1 * lobatto_polynomials[j + 2].value(
1970 quadrature_points[q_point][0]);
1971 const double l_j_2 =
1972 L_i_2 * lobatto_polynomials[j + 2].value(
1973 quadrature_points[q_point][0]);
1974
1975 for (unsigned int k = 0; k < deg; ++k)
1976 {
1977 const double l_k_0 =
1978 l_j_0 * lobatto_polynomials[k + 2].value(
1979 quadrature_points[q_point][2]);
1980 const double l_k_1 =
1981 l_j_1 * lobatto_polynomials[k + 2].value(
1982 quadrature_points[q_point][2]);
1983 const double l_k_2 =
1984 l_j_2 * lobatto_polynomials[k + 2].value(
1985 quadrature_points[q_point][1]);
1986
1987 for (unsigned int l = 0; l < 8; ++l)
1988 {
1989 system_rhs((i * deg + j) * deg + k,
1990 3 * l) += tmp(3 * l) * l_k_0;
1991 system_rhs((i * deg + j) * deg + k,
1992 3 * l + 1) +=
1993 tmp(3 * l + 1) * l_k_1;
1994 system_rhs((i * deg + j) * deg + k,
1995 3 * l + 2) +=
1996 tmp(3 * l + 2) * l_k_2;
1997 }
1998 }
1999 }
2000 }
2001 }
2002
2003 system_matrix_inv.mmult(solution, system_rhs);
2004
2005 for (unsigned int i = 0; i < 2; ++i)
2006 for (unsigned int j = 0; j < 2; ++j)
2007 for (unsigned int k = 0; k < 2; ++k)
2008 for (unsigned int l = 0; l <= deg; ++l)
2009 for (unsigned int m = 0; m < deg; ++m)
2010 for (unsigned int n = 0; n < deg; ++n)
2011 {
2012 if (std::abs(
2013 solution((l * deg + m) * deg + n,
2014 3 * (i + 2 * (j + 2 * k)))) >
2015 1e-14)
2016 this->restriction[index][2 * (2 * i + j) + k](
2017 (l * deg + m) * deg + n + n_boundary_dofs,
2018 dof) = solution((l * deg + m) * deg + n,
2019 3 * (i + 2 * (j + 2 * k)));
2020
2021 if (std::abs(
2022 solution((l * deg + m) * deg + n,
2023 3 * (i + 2 * (j + 2 * k)) + 1)) >
2024 1e-14)
2025 this->restriction[index][2 * (2 * i + j) + k](
2026 (l + (m + deg) * this->degree) * deg + n +
2027 n_boundary_dofs,
2028 dof) =
2029 solution((l * deg + m) * deg + n,
2030 3 * (i + 2 * (j + 2 * k)) + 1);
2031
2032 if (std::abs(
2033 solution((l * deg + m) * deg + n,
2034 3 * (i + 2 * (j + 2 * k)) + 2)) >
2035 1e-14)
2036 this->restriction[index][2 * (2 * i + j) + k](
2037 l +
2038 ((m + 2 * deg) * deg + n) * this->degree +
2039 n_boundary_dofs,
2040 dof) =
2041 solution((l * deg + m) * deg + n,
2042 3 * (i + 2 * (j + 2 * k)) + 2);
2043 }
2044 }
2045 }
2046
2047 break;
2048 }
2049
2050 default:
2052 }
2053}
2054
2055
2056
2057template <int dim>
2058std::vector<unsigned int>
2059FE_Nedelec<dim>::get_dpo_vector(const unsigned int degree, bool dg)
2060{
2061 std::vector<unsigned int> dpo;
2062
2063 if (dg)
2064 {
2065 dpo.resize(dim + 1);
2066 dpo[dim] = PolynomialsNedelec<dim>::n_polynomials(degree);
2067 }
2068 else
2069 {
2070 dpo.push_back(0);
2071 dpo.push_back(degree + 1);
2072 if (dim > 1)
2073 dpo.push_back(2 * degree * (degree + 1));
2074 if (dim > 2)
2075 dpo.push_back(3 * degree * degree * (degree + 1));
2076 }
2077
2078 return dpo;
2079}
2080
2081//---------------------------------------------------------------------------
2082// Data field initialization
2083//---------------------------------------------------------------------------
2084
2085// Check whether a given shape
2086// function has support on a
2087// given face.
2088
2089// We just switch through the
2090// faces of the cell and return
2091// true, if the shape function
2092// has support on the face
2093// and false otherwise.
2094template <int dim>
2095bool
2096FE_Nedelec<dim>::has_support_on_face(const unsigned int shape_index,
2097 const unsigned int face_index) const
2098{
2099 AssertIndexRange(shape_index, this->n_dofs_per_cell());
2101
2102 const unsigned int deg = this->degree - 1;
2103 switch (dim)
2104 {
2105 case 2:
2106 switch (face_index)
2107 {
2108 case 0:
2109 if (!((shape_index > deg) && (shape_index < 2 * this->degree)))
2110 return true;
2111
2112 else
2113 return false;
2114
2115 case 1:
2116 if ((shape_index > deg) &&
2117 (shape_index <
2118 GeometryInfo<2>::lines_per_cell * this->degree))
2119 return true;
2120
2121 else
2122 return false;
2123
2124 case 2:
2125 if (shape_index < 3 * this->degree)
2126 return true;
2127
2128 else
2129 return false;
2130
2131 case 3:
2132 if (!((shape_index >= 2 * this->degree) &&
2133 (shape_index < 3 * this->degree)))
2134 return true;
2135
2136 else
2137 return false;
2138
2139 default:
2140 {
2142 return false;
2143 }
2144 }
2145
2146 case 3:
2147 switch (face_index)
2148 {
2149 case 0:
2150 if (((shape_index > deg) && (shape_index < 2 * this->degree)) ||
2151 ((shape_index >= 5 * this->degree) &&
2152 (shape_index < 6 * this->degree)) ||
2153 ((shape_index >= 9 * this->degree) &&
2154 (shape_index < 10 * this->degree)) ||
2155 ((shape_index >= 11 * this->degree) &&
2156 (shape_index <
2157 GeometryInfo<3>::lines_per_cell * this->degree)) ||
2158 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2159 this->degree) &&
2160 (shape_index < (GeometryInfo<3>::lines_per_cell + 4 * deg) *
2161 this->degree)) ||
2162 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2163 this->degree) &&
2164 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2165 this->degree)) ||
2166 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2167 this->degree) &&
2168 (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2169 this->degree)) ||
2170 ((shape_index >=
2171 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2172 this->degree) &&
2173 (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2174 this->degree)))
2175 return false;
2176
2177 else
2178 return true;
2179
2180 case 1:
2181 if (((shape_index > deg) && (shape_index < 4 * this->degree)) ||
2182 ((shape_index >= 5 * this->degree) &&
2183 (shape_index < 8 * this->degree)) ||
2184 ((shape_index >= 9 * this->degree) &&
2185 (shape_index < 10 * this->degree)) ||
2186 ((shape_index >= 11 * this->degree) &&
2187 (shape_index <
2188 GeometryInfo<3>::lines_per_cell * this->degree)) ||
2189 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2190 this->degree) &&
2191 (shape_index < (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2192 this->degree)) ||
2193 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2194 this->degree) &&
2195 (shape_index < (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2196 this->degree)) ||
2197 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2198 this->degree) &&
2199 (shape_index < (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2200 this->degree)) ||
2201 ((shape_index >=
2202 (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2203 this->degree) &&
2204 (shape_index < (GeometryInfo<3>::lines_per_cell + 12 * deg) *
2205 this->degree)))
2206 return true;
2207
2208 else
2209 return false;
2210
2211 case 2:
2212 if ((shape_index < 3 * this->degree) ||
2213 ((shape_index >= 4 * this->degree) &&
2214 (shape_index < 7 * this->degree)) ||
2215 ((shape_index >= 8 * this->degree) &&
2216 (shape_index < 10 * this->degree)) ||
2217 ((shape_index >=
2218 (GeometryInfo<3>::lines_per_cell + deg) * this->degree) &&
2219 (shape_index < (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2220 this->degree)) ||
2221 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2222 this->degree) &&
2223 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2224 this->degree)) ||
2225 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 8 * deg) *
2226 this->degree) &&
2227 (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2228 this->degree)) ||
2229 ((shape_index >=
2230 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2231 this->degree) &&
2232 (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2233 this->degree)))
2234 return true;
2235
2236 else
2237 return false;
2238
2239 case 3:
2240 if ((shape_index < 2 * this->degree) ||
2241 ((shape_index >= 3 * this->degree) &&
2242 (shape_index < 6 * this->degree)) ||
2243 ((shape_index >= 7 * this->degree) &&
2244 (shape_index < 8 * this->degree)) ||
2245 ((shape_index >= 10 * this->degree) &&
2246 (shape_index <
2247 GeometryInfo<3>::lines_per_cell * this->degree)) ||
2248 ((shape_index >=
2249 (GeometryInfo<3>::lines_per_cell + deg) * this->degree) &&
2250 (shape_index < (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2251 this->degree)) ||
2252 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2253 this->degree) &&
2254 (shape_index < (GeometryInfo<3>::lines_per_cell + 4 * deg) *
2255 this->degree)) ||
2256 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2257 this->degree) &&
2258 (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2259 this->degree)) ||
2260 ((shape_index >=
2261 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2262 this->degree) &&
2263 (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2264 this->degree)))
2265 return true;
2266
2267 else
2268 return false;
2269
2270 case 4:
2271 if ((shape_index < 4 * this->degree) ||
2272 ((shape_index >= 8 * this->degree) &&
2273 (shape_index <
2274 (GeometryInfo<3>::lines_per_cell + deg) * this->degree)) ||
2275 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2276 this->degree) &&
2277 (shape_index < (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2278 this->degree)) ||
2279 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2280 this->degree) &&
2281 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2282 this->degree)) ||
2283 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2284 this->degree) &&
2285 (shape_index < (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2286 this->degree)))
2287 return true;
2288
2289 else
2290 return false;
2291
2292 case 5:
2293 if (((shape_index >= 4 * this->degree) &&
2294 (shape_index <
2295 (GeometryInfo<3>::lines_per_cell + deg) * this->degree)) ||
2296 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2297 this->degree) &&
2298 (shape_index < (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2299 this->degree)) ||
2300 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2301 this->degree) &&
2302 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2303 this->degree)) ||
2304 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2305 this->degree) &&
2306 (shape_index < (GeometryInfo<3>::lines_per_cell + 8 * deg) *
2307 this->degree)) ||
2308 ((shape_index >=
2309 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2310 this->degree) &&
2311 (shape_index < (GeometryInfo<3>::lines_per_cell + 12 * deg) *
2312 this->degree)))
2313 return true;
2314
2315 else
2316 return false;
2317
2318 default:
2319 {
2321 return false;
2322 }
2323 }
2324
2325 default:
2326 {
2328 return false;
2329 }
2330 }
2331}
2332
2333template <int dim>
2336 const unsigned int codim) const
2337{
2338 Assert(codim <= dim, ExcImpossibleInDim(dim));
2339 (void)codim;
2340
2341 // vertex/line/face/cell domination
2342 // --------------------------------
2343 if (const FE_Nedelec<dim> *fe_nedelec_other =
2344 dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2345 {
2346 if (this->degree < fe_nedelec_other->degree)
2348 else if (this->degree == fe_nedelec_other->degree)
2350 else
2352 }
2353 else if (const FE_Nothing<dim> *fe_nothing =
2354 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
2355 {
2356 if (fe_nothing->is_dominating())
2358 else
2359 // the FE_Nothing has no degrees of freedom and it is typically used
2360 // in a context where we don't require any continuity along the
2361 // interface
2363 }
2364
2367}
2368
2369template <int dim>
2370bool
2372{
2373 return true;
2374}
2375
2376template <int dim>
2377std::vector<std::pair<unsigned int, unsigned int>>
2379{
2380 // Nedelec elements do not have any dofs
2381 // on vertices, hence return an empty vector.
2382 return std::vector<std::pair<unsigned int, unsigned int>>();
2383}
2384
2385template <int dim>
2386std::vector<std::pair<unsigned int, unsigned int>>
2388 const FiniteElement<dim> &fe_other) const
2389{
2390 // we can presently only compute these
2391 // identities if both FEs are
2392 // FE_Nedelec or if the other one is an
2393 // FE_Nothing
2394 if (const FE_Nedelec<dim> *fe_nedelec_other =
2395 dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2396 {
2397 // dofs are located on lines, so
2398 // two dofs are identical, if their
2399 // edge shape functions have the
2400 // same polynomial degree.
2401 std::vector<std::pair<unsigned int, unsigned int>> identities;
2402
2403 for (unsigned int i = 0;
2404 i < std::min(fe_nedelec_other->degree, this->degree);
2405 ++i)
2406 identities.emplace_back(i, i);
2407
2408 return identities;
2409 }
2410
2411 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2412 {
2413 // the FE_Nothing has no
2414 // degrees of freedom, so there
2415 // are no equivalencies to be
2416 // recorded
2417 return std::vector<std::pair<unsigned int, unsigned int>>();
2418 }
2419
2420 else
2421 {
2423 return std::vector<std::pair<unsigned int, unsigned int>>();
2424 }
2425}
2426
2427template <int dim>
2428std::vector<std::pair<unsigned int, unsigned int>>
2430 const unsigned int) const
2431{
2432 // we can presently only compute
2433 // these identities if both FEs are
2434 // FE_Nedelec or if the other one is an
2435 // FE_Nothing
2436 if (const FE_Nedelec<dim> *fe_nedelec_other =
2437 dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2438 {
2439 // dofs are located on the interior
2440 // of faces, so two dofs are identical,
2441 // if their face shape functions have
2442 // the same polynomial degree.
2443 const unsigned int p = fe_nedelec_other->degree;
2444 const unsigned int q = this->degree;
2445 const unsigned int p_min = std::min(p, q);
2446 std::vector<std::pair<unsigned int, unsigned int>> identities;
2447
2448 for (unsigned int i = 0; i < p_min; ++i)
2449 for (unsigned int j = 0; j < p_min - 1; ++j)
2450 {
2451 identities.emplace_back(i * (q - 1) + j, i * (p - 1) + j);
2452 identities.emplace_back(i + (j + q - 1) * q, i + (j + p - 1) * p);
2453 }
2454
2455 return identities;
2456 }
2457
2458 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2459 {
2460 // the FE_Nothing has no
2461 // degrees of freedom, so there
2462 // are no equivalencies to be
2463 // recorded
2464 return std::vector<std::pair<unsigned int, unsigned int>>();
2465 }
2466
2467 else
2468 {
2470 return std::vector<std::pair<unsigned int, unsigned int>>();
2471 }
2472}
2473
2474// In this function we compute the face
2475// interpolation matrix. This is usually
2476// done by projection-based interpolation,
2477// but, since one can compute the entries
2478// easy per hand, we save some computation
2479// time at this point and just fill in the
2480// correct values.
2481template <int dim>
2482void
2484 const FiniteElement<dim> &source,
2485 FullMatrix<double> &interpolation_matrix,
2486 const unsigned int face_no) const
2487{
2488 (void)face_no;
2489 // this is only implemented, if the
2490 // source FE is also a
2491 // Nedelec element
2492 AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2493 (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2495 Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
2496 ExcDimensionMismatch(interpolation_matrix.m(),
2497 source.n_dofs_per_face(face_no)));
2498 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
2499 ExcDimensionMismatch(interpolation_matrix.n(),
2500 this->n_dofs_per_face(face_no)));
2501
2502 // ok, source is a Nedelec element, so
2503 // we will be able to do the work
2504 const FE_Nedelec<dim> &source_fe =
2505 dynamic_cast<const FE_Nedelec<dim> &>(source);
2506
2507 // Make sure, that the element,
2508 // for which the DoFs should be
2509 // constrained is the one with
2510 // the higher polynomial degree.
2511 // Actually the procedure will work
2512 // also if this assertion is not
2513 // satisfied. But the matrices
2514 // produced in that case might
2515 // lead to problems in the
2516 // hp-procedures, which use this
2517 // method.
2518 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
2520 interpolation_matrix = 0;
2521
2522 // On lines we can just identify
2523 // all degrees of freedom.
2524 for (unsigned int i = 0; i < this->degree; ++i)
2525 interpolation_matrix(i, i) = 1.0;
2526
2527 // In 3d we have some lines more
2528 // and a face. The procedure stays
2529 // the same as above, but we have
2530 // to take a bit more care of the
2531 // indices of the degrees of
2532 // freedom.
2533 if (dim == 3)
2534 {
2535 const unsigned int p = source_fe.degree;
2536 const unsigned int q = this->degree;
2537
2538 for (unsigned int i = 0; i < q; ++i)
2539 {
2540 for (unsigned int j = 1; j < GeometryInfo<dim>::lines_per_face; ++j)
2541 interpolation_matrix(j * p + i, j * q + i) = 1.0;
2542
2543 for (unsigned int j = 0; j < q - 1; ++j)
2544 {
2545 interpolation_matrix(GeometryInfo<dim>::lines_per_face * p +
2546 i * (p - 1) + j,
2548 i * (q - 1) + j) = 1.0;
2549 interpolation_matrix(GeometryInfo<dim>::lines_per_face * p + i +
2550 (j + p - 1) * p,
2552 (j + q - 1) * q) = 1.0;
2553 }
2554 }
2555 }
2556}
2557
2558
2559
2560template <>
2561void
2563 const unsigned int,
2565 const unsigned int) const
2566{
2568}
2569
2570
2571
2572// In this function we compute the
2573// subface interpolation matrix.
2574// This is done by a projection-
2575// based interpolation. Therefore
2576// we first interpolate the
2577// shape functions of the higher
2578// order element on the lowest
2579// order edge shape functions.
2580// Then the remaining part of
2581// the interpolated shape
2582// functions is projected on the
2583// higher order edge shape
2584// functions, the face shape
2585// functions and the interior
2586// shape functions (if they all
2587// exist).
2588template <int dim>
2589void
2591 const FiniteElement<dim> &source,
2592 const unsigned int subface,
2593 FullMatrix<double> &interpolation_matrix,
2594 const unsigned int face_no) const
2595{
2596 // this is only implemented, if the
2597 // source FE is also a
2598 // Nedelec element
2599 AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2600 (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2602 Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
2603 ExcDimensionMismatch(interpolation_matrix.m(),
2604 source.n_dofs_per_face(face_no)));
2605 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
2606 ExcDimensionMismatch(interpolation_matrix.n(),
2607 this->n_dofs_per_face(face_no)));
2608
2609 // ok, source is a Nedelec element, so
2610 // we will be able to do the work
2611 const FE_Nedelec<dim> &source_fe =
2612 dynamic_cast<const FE_Nedelec<dim> &>(source);
2613
2614 // Make sure, that the element,
2615 // for which the DoFs should be
2616 // constrained is the one with
2617 // the higher polynomial degree.
2618 // Actually the procedure will work
2619 // also if this assertion is not
2620 // satisfied. But the matrices
2621 // produced in that case might
2622 // lead to problems in the
2623 // hp-procedures, which use this
2624 // method.
2625 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
2627 interpolation_matrix = 0.0;
2628 // Perform projection-based interpolation
2629 // as usual.
2630 const QGauss<1> edge_quadrature(source_fe.degree);
2631 const std::vector<Point<1>> &edge_quadrature_points =
2632 edge_quadrature.get_points();
2633 const unsigned int n_edge_quadrature_points = edge_quadrature.size();
2634
2635 switch (dim)
2636 {
2637 case 2:
2638 {
2639 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
2640 ++dof)
2641 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2642 ++q_point)
2643 {
2644 const Point<dim> quadrature_point(
2645 0.0, 0.5 * (edge_quadrature_points[q_point][0] + subface));
2646
2647 interpolation_matrix(0, dof) +=
2648 0.5 * edge_quadrature.weight(q_point) *
2649 this->shape_value_component(dof, quadrature_point, 1);
2650 }
2651
2652 if (source_fe.degree > 1)
2653 {
2654 const std::vector<Polynomials::Polynomial<double>>
2655 &legendre_polynomials =
2657 source_fe.degree - 1);
2658 FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
2659 source_fe.degree - 1);
2660
2661 {
2662 FullMatrix<double> assembling_matrix(source_fe.degree - 1,
2663 n_edge_quadrature_points);
2664
2665 for (unsigned int q_point = 0;
2666 q_point < n_edge_quadrature_points;
2667 ++q_point)
2668 {
2669 const double weight =
2670 std::sqrt(edge_quadrature.weight(q_point));
2671
2672 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2673 assembling_matrix(i, q_point) =
2674 weight * legendre_polynomials[i + 1].value(
2675 edge_quadrature_points[q_point][0]);
2676 }
2677
2678 FullMatrix<double> system_matrix(source_fe.degree - 1,
2679 source_fe.degree - 1);
2680
2681 assembling_matrix.mTmult(system_matrix, assembling_matrix);
2682 system_matrix_inv.invert(system_matrix);
2683 }
2684
2685 Vector<double> solution(source_fe.degree - 1);
2686 Vector<double> system_rhs(source_fe.degree - 1);
2687
2688 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
2689 ++dof)
2690 {
2691 system_rhs = 0.0;
2692
2693 for (unsigned int q_point = 0;
2694 q_point < n_edge_quadrature_points;
2695 ++q_point)
2696 {
2697 const Point<dim> quadrature_point_0(
2698 0.0,
2699 0.5 * (edge_quadrature_points[q_point][0] + subface));
2700 const Point<dim> quadrature_point_1(
2701 0.0, edge_quadrature_points[q_point][0]);
2702 const double tmp =
2703 edge_quadrature.weight(q_point) *
2704 (0.5 * this->shape_value_component(dof,
2705 quadrature_point_0,
2706 1) -
2707 interpolation_matrix(0, dof) *
2708 source_fe.shape_value_component(0,
2709 quadrature_point_1,
2710 1));
2711
2712 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2713 system_rhs(i) +=
2714 tmp * legendre_polynomials[i + 1].value(
2715 edge_quadrature_points[q_point][0]);
2716 }
2717
2718 system_matrix_inv.vmult(solution, system_rhs);
2719
2720 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2721 if (std::abs(solution(i)) > 1e-14)
2722 interpolation_matrix(i + 1, dof) = solution(i);
2723 }
2724 }
2725
2726 break;
2727 }
2728
2729 case 3:
2730 {
2731 const double shifts[4][2] = {{0.0, 0.0},
2732 {1.0, 0.0},
2733 {0.0, 1.0},
2734 {1.0, 1.0}};
2735
2736 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
2737 ++dof)
2738 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2739 ++q_point)
2740 {
2741 const double weight = 0.5 * edge_quadrature.weight(q_point);
2742
2743 for (unsigned int i = 0; i < 2; ++i)
2744 {
2745 Point<dim> quadrature_point(
2746 0.5 * (i + shifts[subface][0]),
2747 0.5 * (edge_quadrature_points[q_point][0] +
2748 shifts[subface][1]),
2749 0.0);
2750
2751 interpolation_matrix(i * source_fe.degree, dof) +=
2752 weight *
2753 this->shape_value_component(
2754 this->face_to_cell_index(dof, 4), quadrature_point, 1);
2755 quadrature_point =
2756 Point<dim>(0.5 * (edge_quadrature_points[q_point][0] +
2757 shifts[subface][0]),
2758 0.5 * (i + shifts[subface][1]),
2759 0.0);
2760 interpolation_matrix((i + 2) * source_fe.degree, dof) +=
2761 weight *
2762 this->shape_value_component(
2763 this->face_to_cell_index(dof, 4), quadrature_point, 0);
2764 }
2765 }
2766
2767 if (source_fe.degree > 1)
2768 {
2769 const std::vector<Polynomials::Polynomial<double>>
2770 &legendre_polynomials =
2772 source_fe.degree - 1);
2773 FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
2774 source_fe.degree - 1);
2775
2776 {
2777 FullMatrix<double> assembling_matrix(source_fe.degree - 1,
2778 n_edge_quadrature_points);
2779
2780 for (unsigned int q_point = 0;
2781 q_point < n_edge_quadrature_points;
2782 ++q_point)
2783 {
2784 const double weight =
2785 std::sqrt(edge_quadrature.weight(q_point));
2786
2787 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2788 assembling_matrix(i, q_point) =
2789 weight * legendre_polynomials[i + 1].value(
2790 edge_quadrature_points[q_point][0]);
2791 }
2792
2793 FullMatrix<double> system_matrix(source_fe.degree - 1,
2794 source_fe.degree - 1);
2795
2796 assembling_matrix.mTmult(system_matrix, assembling_matrix);
2797 system_matrix_inv.invert(system_matrix);
2798 }
2799
2800 FullMatrix<double> solution(source_fe.degree - 1,
2802 FullMatrix<double> system_rhs(source_fe.degree - 1,
2805
2806 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
2807 ++dof)
2808 {
2809 system_rhs = 0.0;
2810
2811 for (unsigned int q_point = 0;
2812 q_point < n_edge_quadrature_points;
2813 ++q_point)
2814 {
2815 const double weight = edge_quadrature.weight(q_point);
2816
2817 for (unsigned int i = 0; i < 2; ++i)
2818 {
2819 Point<dim> quadrature_point_0(
2820 0.5 * (i + shifts[subface][0]),
2821 0.5 * (edge_quadrature_points[q_point][0] +
2822 shifts[subface][1]),
2823 0.0);
2824 Point<dim> quadrature_point_1(
2825 i, edge_quadrature_points[q_point][0], 0.0);
2826
2827 tmp(i) =
2828 weight *
2829 (0.5 * this->shape_value_component(
2830 this->face_to_cell_index(dof, 4),
2831 quadrature_point_0,
2832 1) -
2833 interpolation_matrix(i * source_fe.degree, dof) *
2834 source_fe.shape_value_component(
2835 i * source_fe.degree, quadrature_point_1, 1));
2836 quadrature_point_0 =
2837 Point<dim>(0.5 *
2838 (edge_quadrature_points[q_point][0] +
2839 shifts[subface][0]),
2840 0.5 * (i + shifts[subface][1]),
2841 0.0);
2842 quadrature_point_1 =
2843 Point<dim>(edge_quadrature_points[q_point][0],
2844 i,
2845 0.0);
2846 tmp(i + 2) =
2847 weight *
2848 (0.5 * this->shape_value_component(
2849 this->face_to_cell_index(dof, 4),
2850 quadrature_point_0,
2851 0) -
2852 interpolation_matrix((i + 2) * source_fe.degree,
2853 dof) *
2854 source_fe.shape_value_component(
2855 (i + 2) * source_fe.degree,
2856 quadrature_point_1,
2857 0));
2858 }
2859
2860 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2861 {
2862 const double L_i = legendre_polynomials[i + 1].value(
2863 edge_quadrature_points[q_point][0]);
2864
2865 for (unsigned int j = 0;
2866 j < GeometryInfo<dim>::lines_per_face;
2867 ++j)
2868 system_rhs(i, j) += tmp(j) * L_i;
2869 }
2870 }
2871
2872 system_matrix_inv.mmult(solution, system_rhs);
2873
2874 for (unsigned int i = 0;
2875 i < GeometryInfo<dim>::lines_per_face;
2876 ++i)
2877 for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2878 if (std::abs(solution(j, i)) > 1e-14)
2879 interpolation_matrix(i * source_fe.degree + j + 1,
2880 dof) = solution(j, i);
2881 }
2882
2883 const QGauss<2> quadrature(source_fe.degree);
2884 const std::vector<Point<2>> &quadrature_points =
2885 quadrature.get_points();
2886 const std::vector<Polynomials::Polynomial<double>>
2887 &lobatto_polynomials =
2889 source_fe.degree);
2890 const unsigned int n_boundary_dofs =
2892 const unsigned int n_quadrature_points = quadrature.size();
2893
2894 {
2895 FullMatrix<double> assembling_matrix(source_fe.degree *
2896 (source_fe.degree - 1),
2897 n_quadrature_points);
2898
2899 for (unsigned int q_point = 0; q_point < n_quadrature_points;
2900 ++q_point)
2901 {
2902 const double weight = std::sqrt(quadrature.weight(q_point));
2903
2904 for (unsigned int i = 0; i < source_fe.degree; ++i)
2905 {
2906 const double L_i =
2907 weight * legendre_polynomials[i].value(
2908 quadrature_points[q_point][0]);
2909
2910 for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2911 assembling_matrix(i * (source_fe.degree - 1) + j,
2912 q_point) =
2913 L_i * lobatto_polynomials[j + 2].value(
2914 quadrature_points[q_point][1]);
2915 }
2916 }
2917
2918 FullMatrix<double> system_matrix(assembling_matrix.m(),
2919 assembling_matrix.m());
2920
2921 assembling_matrix.mTmult(system_matrix, assembling_matrix);
2922 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
2923 system_matrix_inv.invert(system_matrix);
2924 }
2925
2926 solution.reinit(system_matrix_inv.m(), 2);
2927 system_rhs.reinit(system_matrix_inv.m(), 2);
2928 tmp.reinit(2);
2929
2930 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
2931 ++dof)
2932 {
2933 system_rhs = 0.0;
2934
2935 for (unsigned int q_point = 0; q_point < n_quadrature_points;
2936 ++q_point)
2937 {
2938 Point<dim> quadrature_point(
2939 0.5 *
2940 (quadrature_points[q_point][0] + shifts[subface][0]),
2941 0.5 *
2942 (quadrature_points[q_point][1] + shifts[subface][1]),
2943 0.0);
2944 tmp(0) = 0.5 * this->shape_value_component(
2945 this->face_to_cell_index(dof, 4),
2946 quadrature_point,
2947 0);
2948 tmp(1) = 0.5 * this->shape_value_component(
2949 this->face_to_cell_index(dof, 4),
2950 quadrature_point,
2951 1);
2952 quadrature_point =
2953 Point<dim>(quadrature_points[q_point][0],
2954 quadrature_points[q_point][1],
2955 0.0);
2956
2957 for (unsigned int i = 0; i < 2; ++i)
2958 for (unsigned int j = 0; j < source_fe.degree; ++j)
2959 {
2960 tmp(0) -= interpolation_matrix(
2961 (i + 2) * source_fe.degree + j, dof) *
2962 source_fe.shape_value_component(
2963 (i + 2) * source_fe.degree + j,
2964 quadrature_point,
2965 0);
2966 tmp(1) -=
2967 interpolation_matrix(i * source_fe.degree + j,
2968 dof) *
2969 source_fe.shape_value_component(
2970 i * source_fe.degree + j, quadrature_point, 1);
2971 }
2972
2973 tmp *= quadrature.weight(q_point);
2974
2975 for (unsigned int i = 0; i < source_fe.degree; ++i)
2976 {
2977 const double L_i_0 = legendre_polynomials[i].value(
2978 quadrature_points[q_point][0]);
2979 const double L_i_1 = legendre_polynomials[i].value(
2980 quadrature_points[q_point][1]);
2981
2982 for (unsigned int j = 0; j < source_fe.degree - 1;
2983 ++j)
2984 {
2985 system_rhs(i * (source_fe.degree - 1) + j, 0) +=
2986 tmp(0) * L_i_0 *
2987 lobatto_polynomials[j + 2].value(
2988 quadrature_points[q_point][1]);
2989 system_rhs(i * (source_fe.degree - 1) + j, 1) +=
2990 tmp(1) * L_i_1 *
2991 lobatto_polynomials[j + 2].value(
2992 quadrature_points[q_point][0]);
2993 }
2994 }
2995 }
2996
2997 system_matrix_inv.mmult(solution, system_rhs);
2998
2999 for (unsigned int i = 0; i < source_fe.degree; ++i)
3000 for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
3001 {
3002 if (std::abs(solution(i * (source_fe.degree - 1) + j,
3003 0)) > 1e-14)
3004 interpolation_matrix(i * (source_fe.degree - 1) + j +
3005 n_boundary_dofs,
3006 dof) =
3007 solution(i * (source_fe.degree - 1) + j, 0);
3008
3009 if (std::abs(solution(i * (source_fe.degree - 1) + j,
3010 1)) > 1e-14)
3011 interpolation_matrix(
3012 i + (j + source_fe.degree - 1) * source_fe.degree +
3013 n_boundary_dofs,
3014 dof) = solution(i * (source_fe.degree - 1) + j, 1);
3015 }
3016 }
3017 }
3018
3019 break;
3020 }
3021
3022 default:
3024 }
3025}
3026
3027template <int dim>
3028const FullMatrix<double> &
3030 const unsigned int child,
3031 const RefinementCase<dim> &refinement_case) const
3032{
3033 AssertIndexRange(refinement_case,
3035 Assert(refinement_case != RefinementCase<dim>::no_refinement,
3036 ExcMessage(
3037 "Prolongation matrices are only available for refined cells!"));
3038 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
3039
3040 // initialization upon first request
3041 if (this->prolongation[refinement_case - 1][child].n() == 0)
3042 {
3043 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
3044
3045 // if matrix got updated while waiting for the lock
3046 if (this->prolongation[refinement_case - 1][child].n() ==
3047 this->n_dofs_per_cell())
3048 return this->prolongation[refinement_case - 1][child];
3049
3050 // now do the work. need to get a non-const version of data in order to
3051 // be able to modify them inside a const function
3052 FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
3053
3054 // Reinit the vectors of
3055 // restriction and prolongation
3056 // matrices to the right sizes.
3057 // Restriction only for isotropic
3058 // refinement
3059#ifdef DEBUG_NEDELEC
3060 deallog << "Embedding" << std::endl;
3061#endif
3063 // Fill prolongation matrices with embedding operators
3065 this_nonconst,
3066 this_nonconst.prolongation,
3067 true,
3068 internal::FE_Nedelec::get_embedding_computation_tolerance(
3069 this->degree));
3070#ifdef DEBUG_NEDELEC
3071 deallog << "Restriction" << std::endl;
3072#endif
3073 this_nonconst.initialize_restriction();
3074 }
3075
3076 // we use refinement_case-1 here. the -1 takes care of the origin of the
3077 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3078 // available and so the vector indices are shifted
3079 return this->prolongation[refinement_case - 1][child];
3080}
3081
3082template <int dim>
3083const FullMatrix<double> &
3085 const unsigned int child,
3086 const RefinementCase<dim> &refinement_case) const
3087{
3088 AssertIndexRange(refinement_case,
3090 Assert(refinement_case != RefinementCase<dim>::no_refinement,
3091 ExcMessage(
3092 "Restriction matrices are only available for refined cells!"));
3094 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
3095
3096 // initialization upon first request
3097 if (this->restriction[refinement_case - 1][child].n() == 0)
3098 {
3099 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
3100
3101 // if matrix got updated while waiting for the lock...
3102 if (this->restriction[refinement_case - 1][child].n() ==
3103 this->n_dofs_per_cell())
3104 return this->restriction[refinement_case - 1][child];
3105
3106 // now do the work. need to get a non-const version of data in order to
3107 // be able to modify them inside a const function
3108 FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
3109
3110 // Reinit the vectors of
3111 // restriction and prolongation
3112 // matrices to the right sizes.
3113 // Restriction only for isotropic
3114 // refinement
3115#ifdef DEBUG_NEDELEC
3116 deallog << "Embedding" << std::endl;
3117#endif
3119 // Fill prolongation matrices with embedding operators
3121 this_nonconst,
3122 this_nonconst.prolongation,
3123 true,
3124 internal::FE_Nedelec::get_embedding_computation_tolerance(
3125 this->degree));
3126#ifdef DEBUG_NEDELEC
3127 deallog << "Restriction" << std::endl;
3128#endif
3129 this_nonconst.initialize_restriction();
3130 }
3131
3132 // we use refinement_case-1 here. the -1 takes care of the origin of the
3133 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3134 // available and so the vector indices are shifted
3135 return this->restriction[refinement_case - 1][child];
3136}
3137
3138
3139// Interpolate a function, which is given by
3140// its values at the generalized support
3141// points in the finite element space on the
3142// reference cell.
3143// This is done as usual by projection-based
3144// interpolation.
3145template <int dim>
3146void
3148 const std::vector<Vector<double>> &support_point_values,
3149 std::vector<double> &nodal_values) const
3150{
3151 // TODO: the implementation makes the assumption that all faces have the
3152 // same number of dofs
3153 AssertDimension(this->n_unique_faces(), 1);
3154 const unsigned int face_no = 0;
3155
3156 const unsigned int deg = this->degree - 1;
3157 Assert(support_point_values.size() == this->generalized_support_points.size(),
3158 ExcDimensionMismatch(support_point_values.size(),
3159 this->generalized_support_points.size()));
3160 Assert(support_point_values[0].size() == this->n_components(),
3161 ExcDimensionMismatch(support_point_values[0].size(),
3162 this->n_components()));
3163 Assert(nodal_values.size() == this->n_dofs_per_cell(),
3164 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
3165 std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3166
3167 switch (dim)
3168 {
3169 case 2:
3170 {
3171 // Let us begin with the
3172 // interpolation part.
3173 const QGauss<1> reference_edge_quadrature(this->degree);
3174 const unsigned int n_edge_points = reference_edge_quadrature.size();
3175
3176 for (unsigned int i = 0; i < 2; ++i)
3177 for (unsigned int j = 0; j < 2; ++j)
3178 {
3179 for (unsigned int q_point = 0; q_point < n_edge_points;
3180 ++q_point)
3181 nodal_values[(i + 2 * j) * this->degree] +=
3182 reference_edge_quadrature.weight(q_point) *
3183 support_point_values[q_point + (i + 2 * j) * n_edge_points]
3184 [1 - j];
3185
3186 // Add the computed support_point_values to the resulting vector
3187 // only, if they are not
3188 // too small.
3189 if (std::abs(nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3190 nodal_values[(i + 2 * j) * this->degree] = 0.0;
3191 }
3192
3193 // If the Nedelec element degree is greater
3194 // than 0 (i.e., the polynomial degree is greater than 1),
3195 // then we have still some higher order edge
3196 // shape functions to consider.
3197 // Note that this->degree returns the polynomial
3198 // degree.
3199 // Here the projection part starts.
3200 // The dof support_point_values are obtained by solving
3201 // a linear system of
3202 // equations.
3203 if (this->degree > 1)
3204 {
3205 // We start with projection
3206 // on the higher order edge
3207 // shape function.
3208 const std::vector<Polynomials::Polynomial<double>>
3209 &lobatto_polynomials =
3211 FullMatrix<double> system_matrix(this->degree - 1,
3212 this->degree - 1);
3213 std::vector<Polynomials::Polynomial<double>>
3214 lobatto_polynomials_grad(this->degree);
3215
3216 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3217 lobatto_polynomials_grad[i] =
3218 lobatto_polynomials[i + 1].derivative();
3219
3220 // Set up the system matrix.
3221 // This can be used for all
3222 // edges.
3223 for (unsigned int i = 0; i < system_matrix.m(); ++i)
3224 for (unsigned int j = 0; j < system_matrix.n(); ++j)
3225 for (unsigned int q_point = 0; q_point < n_edge_points;
3226 ++q_point)
3227 system_matrix(i, j) +=
3228 boundary_weights(q_point, j) *
3229 lobatto_polynomials_grad[i + 1].value(
3230 this->generalized_face_support_points[face_no][q_point]
3231 [0]);
3232
3233 FullMatrix<double> system_matrix_inv(this->degree - 1,
3234 this->degree - 1);
3235
3236 system_matrix_inv.invert(system_matrix);
3237
3238 const unsigned int
3239 line_coordinate[GeometryInfo<2>::lines_per_cell] = {1, 1, 0, 0};
3240 Vector<double> system_rhs(system_matrix.m());
3241 Vector<double> solution(system_rhs.size());
3242
3243 for (unsigned int line = 0;
3244 line < GeometryInfo<dim>::lines_per_cell;
3245 ++line)
3246 {
3247 // Set up the right hand side.
3248 system_rhs = 0;
3249
3250 for (unsigned int q_point = 0; q_point < n_edge_points;
3251 ++q_point)
3252 {
3253 const double tmp =
3254 support_point_values[line * n_edge_points + q_point]
3255 [line_coordinate[line]] -
3256 nodal_values[line * this->degree] *
3257 this->shape_value_component(
3258 line * this->degree,
3259 this->generalized_support_points[line *
3260 n_edge_points +
3261 q_point],
3262 line_coordinate[line]);
3263
3264 for (unsigned int i = 0; i < system_rhs.size(); ++i)
3265 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3266 }
3267
3268 system_matrix_inv.vmult(solution, system_rhs);
3269
3270 // Add the computed support_point_values
3271 // to the resulting vector
3272 // only, if they are not
3273 // too small.
3274 for (unsigned int i = 0; i < solution.size(); ++i)
3275 if (std::abs(solution(i)) > 1e-14)
3276 nodal_values[line * this->degree + i + 1] = solution(i);
3277 }
3278
3279 // Then we go on to the
3280 // interior shape
3281 // functions. Again we
3282 // set up the system
3283 // matrix and use it
3284 // for both, the
3285 // horizontal and the
3286 // vertical, interior
3287 // shape functions.
3288 const QGauss<dim> reference_quadrature(this->degree);
3289 const unsigned int n_interior_points =
3290 reference_quadrature.size();
3291 const std::vector<Polynomials::Polynomial<double>>
3292 &legendre_polynomials =
3294 1);
3295
3296 system_matrix.reinit((this->degree - 1) * this->degree,
3297 (this->degree - 1) * this->degree);
3298 system_matrix = 0;
3299
3300 for (unsigned int i = 0; i < this->degree; ++i)
3301 for (unsigned int j = 0; j < this->degree - 1; ++j)
3302 for (unsigned int k = 0; k < this->degree; ++k)
3303 for (unsigned int l = 0; l < this->degree - 1; ++l)
3304 for (unsigned int q_point = 0;
3305 q_point < n_interior_points;
3306 ++q_point)
3307 system_matrix(i * (this->degree - 1) + j,
3308 k * (this->degree - 1) + l) +=
3309 reference_quadrature.weight(q_point) *
3310 legendre_polynomials[i].value(
3311 this->generalized_support_points
3313 n_edge_points][0]) *
3314 lobatto_polynomials[j + 2].value(
3315 this->generalized_support_points
3317 n_edge_points][1]) *
3318 lobatto_polynomials_grad[k].value(
3319 this->generalized_support_points
3321 n_edge_points][0]) *
3322 lobatto_polynomials[l + 2].value(
3323 this->generalized_support_points
3325 n_edge_points][1]);
3326
3327 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3328 system_matrix_inv.invert(system_matrix);
3329 // Set up the right hand side
3330 // for the horizontal shape
3331 // functions.
3332 system_rhs.reinit(system_matrix_inv.m());
3333 system_rhs = 0;
3334
3335 for (unsigned int q_point = 0; q_point < n_interior_points;
3336 ++q_point)
3337 {
3338 double tmp =
3339 support_point_values[q_point +
3341 n_edge_points][0];
3342
3343 for (unsigned int i = 0; i < 2; ++i)
3344 for (unsigned int j = 0; j <= deg; ++j)
3345 tmp -= nodal_values[(i + 2) * this->degree + j] *
3346 this->shape_value_component(
3347 (i + 2) * this->degree + j,
3348 this->generalized_support_points
3350 n_edge_points],
3351 0);
3352
3353 for (unsigned int i = 0; i <= deg; ++i)
3354 for (unsigned int j = 0; j < deg; ++j)
3355 system_rhs(i * deg + j) +=
3356 reference_quadrature.weight(q_point) * tmp *
3357 lobatto_polynomials_grad[i].value(
3358 this->generalized_support_points
3360 n_edge_points][0]) *
3361 lobatto_polynomials[j + 2].value(
3362 this->generalized_support_points
3364 n_edge_points][1]);
3365 }
3366
3367 solution.reinit(system_matrix.m());
3368 system_matrix_inv.vmult(solution, system_rhs);
3369
3370 // Add the computed support_point_values
3371 // to the resulting vector
3372 // only, if they are not
3373 // too small.
3374 for (unsigned int i = 0; i <= deg; ++i)
3375 for (unsigned int j = 0; j < deg; ++j)
3376 if (std::abs(solution(i * deg + j)) > 1e-14)
3377 nodal_values[(i + GeometryInfo<dim>::lines_per_cell) * deg +
3379 solution(i * deg + j);
3380
3381 system_rhs = 0;
3382 // Set up the right hand side
3383 // for the vertical shape
3384 // functions.
3385
3386 for (unsigned int q_point = 0; q_point < n_interior_points;
3387 ++q_point)
3388 {
3389 double tmp =
3390 support_point_values[q_point +
3392 n_edge_points][1];
3393
3394 for (unsigned int i = 0; i < 2; ++i)
3395 for (unsigned int j = 0; j <= deg; ++j)
3396 tmp -= nodal_values[i * this->degree + j] *
3397 this->shape_value_component(
3398 i * this->degree + j,
3399 this->generalized_support_points
3401 n_edge_points],
3402 1);
3403
3404 for (unsigned int i = 0; i <= deg; ++i)
3405 for (unsigned int j = 0; j < deg; ++j)
3406 system_rhs(i * deg + j) +=
3407 reference_quadrature.weight(q_point) * tmp *
3408 lobatto_polynomials_grad[i].value(
3409 this->generalized_support_points
3411 n_edge_points][1]) *
3412 lobatto_polynomials[j + 2].value(
3413 this->generalized_support_points
3415 n_edge_points][0]);
3416 }
3417
3418 system_matrix_inv.vmult(solution, system_rhs);
3419
3420 // Add the computed support_point_values
3421 // to the resulting vector
3422 // only, if they are not
3423 // too small.
3424 for (unsigned int i = 0; i <= deg; ++i)
3425 for (unsigned int j = 0; j < deg; ++j)
3426 if (std::abs(solution(i * deg + j)) > 1e-14)
3427 nodal_values[i +
3429 this->degree] = solution(i * deg + j);
3430 }
3431
3432 break;
3433 }
3434
3435 case 3:
3436 {
3437 // Let us begin with the
3438 // interpolation part.
3439 const QGauss<1> reference_edge_quadrature(this->degree);
3440 const unsigned int n_edge_points = reference_edge_quadrature.size();
3441
3442 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3443 {
3444 for (unsigned int i = 0; i < 4; ++i)
3445 nodal_values[(i + 8) * this->degree] +=
3446 reference_edge_quadrature.weight(q_point) *
3447 support_point_values[q_point + (i + 8) * n_edge_points][2];
3448
3449 for (unsigned int i = 0; i < 2; ++i)
3450 for (unsigned int j = 0; j < 2; ++j)
3451 for (unsigned int k = 0; k < 2; ++k)
3452 nodal_values[(i + 2 * (2 * j + k)) * this->degree] +=
3453 reference_edge_quadrature.weight(q_point) *
3454 support_point_values[q_point + (i + 2 * (2 * j + k)) *
3455 n_edge_points][1 - k];
3456 }
3457
3458 // Add the computed support_point_values
3459 // to the resulting vector
3460 // only, if they are not
3461 // too small.
3462 for (unsigned int i = 0; i < 4; ++i)
3463 if (std::abs(nodal_values[(i + 8) * this->degree]) < 1e-14)
3464 nodal_values[(i + 8) * this->degree] = 0.0;
3465
3466 for (unsigned int i = 0; i < 2; ++i)
3467 for (unsigned int j = 0; j < 2; ++j)
3468 for (unsigned int k = 0; k < 2; ++k)
3469 if (std::abs(
3470 nodal_values[(i + 2 * (2 * j + k)) * this->degree]) <
3471 1e-14)
3472 nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3473
3474 // If the degree is greater
3475 // than 0, then we have still
3476 // some higher order shape
3477 // functions to consider.
3478 // Here the projection part
3479 // starts. The dof support_point_values
3480 // are obtained by solving
3481 // a linear system of
3482 // equations.
3483 if (this->degree > 1)
3484 {
3485 // We start with projection
3486 // on the higher order edge
3487 // shape function.
3488 const std::vector<Polynomials::Polynomial<double>>
3489 &lobatto_polynomials =
3491 FullMatrix<double> system_matrix(this->degree - 1,
3492 this->degree - 1);
3493 std::vector<Polynomials::Polynomial<double>>
3494 lobatto_polynomials_grad(this->degree);
3495
3496 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3497 lobatto_polynomials_grad[i] =
3498 lobatto_polynomials[i + 1].derivative();
3499
3500 // Set up the system matrix.
3501 // This can be used for all
3502 // edges.
3503 for (unsigned int i = 0; i < system_matrix.m(); ++i)
3504 for (unsigned int j = 0; j < system_matrix.n(); ++j)
3505 for (unsigned int q_point = 0; q_point < n_edge_points;
3506 ++q_point)
3507 system_matrix(i, j) +=
3508 boundary_weights(q_point, j) *
3509 lobatto_polynomials_grad[i + 1].value(
3510 this->generalized_face_support_points[face_no][q_point]
3511 [1]);
3512
3513 FullMatrix<double> system_matrix_inv(this->degree - 1,
3514 this->degree - 1);
3515
3516 system_matrix_inv.invert(system_matrix);
3517
3518 const unsigned int
3519 line_coordinate[GeometryInfo<3>::lines_per_cell] = {
3520 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3521 Vector<double> system_rhs(system_matrix.m());
3522 Vector<double> solution(system_rhs.size());
3523
3524 for (unsigned int line = 0;
3525 line < GeometryInfo<dim>::lines_per_cell;
3526 ++line)
3527 {
3528 // Set up the right hand side.
3529 system_rhs = 0;
3530
3531 for (unsigned int q_point = 0; q_point < this->degree;
3532 ++q_point)
3533 {
3534 const double tmp =
3535 support_point_values[line * this->degree + q_point]
3536 [line_coordinate[line]] -
3537 nodal_values[line * this->degree] *
3538 this->shape_value_component(
3539 line * this->degree,
3540 this
3541 ->generalized_support_points[line * this->degree +
3542 q_point],
3543 line_coordinate[line]);
3544
3545 for (unsigned int i = 0; i < system_rhs.size(); ++i)
3546 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3547 }
3548
3549 system_matrix_inv.vmult(solution, system_rhs);
3550
3551 // Add the computed values
3552 // to the resulting vector
3553 // only, if they are not
3554 // too small.
3555 for (unsigned int i = 0; i < solution.size(); ++i)
3556 if (std::abs(solution(i)) > 1e-14)
3557 nodal_values[line * this->degree + i + 1] = solution(i);
3558 }
3559
3560 // Then we go on to the
3561 // face shape functions.
3562 // Again we set up the
3563 // system matrix and
3564 // use it for both, the
3565 // horizontal and the
3566 // vertical, shape
3567 // functions.
3568 const std::vector<Polynomials::Polynomial<double>>
3569 &legendre_polynomials =
3571 1);
3572 const unsigned int n_face_points = n_edge_points * n_edge_points;
3573
3574 system_matrix.reinit((this->degree - 1) * this->degree,
3575 (this->degree - 1) * this->degree);
3576 system_matrix = 0;
3577
3578 for (unsigned int i = 0; i < this->degree; ++i)
3579 for (unsigned int j = 0; j < this->degree - 1; ++j)
3580 for (unsigned int k = 0; k < this->degree; ++k)
3581 for (unsigned int l = 0; l < this->degree - 1; ++l)
3582 for (unsigned int q_point = 0; q_point < n_face_points;
3583 ++q_point)
3584 system_matrix(i * (this->degree - 1) + j,
3585 k * (this->degree - 1) + l) +=
3586 boundary_weights(q_point + n_edge_points,
3587 2 * (k * (this->degree - 1) + l)) *
3588 legendre_polynomials[i].value(
3589 this->generalized_face_support_points
3590 [face_no][q_point + 4 * n_edge_points][0]) *
3591 lobatto_polynomials[j + 2].value(
3592 this->generalized_face_support_points
3593 [face_no][q_point + 4 * n_edge_points][1]);
3594
3595 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3596 system_matrix_inv.invert(system_matrix);
3597 solution.reinit(system_matrix.m());
3598 system_rhs.reinit(system_matrix.m());
3599
3600 const unsigned int
3601 face_coordinates[GeometryInfo<3>::faces_per_cell][2] = {
3602 {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3605 {{0, 4, 8, 10},
3606 {1, 5, 9, 11},
3607 {8, 9, 2, 6},
3608 {10, 11, 3, 7},
3609 {2, 3, 0, 1},
3610 {6, 7, 4, 5}};
3611
3612 for (const unsigned int face : GeometryInfo<dim>::face_indices())
3613 {
3614 // Set up the right hand side
3615 // for the horizontal shape
3616 // functions.
3617 system_rhs = 0;
3618
3619 for (unsigned int q_point = 0; q_point < n_face_points;
3620 ++q_point)
3621 {
3622 double tmp =
3623 support_point_values[q_point +
3625 n_edge_points]
3626 [face_coordinates[face][0]];
3627
3628 for (unsigned int i = 0; i < 2; ++i)
3629 for (unsigned int j = 0; j <= deg; ++j)
3630 tmp -=
3631 nodal_values[edge_indices[face][i] * this->degree +
3632 j] *
3633 this->shape_value_component(
3634 edge_indices[face][i] * this->degree + j,
3635 this->generalized_support_points
3637 n_edge_points],
3638 face_coordinates[face][0]);
3639
3640 for (unsigned int i = 0; i <= deg; ++i)
3641 for (unsigned int j = 0; j < deg; ++j)
3642 system_rhs(i * deg + j) +=
3643 boundary_weights(q_point + n_edge_points,
3644 2 * (i * deg + j)) *
3645 tmp;
3646 }
3647
3648 system_matrix_inv.vmult(solution, system_rhs);
3649
3650 // Add the computed support_point_values
3651 // to the resulting vector
3652 // only, if they are not
3653 // too small.
3654 for (unsigned int i = 0; i <= deg; ++i)
3655 for (unsigned int j = 0; j < deg; ++j)
3656 if (std::abs(solution(i * deg + j)) > 1e-14)
3657 nodal_values[(2 * face * this->degree + i +
3659 deg +
3661 solution(i * deg + j);
3662
3663 // Set up the right hand side
3664 // for the vertical shape
3665 // functions.
3666 system_rhs = 0;
3667
3668 for (unsigned int q_point = 0; q_point < n_face_points;
3669 ++q_point)
3670 {
3671 double tmp =
3672 support_point_values[q_point +
3674 n_edge_points]
3675 [face_coordinates[face][1]];
3676
3677 for (unsigned int i = 2;
3678 i < GeometryInfo<dim>::lines_per_face;
3679 ++i)
3680 for (unsigned int j = 0; j <= deg; ++j)
3681 tmp -=
3682 nodal_values[edge_indices[face][i] * this->degree +
3683 j] *
3684 this->shape_value_component(
3685 edge_indices[face][i] * this->degree + j,
3686 this->generalized_support_points
3688 n_edge_points],
3689 face_coordinates[face][1]);
3690
3691 for (unsigned int i = 0; i <= deg; ++i)
3692 for (unsigned int j = 0; j < deg; ++j)
3693 system_rhs(i * deg + j) +=
3694 boundary_weights(q_point + n_edge_points,
3695 2 * (i * deg + j) + 1) *
3696 tmp;
3697 }
3698
3699 system_matrix_inv.vmult(solution, system_rhs);
3700
3701 // Add the computed support_point_values
3702 // to the resulting vector
3703 // only, if they are not
3704 // too small.
3705 for (unsigned int i = 0; i <= deg; ++i)
3706 for (unsigned int j = 0; j < deg; ++j)
3707 if (std::abs(solution(i * deg + j)) > 1e-14)
3708 nodal_values[((2 * face + 1) * deg + j +
3710 this->degree +
3711 i] = solution(i * deg + j);
3712 }
3713
3714 // Finally we project
3715 // the remaining parts
3716 // of the function on
3717 // the interior shape
3718 // functions.
3719 const QGauss<dim> reference_quadrature(this->degree);
3720 const unsigned int n_interior_points =
3721 reference_quadrature.size();
3722
3723 // We create the
3724 // system matrix.
3725 system_matrix.reinit(this->degree * deg * deg,
3726 this->degree * deg * deg);
3727 system_matrix = 0;
3728
3729 for (unsigned int i = 0; i <= deg; ++i)
3730 for (unsigned int j = 0; j < deg; ++j)
3731 for (unsigned int k = 0; k < deg; ++k)
3732 for (unsigned int l = 0; l <= deg; ++l)
3733 for (unsigned int m = 0; m < deg; ++m)
3734 for (unsigned int n = 0; n < deg; ++n)
3735 for (unsigned int q_point = 0;
3736 q_point < n_interior_points;
3737 ++q_point)
3738 system_matrix((i * deg + j) * deg + k,
3739 (l * deg + m) * deg + n) +=
3740 reference_quadrature.weight(q_point) *
3741 legendre_polynomials[i].value(
3742 this->generalized_support_points
3743 [q_point +
3745 n_edge_points +
3747 n_face_points][0]) *
3748 lobatto_polynomials[j + 2].value(
3749 this->generalized_support_points
3750 [q_point +
3752 n_edge_points +
3754 n_face_points][1]) *
3755 lobatto_polynomials[k + 2].value(
3756 this->generalized_support_points
3757 [q_point +
3759 n_edge_points +
3761 n_face_points][2]) *
3762 lobatto_polynomials_grad[l].value(
3763 this->generalized_support_points
3764 [q_point +
3766 n_edge_points +
3768 n_face_points][0]) *
3769 lobatto_polynomials[m + 2].value(
3770 this->generalized_support_points
3771 [q_point +
3773 n_edge_points +
3775 n_face_points][1]) *
3776 lobatto_polynomials[n + 2].value(
3777 this->generalized_support_points
3778 [q_point +
3780 n_edge_points +
3782 n_face_points][2]);
3783
3784 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3785 system_matrix_inv.invert(system_matrix);
3786 // Set up the right hand side.
3787 system_rhs.reinit(system_matrix.m());
3788 system_rhs = 0;
3789
3790 for (unsigned int q_point = 0; q_point < n_interior_points;
3791 ++q_point)
3792 {
3793 double tmp =
3794 support_point_values[q_point +
3796 n_edge_points +
3798 n_face_points][0];
3799
3800 for (unsigned int i = 0; i <= deg; ++i)
3801 {
3802 for (unsigned int j = 0; j < 2; ++j)
3803 for (unsigned int k = 0; k < 2; ++k)
3804 tmp -=
3805 nodal_values[i + (j + 4 * k + 2) * this->degree] *
3806 this->shape_value_component(
3807 i + (j + 4 * k + 2) * this->degree,
3808 this->generalized_support_points
3809 [q_point +
3811 n_edge_points +
3813 n_face_points],
3814 0);
3815
3816 for (unsigned int j = 0; j < deg; ++j)
3817 for (unsigned int k = 0; k < 4; ++k)
3818 tmp -=
3819 nodal_values[(i + 2 * (k + 2) * this->degree +
3821 deg +
3822 j +
3824 this->shape_value_component(
3825 (i + 2 * (k + 2) * this->degree +
3827 deg +
3829 this->generalized_support_points
3830 [q_point +
3832 n_edge_points +
3834 n_face_points],
3835 0);
3836 }
3837
3838 for (unsigned int i = 0; i <= deg; ++i)
3839 for (unsigned int j = 0; j < deg; ++j)
3840 for (unsigned int k = 0; k < deg; ++k)
3841 system_rhs((i * deg + j) * deg + k) +=
3842 reference_quadrature.weight(q_point) * tmp *
3843 lobatto_polynomials_grad[i].value(
3844 this->generalized_support_points
3845 [q_point +
3847 n_edge_points +
3849 n_face_points][0]) *
3850 lobatto_polynomials[j + 2].value(
3851 this->generalized_support_points
3852 [q_point +
3854 n_edge_points +
3856 n_face_points][1]) *
3857 lobatto_polynomials[k + 2].value(
3858 this->generalized_support_points
3859 [q_point +
3861 n_edge_points +
3863 n_face_points][2]);
3864 }
3865
3866 solution.reinit(system_rhs.size());
3867 system_matrix_inv.vmult(solution, system_rhs);
3868
3869 // Add the computed values
3870 // to the resulting vector
3871 // only, if they are not
3872 // too small.
3873 for (unsigned int i = 0; i <= deg; ++i)
3874 for (unsigned int j = 0; j < deg; ++j)
3875 for (unsigned int k = 0; k < deg; ++k)
3876 if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3877 nodal_values
3878 [((i + 2 * GeometryInfo<dim>::faces_per_cell) * deg +
3881 deg +
3883 solution((i * deg + j) * deg + k);
3884
3885 // Set up the right hand side.
3886 system_rhs = 0;
3887
3888 for (unsigned int q_point = 0; q_point < n_interior_points;
3889 ++q_point)
3890 {
3891 double tmp =
3892 support_point_values[q_point +
3894 n_edge_points +
3896 n_face_points][1];
3897
3898 for (unsigned int i = 0; i <= deg; ++i)
3899 for (unsigned int j = 0; j < 2; ++j)
3900 {
3901 for (unsigned int k = 0; k < 2; ++k)
3902 tmp -= nodal_values[i + (4 * j + k) * this->degree] *
3903 this->shape_value_component(
3904 i + (4 * j + k) * this->degree,
3905 this->generalized_support_points
3906 [q_point +
3908 n_edge_points +
3910 n_face_points],
3911 1);
3912
3913 for (unsigned int k = 0; k < deg; ++k)
3914 tmp -=
3915 nodal_values[(i + 2 * j * this->degree +
3917 deg +
3918 k +
3920 this->shape_value_component(
3921 (i + 2 * j * this->degree +
3923 deg +
3925 this->generalized_support_points
3926 [q_point +
3928 n_edge_points +
3930 n_face_points],
3931 1) +
3932 nodal_values[i +
3933 ((2 * j + 9) * deg + k +
3935 this->degree] *
3936 this->shape_value_component(
3937 i + ((2 * j + 9) * deg + k +
3939 this->degree,
3940 this->generalized_support_points
3941 [q_point +
3943 n_edge_points +
3945 n_face_points],
3946 1);
3947 }
3948
3949 for (unsigned int i = 0; i <= deg; ++i)
3950 for (unsigned int j = 0; j < deg; ++j)
3951 for (unsigned int k = 0; k < deg; ++k)
3952 system_rhs((i * deg + j) * deg + k) +=
3953 reference_quadrature.weight(q_point) * tmp *
3954 lobatto_polynomials_grad[i].value(
3955 this->generalized_support_points
3956 [q_point +
3958 n_edge_points +
3960 n_face_points][1]) *
3961 lobatto_polynomials[j + 2].value(
3962 this->generalized_support_points
3963 [q_point +
3965 n_edge_points +
3967 n_face_points][0]) *
3968 lobatto_polynomials[k + 2].value(
3969 this->generalized_support_points
3970 [q_point +
3972 n_edge_points +
3974 n_face_points][2]);
3975 }
3976
3977 system_matrix_inv.vmult(solution, system_rhs);
3978
3979 // Add the computed support_point_values
3980 // to the resulting vector
3981 // only, if they are not
3982 // too small.
3983 for (unsigned int i = 0; i <= deg; ++i)
3984 for (unsigned int j = 0; j < deg; ++j)
3985 for (unsigned int k = 0; k < deg; ++k)
3986 if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3987 nodal_values[((i + this->degree +
3989 deg +
3992 deg +
3994 solution((i * deg + j) * deg + k);
3995
3996 // Set up the right hand side.
3997 system_rhs = 0;
3998
3999 for (unsigned int q_point = 0; q_point < n_interior_points;
4000 ++q_point)
4001 {
4002 double tmp =
4003 support_point_values[q_point +
4005 n_edge_points +
4007 n_face_points][2];
4008
4009 for (unsigned int i = 0; i <= deg; ++i)
4010 for (unsigned int j = 0; j < 4; ++j)
4011 {
4012 tmp -= nodal_values[i + (j + 8) * this->degree] *
4013 this->shape_value_component(
4014 i + (j + 8) * this->degree,
4015 this->generalized_support_points
4016 [q_point +
4018 n_edge_points +
4020 n_face_points],
4021 2);
4022
4023 for (unsigned int k = 0; k < deg; ++k)
4024 tmp -=
4025 nodal_values[i +
4026 ((2 * j + 1) * deg + k +
4028 this->degree] *
4029 this->shape_value_component(
4030 i + ((2 * j + 1) * deg + k +
4032 this->degree,
4033 this->generalized_support_points
4034 [q_point +
4036 n_edge_points +
4038 n_face_points],
4039 2);
4040 }
4041
4042 for (unsigned int i = 0; i <= deg; ++i)
4043 for (unsigned int j = 0; j < deg; ++j)
4044 for (unsigned int k = 0; k < deg; ++k)
4045 system_rhs((i * deg + j) * deg + k) +=
4046 reference_quadrature.weight(q_point) * tmp *
4047 lobatto_polynomials_grad[i].value(
4048 this->generalized_support_points
4049 [q_point +
4051 n_edge_points +
4053 n_face_points][2]) *
4054 lobatto_polynomials[j + 2].value(
4055 this->generalized_support_points
4056 [q_point +
4058 n_edge_points +
4060 n_face_points][0]) *
4061 lobatto_polynomials[k + 2].value(
4062 this->generalized_support_points
4063 [q_point +
4065 n_edge_points +
4067 n_face_points][1]);
4068 }
4069
4070 system_matrix_inv.vmult(solution, system_rhs);
4071
4072 // Add the computed support_point_values
4073 // to the resulting vector
4074 // only, if they are not
4075 // too small.
4076 for (unsigned int i = 0; i <= deg; ++i)
4077 for (unsigned int j = 0; j < deg; ++j)
4078 for (unsigned int k = 0; k < deg; ++k)
4079 if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4080 nodal_values
4081 [i +
4082 ((j + 2 * (deg + GeometryInfo<dim>::faces_per_cell)) *
4083 deg +
4085 this->degree] = solution((i * deg + j) * deg + k);
4086 }
4087
4088 break;
4089 }
4090
4091 default:
4093 }
4094}
4095
4096
4097
4098template <int dim>
4099std::pair<Table<2, bool>, std::vector<unsigned int>>
4101{
4102 Table<2, bool> constant_modes(dim, this->n_dofs_per_cell());
4103 for (unsigned int d = 0; d < dim; ++d)
4104 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
4105 constant_modes(d, i) = true;
4106 std::vector<unsigned int> components;
4107 for (unsigned int d = 0; d < dim; ++d)
4108 components.push_back(d);
4109 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
4110 components);
4111}
4112
4113
4114template <int dim>
4115std::size_t
4117{
4119 return 0;
4120}
4121
4122template <int dim>
4123std::vector<unsigned int>
4124FE_Nedelec<dim>::get_embedding_dofs(const unsigned int sub_degree) const
4125{
4126 Assert((sub_degree > 0) && (sub_degree <= this->degree),
4127 ExcIndexRange(sub_degree, 1, this->degree));
4128
4129 switch (dim)
4130 {
4131 case 2:
4132 {
4133 // The Nedelec cell has only Face (Line) and Cell DoFs...
4134 const unsigned int n_face_dofs_sub =
4136 const unsigned int n_cell_dofs_sub =
4137 2 * (sub_degree - 1) * sub_degree;
4138
4139 std::vector<unsigned int> embedding_dofs(n_face_dofs_sub +
4140 n_cell_dofs_sub);
4141
4142 unsigned int i = 0;
4143
4144 // Identify the Face/Line DoFs
4145 while (i < n_face_dofs_sub)
4146 {
4147 const unsigned int face_index = i / sub_degree;
4148 embedding_dofs[i] = i % sub_degree + face_index * this->degree;
4149 ++i;
4150 }
4151
4152 // Identify the Cell DoFs
4153 if (sub_degree >= 2)
4154 {
4155 const unsigned int n_face_dofs =
4156 GeometryInfo<dim>::lines_per_cell * this->degree;
4157
4158 // For the first component
4159 for (unsigned ku = 0; ku < sub_degree; ++ku)
4160 for (unsigned kv = 2; kv <= sub_degree; ++kv)
4161 embedding_dofs[i++] =
4162 n_face_dofs + ku * (this->degree - 1) + (kv - 2);
4163
4164 // For the second component
4165 for (unsigned ku = 2; ku <= sub_degree; ++ku)
4166 for (unsigned kv = 0; kv < sub_degree; ++kv)
4167 embedding_dofs[i++] = n_face_dofs +
4168 this->degree * (this->degree - 1) +
4169 (ku - 2) * (this->degree) + kv;
4170 }
4171 Assert(i == (n_face_dofs_sub + n_cell_dofs_sub), ExcInternalError());
4172 return embedding_dofs;
4173 }
4174 default:
4176 return std::vector<unsigned int>();
4177 }
4178}
4179//----------------------------------------------------------------------//
4180
4181
4182// explicit instantiations
4183#include "fe_nedelec.inst"
4184
4185
void initialize_quad_dof_index_permutation_and_sign_change()
std::vector< unsigned int > get_embedding_dofs(const unsigned int sub_degree) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual bool hp_constraints_are_implemented() const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
void initialize_restriction()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree, bool dg=false)
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual std::size_t memory_consumption() const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual std::string get_name() const override
friend class FE_Nedelec
Definition fe_nedelec.h:656
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
void initialize_support_points(const unsigned int order)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
FullMatrix< double > inverse_node_matrix
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
std::vector< MappingKind > mapping_kind
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
virtual std::string get_name() const =0
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
Definition fe.h:2423
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2411
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
size_type n() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void invert(const FullMatrix< number2 > &M)
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
size_type m() const
Definition point.h:111
static unsigned int n_polynomials(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
unsigned int edge_indices[GeometryInfo< dim >::lines_per_cell]
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ mapping_nedelec
Definition mapping.h:127
#define DEAL_II_NOT_IMPLEMENTED()
LogStream deallog
Definition logstream.cc:36
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()