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