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