Reference documentation for deal.II version Git a39983f358 2020-01-20 20:47:10 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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<dim>(
73  PolynomialsNedelec<dim>(order),
74  FiniteElementData<dim>(get_dpo_vector(order),
75  dim,
76  order + 1,
77  FiniteElementData<dim>::Hcurl),
78  std::vector<bool>(PolynomialsNedelec<dim>::n_polynomials(order), true),
79  std::vector<ComponentMask>(PolynomialsNedelec<dim>::n_polynomials(order),
80  std::vector<bool>(dim, true)))
81 {
82 #ifdef DEBUG_NEDELEC
83  deallog << get_name() << std::endl;
84 #endif
85 
86  Assert(dim >= 2, ExcImpossibleInDim(dim));
87 
88  const unsigned int n_dofs = this->dofs_per_cell;
89 
90  this->mapping_kind = {mapping_nedelec};
91  // First, initialize the
92  // generalized support points and
93  // quadrature weights, since they
94  // are required for interpolation.
96  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
98  // From now on, the shape functions
99  // will be the correct ones, not
100  // the raw shape functions anymore.
101 
102  // do not initialize embedding and restriction here. these matrices are
103  // initialized on demand in get_restriction_matrix and
104  // get_prolongation_matrix
105 
106 #ifdef DEBUG_NEDELEC
107  deallog << "Face Embedding" << std::endl;
108 #endif
110 
111  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  {
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
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>>
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
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 =
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] =
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] =
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 =
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(
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 +
451  n_edge_points] =
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(
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
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
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 =
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 =
765  const unsigned int n_boundary_dofs =
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 =
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 =
1262  const unsigned int n_edge_dofs =
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>::n_polynomials(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  AssertIndexRange(shape_index, this->dofs_per_cell);
2041 
2042  const unsigned int deg = this->degree - 1;
2043  switch (dim)
2044  {
2045  case 2:
2046  switch (face_index)
2047  {
2048  case 0:
2049  if (!((shape_index > deg) && (shape_index < 2 * this->degree)))
2050  return true;
2051 
2052  else
2053  return false;
2054 
2055  case 1:
2056  if ((shape_index > deg) &&
2057  (shape_index <
2058  GeometryInfo<2>::lines_per_cell * this->degree))
2059  return true;
2060 
2061  else
2062  return false;
2063 
2064  case 2:
2065  if (shape_index < 3 * this->degree)
2066  return true;
2067 
2068  else
2069  return false;
2070 
2071  case 3:
2072  if (!((shape_index >= 2 * this->degree) &&
2073  (shape_index < 3 * this->degree)))
2074  return true;
2075 
2076  else
2077  return false;
2078 
2079  default:
2080  {
2081  Assert(false, ExcNotImplemented());
2082  return false;
2083  }
2084  }
2085 
2086  case 3:
2087  switch (face_index)
2088  {
2089  case 0:
2090  if (((shape_index > deg) && (shape_index < 2 * this->degree)) ||
2091  ((shape_index >= 5 * this->degree) &&
2092  (shape_index < 6 * this->degree)) ||
2093  ((shape_index >= 9 * this->degree) &&
2094  (shape_index < 10 * this->degree)) ||
2095  ((shape_index >= 11 * this->degree) &&
2096  (shape_index <
2097  GeometryInfo<3>::lines_per_cell * this->degree)) ||
2098  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2099  this->degree) &&
2100  (shape_index < (GeometryInfo<3>::lines_per_cell + 4 * deg) *
2101  this->degree)) ||
2102  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2103  this->degree) &&
2104  (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2105  this->degree)) ||
2106  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2107  this->degree) &&
2108  (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2109  this->degree)) ||
2110  ((shape_index >=
2111  (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2112  this->degree) &&
2113  (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2114  this->degree)))
2115  return false;
2116 
2117  else
2118  return true;
2119 
2120  case 1:
2121  if (((shape_index > deg) && (shape_index < 4 * this->degree)) ||
2122  ((shape_index >= 5 * this->degree) &&
2123  (shape_index < 8 * this->degree)) ||
2124  ((shape_index >= 9 * this->degree) &&
2125  (shape_index < 10 * this->degree)) ||
2126  ((shape_index >= 11 * this->degree) &&
2127  (shape_index <
2128  GeometryInfo<3>::lines_per_cell * this->degree)) ||
2129  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2130  this->degree) &&
2131  (shape_index < (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2132  this->degree)) ||
2133  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2134  this->degree) &&
2135  (shape_index < (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2136  this->degree)) ||
2137  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2138  this->degree) &&
2139  (shape_index < (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2140  this->degree)) ||
2141  ((shape_index >=
2142  (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2143  this->degree) &&
2144  (shape_index < (GeometryInfo<3>::lines_per_cell + 12 * deg) *
2145  this->degree)))
2146  return true;
2147 
2148  else
2149  return false;
2150 
2151  case 2:
2152  if ((shape_index < 3 * this->degree) ||
2153  ((shape_index >= 4 * this->degree) &&
2154  (shape_index < 7 * this->degree)) ||
2155  ((shape_index >= 8 * this->degree) &&
2156  (shape_index < 10 * this->degree)) ||
2157  ((shape_index >=
2158  (GeometryInfo<3>::lines_per_cell + deg) * this->degree) &&
2159  (shape_index < (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2160  this->degree)) ||
2161  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2162  this->degree) &&
2163  (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2164  this->degree)) ||
2165  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 8 * deg) *
2166  this->degree) &&
2167  (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2168  this->degree)) ||
2169  ((shape_index >=
2170  (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2171  this->degree) &&
2172  (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2173  this->degree)))
2174  return true;
2175 
2176  else
2177  return false;
2178 
2179  case 3:
2180  if ((shape_index < 2 * this->degree) ||
2181  ((shape_index >= 3 * this->degree) &&
2182  (shape_index < 6 * this->degree)) ||
2183  ((shape_index >= 7 * this->degree) &&
2184  (shape_index < 8 * this->degree)) ||
2185  ((shape_index >= 10 * this->degree) &&
2186  (shape_index <
2187  GeometryInfo<3>::lines_per_cell * this->degree)) ||
2188  ((shape_index >=
2189  (GeometryInfo<3>::lines_per_cell + deg) * this->degree) &&
2190  (shape_index < (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2191  this->degree)) ||
2192  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2193  this->degree) &&
2194  (shape_index < (GeometryInfo<3>::lines_per_cell + 4 * deg) *
2195  this->degree)) ||
2196  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2197  this->degree) &&
2198  (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2199  this->degree)) ||
2200  ((shape_index >=
2201  (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2202  this->degree) &&
2203  (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2204  this->degree)))
2205  return true;
2206 
2207  else
2208  return false;
2209 
2210  case 4:
2211  if ((shape_index < 4 * this->degree) ||
2212  ((shape_index >= 8 * this->degree) &&
2213  (shape_index <
2214  (GeometryInfo<3>::lines_per_cell + deg) * this->degree)) ||
2215  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2216  this->degree) &&
2217  (shape_index < (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2218  this->degree)) ||
2219  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2220  this->degree) &&
2221  (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2222  this->degree)) ||
2223  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2224  this->degree) &&
2225  (shape_index < (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2226  this->degree)))
2227  return true;
2228 
2229  else
2230  return false;
2231 
2232  case 5:
2233  if (((shape_index >= 4 * this->degree) &&
2234  (shape_index <
2235  (GeometryInfo<3>::lines_per_cell + deg) * this->degree)) ||
2236  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2237  this->degree) &&
2238  (shape_index < (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2239  this->degree)) ||
2240  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2241  this->degree) &&
2242  (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2243  this->degree)) ||
2244  ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2245  this->degree) &&
2246  (shape_index < (GeometryInfo<3>::lines_per_cell + 8 * deg) *
2247  this->degree)) ||
2248  ((shape_index >=
2249  (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2250  this->degree) &&
2251  (shape_index < (GeometryInfo<3>::lines_per_cell + 12 * deg) *
2252  this->degree)))
2253  return true;
2254 
2255  else
2256  return false;
2257 
2258  default:
2259  {
2260  Assert(false, ExcNotImplemented());
2261  return false;
2262  }
2263  }
2264 
2265  default:
2266  {
2267  Assert(false, ExcNotImplemented());
2268  return false;
2269  }
2270  }
2271 }
2272 
2273 template <int dim>
2276  const unsigned int codim) const
2277 {
2278  Assert(codim <= dim, ExcImpossibleInDim(dim));
2279  (void)codim;
2280 
2281  // vertex/line/face/cell domination
2282  // --------------------------------
2283  if (const FE_Nedelec<dim> *fe_nedelec_other =
2284  dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2285  {
2286  if (this->degree < fe_nedelec_other->degree)
2288  else if (this->degree == fe_nedelec_other->degree)
2290  else
2292  }
2293  else if (const FE_Nothing<dim> *fe_nothing =
2294  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
2295  {
2296  if (fe_nothing->is_dominating())
2298  else
2299  // the FE_Nothing has no degrees of freedom and it is typically used
2300  // in a context where we don't require any continuity along the
2301  // interface
2303  }
2304 
2305  Assert(false, ExcNotImplemented());
2307 }
2308 
2309 template <int dim>
2310 bool
2312 {
2313  return true;
2314 }
2315 
2316 template <int dim>
2317 std::vector<std::pair<unsigned int, unsigned int>>
2319 {
2320  // Nedelec elements do not have any dofs
2321  // on vertices, hence return an empty vector.
2322  return std::vector<std::pair<unsigned int, unsigned int>>();
2323 }
2324 
2325 template <int dim>
2326 std::vector<std::pair<unsigned int, unsigned int>>
2328  const FiniteElement<dim> &fe_other) const
2329 {
2330  // we can presently only compute these
2331  // identities if both FEs are
2332  // FE_Nedelec or if the other one is an
2333  // FE_Nothing
2334  if (const FE_Nedelec<dim> *fe_nedelec_other =
2335  dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2336  {
2337  // dofs are located on lines, so
2338  // two dofs are identical, if their
2339  // edge shape functions have the
2340  // same polynomial degree.
2341  std::vector<std::pair<unsigned int, unsigned int>> identities;
2342 
2343  for (unsigned int i = 0;
2344  i < std::min(fe_nedelec_other->degree, this->degree);
2345  ++i)
2346  identities.emplace_back(i, i);
2347 
2348  return identities;
2349  }
2350 
2351  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2352  {
2353  // the FE_Nothing has no
2354  // degrees of freedom, so there
2355  // are no equivalencies to be
2356  // recorded
2357  return std::vector<std::pair<unsigned int, unsigned int>>();
2358  }
2359 
2360  else
2361  {
2362  Assert(false, ExcNotImplemented());
2363  return std::vector<std::pair<unsigned int, unsigned int>>();
2364  }
2365 }
2366 
2367 template <int dim>
2368 std::vector<std::pair<unsigned int, unsigned int>>
2370  const FiniteElement<dim> &fe_other) const
2371 {
2372  // we can presently only compute
2373  // these identities if both FEs are
2374  // FE_Nedelec or if the other one is an
2375  // FE_Nothing
2376  if (const FE_Nedelec<dim> *fe_nedelec_other =
2377  dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2378  {
2379  // dofs are located on the interior
2380  // of faces, so two dofs are identical,
2381  // if their face shape functions have
2382  // the same polynomial degree.
2383  const unsigned int p = fe_nedelec_other->degree;
2384  const unsigned int q = this->degree;
2385  const unsigned int p_min = std::min(p, q);
2386  std::vector<std::pair<unsigned int, unsigned int>> identities;
2387 
2388  for (unsigned int i = 0; i < p_min; ++i)
2389  for (unsigned int j = 0; j < p_min - 1; ++j)
2390  {
2391  identities.emplace_back(i * (q - 1) + j, i * (p - 1) + j);
2392  identities.emplace_back(i + (j + q - 1) * q, i + (j + p - 1) * p);
2393  }
2394 
2395  return identities;
2396  }
2397 
2398  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2399  {
2400  // the FE_Nothing has no
2401  // degrees of freedom, so there
2402  // are no equivalencies to be
2403  // recorded
2404  return std::vector<std::pair<unsigned int, unsigned int>>();
2405  }
2406 
2407  else
2408  {
2409  Assert(false, ExcNotImplemented());
2410  return std::vector<std::pair<unsigned int, unsigned int>>();
2411  }
2412 }
2413 
2414 // In this function we compute the face
2415 // interpolation matrix. This is usually
2416 // done by projection-based interpolation,
2417 // but, since one can compute the entries
2418 // easy per hand, we save some computation
2419 // time at this point and just fill in the
2420 // correct values.
2421 template <int dim>
2422 void
2424  const FiniteElement<dim> &source,
2425  FullMatrix<double> & interpolation_matrix) const
2426 {
2427  // this is only implemented, if the
2428  // source FE is also a
2429  // Nedelec element
2430  AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2431  (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2433  Assert(interpolation_matrix.m() == source.dofs_per_face,
2434  ExcDimensionMismatch(interpolation_matrix.m(), source.dofs_per_face));
2435  Assert(interpolation_matrix.n() == this->dofs_per_face,
2436  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
2437 
2438  // ok, source is a Nedelec element, so
2439  // we will be able to do the work
2440  const FE_Nedelec<dim> &source_fe =
2441  dynamic_cast<const FE_Nedelec<dim> &>(source);
2442 
2443  // Make sure, that the element,
2444  // for which the DoFs should be
2445  // constrained is the one with
2446  // the higher polynomial degree.
2447  // Actually the procedure will work
2448  // also if this assertion is not
2449  // satisfied. But the matrices
2450  // produced in that case might
2451  // lead to problems in the
2452  // hp procedures, which use this
2453  // method.
2454  Assert(this->dofs_per_face <= source_fe.dofs_per_face,
2456  interpolation_matrix = 0;
2457 
2458  // On lines we can just identify
2459  // all degrees of freedom.
2460  for (unsigned int i = 0; i < this->degree; ++i)
2461  interpolation_matrix(i, i) = 1.0;
2462 
2463  // In 3d we have some lines more
2464  // and a face. The procedure stays
2465  // the same as above, but we have
2466  // to take a bit more care of the
2467  // indices of the degrees of
2468  // freedom.
2469  if (dim == 3)
2470  {
2471  const unsigned int p = source_fe.degree;
2472  const unsigned int q = this->degree;
2473 
2474  for (unsigned int i = 0; i < q; ++i)
2475  {
2476  for (unsigned int j = 1; j < GeometryInfo<dim>::lines_per_face; ++j)
2477  interpolation_matrix(j * p + i, j * q + i) = 1.0;
2478 
2479  for (unsigned int j = 0; j < q - 1; ++j)
2480  {
2481  interpolation_matrix(GeometryInfo<dim>::lines_per_face * p +
2482  i * (p - 1) + j,
2484  i * (q - 1) + j) = 1.0;
2485  interpolation_matrix(GeometryInfo<dim>::lines_per_face * p + i +
2486  (j + p - 1) * p,
2488  (j + q - 1) * q) = 1.0;
2489  }
2490  }
2491  }
2492 }
2493 
2494 
2495 
2496 template <>
2497 void
2499  const unsigned int,
2500  FullMatrix<double> &) const
2501 {
2502  Assert(false, ExcNotImplemented());
2503 }
2504 
2505 
2506 
2507 // In this function we compute the
2508 // subface interpolation matrix.
2509 // This is done by a projection-
2510 // based interpolation. Therefore
2511 // we first interpolate the
2512 // shape functions of the higher
2513 // order element on the lowest
2514 // order edge shape functions.
2515 // Then the remaining part of
2516 // the interpolated shape
2517 // functions is projected on the
2518 // higher order edge shape
2519 // functions, the face shape
2520 // functions and the interior
2521 // shape functions (if they all
2522 // exist).
2523 template <int dim>
2524 void
2526  const FiniteElement<dim> &source,
2527  const unsigned int subface,
2528  FullMatrix<double> & interpolation_matrix) const
2529 {
2530  // this is only implemented, if the
2531  // source FE is also a
2532  // Nedelec element
2533  AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2534  (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2536  Assert(interpolation_matrix.m() == source.dofs_per_face,
2537  ExcDimensionMismatch(interpolation_matrix.m(), source.dofs_per_face));
2538  Assert(interpolation_matrix.n() == this->dofs_per_face,
2539  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
2540 
2541  // ok, source is a Nedelec element, so
2542  // we will be able to do the work
2543  const FE_Nedelec<dim> &source_fe =
2544  dynamic_cast<const FE_Nedelec<dim> &>(source);
2545 
2546  // Make sure, that the element,
2547  // for which the DoFs should be
2548  // constrained is the one with
2549  // the higher polynomial degree.
2550  // Actually the procedure will work
2551  // also if this assertion is not
2552  // satisfied. But the matrices
2553  // produced in that case might
2554  // lead to problems in the
2555  // hp procedures, which use this
2556  // method.
2557  Assert(this->dofs_per_face <= source_fe.dofs_per_face,
2559  interpolation_matrix = 0.0;
2560  // Perform projection-based interpolation
2561  // as usual.
2562  const QGauss<1> edge_quadrature(source_fe.degree);
2563  const std::vector<Point<1>> &edge_quadrature_points =
2564  edge_quadrature.get_points();
2565  const unsigned int n_edge_quadrature_points = edge_quadrature.size();
2566 
2567  switch (dim)
2568  {
2569  case 2:
2570  {
2571  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2572  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2573  ++q_point)
2574  {
2575  const Point<dim> quadrature_point(
2576  0.0, 0.5 * (edge_quadrature_points[q_point](0) + subface));
2577 
2578  interpolation_matrix(0, dof) +=
2579  0.5 * edge_quadrature.weight(q_point) *
2580  this->shape_value_component(dof, quadrature_point, 1);
2581  }
2582 
2583  if (source_fe.degree > 1)
2584  {
2585  const std::vector<Polynomials::Polynomial<double>>
2586  &legendre_polynomials =
2588  source_fe.degree - 1);
2589  FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
2590  source_fe.degree - 1);
2591 
2592  {
2593  FullMatrix<double> assembling_matrix(source_fe.degree - 1,
2594  n_edge_quadrature_points);
2595 
2596  for (unsigned int q_point = 0;
2597  q_point < n_edge_quadrature_points;
2598  ++q_point)
2599  {
2600  const double weight =
2601  std::sqrt(edge_quadrature.weight(q_point));
2602 
2603  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2604  assembling_matrix(i, q_point) =
2605  weight * legendre_polynomials[i + 1].value(
2606  edge_quadrature_points[q_point](0));
2607  }
2608 
2609  FullMatrix<double> system_matrix(source_fe.degree - 1,
2610  source_fe.degree - 1);
2611 
2612  assembling_matrix.mTmult(system_matrix, assembling_matrix);
2613  system_matrix_inv.invert(system_matrix);
2614  }
2615 
2616  Vector<double> solution(source_fe.degree - 1);
2617  Vector<double> system_rhs(source_fe.degree - 1);
2618 
2619  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2620  {
2621  system_rhs = 0.0;
2622 
2623  for (unsigned int q_point = 0;
2624  q_point < n_edge_quadrature_points;
2625  ++q_point)
2626  {
2627  const Point<dim> quadrature_point_0(
2628  0.0,
2629  0.5 * (edge_quadrature_points[q_point](0) + subface));
2630  const Point<dim> quadrature_point_1(
2631  0.0, edge_quadrature_points[q_point](0));
2632  const double tmp =
2633  edge_quadrature.weight(q_point) *
2634  (0.5 * this->shape_value_component(dof,
2635  quadrature_point_0,
2636  1) -
2637  interpolation_matrix(0, dof) *
2638  source_fe.shape_value_component(0,
2639  quadrature_point_1,
2640  1));
2641 
2642  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2643  system_rhs(i) +=
2644  tmp * legendre_polynomials[i + 1].value(
2645  edge_quadrature_points[q_point](0));
2646  }
2647 
2648  system_matrix_inv.vmult(solution, system_rhs);
2649 
2650  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2651  if (std::abs(solution(i)) > 1e-14)
2652  interpolation_matrix(i + 1, dof) = solution(i);
2653  }
2654  }
2655 
2656  break;
2657  }
2658 
2659  case 3:
2660  {
2661  const double shifts[4][2] = {{0.0, 0.0},
2662  {1.0, 0.0},
2663  {0.0, 1.0},
2664  {1.0, 1.0}};
2665 
2666  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2667  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2668  ++q_point)
2669  {
2670  const double weight = 0.5 * edge_quadrature.weight(q_point);
2671 
2672  for (unsigned int i = 0; i < 2; ++i)
2673  {
2674  Point<dim> quadrature_point(
2675  0.5 * (i + shifts[subface][0]),
2676  0.5 * (edge_quadrature_points[q_point](0) +
2677  shifts[subface][1]),
2678  0.0);
2679 
2680  interpolation_matrix(i * source_fe.degree, dof) +=
2681  weight *
2682  this->shape_value_component(
2683  this->face_to_cell_index(dof, 4), quadrature_point, 1);
2684  quadrature_point =
2685  Point<dim>(0.5 * (edge_quadrature_points[q_point](0) +
2686  shifts[subface][0]),
2687  0.5 * (i + shifts[subface][1]),
2688  0.0);
2689  interpolation_matrix((i + 2) * source_fe.degree, dof) +=
2690  weight *
2691  this->shape_value_component(
2692  this->face_to_cell_index(dof, 4), quadrature_point, 0);
2693  }
2694  }
2695 
2696  if (source_fe.degree > 1)
2697  {
2698  const std::vector<Polynomials::Polynomial<double>>
2699  &legendre_polynomials =
2701  source_fe.degree - 1);
2702  FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
2703  source_fe.degree - 1);
2704 
2705  {
2706  FullMatrix<double> assembling_matrix(source_fe.degree - 1,
2707  n_edge_quadrature_points);
2708 
2709  for (unsigned int q_point = 0;
2710  q_point < n_edge_quadrature_points;
2711  ++q_point)
2712  {
2713  const double weight =
2714  std::sqrt(edge_quadrature.weight(q_point));
2715 
2716  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2717  assembling_matrix(i, q_point) =
2718  weight * legendre_polynomials[i + 1].value(
2719  edge_quadrature_points[q_point](0));
2720  }
2721 
2722  FullMatrix<double> system_matrix(source_fe.degree - 1,
2723  source_fe.degree - 1);
2724 
2725  assembling_matrix.mTmult(system_matrix, assembling_matrix);
2726  system_matrix_inv.invert(system_matrix);
2727  }
2728 
2729  FullMatrix<double> solution(source_fe.degree - 1,
2731  FullMatrix<double> system_rhs(source_fe.degree - 1,
2734 
2735  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2736  {
2737  system_rhs = 0.0;
2738 
2739  for (unsigned int q_point = 0;
2740  q_point < n_edge_quadrature_points;
2741  ++q_point)
2742  {
2743  const double weight = edge_quadrature.weight(q_point);
2744 
2745  for (unsigned int i = 0; i < 2; ++i)
2746  {
2747  Point<dim> quadrature_point_0(
2748  0.5 * (i + shifts[subface][0]),
2749  0.5 * (edge_quadrature_points[q_point](0) +
2750  shifts[subface][1]),
2751  0.0);
2752  Point<dim> quadrature_point_1(
2753  i, edge_quadrature_points[q_point](0), 0.0);
2754 
2755  tmp(i) =
2756  weight *
2757  (0.5 * this->shape_value_component(
2758  this->face_to_cell_index(dof, 4),
2759  quadrature_point_0,
2760  1) -
2761  interpolation_matrix(i * source_fe.degree, dof) *
2762  source_fe.shape_value_component(
2763  i * source_fe.degree, quadrature_point_1, 1));
2764  quadrature_point_0 =
2765  Point<dim>(0.5 *
2766  (edge_quadrature_points[q_point](0) +
2767  shifts[subface][0]),
2768  0.5 * (i + shifts[subface][1]),
2769  0.0);
2770  quadrature_point_1 =
2771  Point<dim>(edge_quadrature_points[q_point](0),
2772  i,
2773  0.0);
2774  tmp(i + 2) =
2775  weight *
2776  (0.5 * this->shape_value_component(
2777  this->face_to_cell_index(dof, 4),
2778  quadrature_point_0,
2779  0) -
2780  interpolation_matrix((i + 2) * source_fe.degree,
2781  dof) *
2782  source_fe.shape_value_component(
2783  (i + 2) * source_fe.degree,
2784  quadrature_point_1,
2785  0));
2786  }
2787 
2788  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2789  {
2790  const double L_i = legendre_polynomials[i + 1].value(
2791  edge_quadrature_points[q_point](0));
2792 
2793  for (unsigned int j = 0;
2794  j < GeometryInfo<dim>::lines_per_face;
2795  ++j)
2796  system_rhs(i, j) += tmp(j) * L_i;
2797  }
2798  }
2799 
2800  system_matrix_inv.mmult(solution, system_rhs);
2801 
2802  for (unsigned int i = 0;
2803  i < GeometryInfo<dim>::lines_per_face;
2804  ++i)
2805  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2806  if (std::abs(solution(j, i)) > 1e-14)
2807  interpolation_matrix(i * source_fe.degree + j + 1,
2808  dof) = solution(j, i);
2809  }
2810 
2811  const QGauss<2> quadrature(source_fe.degree);
2812  const std::vector<Point<2>> &quadrature_points =
2813  quadrature.get_points();
2814  const std::vector<Polynomials::Polynomial<double>>
2815  &lobatto_polynomials =
2817  source_fe.degree);
2818  const unsigned int n_boundary_dofs =
2820  const unsigned int n_quadrature_points = quadrature.size();
2821 
2822  {
2823  FullMatrix<double> assembling_matrix(source_fe.degree *
2824  (source_fe.degree - 1),
2825  n_quadrature_points);
2826 
2827  for (unsigned int q_point = 0; q_point < n_quadrature_points;
2828  ++q_point)
2829  {
2830  const double weight = std::sqrt(quadrature.weight(q_point));
2831 
2832  for (unsigned int i = 0; i < source_fe.degree; ++i)
2833  {
2834  const double L_i =
2835  weight * legendre_polynomials[i].value(
2836  quadrature_points[q_point](0));
2837 
2838  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2839  assembling_matrix(i * (source_fe.degree - 1) + j,
2840  q_point) =
2841  L_i * lobatto_polynomials[j + 2].value(
2842  quadrature_points[q_point](1));
2843  }
2844  }
2845 
2846  FullMatrix<double> system_matrix(assembling_matrix.m(),
2847  assembling_matrix.m());
2848 
2849  assembling_matrix.mTmult(system_matrix, assembling_matrix);
2850  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
2851  system_matrix_inv.invert(system_matrix);
2852  }
2853 
2854  solution.reinit(system_matrix_inv.m(), 2);
2855  system_rhs.reinit(system_matrix_inv.m(), 2);
2856  tmp.reinit(2);
2857 
2858  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2859  {
2860  system_rhs = 0.0;
2861 
2862  for (unsigned int q_point = 0; q_point < n_quadrature_points;
2863  ++q_point)
2864  {
2865  Point<dim> quadrature_point(
2866  0.5 *
2867  (quadrature_points[q_point](0) + shifts[subface][0]),
2868  0.5 *
2869  (quadrature_points[q_point](1) + shifts[subface][1]),
2870  0.0);
2871  tmp(0) = 0.5 * this->shape_value_component(
2872  this->face_to_cell_index(dof, 4),
2873  quadrature_point,
2874  0);
2875  tmp(1) = 0.5 * this->shape_value_component(
2876  this->face_to_cell_index(dof, 4),
2877  quadrature_point,
2878  1);
2879  quadrature_point =
2880  Point<dim>(quadrature_points[q_point](0),
2881  quadrature_points[q_point](1),
2882  0.0);
2883 
2884  for (unsigned int i = 0; i < 2; ++i)
2885  for (unsigned int j = 0; j < source_fe.degree; ++j)
2886  {
2887  tmp(0) -= interpolation_matrix(
2888  (i + 2) * source_fe.degree + j, dof) *
2889  source_fe.shape_value_component(
2890  (i + 2) * source_fe.degree + j,
2891  quadrature_point,
2892  0);
2893  tmp(1) -=
2894  interpolation_matrix(i * source_fe.degree + j,
2895  dof) *
2896  source_fe.shape_value_component(
2897  i * source_fe.degree + j, quadrature_point, 1);
2898  }
2899 
2900  tmp *= quadrature.weight(q_point);
2901 
2902  for (unsigned int i = 0; i < source_fe.degree; ++i)
2903  {
2904  const double L_i_0 = legendre_polynomials[i].value(
2905  quadrature_points[q_point](0));
2906  const double L_i_1 = legendre_polynomials[i].value(
2907  quadrature_points[q_point](1));
2908 
2909  for (unsigned int j = 0; j < source_fe.degree - 1;
2910  ++j)
2911  {
2912  system_rhs(i * (source_fe.degree - 1) + j, 0) +=
2913  tmp(0) * L_i_0 *
2914  lobatto_polynomials[j + 2].value(
2915  quadrature_points[q_point](1));
2916  system_rhs(i * (source_fe.degree - 1) + j, 1) +=
2917  tmp(1) * L_i_1 *
2918  lobatto_polynomials[j + 2].value(
2919  quadrature_points[q_point](0));
2920  }
2921  }
2922  }
2923 
2924  system_matrix_inv.mmult(solution, system_rhs);
2925 
2926  for (unsigned int i = 0; i < source_fe.degree; ++i)
2927  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2928  {
2929  if (std::abs(solution(i * (source_fe.degree - 1) + j,
2930  0)) > 1e-14)
2931  interpolation_matrix(i * (source_fe.degree - 1) + j +
2932  n_boundary_dofs,
2933  dof) =
2934  solution(i * (source_fe.degree - 1) + j, 0);
2935 
2936  if (std::abs(solution(i * (source_fe.degree - 1) + j,
2937  1)) > 1e-14)
2938  interpolation_matrix(
2939  i + (j + source_fe.degree - 1) * source_fe.degree +
2940  n_boundary_dofs,
2941  dof) = solution(i * (source_fe.degree - 1) + j, 1);
2942  }
2943  }
2944  }
2945 
2946  break;
2947  }
2948 
2949  default:
2950  Assert(false, ExcNotImplemented());
2951  }
2952 }
2953 
2954 template <int dim>
2955 const FullMatrix<double> &
2957  const unsigned int child,
2958  const RefinementCase<dim> &refinement_case) const
2959 {
2960  AssertIndexRange(refinement_case,
2962  Assert(refinement_case != RefinementCase<dim>::no_refinement,
2963  ExcMessage(
2964  "Prolongation matrices are only available for refined cells!"));
2965  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
2966 
2967  // initialization upon first request
2968  if (this->prolongation[refinement_case - 1][child].n() == 0)
2969  {
2970  std::lock_guard<std::mutex> lock(this->mutex);
2971 
2972  // if matrix got updated while waiting for the lock
2973  if (this->prolongation[refinement_case - 1][child].n() ==
2974  this->dofs_per_cell)
2975  return this->prolongation[refinement_case - 1][child];
2976 
2977  // now do the work. need to get a non-const version of data in order to
2978  // be able to modify them inside a const function
2979  FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
2980 
2981  // Reinit the vectors of
2982  // restriction and prolongation
2983  // matrices to the right sizes.
2984  // Restriction only for isotropic
2985  // refinement
2986 #ifdef DEBUG_NEDELEC
2987  deallog << "Embedding" << std::endl;
2988 #endif
2990  // Fill prolongation matrices with embedding operators
2992  this_nonconst,
2993  this_nonconst.prolongation,
2994  true,
2995  internal::FE_Nedelec::get_embedding_computation_tolerance(
2996  this->degree));
2997 #ifdef DEBUG_NEDELEC
2998  deallog << "Restriction" << std::endl;
2999 #endif
3000  this_nonconst.initialize_restriction();
3001  }
3002 
3003  // we use refinement_case-1 here. the -1 takes care of the origin of the
3004  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3005  // available and so the vector indices are shifted
3006  return this->prolongation[refinement_case - 1][child];
3007 }
3008 
3009 template <int dim>
3010 const FullMatrix<double> &
3012  const unsigned int child,
3013  const RefinementCase<dim> &refinement_case) const
3014 {
3015  AssertIndexRange(refinement_case,
3017  Assert(refinement_case != RefinementCase<dim>::no_refinement,
3018  ExcMessage(
3019  "Restriction matrices are only available for refined cells!"));
3021  child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
3022 
3023  // initialization upon first request
3024  if (this->restriction[refinement_case - 1][child].n() == 0)
3025  {
3026  std::lock_guard<std::mutex> lock(this->mutex);
3027 
3028  // if matrix got updated while waiting for the lock...
3029  if (this->restriction[refinement_case - 1][child].n() ==
3030  this->dofs_per_cell)
3031  return this->restriction[refinement_case - 1][child];
3032 
3033  // now do the work. need to get a non-const version of data in order to
3034  // be able to modify them inside a const function
3035  FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
3036 
3037  // Reinit the vectors of
3038  // restriction and prolongation
3039  // matrices to the right sizes.
3040  // Restriction only for isotropic
3041  // refinement
3042 #ifdef DEBUG_NEDELEC
3043  deallog << "Embedding" << std::endl;
3044 #endif
3046  // Fill prolongation matrices with embedding operators
3048  this_nonconst,
3049  this_nonconst.prolongation,
3050  true,
3051  internal::FE_Nedelec::get_embedding_computation_tolerance(
3052  this->degree));
3053 #ifdef DEBUG_NEDELEC
3054  deallog << "Restriction" << std::endl;
3055 #endif
3056  this_nonconst.initialize_restriction();
3057  }
3058 
3059  // we use refinement_case-1 here. the -1 takes care of the origin of the
3060  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3061  // available and so the vector indices are shifted
3062  return this->restriction[refinement_case - 1][child];
3063 }
3064 
3065 
3066 // Interpolate a function, which is given by
3067 // its values at the generalized support
3068 // points in the finite element space on the
3069 // reference cell.
3070 // This is done as usual by projection-based
3071 // interpolation.
3072 template <int dim>
3073 void
3075  const std::vector<Vector<double>> &support_point_values,
3076  std::vector<double> & nodal_values) const
3077 {
3078  const unsigned int deg = this->degree - 1;
3079  Assert(support_point_values.size() == this->generalized_support_points.size(),
3080  ExcDimensionMismatch(support_point_values.size(),
3081  this->generalized_support_points.size()));
3082  Assert(support_point_values[0].size() == this->n_components(),
3083  ExcDimensionMismatch(support_point_values[0].size(),
3084  this->n_components()));
3085  Assert(nodal_values.size() == this->dofs_per_cell,
3086  ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
3087  std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3088 
3089  switch (dim)
3090  {
3091  case 2:
3092  {
3093  // Let us begin with the
3094  // interpolation part.
3095  const QGauss<dim - 1> reference_edge_quadrature(this->degree);
3096  const unsigned int n_edge_points = reference_edge_quadrature.size();
3097 
3098  for (unsigned int i = 0; i < 2; ++i)
3099  for (unsigned int j = 0; j < 2; ++j)
3100  {
3101  for (unsigned int q_point = 0; q_point < n_edge_points;
3102  ++q_point)
3103  nodal_values[(i + 2 * j) * this->degree] +=
3104  reference_edge_quadrature.weight(q_point) *
3105  support_point_values[q_point + (i + 2 * j) * n_edge_points]
3106  [1 - j];
3107 
3108  // Add the computed support_point_values to the resulting vector
3109  // only, if they are not
3110  // too small.
3111  if (std::abs(nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3112  nodal_values[(i + 2 * j) * this->degree] = 0.0;
3113  }
3114 
3115  // If the degree is greater
3116  // than 0, then we have still
3117  // some higher order edge
3118  // shape functions to
3119  // consider.
3120  // Here the projection part
3121  // starts. The dof support_point_values
3122  // are obtained by solving
3123  // a linear system of
3124  // equations.
3125  if (this->degree - 1 > 1)
3126  {
3127  // We start with projection
3128  // on the higher order edge
3129  // shape function.
3130  const std::vector<Polynomials::Polynomial<double>>
3131  &lobatto_polynomials =
3133  FullMatrix<double> system_matrix(this->degree - 1,
3134  this->degree - 1);
3135  std::vector<Polynomials::Polynomial<double>>
3136  lobatto_polynomials_grad(this->degree);
3137 
3138  for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3139  lobatto_polynomials_grad[i] =
3140  lobatto_polynomials[i + 1].derivative();
3141 
3142  // Set up the system matrix.
3143  // This can be used for all
3144  // edges.
3145  for (unsigned int i = 0; i < system_matrix.m(); ++i)
3146  for (unsigned int j = 0; j < system_matrix.n(); ++j)
3147  for (unsigned int q_point = 0; q_point < n_edge_points;
3148  ++q_point)
3149  system_matrix(i, j) +=
3150  boundary_weights(q_point, j) *
3151  lobatto_polynomials_grad[i + 1].value(
3152  this->generalized_face_support_points[q_point](0));
3153 
3154  FullMatrix<double> system_matrix_inv(this->degree - 1,
3155  this->degree - 1);
3156 
3157  system_matrix_inv.invert(system_matrix);
3158 
3159  const unsigned int
3160  line_coordinate[GeometryInfo<2>::lines_per_cell] = {1, 1, 0, 0};
3161  Vector<double> system_rhs(system_matrix.m());
3162  Vector<double> solution(system_rhs.size());
3163 
3164  for (unsigned int line = 0;
3165  line < GeometryInfo<dim>::lines_per_cell;
3166  ++line)
3167  {
3168  // Set up the right hand side.
3169  system_rhs = 0;
3170 
3171  for (unsigned int q_point = 0; q_point < n_edge_points;
3172  ++q_point)
3173  {
3174  const double tmp =
3175  support_point_values[line * n_edge_points + q_point]
3176  [line_coordinate[line]] -
3177  nodal_values[line * this->degree] *
3178  this->shape_value_component(
3179  line * this->degree,
3180  this->generalized_support_points[line *
3181  n_edge_points +
3182  q_point],
3183  line_coordinate[line]);
3184 
3185  for (unsigned int i = 0; i < system_rhs.size(); ++i)
3186  system_rhs(i) += boundary_weights(q_point, i) * tmp;
3187  }
3188 
3189  system_matrix_inv.vmult(solution, system_rhs);
3190 
3191  // Add the computed support_point_values
3192  // to the resulting vector
3193  // only, if they are not
3194  // too small.
3195  for (unsigned int i = 0; i < solution.size(); ++i)
3196  if (std::abs(solution(i)) > 1e-14)
3197  nodal_values[line * this->degree + i + 1] = solution(i);
3198  }
3199 
3200  // Then we go on to the
3201  // interior shape
3202  // functions. Again we
3203  // set up the system
3204  // matrix and use it
3205  // for both, the
3206  // horizontal and the
3207  // vertical, interior
3208  // shape functions.
3209  const QGauss<dim> reference_quadrature(this->degree);
3210  const unsigned int n_interior_points =
3211  reference_quadrature.size();
3212  const std::vector<Polynomials::Polynomial<double>>
3213  &legendre_polynomials =
3215  1);
3216 
3217  system_matrix.reinit((this->degree - 1) * this->degree,
3218  (this->degree - 1) * this->degree);
3219  system_matrix = 0;
3220 
3221  for (unsigned int i = 0; i < this->degree; ++i)
3222  for (unsigned int j = 0; j < this->degree - 1; ++j)
3223  for (unsigned int k = 0; k < this->degree; ++k)
3224  for (unsigned int l = 0; l < this->degree - 1; ++l)
3225  for (unsigned int q_point = 0;
3226  q_point < n_interior_points;
3227  ++q_point)
3228  system_matrix(i * (this->degree - 1) + j,
3229  k * (this->degree - 1) + l) +=
3230  reference_quadrature.weight(q_point) *
3231  legendre_polynomials[i].value(
3234  n_edge_points](0)) *
3235  lobatto_polynomials[j + 2].value(
3238  n_edge_points](1)) *
3239  lobatto_polynomials_grad[k].value(
3242  n_edge_points](0)) *
3243  lobatto_polynomials[l + 2].value(
3246  n_edge_points](1));
3247 
3248  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3249  system_matrix_inv.invert(system_matrix);
3250  // Set up the right hand side
3251  // for the horizontal shape
3252  // functions.
3253  system_rhs.reinit(system_matrix_inv.m());
3254  system_rhs = 0;
3255 
3256  for (unsigned int q_point = 0; q_point < n_interior_points;
3257  ++q_point)
3258  {
3259  double tmp =
3260  support_point_values[q_point +
3262  n_edge_points][0];
3263 
3264  for (unsigned int i = 0; i < 2; ++i)
3265  for (unsigned int j = 0; j <= deg; ++j)
3266  tmp -= nodal_values[(i + 2) * this->degree + j] *
3267  this->shape_value_component(
3268  (i + 2) * this->degree + j,
3271  n_edge_points],
3272  0);
3273 
3274  for (unsigned int i = 0; i <= deg; ++i)
3275  for (unsigned int j = 0; j < deg; ++j)
3276  system_rhs(i * deg + j) +=
3277  reference_quadrature.weight(q_point) * tmp *
3278  lobatto_polynomials_grad[i].value(
3281  n_edge_points](0)) *
3282  lobatto_polynomials[j + 2].value(
3285  n_edge_points](1));
3286  }
3287 
3288  solution.reinit(system_matrix.m());
3289  system_matrix_inv.vmult(solution, system_rhs);
3290 
3291  // Add the computed support_point_values
3292  // to the resulting vector
3293  // only, if they are not
3294  // too small.
3295  for (unsigned int i = 0; i <= deg; ++i)
3296  for (unsigned int j = 0; j < deg; ++j)
3297  if (std::abs(solution(i * deg + j)) > 1e-14)
3298  nodal_values[(i + GeometryInfo<dim>::lines_per_cell) * deg +
3300  solution(i * deg + j);
3301 
3302  system_rhs = 0;
3303  // Set up the right hand side
3304  // for the vertical shape
3305  // functions.
3306 
3307  for (unsigned int q_point = 0; q_point < n_interior_points;
3308  ++q_point)
3309  {
3310  double tmp =
3311  support_point_values[q_point +
3313  n_edge_points][1];
3314 
3315  for (unsigned int i = 0; i < 2; ++i)
3316  for (unsigned int j = 0; j <= deg; ++j)
3317  tmp -= nodal_values[i * this->degree + j] *
3318  this->shape_value_component(
3319  i * this->degree + j,
3322  n_edge_points],
3323  1);
3324 
3325  for (unsigned int i = 0; i <= deg; ++i)
3326  for (unsigned int j = 0; j < deg; ++j)
3327  system_rhs(i * deg + j) +=
3328  reference_quadrature.weight(q_point) * tmp *
3329  lobatto_polynomials_grad[i].value(
3332  n_edge_points](1)) *
3333  lobatto_polynomials[j + 2].value(
3336  n_edge_points](0));
3337  }
3338 
3339  system_matrix_inv.vmult(solution, system_rhs);
3340 
3341  // Add the computed support_point_values
3342  // to the resulting vector
3343  // only, if they are not
3344  // too small.
3345  for (unsigned int i = 0; i <= deg; ++i)
3346  for (unsigned int j = 0; j < deg; ++j)
3347  if (std::abs(solution(i * deg + j)) > 1e-14)
3348  nodal_values[i +
3349  (j + GeometryInfo<dim>::lines_per_cell + deg) *
3350  this->degree] = solution(i * deg + j);
3351  }
3352 
3353  break;
3354  }
3355 
3356  case 3:
3357  {
3358  // Let us begin with the
3359  // interpolation part.
3360  const QGauss<1> reference_edge_quadrature(this->degree);
3361  const unsigned int n_edge_points = reference_edge_quadrature.size();
3362 
3363  for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3364  {
3365  for (unsigned int i = 0; i < 4; ++i)
3366  nodal_values[(i + 8) * this->degree] +=
3367  reference_edge_quadrature.weight(q_point) *
3368  support_point_values[q_point + (i + 8) * n_edge_points][2];
3369 
3370  for (unsigned int i = 0; i < 2; ++i)
3371  for (unsigned int j = 0; j < 2; ++j)
3372  for (unsigned int k = 0; k < 2; ++k)
3373  nodal_values[(i + 2 * (2 * j + k)) * this->degree] +=
3374  reference_edge_quadrature.weight(q_point) *
3375  support_point_values[q_point + (i + 2 * (2 * j + k)) *
3376  n_edge_points][1 - k];
3377  }
3378 
3379  // Add the computed support_point_values
3380  // to the resulting vector
3381  // only, if they are not
3382  // too small.
3383  for (unsigned int i = 0; i < 4; ++i)
3384  if (std::abs(nodal_values[(i + 8) * this->degree]) < 1e-14)
3385  nodal_values[(i + 8) * this->degree] = 0.0;
3386 
3387  for (unsigned int i = 0; i < 2; ++i)
3388  for (unsigned int j = 0; j < 2; ++j)
3389  for (unsigned int k = 0; k < 2; ++k)
3390  if (std::abs(
3391  nodal_values[(i + 2 * (2 * j + k)) * this->degree]) <
3392  1e-14)
3393  nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3394 
3395  // If the degree is greater
3396  // than 0, then we have still
3397  // some higher order shape
3398  // functions to consider.
3399  // Here the projection part
3400  // starts. The dof support_point_values
3401  // are obtained by solving
3402  // a linear system of
3403  // equations.
3404  if (this->degree > 1)
3405  {
3406  // We start with projection
3407  // on the higher order edge
3408  // shape function.
3409  const std::vector<Polynomials::Polynomial<double>>
3410  &lobatto_polynomials =
3412  FullMatrix<double> system_matrix(this->degree - 1,
3413  this->degree - 1);
3414  std::vector<Polynomials::Polynomial<double>>
3415  lobatto_polynomials_grad(this->degree);
3416 
3417  for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3418  lobatto_polynomials_grad[i] =
3419  lobatto_polynomials[i + 1].derivative();
3420 
3421  // Set up the system matrix.
3422  // This can be used for all
3423  // edges.
3424  for (unsigned int i = 0; i < system_matrix.m(); ++i)
3425  for (unsigned int j = 0; j < system_matrix.n(); ++j)
3426  for (unsigned int q_point = 0; q_point < n_edge_points;
3427  ++q_point)
3428  system_matrix(i, j) +=
3429  boundary_weights(q_point, j) *
3430  lobatto_polynomials_grad[i + 1].value(
3431  this->generalized_face_support_points[q_point](1));
3432 
3433  FullMatrix<double> system_matrix_inv(this->degree - 1,
3434  this->degree - 1);
3435 
3436  system_matrix_inv.invert(system_matrix);
3437 
3438  const unsigned int
3439  line_coordinate[GeometryInfo<3>::lines_per_cell] = {
3440  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3441  Vector<double> system_rhs(system_matrix.m());
3442  Vector<double> solution(system_rhs.size());
3443 
3444  for (unsigned int line = 0;
3445  line < GeometryInfo<dim>::lines_per_cell;
3446  ++line)
3447  {
3448  // Set up the right hand side.
3449  system_rhs = 0;
3450 
3451  for (unsigned int q_point = 0; q_point < this->degree;
3452  ++q_point)
3453  {
3454  const double tmp =
3455  support_point_values[line * this->degree + q_point]
3456  [line_coordinate[line]] -
3457  nodal_values[line * this->degree] *
3458  this->shape_value_component(
3459  line * this->degree,
3460  this
3461  ->generalized_support_points[line * this->degree +
3462  q_point],
3463  line_coordinate[line]);
3464 
3465  for (unsigned int i = 0; i < system_rhs.size(); ++i)
3466  system_rhs(i) += boundary_weights(q_point, i) * tmp;
3467  }
3468 
3469  system_matrix_inv.vmult(solution, system_rhs);
3470 
3471  // Add the computed values
3472  // to the resulting vector
3473  // only, if they are not
3474  // too small.
3475  for (unsigned int i = 0; i < solution.size(); ++i)
3476  if (std::abs(solution(i)) > 1e-14)
3477  nodal_values[line * this->degree + i + 1] = solution(i);
3478  }
3479 
3480  // Then we go on to the
3481  // face shape functions.
3482  // Again we set up the
3483  // system matrix and
3484  // use it for both, the
3485  // horizontal and the
3486  // vertical, shape
3487  // functions.
3488  const std::vector<Polynomials::Polynomial<double>>
3489  &legendre_polynomials =
3491  1);
3492  const unsigned int n_face_points = n_edge_points * n_edge_points;
3493 
3494  system_matrix.reinit((this->degree - 1) * this->degree,
3495  (this->degree - 1) * this->degree);
3496  system_matrix = 0;
3497 
3498  for (unsigned int i = 0; i < this->degree; ++i)
3499  for (unsigned int j = 0; j < this->degree - 1; ++j)
3500  for (unsigned int k = 0; k < this->degree; ++k)
3501  for (unsigned int l = 0; l < this->degree - 1; ++l)
3502  for (unsigned int q_point = 0; q_point < n_face_points;
3503  ++q_point)
3504  system_matrix(i * (this->degree - 1) + j,
3505  k * (this->degree - 1) + l) +=
3506  boundary_weights(q_point + n_edge_points,
3507  2 * (k * (this->degree - 1) + l)) *
3508  legendre_polynomials[i].value(
3510  [q_point + 4 * n_edge_points](0)) *
3511  lobatto_polynomials[j + 2].value(
3513  [q_point + 4 * n_edge_points](1));
3514 
3515  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3516  system_matrix_inv.invert(system_matrix);
3517  solution.reinit(system_matrix.m());
3518  system_rhs.reinit(system_matrix.m());
3519 
3520  const unsigned int
3521  face_coordinates[GeometryInfo<3>::faces_per_cell][2] = {
3522  {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3523  const unsigned int edge_indices[GeometryInfo<3>::faces_per_cell]
3525  {{0, 4, 8, 10},
3526  {1, 5, 9, 11},
3527  {8, 9, 2, 6},
3528  {10, 11, 3, 7},
3529  {2, 3, 0, 1},
3530  {6, 7, 4, 5}};
3531 
3532  for (unsigned int face = 0;
3533  face < GeometryInfo<dim>::faces_per_cell;
3534  ++face)
3535  {
3536  // Set up the right hand side
3537  // for the horizontal shape
3538  // functions.
3539  system_rhs = 0;
3540 
3541  for (unsigned int q_point = 0; q_point < n_face_points;
3542  ++q_point)
3543  {
3544  double tmp =
3545  support_point_values[q_point +
3547  n_edge_points]
3548  [face_coordinates[face][0]];
3549 
3550  for (unsigned int i = 0; i < 2; ++i)
3551  for (unsigned int j = 0; j <= deg; ++j)
3552  tmp -=
3553  nodal_values[edge_indices[face][i] * this->degree +
3554  j] *
3555  this->shape_value_component(
3556  edge_indices[face][i] * this->degree + j,
3559  n_edge_points],
3560  face_coordinates[face][0]);
3561 
3562  for (unsigned int i = 0; i <= deg; ++i)
3563  for (unsigned int j = 0; j < deg; ++j)
3564  system_rhs(i * deg + j) +=
3565  boundary_weights(q_point + n_edge_points,
3566  2 * (i * deg + j)) *
3567  tmp;
3568  }
3569 
3570  system_matrix_inv.vmult(solution, system_rhs);
3571 
3572  // Add the computed support_point_values
3573  // to the resulting vector
3574  // only, if they are not
3575  // too small.
3576  for (unsigned int i = 0; i <= deg; ++i)
3577  for (unsigned int j = 0; j < deg; ++j)
3578  if (std::abs(solution(i * deg + j)) > 1e-14)
3579  nodal_values[(2 * face * this->degree + i +
3581  deg +
3583  solution(i * deg + j);
3584 
3585  // Set up the right hand side
3586  // for the vertical shape
3587  // functions.
3588  system_rhs = 0;
3589 
3590  for (unsigned int q_point = 0; q_point < n_face_points;
3591  ++q_point)
3592  {
3593  double tmp =
3594  support_point_values[q_point +
3596  n_edge_points]
3597  [face_coordinates[face][1]];
3598 
3599  for (unsigned int i = 2;
3600  i < GeometryInfo<dim>::lines_per_face;
3601  ++i)
3602  for (unsigned int j = 0; j <= deg; ++j)
3603  tmp -=
3604  nodal_values[edge_indices[face][i] * this->degree +
3605  j] *
3606  this->shape_value_component(
3607  edge_indices[face][i] * this->degree + j,
3610  n_edge_points],
3611  face_coordinates[face][1]);
3612 
3613  for (unsigned int i = 0; i <= deg; ++i)
3614  for (unsigned int j = 0; j < deg; ++j)
3615  system_rhs(i * deg + j) +=
3616  boundary_weights(q_point + n_edge_points,
3617  2 * (i * deg + j) + 1) *
3618  tmp;
3619  }
3620 
3621  system_matrix_inv.vmult(solution, system_rhs);
3622 
3623  // Add the computed support_point_values
3624  // to the resulting vector
3625  // only, if they are not
3626  // too small.
3627  for (unsigned int i = 0; i <= deg; ++i)
3628  for (unsigned int j = 0; j < deg; ++j)
3629  if (std::abs(solution(i * deg + j)) > 1e-14)
3630  nodal_values[((2 * face + 1) * deg + j +
3632  this->degree +
3633  i] = solution(i * deg + j);
3634  }
3635 
3636  // Finally we project
3637  // the remaining parts
3638  // of the function on
3639  // the interior shape
3640  // functions.
3641  const QGauss<dim> reference_quadrature(this->degree);
3642  const unsigned int n_interior_points =
3643  reference_quadrature.size();
3644 
3645  // We create the
3646  // system matrix.
3647  system_matrix.reinit(this->degree * deg * deg,
3648  this->degree * deg * deg);
3649  system_matrix = 0;
3650 
3651  for (unsigned int i = 0; i <= deg; ++i)
3652  for (unsigned int j = 0; j < deg; ++j)
3653  for (unsigned int k = 0; k < deg; ++k)
3654  for (unsigned int l = 0; l <= deg; ++l)
3655  for (unsigned int m = 0; m < deg; ++m)
3656  for (unsigned int n = 0; n < deg; ++n)
3657  for (unsigned int q_point = 0;
3658  q_point < n_interior_points;
3659  ++q_point)
3660  system_matrix((i * deg + j) * deg + k,
3661  (l * deg + m) * deg + n) +=
3662  reference_quadrature.weight(q_point) *
3663  legendre_polynomials[i].value(
3665  [q_point +
3667  n_edge_points +
3669  n_face_points](0)) *
3670  lobatto_polynomials[j + 2].value(
3672  [q_point +
3674  n_edge_points +
3676  n_face_points](1)) *
3677  lobatto_polynomials[k + 2].value(
3679  [q_point +
3681  n_edge_points +
3683  n_face_points](2)) *
3684  lobatto_polynomials_grad[l].value(
3686  [q_point +
3688  n_edge_points +
3690  n_face_points](0)) *
3691  lobatto_polynomials[m + 2].value(
3693  [q_point +
3695  n_edge_points +
3697  n_face_points](1)) *
3698  lobatto_polynomials[n + 2].value(
3700  [q_point +
3702  n_edge_points +
3704  n_face_points](2));
3705 
3706  system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3707  system_matrix_inv.invert(system_matrix);
3708  // Set up the right hand side.
3709  system_rhs.reinit(system_matrix.m());
3710  system_rhs = 0;
3711 
3712  for (unsigned int q_point = 0; q_point < n_interior_points;
3713  ++q_point)
3714  {
3715  double tmp =
3716  support_point_values[q_point +
3718  n_edge_points +
3720  n_face_points][0];
3721 
3722  for (unsigned int i = 0; i <= deg; ++i)
3723  {
3724  for (unsigned int j = 0; j < 2; ++j)
3725  for (unsigned int k = 0; k < 2; ++k)
3726  tmp -=
3727  nodal_values[i + (j + 4 * k + 2) * this->degree] *
3728  this->shape_value_component(
3729  i + (j + 4 * k + 2) * this->degree,
3731  [q_point +
3733  n_edge_points +
3735  n_face_points],
3736  0);
3737 
3738  for (unsigned int j = 0; j < deg; ++j)
3739  for (unsigned int k = 0; k < 4; ++k)
3740  tmp -=
3741  nodal_values[(i + 2 * (k + 2) * this->degree +
3743  deg +
3744  j +
3746  this->shape_value_component(
3747  (i + 2 * (k + 2) * this->degree +
3749  deg +
3752  [q_point +
3754  n_edge_points +
3756  n_face_points],
3757  0);
3758  }
3759 
3760  for (unsigned int i = 0; i <= deg; ++i)
3761  for (unsigned int j = 0; j < deg; ++j)
3762  for (unsigned int k = 0; k < deg; ++k)
3763  system_rhs((i * deg + j) * deg + k) +=
3764  reference_quadrature.weight(q_point) * tmp *
3765  lobatto_polynomials_grad[i].value(
3767  [q_point +
3769  n_edge_points +
3771  n_face_points](0)) *
3772  lobatto_polynomials[j + 2].value(
3774  [q_point +
3776  n_edge_points +
3778  n_face_points](1)) *
3779  lobatto_polynomials[k + 2].value(
3781  [q_point +
3783  n_edge_points +
3785  n_face_points](2));
3786  }
3787 
3788  solution.reinit(system_rhs.size());
3789  system_matrix_inv.vmult(solution, system_rhs);
3790 
3791  // Add the computed values
3792  // to the resulting vector
3793  // only, if they are not
3794  // too small.
3795  for (unsigned int i = 0; i <= deg; ++i)
3796  for (unsigned int j = 0; j < deg; ++j)
3797  for (unsigned int k = 0; k < deg; ++k)
3798  if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3799  nodal_values
3800  [((i + 2 * GeometryInfo<dim>::faces_per_cell) * deg +
3803  deg +
3805  solution((i * deg + j) * deg + k);
3806 
3807  // Set up the right hand side.
3808  system_rhs = 0;
3809 
3810  for (unsigned int q_point = 0; q_point < n_interior_points;
3811  ++q_point)
3812  {
3813  double tmp =
3814  support_point_values[q_point +
3816  n_edge_points +
3818  n_face_points][1];
3819 
3820  for (unsigned int i = 0; i <= deg; ++i)
3821  for (unsigned int j = 0; j < 2; ++j)
3822  {
3823  for (unsigned int k = 0; k < 2; ++k)
3824  tmp -= nodal_values[i + (4 * j + k) * this->degree] *
3825  this->shape_value_component(
3826  i + (4 * j + k) * this->degree,
3828  [q_point +
3830  n_edge_points +
3832  n_face_points],
3833  1);
3834 
3835  for (unsigned int k = 0; k < deg; ++k)
3836  tmp -=
3837  nodal_values[(i + 2 * j * this->degree +
3839  deg +
3840  k +
3842  this->shape_value_component(
3843  (i + 2 * j * this->degree +
3845  deg +
3848  [q_point +
3850  n_edge_points +
3852  n_face_points],
3853  1) +
3854  nodal_values[i +
3855  ((2 * j + 9) * deg + k +
3857  this->degree] *
3858  this->shape_value_component(
3859  i + ((2 * j + 9) * deg + k +
3861  this->degree,
3863  [q_point +
3865  n_edge_points +
3867  n_face_points],
3868  1);
3869  }
3870 
3871  for (unsigned int i = 0; i <= deg; ++i)
3872  for (unsigned int j = 0; j < deg; ++j)
3873  for (unsigned int k = 0; k < deg; ++k)
3874  system_rhs((i * deg + j) * deg + k) +=
3875  reference_quadrature.weight(q_point) * tmp *
3876  lobatto_polynomials_grad[i].value(
3878  [q_point +
3880  n_edge_points +
3882  n_face_points](1)) *
3883  lobatto_polynomials[j + 2].value(
3885  [q_point +
3887  n_edge_points +
3889  n_face_points](0)) *
3890  lobatto_polynomials[k + 2].value(
3892  [q_point +
3894  n_edge_points +
3896  n_face_points](2));
3897  }
3898 
3899  system_matrix_inv.vmult(solution, system_rhs);
3900 
3901  // Add the computed support_point_values
3902  // to the resulting vector
3903  // only, if they are not
3904  // too small.
3905  for (unsigned int i = 0; i <= deg; ++i)
3906  for (unsigned int j = 0; j < deg; ++j)
3907  for (unsigned int k = 0; k < deg; ++k)
3908  if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3909  nodal_values[((i + this->degree +
3911  deg +
3914  deg +
3916  solution((i * deg + j) * deg + k);
3917 
3918  // Set up the right hand side.
3919  system_rhs = 0;
3920 
3921  for (unsigned int q_point = 0; q_point < n_interior_points;
3922  ++q_point)
3923  {
3924  double tmp =
3925  support_point_values[q_point +
3927  n_edge_points +
3929  n_face_points][2];
3930 
3931  for (unsigned int i = 0; i <= deg; ++i)
3932  for (unsigned int j = 0; j < 4; ++j)
3933  {
3934  tmp -= nodal_values[i + (j + 8) * this->degree] *
3935  this->shape_value_component(
3936  i + (j + 8) * this->degree,
3938  [q_point +
3940  n_edge_points +
3942  n_face_points],
3943  2);
3944 
3945  for (unsigned int k = 0; k < deg; ++k)
3946  tmp -=
3947  nodal_values[i +
3948  ((2 * j + 1) * deg + k +
3950  this->degree] *
3951  this->shape_value_component(
3952  i + ((2 * j + 1) * deg + k +
3954  this->degree,
3955  this->generalized_support_points
3956  [q_point +
3958  n_edge_points +
3960  n_face_points],
3961  2);
3962  }
3963 
3964  for (unsigned int i = 0; i <= deg; ++i)
3965  for (unsigned int j = 0; j < deg; ++j)
3966  for (unsigned int k = 0; k < deg; ++k)
3967  system_rhs((i * deg + j) * deg + k) +=
3968  reference_quadrature.weight(q_point) * tmp *
3969  lobatto_polynomials_grad[i].value(
3971  [q_point +
3973  n_edge_points +
3975  n_face_points](2)) *
3976  lobatto_polynomials[j + 2].value(
3978  [q_point +
3980  n_edge_points +
3982  n_face_points](0)) *
3983  lobatto_polynomials[k + 2].value(
3985  [q_point +
3987  n_edge_points +
3989  n_face_points](1));
3990  }
3991 
3992  system_matrix_inv.vmult(solution, system_rhs);
3993 
3994  // Add the computed support_point_values
3995  // to the resulting vector
3996  // only, if they are not
3997  // too small.
3998  for (unsigned int i = 0; i <= deg; ++i)
3999  for (unsigned int j = 0; j < deg; ++j)
4000  for (unsigned int k = 0; k < deg; ++k)
4001  if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4002  nodal_values
4003  [i +
4004  ((j + 2 * (deg + GeometryInfo<dim>::faces_per_cell)) *
4005  deg +
4007  this->degree] = solution((i * deg + j) * deg + k);
4008  }
4009 
4010  break;
4011  }
4012 
4013  default:
4014  Assert(false, ExcNotImplemented());
4015  }
4016 }
4017 
4018 
4019 
4020 template <int dim>
4021 std::pair<Table<2, bool>, std::vector<unsigned int>>
4023 {
4024  Table<2, bool> constant_modes(dim, this->dofs_per_cell);
4025  for (unsigned int d = 0; d < dim; ++d)
4026  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
4027  constant_modes(d, i) = true;
4028  std::vector<unsigned int> components;
4029  for (unsigned int d = 0; d < dim; ++d)
4030  components.push_back(d);
4031  return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
4032  components);
4033 }
4034 
4035 
4036 template <int dim>
4037 std::size_t
4039 {
4040  Assert(false, ExcNotImplemented());
4041  return 0;
4042 }
4043 
4044 
4045 //----------------------------------------------------------------------//
4046 
4047 
4048 // explicit instantiations
4049 #include "fe_nedelec.inst"
4050 
4051 
4052 DEAL_II_NAMESPACE_CLOSE
size_type m() const
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_nedelec.cc:2525
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:3011
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
Definition: fe_nedelec.cc:2423
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2442
virtual bool hp_constraints_are_implemented() const override
Definition: fe_nedelec.cc:2311
const unsigned int components
Definition: fe_base.h:292
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2493
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
FullMatrix< double > interface_constraints
Definition: fe.h:2468
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:982
const std::vector< Point< dim > > & get_points() const
static unsigned int n_polynomials(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:877
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_nedelec.cc:2275
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
const unsigned int degree
Definition: fe_base.h:298
const Point< dim > & point(const unsigned int i) const
void invert(const FullMatrix< number2 > &M)
STL namespace.
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_nedelec.cc:2956
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
FullMatrix< double > inverse_node_matrix
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_nedelec.cc:4022
std::vector< Point< dim - 1 > > generalized_face_support_points
Definition: fe.h:2499
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double >> &support_point_values, std::vector< double > &nodal_values) const override
Definition: fe_nedelec.cc:3074
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree, bool dg=false)
Definition: fe_nedelec.cc:2000
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
Definition: fe_nedelec.cc:2327
void initialize_support_points(const unsigned int order)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2456
Table< 2, double > boundary_weights
Definition: fe_nedelec.h:364
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1411
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
Definition: fe_nedelec.cc:2318
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
virtual std::string get_name() const =0
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other) const override
Definition: fe_nedelec.cc:2369
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Definition: fe.cc:537
Threads::Mutex mutex
Definition: fe_nedelec.h:369
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
unsigned int n_components() const
virtual std::size_t memory_consumption() const override
Definition: fe_nedelec.cc:4038
virtual std::string get_name() const override
Definition: fe_nedelec.cc:207
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
const unsigned int dofs_per_face
Definition: fe_base.h:275
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
Definition: fe.cc:275
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:699
void initialize_restriction()
Definition: fe_nedelec.cc:516
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
Definition: fe_nedelec.cc:225
void fill(const FullMatrix< number2 > &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0)
double weight(const unsigned int i) const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_nedelec.cc:2036
std::vector< MappingKind > mapping_kind