Reference documentation for deal.II version 9.1.0-pre
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_face_domination(
2278  const FiniteElement<dim> &fe_other) const
2279 {
2280  if (const FE_Nedelec<dim> *fe_nedelec_other =
2281  dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2282  {
2283  if (this->degree < fe_nedelec_other->degree)
2284  return FiniteElementDomination::this_element_dominates;
2285  else if (this->degree == fe_nedelec_other->degree)
2286  return FiniteElementDomination::either_element_can_dominate;
2287  else
2288  return FiniteElementDomination::other_element_dominates;
2289  }
2290  else if (const FE_Nothing<dim> *fe_nothing =
2291  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
2292  {
2293  // TODO: ???
2294  // the FE_Nothing has no
2295  // degrees of
2296  // freedom. nevertheless, we
2297  // say that the FE_Q element
2298  // dominates so that we don't
2299  // have to force the FE_Q side
2300  // to become a zero function
2301  // and rather allow the
2302  // function to be discontinuous
2303  // along the interface
2304  // return FiniteElementDomination::other_element_dominates;
2305  if (fe_nothing->is_dominating())
2306  {
2307  return FiniteElementDomination::other_element_dominates;
2308  }
2309  else
2310  {
2311  // the FE_Nothing has no degrees of freedom and it is typically used
2312  // in a context where we don't require any continuity along the
2313  // interface
2314  return FiniteElementDomination::no_requirements;
2315  }
2316  }
2317 
2318  Assert(false, ExcNotImplemented());
2319  return FiniteElementDomination::neither_element_dominates;
2320 }
2321 
2322 template <int dim>
2323 bool
2324 FE_Nedelec<dim>::hp_constraints_are_implemented() const
2325 {
2326  return true;
2327 }
2328 
2329 template <int dim>
2330 std::vector<std::pair<unsigned int, unsigned int>>
2331 FE_Nedelec<dim>::hp_vertex_dof_identities(const FiniteElement<dim> &) const
2332 {
2333  // Nedelec elements do not have any dofs
2334  // on vertices, hence return an empty vector.
2335  return std::vector<std::pair<unsigned int, unsigned int>>();
2336 }
2337 
2338 template <int dim>
2339 std::vector<std::pair<unsigned int, unsigned int>>
2340 FE_Nedelec<dim>::hp_line_dof_identities(
2341  const FiniteElement<dim> &fe_other) const
2342 {
2343  // we can presently only compute these
2344  // identities if both FEs are
2345  // FE_Nedelec or if the other one is an
2346  // FE_Nothing
2347  if (const FE_Nedelec<dim> *fe_nedelec_other =
2348  dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2349  {
2350  // dofs are located on lines, so
2351  // two dofs are identical, if their
2352  // edge shape functions have the
2353  // same polynomial degree.
2354  std::vector<std::pair<unsigned int, unsigned int>> identities;
2355 
2356  for (unsigned int i = 0;
2357  i < std::min(fe_nedelec_other->degree, this->degree);
2358  ++i)
2359  identities.emplace_back(i, i);
2360 
2361  return identities;
2362  }
2363 
2364  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2365  {
2366  // the FE_Nothing has no
2367  // degrees of freedom, so there
2368  // are no equivalencies to be
2369  // recorded
2370  return std::vector<std::pair<unsigned int, unsigned int>>();
2371  }
2372 
2373  else
2374  {
2375  Assert(false, ExcNotImplemented());
2376  return std::vector<std::pair<unsigned int, unsigned int>>();
2377  }
2378 }
2379 
2380 template <int dim>
2381 std::vector<std::pair<unsigned int, unsigned int>>
2382 FE_Nedelec<dim>::hp_quad_dof_identities(
2383  const FiniteElement<dim> &fe_other) const
2384 {
2385  // we can presently only compute
2386  // these identities if both FEs are
2387  // FE_Nedelec or if the other one is an
2388  // FE_Nothing
2389  if (const FE_Nedelec<dim> *fe_nedelec_other =
2390  dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2391  {
2392  // dofs are located on the interior
2393  // of faces, so two dofs are identical,
2394  // if their face shape functions have
2395  // the same polynomial degree.
2396  const unsigned int p = fe_nedelec_other->degree;
2397  const unsigned int q = this->degree;
2398  const unsigned int p_min = std::min(p, q);
2399  std::vector<std::pair<unsigned int, unsigned int>> identities;
2400 
2401  for (unsigned int i = 0; i < p_min; ++i)
2402  for (unsigned int j = 0; j < p_min - 1; ++j)
2403  {
2404  identities.emplace_back(i * (q - 1) + j, i * (p - 1) + j);
2405  identities.emplace_back(i + (j + q - 1) * q, i + (j + p - 1) * p);
2406  }
2407 
2408  return identities;
2409  }
2410 
2411  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2412  {
2413  // the FE_Nothing has no
2414  // degrees of freedom, so there
2415  // are no equivalencies to be
2416  // recorded
2417  return std::vector<std::pair<unsigned int, unsigned int>>();
2418  }
2419 
2420  else
2421  {
2422  Assert(false, ExcNotImplemented());
2423  return std::vector<std::pair<unsigned int, unsigned int>>();
2424  }
2425 }
2426 
2427 // In this function we compute the face
2428 // interpolation matrix. This is usually
2429 // done by projection-based interpolation,
2430 // but, since one can compute the entries
2431 // easy per hand, we save some computation
2432 // time at this point and just fill in the
2433 // correct values.
2434 template <int dim>
2435 void
2436 FE_Nedelec<dim>::get_face_interpolation_matrix(
2437  const FiniteElement<dim> &source,
2438  FullMatrix<double> & interpolation_matrix) const
2439 {
2440  // this is only implemented, if the
2441  // source FE is also a
2442  // Nedelec element
2443  AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2444  (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2445  (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
2446  Assert(interpolation_matrix.m() == source.dofs_per_face,
2447  ExcDimensionMismatch(interpolation_matrix.m(), source.dofs_per_face));
2448  Assert(interpolation_matrix.n() == this->dofs_per_face,
2449  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
2450 
2451  // ok, source is a Nedelec element, so
2452  // we will be able to do the work
2453  const FE_Nedelec<dim> &source_fe =
2454  dynamic_cast<const FE_Nedelec<dim> &>(source);
2455 
2456  // Make sure, that the element,
2457  // for which the DoFs should be
2458  // constrained is the one with
2459  // the higher polynomial degree.
2460  // Actually the procedure will work
2461  // also if this assertion is not
2462  // satisfied. But the matrices
2463  // produced in that case might
2464  // lead to problems in the
2465  // hp procedures, which use this
2466  // method.
2467  Assert(this->dofs_per_face <= source_fe.dofs_per_face,
2468  (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
2469  interpolation_matrix = 0;
2470 
2471  // On lines we can just identify
2472  // all degrees of freedom.
2473  for (unsigned int i = 0; i < this->degree; ++i)
2474  interpolation_matrix(i, i) = 1.0;
2475 
2476  // In 3d we have some lines more
2477  // and a face. The procedure stays
2478  // the same as above, but we have
2479  // to take a bit more care of the
2480  // indices of the degrees of
2481  // freedom.
2482  if (dim == 3)
2483  {
2484  const unsigned int p = source_fe.degree;
2485  const unsigned int q = this->degree;
2486 
2487  for (unsigned int i = 0; i < q; ++i)
2488  {
2489  for (int j = 1; j < (int)GeometryInfo<dim>::lines_per_face; ++j)
2490  interpolation_matrix(j * p + i, j * q + i) = 1.0;
2491 
2492  for (unsigned int j = 0; j < q - 1; ++j)
2493  {
2494  interpolation_matrix(GeometryInfo<dim>::lines_per_face * p +
2495  i * (p - 1) + j,
2496  GeometryInfo<dim>::lines_per_face * q +
2497  i * (q - 1) + j) = 1.0;
2498  interpolation_matrix(GeometryInfo<dim>::lines_per_face * p + i +
2499  (j + p - 1) * p,
2500  GeometryInfo<dim>::lines_per_face * q + i +
2501  (j + q - 1) * q) = 1.0;
2502  }
2503  }
2504  }
2505 }
2506 
2507 
2508 
2509 template <>
2510 void
2511 FE_Nedelec<1>::get_subface_interpolation_matrix(const FiniteElement<1, 1> &,
2512  const unsigned int,
2513  FullMatrix<double> &) const
2514 {
2515  Assert(false, ExcNotImplemented());
2516 }
2517 
2518 
2519 
2520 // In this function we compute the
2521 // subface interpolation matrix.
2522 // This is done by a projection-
2523 // based interpolation. Therefore
2524 // we first interpolate the
2525 // shape functions of the higher
2526 // order element on the lowest
2527 // order edge shape functions.
2528 // Then the remaining part of
2529 // the interpolated shape
2530 // functions is projected on the
2531 // higher order edge shape
2532 // functions, the face shape
2533 // functions and the interior
2534 // shape functions (if they all
2535 // exist).
2536 template <int dim>
2537 void
2538 FE_Nedelec<dim>::get_subface_interpolation_matrix(
2539  const FiniteElement<dim> &source,
2540  const unsigned int subface,
2541  FullMatrix<double> & interpolation_matrix) const
2542 {
2543  // this is only implemented, if the
2544  // source FE is also a
2545  // Nedelec element
2546  AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2547  (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2548  typename FiniteElement<dim>::ExcInterpolationNotImplemented());
2549  Assert(interpolation_matrix.m() == source.dofs_per_face,
2550  ExcDimensionMismatch(interpolation_matrix.m(), source.dofs_per_face));
2551  Assert(interpolation_matrix.n() == this->dofs_per_face,
2552  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
2553 
2554  // ok, source is a Nedelec element, so
2555  // we will be able to do the work
2556  const FE_Nedelec<dim> &source_fe =
2557  dynamic_cast<const FE_Nedelec<dim> &>(source);
2558 
2559  // Make sure, that the element,
2560  // for which the DoFs should be
2561  // constrained is the one with
2562  // the higher polynomial degree.
2563  // Actually the procedure will work
2564  // also if this assertion is not
2565  // satisfied. But the matrices
2566  // produced in that case might
2567  // lead to problems in the
2568  // hp procedures, which use this
2569  // method.
2570  Assert(this->dofs_per_face <= source_fe.dofs_per_face,
2571  (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
2572  interpolation_matrix = 0.0;
2573  // Perform projection-based interpolation
2574  // as usual.
2575  const QGauss<1> edge_quadrature(source_fe.degree);
2576  const std::vector<Point<1>> &edge_quadrature_points =
2577  edge_quadrature.get_points();
2578  const unsigned int &n_edge_quadrature_points = edge_quadrature.size();
2579 
2580  switch (dim)
2581  {
2582  case 2:
2583  {
2584  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2585  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2586  ++q_point)
2587  {
2588  const Point<dim> quadrature_point(
2589  0.0, 0.5 * (edge_quadrature_points[q_point](0) + subface));
2590 
2591  interpolation_matrix(0, dof) +=
2592  0.5 * edge_quadrature.weight(q_point) *
2593  this->shape_value_component(dof, quadrature_point, 1);
2594  }
2595 
2596  if (source_fe.degree > 1)
2597  {
2598  const std::vector<Polynomials::Polynomial<double>>
2599  &legendre_polynomials =
2600  Polynomials::Legendre::generate_complete_basis(
2601  source_fe.degree - 1);
2602  FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
2603  source_fe.degree - 1);
2604 
2605  {
2606  FullMatrix<double> assembling_matrix(source_fe.degree - 1,
2607  n_edge_quadrature_points);
2608 
2609  for (unsigned int q_point = 0;
2610  q_point < n_edge_quadrature_points;
2611  ++q_point)
2612  {
2613  const double weight =
2614  std::sqrt(edge_quadrature.weight(q_point));
2615 
2616  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2617  assembling_matrix(i, q_point) =
2618  weight * legendre_polynomials[i + 1].value(
2619  edge_quadrature_points[q_point](0));
2620  }
2621 
2622  FullMatrix<double> system_matrix(source_fe.degree - 1,
2623  source_fe.degree - 1);
2624 
2625  assembling_matrix.mTmult(system_matrix, assembling_matrix);
2626  system_matrix_inv.invert(system_matrix);
2627  }
2628 
2629  Vector<double> solution(source_fe.degree - 1);
2630  Vector<double> system_rhs(source_fe.degree - 1);
2631 
2632  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2633  {
2634  system_rhs = 0.0;
2635 
2636  for (unsigned int q_point = 0;
2637  q_point < n_edge_quadrature_points;
2638  ++q_point)
2639  {
2640  const Point<dim> quadrature_point_0(
2641  0.0,
2642  0.5 * (edge_quadrature_points[q_point](0) + subface));
2643  const Point<dim> quadrature_point_1(
2644  0.0, edge_quadrature_points[q_point](0));
2645  const double tmp =
2646  edge_quadrature.weight(q_point) *
2647  (0.5 * this->shape_value_component(dof,
2648  quadrature_point_0,
2649  1) -
2650  interpolation_matrix(0, dof) *
2651  source_fe.shape_value_component(0,
2652  quadrature_point_1,
2653  1));
2654 
2655  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2656  system_rhs(i) +=
2657  tmp * legendre_polynomials[i + 1].value(
2658  edge_quadrature_points[q_point](0));
2659  }
2660 
2661  system_matrix_inv.vmult(solution, system_rhs);
2662 
2663  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2664  if (std::abs(solution(i)) > 1e-14)
2665  interpolation_matrix(i + 1, dof) = solution(i);
2666  }
2667  }
2668 
2669  break;
2670  }
2671 
2672  case 3:
2673  {
2674  const double shifts[4][2] = {{0.0, 0.0},
2675  {1.0, 0.0},
2676  {0.0, 1.0},
2677  {1.0, 1.0}};
2678 
2679  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2680  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2681  ++q_point)
2682  {
2683  const double weight = 0.5 * edge_quadrature.weight(q_point);
2684 
2685  for (unsigned int i = 0; i < 2; ++i)
2686  {
2687  Point<dim> quadrature_point(
2688  0.5 * (i + shifts[subface][0]),
2689  0.5 * (edge_quadrature_points[q_point](0) +
2690  shifts[subface][1]),
2691  0.0);
2692 
2693  interpolation_matrix(i * source_fe.degree, dof) +=
2694  weight *
2695  this->shape_value_component(
2696  this->face_to_cell_index(dof, 4), quadrature_point, 1);
2697  quadrature_point =
2698  Point<dim>(0.5 * (edge_quadrature_points[q_point](0) +
2699  shifts[subface][0]),
2700  0.5 * (i + shifts[subface][1]),
2701  0.0);
2702  interpolation_matrix((i + 2) * source_fe.degree, dof) +=
2703  weight *
2704  this->shape_value_component(
2705  this->face_to_cell_index(dof, 4), quadrature_point, 0);
2706  }
2707  }
2708 
2709  if (source_fe.degree > 1)
2710  {
2711  const std::vector<Polynomials::Polynomial<double>>
2712  &legendre_polynomials =
2713  Polynomials::Legendre::generate_complete_basis(
2714  source_fe.degree - 1);
2715  FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
2716  source_fe.degree - 1);
2717 
2718  {
2719  FullMatrix<double> assembling_matrix(source_fe.degree - 1,
2720  n_edge_quadrature_points);
2721 
2722  for (unsigned int q_point = 0;
2723  q_point < n_edge_quadrature_points;
2724  ++q_point)
2725  {
2726  const double weight =
2727  std::sqrt(edge_quadrature.weight(q_point));
2728 
2729  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2730  assembling_matrix(i, q_point) =
2731  weight * legendre_polynomials[i + 1].value(
2732  edge_quadrature_points[q_point](0));
2733  }
2734 
2735  FullMatrix<double> system_matrix(source_fe.degree - 1,
2736  source_fe.degree - 1);
2737 
2738  assembling_matrix.mTmult(system_matrix, assembling_matrix);
2739  system_matrix_inv.invert(system_matrix);
2740  }
2741 
2742  FullMatrix<double> solution(source_fe.degree - 1,
2743  GeometryInfo<dim>::lines_per_face);
2744  FullMatrix<double> system_rhs(source_fe.degree - 1,
2745  GeometryInfo<dim>::lines_per_face);
2746  Vector<double> tmp(GeometryInfo<dim>::lines_per_face);
2747 
2748  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2749  {
2750  system_rhs = 0.0;
2751 
2752  for (unsigned int q_point = 0;
2753  q_point < n_edge_quadrature_points;
2754  ++q_point)
2755  {
2756  const double weight = edge_quadrature.weight(q_point);
2757 
2758  for (unsigned int i = 0; i < 2; ++i)
2759  {
2760  Point<dim> quadrature_point_0(
2761  0.5 * (i + shifts[subface][0]),
2762  0.5 * (edge_quadrature_points[q_point](0) +
2763  shifts[subface][1]),
2764  0.0);
2765  Point<dim> quadrature_point_1(
2766  i, edge_quadrature_points[q_point](0), 0.0);
2767 
2768  tmp(i) =
2769  weight *
2770  (0.5 * this->shape_value_component(
2771  this->face_to_cell_index(dof, 4),
2772  quadrature_point_0,
2773  1) -
2774  interpolation_matrix(i * source_fe.degree, dof) *
2775  source_fe.shape_value_component(
2776  i * source_fe.degree, quadrature_point_1, 1));
2777  quadrature_point_0 =
2778  Point<dim>(0.5 *
2779  (edge_quadrature_points[q_point](0) +
2780  shifts[subface][0]),
2781  0.5 * (i + shifts[subface][1]),
2782  0.0);
2783  quadrature_point_1 =
2784  Point<dim>(edge_quadrature_points[q_point](0),
2785  i,
2786  0.0);
2787  tmp(i + 2) =
2788  weight *
2789  (0.5 * this->shape_value_component(
2790  this->face_to_cell_index(dof, 4),
2791  quadrature_point_0,
2792  0) -
2793  interpolation_matrix((i + 2) * source_fe.degree,
2794  dof) *
2795  source_fe.shape_value_component(
2796  (i + 2) * source_fe.degree,
2797  quadrature_point_1,
2798  0));
2799  }
2800 
2801  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2802  {
2803  const double L_i = legendre_polynomials[i + 1].value(
2804  edge_quadrature_points[q_point](0));
2805 
2806  for (unsigned int j = 0;
2807  j < GeometryInfo<dim>::lines_per_face;
2808  ++j)
2809  system_rhs(i, j) += tmp(j) * L_i;
2810  }
2811  }
2812 
2813  system_matrix_inv.mmult(solution, system_rhs);
2814 
2815  for (unsigned int i = 0;
2816  i < GeometryInfo<dim>::lines_per_face;
2817  ++i)
2818  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2819  if (std::abs(solution(j, i)) > 1e-14)
2820  interpolation_matrix(i * source_fe.degree + j + 1,
2821  dof) = solution(j, i);
2822  }
2823 
2824  const QGauss<2> quadrature(source_fe.degree);
2825  const std::vector<Point<2>> &quadrature_points =
2826  quadrature.get_points();
2827  const std::vector<Polynomials::Polynomial<double>>
2828  &lobatto_polynomials =
2829  Polynomials::Lobatto::generate_complete_basis(
2830  source_fe.degree);
2831  const unsigned int n_boundary_dofs =
2832  GeometryInfo<dim>::lines_per_face * source_fe.degree;
2833  const unsigned int &n_quadrature_points = quadrature.size();
2834 
2835  {
2836  FullMatrix<double> assembling_matrix(source_fe.degree *
2837  (source_fe.degree - 1),
2838  n_quadrature_points);
2839 
2840  for (unsigned int q_point = 0; q_point < n_quadrature_points;
2841  ++q_point)
2842  {
2843  const double weight = std::sqrt(quadrature.weight(q_point));
2844 
2845  for (unsigned int i = 0; i < source_fe.degree; ++i)
2846  {
2847  const double L_i =
2848  weight * legendre_polynomials[i].value(
2849  quadrature_points[q_point](0));
2850 
2851  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2852  assembling_matrix(i * (source_fe.degree - 1) + j,
2853  q_point) =
2854  L_i * lobatto_polynomials[j + 2].value(
2855  quadrature_points[q_point](1));
2856  }
2857  }
2858 
2859  FullMatrix<double> system_matrix(assembling_matrix.m(),
2860  assembling_matrix.m());
2861 
2862  assembling_matrix.mTmult(system_matrix, assembling_matrix);
2863  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
2864  system_matrix_inv.invert(system_matrix);
2865  }
2866 
2867  solution.reinit(system_matrix_inv.m(), 2);
2868  system_rhs.reinit(system_matrix_inv.m(), 2);
2869  tmp.reinit(2);
2870 
2871  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2872  {
2873  system_rhs = 0.0;
2874 
2875  for (unsigned int q_point = 0; q_point < n_quadrature_points;
2876  ++q_point)
2877  {
2878  Point<dim> quadrature_point(
2879  0.5 *
2880  (quadrature_points[q_point](0) + shifts[subface][0]),
2881  0.5 *
2882  (quadrature_points[q_point](1) + shifts[subface][1]),
2883  0.0);
2884  tmp(0) = 0.5 * this->shape_value_component(
2885  this->face_to_cell_index(dof, 4),
2886  quadrature_point,
2887  0);
2888  tmp(1) = 0.5 * this->shape_value_component(
2889  this->face_to_cell_index(dof, 4),
2890  quadrature_point,
2891  1);
2892  quadrature_point =
2893  Point<dim>(quadrature_points[q_point](0),
2894  quadrature_points[q_point](1),
2895  0.0);
2896 
2897  for (unsigned int i = 0; i < 2; ++i)
2898  for (unsigned int j = 0; j < source_fe.degree; ++j)
2899  {
2900  tmp(0) -= interpolation_matrix(
2901  (i + 2) * source_fe.degree + j, dof) *
2902  source_fe.shape_value_component(
2903  (i + 2) * source_fe.degree + j,
2904  quadrature_point,
2905  0);
2906  tmp(1) -=
2907  interpolation_matrix(i * source_fe.degree + j,
2908  dof) *
2909  source_fe.shape_value_component(
2910  i * source_fe.degree + j, quadrature_point, 1);
2911  }
2912 
2913  tmp *= quadrature.weight(q_point);
2914 
2915  for (unsigned int i = 0; i < source_fe.degree; ++i)
2916  {
2917  const double L_i_0 = legendre_polynomials[i].value(
2918  quadrature_points[q_point](0));
2919  const double L_i_1 = legendre_polynomials[i].value(
2920  quadrature_points[q_point](1));
2921 
2922  for (unsigned int j = 0; j < source_fe.degree - 1;
2923  ++j)
2924  {
2925  system_rhs(i * (source_fe.degree - 1) + j, 0) +=
2926  tmp(0) * L_i_0 *
2927  lobatto_polynomials[j + 2].value(
2928  quadrature_points[q_point](1));
2929  system_rhs(i * (source_fe.degree - 1) + j, 1) +=
2930  tmp(1) * L_i_1 *
2931  lobatto_polynomials[j + 2].value(
2932  quadrature_points[q_point](0));
2933  }
2934  }
2935  }
2936 
2937  system_matrix_inv.mmult(solution, system_rhs);
2938 
2939  for (unsigned int i = 0; i < source_fe.degree; ++i)
2940  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2941  {
2942  if (std::abs(solution(i * (source_fe.degree - 1) + j,
2943  0)) > 1e-14)
2944  interpolation_matrix(i * (source_fe.degree - 1) + j +
2945  n_boundary_dofs,
2946  dof) =
2947  solution(i * (source_fe.degree - 1) + j, 0);
2948 
2949  if (std::abs(solution(i * (source_fe.degree - 1) + j,
2950  1)) > 1e-14)
2951  interpolation_matrix(
2952  i + (j + source_fe.degree - 1) * source_fe.degree +
2953  n_boundary_dofs,
2954  dof) = solution(i * (source_fe.degree - 1) + j, 1);
2955  }
2956  }
2957  }
2958 
2959  break;
2960  }
2961 
2962  default:
2963  Assert(false, ExcNotImplemented());
2964  }
2965 }
2966 
2967 template <int dim>
2968 const FullMatrix<double> &
2969 FE_Nedelec<dim>::get_prolongation_matrix(
2970  const unsigned int child,
2971  const RefinementCase<dim> &refinement_case) const
2972 {
2973  Assert(refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
2974  ExcIndexRange(refinement_case,
2975  0,
2976  RefinementCase<dim>::isotropic_refinement + 1));
2977  Assert(refinement_case != RefinementCase<dim>::no_refinement,
2978  ExcMessage(
2979  "Prolongation matrices are only available for refined cells!"));
2980  Assert(child < GeometryInfo<dim>::n_children(refinement_case),
2981  ExcIndexRange(child,
2982  0,
2983  GeometryInfo<dim>::n_children(refinement_case)));
2984 
2985  // initialization upon first request
2986  if (this->prolongation[refinement_case - 1][child].n() == 0)
2987  {
2988  Threads::Mutex::ScopedLock lock(this->mutex);
2989 
2990  // if matrix got updated while waiting for the lock
2991  if (this->prolongation[refinement_case - 1][child].n() ==
2992  this->dofs_per_cell)
2993  return this->prolongation[refinement_case - 1][child];
2994 
2995  // now do the work. need to get a non-const version of data in order to
2996  // be able to modify them inside a const function
2997  FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
2998 
2999  // Reinit the vectors of
3000  // restriction and prolongation
3001  // matrices to the right sizes.
3002  // Restriction only for isotropic
3003  // refinement
3004 #ifdef DEBUG_NEDELEC
3005  deallog << "Embedding" << std::endl;
3006 #endif
3007  this_nonconst.reinit_restriction_and_prolongation_matrices();
3008  // Fill prolongation matrices with embedding operators
3009  FETools::compute_embedding_matrices(
3010  this_nonconst,
3011  this_nonconst.prolongation,
3012  true,
3013  internal::FE_Nedelec::get_embedding_computation_tolerance(
3014  this->degree));
3015 #ifdef DEBUG_NEDELEC
3016  deallog << "Restriction" << std::endl;
3017 #endif
3018  this_nonconst.initialize_restriction();
3019  }
3020 
3021  // we use refinement_case-1 here. the -1 takes care of the origin of the
3022  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3023  // available and so the vector indices are shifted
3024  return this->prolongation[refinement_case - 1][child];
3025 }
3026 
3027 template <int dim>
3028 const FullMatrix<double> &
3029 FE_Nedelec<dim>::get_restriction_matrix(
3030  const unsigned int child,
3031  const RefinementCase<dim> &refinement_case) const
3032 {
3033  Assert(refinement_case < RefinementCase<dim>::isotropic_refinement + 1,
3034  ExcIndexRange(refinement_case,
3035  0,
3036  RefinementCase<dim>::isotropic_refinement + 1));
3037  Assert(refinement_case != RefinementCase<dim>::no_refinement,
3038  ExcMessage(
3039  "Restriction matrices are only available for refined cells!"));
3040  Assert(child <
3041  GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
3042  ExcIndexRange(child,
3043  0,
3044  GeometryInfo<dim>::n_children(
3045  RefinementCase<dim>(refinement_case))));
3046 
3047  // initialization upon first request
3048  if (this->restriction[refinement_case - 1][child].n() == 0)
3049  {
3050  Threads::Mutex::ScopedLock lock(this->mutex);
3051 
3052  // if matrix got updated while waiting for the lock...
3053  if (this->restriction[refinement_case - 1][child].n() ==
3054  this->dofs_per_cell)
3055  return this->restriction[refinement_case - 1][child];
3056 
3057  // now do the work. need to get a non-const version of data in order to
3058  // be able to modify them inside a const function
3059  FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
3060 
3061  // Reinit the vectors of
3062  // restriction and prolongation
3063  // matrices to the right sizes.
3064  // Restriction only for isotropic
3065  // refinement
3066 #ifdef DEBUG_NEDELEC
3067  deallog << "Embedding" << std::endl;
3068 #endif
3069  this_nonconst.reinit_restriction_and_prolongation_matrices();
3070  // Fill prolongation matrices with embedding operators
3071  FETools::compute_embedding_matrices(
3072  this_nonconst,
3073  this_nonconst.prolongation,
3074  true,
3075  internal::FE_Nedelec::get_embedding_computation_tolerance(
3076  this->degree));
3077 #ifdef DEBUG_NEDELEC
3078  deallog << "Restriction" << std::endl;
3079 #endif
3080  this_nonconst.initialize_restriction();
3081  }
3082 
3083  // we use refinement_case-1 here. the -1 takes care of the origin of the
3084  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3085  // available and so the vector indices are shifted
3086  return this->restriction[refinement_case - 1][child];
3087 }
3088 
3089 
3090 // Interpolate a function, which is given by
3091 // its values at the generalized support
3092 // points in the finite element space on the
3093 // reference cell.
3094 // This is done as usual by projection-based
3095 // interpolation.
3096 template <int dim>
3097 void
3098 FE_Nedelec<dim>::convert_generalized_support_point_values_to_dof_values(
3099  const std::vector<Vector<double>> &support_point_values,
3100  std::vector<double> & nodal_values) const
3101 {
3102  const unsigned int deg = this->degree - 1;
3103  Assert(support_point_values.size() == this->generalized_support_points.size(),
3104  ExcDimensionMismatch(support_point_values.size(),
3105  this->generalized_support_points.size()));
3106  Assert(support_point_values[0].size() == this->n_components(),
3107  ExcDimensionMismatch(support_point_values[0].size(),
3108  this->n_components()));
3109  Assert(nodal_values.size() == this->dofs_per_cell,
3110  ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
3111  std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3112 
3113  switch (dim)
3114  {
3115  case 2:
3116  {
3117  // Let us begin with the
3118  // interpolation part.
3119  const QGauss<dim - 1> reference_edge_quadrature(this->degree);
3120  const unsigned int &n_edge_points = reference_edge_quadrature.size();
3121 
3122  for (unsigned int i = 0; i < 2; ++i)
3123  for (unsigned int j = 0; j < 2; ++j)
3124  {
3125  for (unsigned int q_point = 0; q_point < n_edge_points;
3126  ++q_point)
3127  nodal_values[(i + 2 * j) * this->degree] +=
3128  reference_edge_quadrature.weight(q_point) *
3129  support_point_values[q_point + (i + 2 * j) * n_edge_points]
3130  [1 - j];
3131 
3132  // Add the computed support_point_values to the resulting vector
3133  // only, if they are not
3134  // too small.
3135  if (std::abs(nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3136  nodal_values[(i + 2 * j) * this->degree] = 0.0;
3137  }
3138 
3139  // If the degree is greater
3140  // than 0, then we have still
3141  // some higher order edge
3142  // shape functions to
3143  // consider.
3144  // Here the projection part
3145  // starts. The dof support_point_values
3146  // are obtained by solving
3147  // a linear system of
3148  // equations.
3149  if (this->degree - 1 > 1)
3150  {
3151  // We start with projection
3152  // on the higher order edge
3153  // shape function.
3154  const std::vector<Polynomials::Polynomial<double>>
3155  &lobatto_polynomials =
3156  Polynomials::Lobatto::generate_complete_basis(this->degree);
3157  FullMatrix<double> system_matrix(this->degree - 1,
3158  this->degree - 1);
3159  std::vector<Polynomials::Polynomial<double>>
3160  lobatto_polynomials_grad(this->degree);
3161 
3162  for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3163  lobatto_polynomials_grad[i] =
3164  lobatto_polynomials[i + 1].derivative();
3165 
3166  // Set up the system matrix.
3167  // This can be used for all
3168  // edges.
3169  for (unsigned int i = 0; i < system_matrix.m(); ++i)
3170  for (unsigned int j = 0; j < system_matrix.n(); ++j)
3171  for (unsigned int q_point = 0; q_point < n_edge_points;
3172  ++q_point)
3173  system_matrix(i, j) +=
3174  boundary_weights(q_point, j) *
3175  lobatto_polynomials_grad[i + 1].value(
3176  this->generalized_face_support_points[q_point](0));
3177 
3178  FullMatrix<double> system_matrix_inv(this->degree - 1,
3179  this->degree - 1);
3180 
3181  system_matrix_inv.invert(system_matrix);
3182 
3183  const unsigned int
3184  line_coordinate[GeometryInfo<2>::lines_per_cell] = {1, 1, 0, 0};
3185  Vector<double> system_rhs(system_matrix.m());
3186  Vector<double> solution(system_rhs.size());
3187 
3188  for (unsigned int line = 0;
3189  line < GeometryInfo<dim>::lines_per_cell;
3190  ++line)
3191  {
3192  // Set up the right hand side.
3193  system_rhs = 0;
3194 
3195  for (unsigned int q_point = 0; q_point < n_edge_points;
3196  ++q_point)
3197  {
3198  const double tmp =
3199  support_point_values[line * n_edge_points + q_point]
3200  [line_coordinate[line]] -
3201  nodal_values[line * this->degree] *
3202  this->shape_value_component(
3203  line * this->degree,
3204  this->generalized_support_points[line *
3205  n_edge_points +
3206  q_point],
3207  line_coordinate[line]);
3208 
3209  for (unsigned int i = 0; i < system_rhs.size(); ++i)
3210  system_rhs(i) += boundary_weights(q_point, i) * tmp;
3211  }
3212 
3213  system_matrix_inv.vmult(solution, system_rhs);
3214 
3215  // Add the computed support_point_values
3216  // to the resulting vector
3217  // only, if they are not
3218  // too small.
3219  for (unsigned int i = 0; i < solution.size(); ++i)
3220  if (std::abs(solution(i)) > 1e-14)
3221  nodal_values[line * this->degree + i + 1] = solution(i);
3222  }
3223 
3224  // Then we go on to the
3225  // interior shape
3226  // functions. Again we
3227  // set up the system
3228  // matrix and use it
3229  // for both, the
3230  // horizontal and the
3231  // vertical, interior
3232  // shape functions.
3233  const QGauss<dim> reference_quadrature(this->degree);
3234  const unsigned int &n_interior_points =
3235  reference_quadrature.size();
3236  const std::vector<Polynomials::Polynomial<double>>
3237  &legendre_polynomials =
3238  Polynomials::Legendre::generate_complete_basis(this->degree -
3239  1);
3240 
3241  system_matrix.reinit((this->degree - 1) * this->degree,
3242  (this->degree - 1) * this->degree);
3243  system_matrix = 0;
3244 
3245  for (unsigned int i = 0; i < this->degree; ++i)
3246  for (unsigned int j = 0; j < this->degree - 1; ++j)
3247  for (unsigned int k = 0; k < this->degree; ++k)
3248  for (unsigned int l = 0; l < this->degree - 1; ++l)
3249  for (unsigned int q_point = 0;
3250  q_point < n_interior_points;
3251  ++q_point)
3252  system_matrix(i * (this->degree - 1) + j,
3253  k * (this->degree - 1) + l) +=
3254  reference_quadrature.weight(q_point) *
3255  legendre_polynomials[i].value(
3256  this->generalized_support_points
3257  [q_point + GeometryInfo<dim>::lines_per_cell *
3258  n_edge_points](0)) *
3259  lobatto_polynomials[j + 2].value(
3260  this->generalized_support_points
3261  [q_point + GeometryInfo<dim>::lines_per_cell *
3262  n_edge_points](1)) *
3263  lobatto_polynomials_grad[k].value(
3264  this->generalized_support_points
3265  [q_point + GeometryInfo<dim>::lines_per_cell *
3266  n_edge_points](0)) *
3267  lobatto_polynomials[l + 2].value(
3268  this->generalized_support_points
3269  [q_point + GeometryInfo<dim>::lines_per_cell *
3270  n_edge_points](1));
3271 
3272  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3273  system_matrix_inv.invert(system_matrix);
3274  // Set up the right hand side
3275  // for the horizontal shape
3276  // functions.
3277  system_rhs.reinit(system_matrix_inv.m());
3278  system_rhs = 0;
3279 
3280  for (unsigned int q_point = 0; q_point < n_interior_points;
3281  ++q_point)
3282  {
3283  double tmp =
3284  support_point_values[q_point +
3285  GeometryInfo<dim>::lines_per_cell *
3286  n_edge_points][0];
3287 
3288  for (unsigned int i = 0; i < 2; ++i)
3289  for (unsigned int j = 0; j <= deg; ++j)
3290  tmp -= nodal_values[(i + 2) * this->degree + j] *
3291  this->shape_value_component(
3292  (i + 2) * this->degree + j,
3293  this->generalized_support_points
3294  [q_point + GeometryInfo<dim>::lines_per_cell *
3295  n_edge_points],
3296  0);
3297 
3298  for (unsigned int i = 0; i <= deg; ++i)
3299  for (unsigned int j = 0; j < deg; ++j)
3300  system_rhs(i * deg + j) +=
3301  reference_quadrature.weight(q_point) * tmp *
3302  lobatto_polynomials_grad[i].value(
3303  this->generalized_support_points
3304  [q_point + GeometryInfo<dim>::lines_per_cell *
3305  n_edge_points](0)) *
3306  lobatto_polynomials[j + 2].value(
3307  this->generalized_support_points
3308  [q_point + GeometryInfo<dim>::lines_per_cell *
3309  n_edge_points](1));
3310  }
3311 
3312  solution.reinit(system_matrix.m());
3313  system_matrix_inv.vmult(solution, system_rhs);
3314 
3315  // Add the computed support_point_values
3316  // to the resulting vector
3317  // only, if they are not
3318  // too small.
3319  for (unsigned int i = 0; i <= deg; ++i)
3320  for (unsigned int j = 0; j < deg; ++j)
3321  if (std::abs(solution(i * deg + j)) > 1e-14)
3322  nodal_values[(i + GeometryInfo<dim>::lines_per_cell) * deg +
3323  j + GeometryInfo<dim>::lines_per_cell] =
3324  solution(i * deg + j);
3325 
3326  system_rhs = 0;
3327  // Set up the right hand side
3328  // for the vertical shape
3329  // functions.
3330 
3331  for (unsigned int q_point = 0; q_point < n_interior_points;
3332  ++q_point)
3333  {
3334  double tmp =
3335  support_point_values[q_point +
3336  GeometryInfo<dim>::lines_per_cell *
3337  n_edge_points][1];
3338 
3339  for (unsigned int i = 0; i < 2; ++i)
3340  for (unsigned int j = 0; j <= deg; ++j)
3341  tmp -= nodal_values[i * this->degree + j] *
3342  this->shape_value_component(
3343  i * this->degree + j,
3344  this->generalized_support_points
3345  [q_point + GeometryInfo<dim>::lines_per_cell *
3346  n_edge_points],
3347  1);
3348 
3349  for (unsigned int i = 0; i <= deg; ++i)
3350  for (unsigned int j = 0; j < deg; ++j)
3351  system_rhs(i * deg + j) +=
3352  reference_quadrature.weight(q_point) * tmp *
3353  lobatto_polynomials_grad[i].value(
3354  this->generalized_support_points
3355  [q_point + GeometryInfo<dim>::lines_per_cell *
3356  n_edge_points](1)) *
3357  lobatto_polynomials[j + 2].value(
3358  this->generalized_support_points
3359  [q_point + GeometryInfo<dim>::lines_per_cell *
3360  n_edge_points](0));
3361  }
3362 
3363  system_matrix_inv.vmult(solution, system_rhs);
3364 
3365  // Add the computed support_point_values
3366  // to the resulting vector
3367  // only, if they are not
3368  // too small.
3369  for (unsigned int i = 0; i <= deg; ++i)
3370  for (unsigned int j = 0; j < deg; ++j)
3371  if (std::abs(solution(i * deg + j)) > 1e-14)
3372  nodal_values[i +
3373  (j + GeometryInfo<dim>::lines_per_cell + deg) *
3374  this->degree] = solution(i * deg + j);
3375  }
3376 
3377  break;
3378  }
3379 
3380  case 3:
3381  {
3382  // Let us begin with the
3383  // interpolation part.
3384  const QGauss<1> reference_edge_quadrature(this->degree);
3385  const unsigned int &n_edge_points = reference_edge_quadrature.size();
3386 
3387  for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3388  {
3389  for (unsigned int i = 0; i < 4; ++i)
3390  nodal_values[(i + 8) * this->degree] +=
3391  reference_edge_quadrature.weight(q_point) *
3392  support_point_values[q_point + (i + 8) * n_edge_points][2];
3393 
3394  for (unsigned int i = 0; i < 2; ++i)
3395  for (unsigned int j = 0; j < 2; ++j)
3396  for (unsigned int k = 0; k < 2; ++k)
3397  nodal_values[(i + 2 * (2 * j + k)) * this->degree] +=
3398  reference_edge_quadrature.weight(q_point) *
3399  support_point_values[q_point + (i + 2 * (2 * j + k)) *
3400  n_edge_points][1 - k];
3401  }
3402 
3403  // Add the computed support_point_values
3404  // to the resulting vector
3405  // only, if they are not
3406  // too small.
3407  for (unsigned int i = 0; i < 4; ++i)
3408  if (std::abs(nodal_values[(i + 8) * this->degree]) < 1e-14)
3409  nodal_values[(i + 8) * this->degree] = 0.0;
3410 
3411  for (unsigned int i = 0; i < 2; ++i)
3412  for (unsigned int j = 0; j < 2; ++j)
3413  for (unsigned int k = 0; k < 2; ++k)
3414  if (std::abs(
3415  nodal_values[(i + 2 * (2 * j + k)) * this->degree]) <
3416  1e-14)
3417  nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3418 
3419  // If the degree is greater
3420  // than 0, then we have still
3421  // some higher order shape
3422  // functions to consider.
3423  // Here the projection part
3424  // starts. The dof support_point_values
3425  // are obtained by solving
3426  // a linear system of
3427  // equations.
3428  if (this->degree > 1)
3429  {
3430  // We start with projection
3431  // on the higher order edge
3432  // shape function.
3433  const std::vector<Polynomials::Polynomial<double>>
3434  &lobatto_polynomials =
3435  Polynomials::Lobatto::generate_complete_basis(this->degree);
3436  FullMatrix<double> system_matrix(this->degree - 1,
3437  this->degree - 1);
3438  std::vector<Polynomials::Polynomial<double>>
3439  lobatto_polynomials_grad(this->degree);
3440 
3441  for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3442  lobatto_polynomials_grad[i] =
3443  lobatto_polynomials[i + 1].derivative();
3444 
3445  // Set up the system matrix.
3446  // This can be used for all
3447  // edges.
3448  for (unsigned int i = 0; i < system_matrix.m(); ++i)
3449  for (unsigned int j = 0; j < system_matrix.n(); ++j)
3450  for (unsigned int q_point = 0; q_point < n_edge_points;
3451  ++q_point)
3452  system_matrix(i, j) +=
3453  boundary_weights(q_point, j) *
3454  lobatto_polynomials_grad[i + 1].value(
3455  this->generalized_face_support_points[q_point](1));
3456 
3457  FullMatrix<double> system_matrix_inv(this->degree - 1,
3458  this->degree - 1);
3459 
3460  system_matrix_inv.invert(system_matrix);
3461 
3462  const unsigned int
3463  line_coordinate[GeometryInfo<3>::lines_per_cell] = {
3464  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3465  Vector<double> system_rhs(system_matrix.m());
3466  Vector<double> solution(system_rhs.size());
3467 
3468  for (unsigned int line = 0;
3469  line < GeometryInfo<dim>::lines_per_cell;
3470  ++line)
3471  {
3472  // Set up the right hand side.
3473  system_rhs = 0;
3474 
3475  for (unsigned int q_point = 0; q_point < this->degree;
3476  ++q_point)
3477  {
3478  const double tmp =
3479  support_point_values[line * this->degree + q_point]
3480  [line_coordinate[line]] -
3481  nodal_values[line * this->degree] *
3482  this->shape_value_component(
3483  line * this->degree,
3484  this
3485  ->generalized_support_points[line * this->degree +
3486  q_point],
3487  line_coordinate[line]);
3488 
3489  for (unsigned int i = 0; i < system_rhs.size(); ++i)
3490  system_rhs(i) += boundary_weights(q_point, i) * tmp;
3491  }
3492 
3493  system_matrix_inv.vmult(solution, system_rhs);
3494 
3495  // Add the computed values
3496  // to the resulting vector
3497  // only, if they are not
3498  // too small.
3499  for (unsigned int i = 0; i < solution.size(); ++i)
3500  if (std::abs(solution(i)) > 1e-14)
3501  nodal_values[line * this->degree + i + 1] = solution(i);
3502  }
3503 
3504  // Then we go on to the
3505  // face shape functions.
3506  // Again we set up the
3507  // system matrix and
3508  // use it for both, the
3509  // horizontal and the
3510  // vertical, shape
3511  // functions.
3512  const std::vector<Polynomials::Polynomial<double>>
3513  &legendre_polynomials =
3514  Polynomials::Legendre::generate_complete_basis(this->degree -
3515  1);
3516  const unsigned int n_face_points = n_edge_points * n_edge_points;
3517 
3518  system_matrix.reinit((this->degree - 1) * this->degree,
3519  (this->degree - 1) * this->degree);
3520  system_matrix = 0;
3521 
3522  for (unsigned int i = 0; i < this->degree; ++i)
3523  for (unsigned int j = 0; j < this->degree - 1; ++j)
3524  for (unsigned int k = 0; k < this->degree; ++k)
3525  for (unsigned int l = 0; l < this->degree - 1; ++l)
3526  for (unsigned int q_point = 0; q_point < n_face_points;
3527  ++q_point)
3528  system_matrix(i * (this->degree - 1) + j,
3529  k * (this->degree - 1) + l) +=
3530  boundary_weights(q_point + n_edge_points,
3531  2 * (k * (this->degree - 1) + l)) *
3532  legendre_polynomials[i].value(
3533  this->generalized_face_support_points
3534  [q_point + 4 * n_edge_points](0)) *
3535  lobatto_polynomials[j + 2].value(
3536  this->generalized_face_support_points
3537  [q_point + 4 * n_edge_points](1));
3538 
3539  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3540  system_matrix_inv.invert(system_matrix);
3541  solution.reinit(system_matrix.m());
3542  system_rhs.reinit(system_matrix.m());
3543 
3544  const unsigned int
3545  face_coordinates[GeometryInfo<3>::faces_per_cell][2] = {
3546  {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3547  const unsigned int edge_indices[GeometryInfo<3>::faces_per_cell]
3548  [GeometryInfo<3>::lines_per_face] =
3549  {{0, 4, 8, 10},
3550  {1, 5, 9, 11},
3551  {8, 9, 2, 6},
3552  {10, 11, 3, 7},
3553  {2, 3, 0, 1},
3554  {6, 7, 4, 5}};
3555 
3556  for (unsigned int face = 0;
3557  face < GeometryInfo<dim>::faces_per_cell;
3558  ++face)
3559  {
3560  // Set up the right hand side
3561  // for the horizontal shape
3562  // functions.
3563  system_rhs = 0;
3564 
3565  for (unsigned int q_point = 0; q_point < n_face_points;
3566  ++q_point)
3567  {
3568  double tmp =
3569  support_point_values[q_point +
3570  GeometryInfo<dim>::lines_per_cell *
3571  n_edge_points]
3572  [face_coordinates[face][0]];
3573 
3574  for (unsigned int i = 0; i < 2; ++i)
3575  for (unsigned int j = 0; j <= deg; ++j)
3576  tmp -=
3577  nodal_values[edge_indices[face][i] * this->degree +
3578  j] *
3579  this->shape_value_component(
3580  edge_indices[face][i] * this->degree + j,
3581  this->generalized_support_points
3582  [q_point + GeometryInfo<dim>::lines_per_cell *
3583  n_edge_points],
3584  face_coordinates[face][0]);
3585 
3586  for (unsigned int i = 0; i <= deg; ++i)
3587  for (unsigned int j = 0; j < deg; ++j)
3588  system_rhs(i * deg + j) +=
3589  boundary_weights(q_point + n_edge_points,
3590  2 * (i * deg + j)) *
3591  tmp;
3592  }
3593 
3594  system_matrix_inv.vmult(solution, system_rhs);
3595 
3596  // Add the computed support_point_values
3597  // to the resulting vector
3598  // only, if they are not
3599  // too small.
3600  for (unsigned int i = 0; i <= deg; ++i)
3601  for (unsigned int j = 0; j < deg; ++j)
3602  if (std::abs(solution(i * deg + j)) > 1e-14)
3603  nodal_values[(2 * face * this->degree + i +
3604  GeometryInfo<dim>::lines_per_cell) *
3605  deg +
3606  j + GeometryInfo<dim>::lines_per_cell] =
3607  solution(i * deg + j);
3608 
3609  // Set up the right hand side
3610  // for the vertical shape
3611  // functions.
3612  system_rhs = 0;
3613 
3614  for (unsigned int q_point = 0; q_point < n_face_points;
3615  ++q_point)
3616  {
3617  double tmp =
3618  support_point_values[q_point +
3619  GeometryInfo<dim>::lines_per_cell *
3620  n_edge_points]
3621  [face_coordinates[face][1]];
3622 
3623  for (int i = 2;
3624  i < (int)GeometryInfo<dim>::lines_per_face;
3625  ++i)
3626  for (unsigned int j = 0; j <= deg; ++j)
3627  tmp -=
3628  nodal_values[edge_indices[face][i] * this->degree +
3629  j] *
3630  this->shape_value_component(
3631  edge_indices[face][i] * this->degree + j,
3632  this->generalized_support_points
3633  [q_point + GeometryInfo<dim>::lines_per_cell *
3634  n_edge_points],
3635  face_coordinates[face][1]);
3636 
3637  for (unsigned int i = 0; i <= deg; ++i)
3638  for (unsigned int j = 0; j < deg; ++j)
3639  system_rhs(i * deg + j) +=
3640  boundary_weights(q_point + n_edge_points,
3641  2 * (i * deg + j) + 1) *
3642  tmp;
3643  }
3644 
3645  system_matrix_inv.vmult(solution, system_rhs);
3646 
3647  // Add the computed support_point_values
3648  // to the resulting vector
3649  // only, if they are not
3650  // too small.
3651  for (unsigned int i = 0; i <= deg; ++i)
3652  for (unsigned int j = 0; j < deg; ++j)
3653  if (std::abs(solution(i * deg + j)) > 1e-14)
3654  nodal_values[((2 * face + 1) * deg + j +
3655  GeometryInfo<dim>::lines_per_cell) *
3656  this->degree +
3657  i] = solution(i * deg + j);
3658  }
3659 
3660  // Finally we project
3661  // the remaining parts
3662  // of the function on
3663  // the interior shape
3664  // functions.
3665  const QGauss<dim> reference_quadrature(this->degree);
3666  const unsigned int n_interior_points =
3667  reference_quadrature.size();
3668 
3669  // We create the
3670  // system matrix.
3671  system_matrix.reinit(this->degree * deg * deg,
3672  this->degree * deg * deg);
3673  system_matrix = 0;
3674 
3675  for (unsigned int i = 0; i <= deg; ++i)
3676  for (unsigned int j = 0; j < deg; ++j)
3677  for (unsigned int k = 0; k < deg; ++k)
3678  for (unsigned int l = 0; l <= deg; ++l)
3679  for (unsigned int m = 0; m < deg; ++m)
3680  for (unsigned int n = 0; n < deg; ++n)
3681  for (unsigned int q_point = 0;
3682  q_point < n_interior_points;
3683  ++q_point)
3684  system_matrix((i * deg + j) * deg + k,
3685  (l * deg + m) * deg + n) +=
3686  reference_quadrature.weight(q_point) *
3687  legendre_polynomials[i].value(
3688  this->generalized_support_points
3689  [q_point +
3690  GeometryInfo<dim>::lines_per_cell *
3691  n_edge_points +
3692  GeometryInfo<dim>::faces_per_cell *
3693  n_face_points](0)) *
3694  lobatto_polynomials[j + 2].value(
3695  this->generalized_support_points
3696  [q_point +
3697  GeometryInfo<dim>::lines_per_cell *
3698  n_edge_points +
3699  GeometryInfo<dim>::faces_per_cell *
3700  n_face_points](1)) *
3701  lobatto_polynomials[k + 2].value(
3702  this->generalized_support_points
3703  [q_point +
3704  GeometryInfo<dim>::lines_per_cell *
3705  n_edge_points +
3706  GeometryInfo<dim>::faces_per_cell *
3707  n_face_points](2)) *
3708  lobatto_polynomials_grad[l].value(
3709  this->generalized_support_points
3710  [q_point +
3711  GeometryInfo<dim>::lines_per_cell *
3712  n_edge_points +
3713  GeometryInfo<dim>::faces_per_cell *
3714  n_face_points](0)) *
3715  lobatto_polynomials[m + 2].value(
3716  this->generalized_support_points
3717  [q_point +
3718  GeometryInfo<dim>::lines_per_cell *
3719  n_edge_points +
3720  GeometryInfo<dim>::faces_per_cell *
3721  n_face_points](1)) *
3722  lobatto_polynomials[n + 2].value(
3723  this->generalized_support_points
3724  [q_point +
3725  GeometryInfo<dim>::lines_per_cell *
3726  n_edge_points +
3727  GeometryInfo<dim>::faces_per_cell *
3728  n_face_points](2));
3729 
3730  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3731  system_matrix_inv.invert(system_matrix);
3732  // Set up the right hand side.
3733  system_rhs.reinit(system_matrix.m());
3734  system_rhs = 0;
3735 
3736  for (unsigned int q_point = 0; q_point < n_interior_points;
3737  ++q_point)
3738  {
3739  double tmp =
3740  support_point_values[q_point +
3741  GeometryInfo<dim>::lines_per_cell *
3742  n_edge_points +
3743  GeometryInfo<dim>::faces_per_cell *
3744  n_face_points][0];
3745 
3746  for (unsigned int i = 0; i <= deg; ++i)
3747  {
3748  for (unsigned int j = 0; j < 2; ++j)
3749  for (unsigned int k = 0; k < 2; ++k)
3750  tmp -=
3751  nodal_values[i + (j + 4 * k + 2) * this->degree] *
3752  this->shape_value_component(
3753  i + (j + 4 * k + 2) * this->degree,
3754  this->generalized_support_points
3755  [q_point +
3756  GeometryInfo<dim>::lines_per_cell *
3757  n_edge_points +
3758  GeometryInfo<dim>::faces_per_cell *
3759  n_face_points],
3760  0);
3761 
3762  for (unsigned int j = 0; j < deg; ++j)
3763  for (unsigned int k = 0; k < 4; ++k)
3764  tmp -=
3765  nodal_values[(i + 2 * (k + 2) * this->degree +
3766  GeometryInfo<dim>::lines_per_cell) *
3767  deg +
3768  j +
3769  GeometryInfo<dim>::lines_per_cell] *
3770  this->shape_value_component(
3771  (i + 2 * (k + 2) * this->degree +
3772  GeometryInfo<dim>::lines_per_cell) *
3773  deg +
3774  j + GeometryInfo<dim>::lines_per_cell,
3775  this->generalized_support_points
3776  [q_point +
3777  GeometryInfo<dim>::lines_per_cell *
3778  n_edge_points +
3779  GeometryInfo<dim>::faces_per_cell *
3780  n_face_points],
3781  0);
3782  }
3783 
3784  for (unsigned int i = 0; i <= deg; ++i)
3785  for (unsigned int j = 0; j < deg; ++j)
3786  for (unsigned int k = 0; k < deg; ++k)
3787  system_rhs((i * deg + j) * deg + k) +=
3788  reference_quadrature.weight(q_point) * tmp *
3789  lobatto_polynomials_grad[i].value(
3790  this->generalized_support_points
3791  [q_point +
3792  GeometryInfo<dim>::lines_per_cell *
3793  n_edge_points +
3794  GeometryInfo<dim>::faces_per_cell *
3795  n_face_points](0)) *
3796  lobatto_polynomials[j + 2].value(
3797  this->generalized_support_points
3798  [q_point +
3799  GeometryInfo<dim>::lines_per_cell *
3800  n_edge_points +
3801  GeometryInfo<dim>::faces_per_cell *
3802  n_face_points](1)) *
3803  lobatto_polynomials[k + 2].value(
3804  this->generalized_support_points
3805  [q_point +
3806  GeometryInfo<dim>::lines_per_cell *
3807  n_edge_points +
3808  GeometryInfo<dim>::faces_per_cell *
3809  n_face_points](2));
3810  }
3811 
3812  solution.reinit(system_rhs.size());
3813  system_matrix_inv.vmult(solution, system_rhs);
3814 
3815  // Add the computed values
3816  // to the resulting vector
3817  // only, if they are not
3818  // too small.
3819  for (unsigned int i = 0; i <= deg; ++i)
3820  for (unsigned int j = 0; j < deg; ++j)
3821  for (unsigned int k = 0; k < deg; ++k)
3822  if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3823  nodal_values
3824  [((i + 2 * GeometryInfo<dim>::faces_per_cell) * deg +
3825  j + GeometryInfo<dim>::lines_per_cell +
3826  2 * GeometryInfo<dim>::faces_per_cell) *
3827  deg +
3828  k + GeometryInfo<dim>::lines_per_cell] =
3829  solution((i * deg + j) * deg + k);
3830 
3831  // Set up the right hand side.
3832  system_rhs = 0;
3833 
3834  for (unsigned int q_point = 0; q_point < n_interior_points;
3835  ++q_point)
3836  {
3837  double tmp =
3838  support_point_values[q_point +
3839  GeometryInfo<dim>::lines_per_cell *
3840  n_edge_points +
3841  GeometryInfo<dim>::faces_per_cell *
3842  n_face_points][1];
3843 
3844  for (unsigned int i = 0; i <= deg; ++i)
3845  for (unsigned int j = 0; j < 2; ++j)
3846  {
3847  for (unsigned int k = 0; k < 2; ++k)
3848  tmp -= nodal_values[i + (4 * j + k) * this->degree] *
3849  this->shape_value_component(
3850  i + (4 * j + k) * this->degree,
3851  this->generalized_support_points
3852  [q_point +
3853  GeometryInfo<dim>::lines_per_cell *
3854  n_edge_points +
3855  GeometryInfo<dim>::faces_per_cell *
3856  n_face_points],
3857  1);
3858 
3859  for (unsigned int k = 0; k < deg; ++k)
3860  tmp -=
3861  nodal_values[(i + 2 * j * this->degree +
3862  GeometryInfo<dim>::lines_per_cell) *
3863  deg +
3864  k +
3865  GeometryInfo<dim>::lines_per_cell] *
3866  this->shape_value_component(
3867  (i + 2 * j * this->degree +
3868  GeometryInfo<dim>::lines_per_cell) *
3869  deg +
3870  k + GeometryInfo<dim>::lines_per_cell,
3871  this->generalized_support_points
3872  [q_point +
3873  GeometryInfo<dim>::lines_per_cell *
3874  n_edge_points +
3875  GeometryInfo<dim>::faces_per_cell *
3876  n_face_points],
3877  1) +
3878  nodal_values[i +
3879  ((2 * j + 9) * deg + k +
3880  GeometryInfo<dim>::lines_per_cell) *
3881  this->degree] *
3882  this->shape_value_component(
3883  i + ((2 * j + 9) * deg + k +
3884  GeometryInfo<dim>::lines_per_cell) *
3885  this->degree,
3886  this->generalized_support_points
3887  [q_point +
3888  GeometryInfo<dim>::lines_per_cell *
3889  n_edge_points +
3890  GeometryInfo<dim>::faces_per_cell *
3891  n_face_points],
3892  1);
3893  }
3894 
3895  for (unsigned int i = 0; i <= deg; ++i)
3896  for (unsigned int j = 0; j < deg; ++j)
3897  for (unsigned int k = 0; k < deg; ++k)
3898  system_rhs((i * deg + j) * deg + k) +=
3899  reference_quadrature.weight(q_point) * tmp *
3900  lobatto_polynomials_grad[i].value(
3901  this->generalized_support_points
3902  [q_point +
3903  GeometryInfo<dim>::lines_per_cell *
3904  n_edge_points +
3905  GeometryInfo<dim>::faces_per_cell *
3906  n_face_points](1)) *
3907  lobatto_polynomials[j + 2].value(
3908  this->generalized_support_points
3909  [q_point +
3910  GeometryInfo<dim>::lines_per_cell *
3911  n_edge_points +
3912  GeometryInfo<dim>::faces_per_cell *
3913  n_face_points](0)) *
3914  lobatto_polynomials[k + 2].value(
3915  this->generalized_support_points
3916  [q_point +
3917  GeometryInfo<dim>::lines_per_cell *
3918  n_edge_points +
3919  GeometryInfo<dim>::faces_per_cell *
3920  n_face_points](2));
3921  }
3922 
3923  system_matrix_inv.vmult(solution, system_rhs);
3924 
3925  // Add the computed support_point_values
3926  // to the resulting vector
3927  // only, if they are not
3928  // too small.
3929  for (unsigned int i = 0; i <= deg; ++i)
3930  for (unsigned int j = 0; j < deg; ++j)
3931  for (unsigned int k = 0; k < deg; ++k)
3932  if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3933  nodal_values[((i + this->degree +
3934  2 * GeometryInfo<dim>::faces_per_cell) *
3935  deg +
3936  j + GeometryInfo<dim>::lines_per_cell +
3937  2 * GeometryInfo<dim>::faces_per_cell) *
3938  deg +
3939  k + GeometryInfo<dim>::lines_per_cell] =
3940  solution((i * deg + j) * deg + k);
3941 
3942  // Set up the right hand side.
3943  system_rhs = 0;
3944 
3945  for (unsigned int q_point = 0; q_point < n_interior_points;
3946  ++q_point)
3947  {
3948  double tmp =
3949  support_point_values[q_point +
3950  GeometryInfo<dim>::lines_per_cell *
3951  n_edge_points +
3952  GeometryInfo<dim>::faces_per_cell *
3953  n_face_points][2];
3954 
3955  for (unsigned int i = 0; i <= deg; ++i)
3956  for (unsigned int j = 0; j < 4; ++j)
3957  {
3958  tmp -= nodal_values[i + (j + 8) * this->degree] *
3959  this->shape_value_component(
3960  i + (j + 8) * this->degree,
3961  this->generalized_support_points
3962  [q_point +
3963  GeometryInfo<dim>::lines_per_cell *
3964  n_edge_points +
3965  GeometryInfo<dim>::faces_per_cell *
3966  n_face_points],
3967  2);
3968 
3969  for (unsigned int k = 0; k < deg; ++k)
3970  tmp -=
3971  nodal_values[i +
3972  ((2 * j + 1) * deg + k +
3973  GeometryInfo<dim>::lines_per_cell) *
3974  this->degree] *
3975  this->shape_value_component(
3976  i + ((2 * j + 1) * deg + k +
3977  GeometryInfo<dim>::lines_per_cell) *
3978  this->degree,
3979  this->generalized_support_points
3980  [q_point +
3981  GeometryInfo<dim>::lines_per_cell *
3982  n_edge_points +
3983  GeometryInfo<dim>::faces_per_cell *
3984  n_face_points],
3985  2);
3986  }
3987 
3988  for (unsigned int i = 0; i <= deg; ++i)
3989  for (unsigned int j = 0; j < deg; ++j)
3990  for (unsigned int k = 0; k < deg; ++k)
3991  system_rhs((i * deg + j) * deg + k) +=
3992  reference_quadrature.weight(q_point) * tmp *
3993  lobatto_polynomials_grad[i].value(
3994  this->generalized_support_points
3995  [q_point +
3996  GeometryInfo<dim>::lines_per_cell *
3997  n_edge_points +
3998  GeometryInfo<dim>::faces_per_cell *
3999  n_face_points](2)) *
4000  lobatto_polynomials[j + 2].value(
4001  this->generalized_support_points
4002  [q_point +
4003  GeometryInfo<dim>::lines_per_cell *
4004  n_edge_points +
4005  GeometryInfo<dim>::faces_per_cell *
4006  n_face_points](0)) *
4007  lobatto_polynomials[k + 2].value(
4008  this->generalized_support_points
4009  [q_point +
4010  GeometryInfo<dim>::lines_per_cell *
4011  n_edge_points +
4012  GeometryInfo<dim>::faces_per_cell *
4013  n_face_points](1));
4014  }
4015 
4016  system_matrix_inv.vmult(solution, system_rhs);
4017 
4018  // Add the computed support_point_values
4019  // to the resulting vector
4020  // only, if they are not
4021  // too small.
4022  for (unsigned int i = 0; i <= deg; ++i)
4023  for (unsigned int j = 0; j < deg; ++j)
4024  for (unsigned int k = 0; k < deg; ++k)
4025  if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4026  nodal_values
4027  [i +
4028  ((j + 2 * (deg + GeometryInfo<dim>::faces_per_cell)) *
4029  deg +
4030  k + GeometryInfo<dim>::lines_per_cell) *
4031  this->degree] = solution((i * deg + j) * deg + k);
4032  }
4033 
4034  break;
4035  }
4036 
4037  default:
4038  Assert(false, ExcNotImplemented());
4039  }
4040 }
4041 
4042 
4043 
4044 template <int dim>
4045 std::pair<Table<2, bool>, std::vector<unsigned int>>
4046 FE_Nedelec<dim>::get_constant_modes() const
4047 {
4048  Table<2, bool> constant_modes(dim, this->dofs_per_cell);
4049  for (unsigned int d = 0; d < dim; ++d)
4050  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
4051  constant_modes(d, i) = true;
4052  std::vector<unsigned int> components;
4053  for (unsigned int d = 0; d < dim; ++d)
4054  components.push_back(d);
4055  return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
4056  components);
4057 }
4058 
4059 
4060 template <int dim>
4061 std::size_t
4062 FE_Nedelec<dim>::memory_consumption() const
4063 {
4064  Assert(false, ExcNotImplemented());
4065  return 0;
4066 }
4067 
4068 
4069 //----------------------------------------------------------------------//
4070 
4071 
4072 // explicit instantiations
4073 #include "fe_nedelec.inst"
4074 
4075 
4076 DEAL_II_NAMESPACE_CLOSE
FE_Nedelec
Definition: fe_nedelec.h:113
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:2538
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:3029
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:2436
FiniteElement< dim, dim >::restriction
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2419
FE_Nedelec::hp_constraints_are_implemented
virtual bool hp_constraints_are_implemented() const override
Definition: fe_nedelec.cc:2324
FiniteElementData::components
const unsigned int components
Definition: fe_base.h:307
FE_Nedelec::FE_Nedelec
friend class FE_Nedelec
Definition: fe_nedelec.h:335
FiniteElement< dim, dim >::generalized_support_points
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2470
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:2445
Polynomials::Lobatto::generate_complete_basis
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:971
Polynomials::Legendre::generate_complete_basis
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:866
Vector::size
size_type size() const
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:313
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:2969
FullMatrix::vmult
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
FE_Nedelec::get_constant_modes
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_nedelec.cc:4046
FiniteElement< dim, dim >::generalized_face_support_points
std::vector< Point< dim-1 > > generalized_face_support_points
Definition: fe.h:2476
StandardExceptions::ExcIndexRange
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Quadrature::point
const Point< dim > & point(const unsigned int i) const
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:3098
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:2340
FE_Nedelec::initialize_support_points
void initialize_support_points(const unsigned int order)
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
Quadrature::size
unsigned int size() const
StandardExceptions::ExcMessage
static::ExceptionBase & ExcMessage(std::string arg1)
FullMatrix::n
size_type n() const
StandardExceptions::ExcImpossibleInDim
static::ExceptionBase & ExcImpossibleInDim(int arg1)
FiniteElement< dim, dim >::prolongation
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2433
FullMatrix::invert
void invert(const FullMatrix< number2 > &M)
FE_Nedelec::boundary_weights
Table< 2, double > boundary_weights
Definition: fe_nedelec.h:324
TableBase::reinit
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1227
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:2331
StandardExceptions::ExcDimensionMismatch
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Quadrature::weight
double weight(const unsigned int i) const
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
FullMatrix::mTmult
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
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:2382
FullMatrix::mmult
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FE_Nothing< dim >
FiniteElementData::n_components
unsigned int n_components() const
FullMatrix::m
size_type m() const
PolynomialsNedelec::compute_n_pols
static unsigned int compute_n_pols(unsigned int degree)
Definition: polynomials_nedelec.cc:1487
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_base.h:297
FE_Nedelec::compare_for_face_domination
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim > &fe_other) const override
Definition: fe_nedelec.cc:2277
QProjector::project_to_all_faces
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
FiniteElement< dim, dim >::face_to_cell_index
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
FE_Nedelec::memory_consumption
virtual std::size_t memory_consumption() const override
Definition: fe_nedelec.cc:4062
FE_Nedelec::get_name
virtual std::string get_name() const override
Definition: fe_nedelec.cc:207
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:290
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)
StandardExceptions::ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
Table
Definition: table.h:37
FE_Nedelec::initialize_restriction
void initialize_restriction()
Definition: fe_nedelec.cc:516
FE_PolyTensor< PolynomialsNedelec< dim >, dim >::inverse_node_matrix
FullMatrix< double > inverse_node_matrix
Definition: fe_poly_tensor.h:454
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)
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.11