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