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