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