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