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