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