Reference documentation for deal.II version 8.5.1
fe_nedelec.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2013 - 2017 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 #include <sstream>
32 #include <iostream>
33 
34 //TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
35 //adjust_line_dof_index_for_line_orientation_table fields, and write tests
36 //similar to bits/face_orientation_and_fe_q_*
37 
38 
39 DEAL_II_NAMESPACE_OPEN
40 
41 //#define DEBUG_NEDELEC
42 
43 namespace internal
44 {
45  namespace
46  {
47  double
48  get_embedding_computation_tolerance(const unsigned int p)
49  {
50  // This heuristic was computed by monitoring the worst residual
51  // resulting from the least squares computation when computing
52  // the face embedding matrices in the FE_Nedelec constructor.
53  // The residual growth is exponential, but is bounded by this
54  // function up to degree 12.
55  return 1.e-15*std::exp(std::pow(p,1.075));
56  }
57  }
58 }
59 
60 
61 template <int dim>
62 FE_Nedelec<dim>::FE_Nedelec (const unsigned int order)
63  :
64  FE_PolyTensor<PolynomialsNedelec<dim>, dim>
65  (order,
66  FiniteElementData<dim> (get_dpo_vector (order), dim, order + 1,
67  FiniteElementData<dim>::Hcurl),
68  std::vector<bool> (PolynomialsNedelec<dim>::compute_n_pols (order), true),
69  std::vector<ComponentMask>
70  (PolynomialsNedelec<dim>::compute_n_pols (order),
71  std::vector<bool> (dim, true)))
72 {
73 #ifdef DEBUG_NEDELEC
74  deallog << get_name() << std::endl;
75 #endif
76 
77  Assert (dim >= 2, ExcImpossibleInDim(dim));
78 
79  const unsigned int n_dofs = this->dofs_per_cell;
80 
81  this->mapping_type = mapping_nedelec;
82  // First, initialize the
83  // generalized support points and
84  // quadrature weights, since they
85  // are required for interpolation.
86  initialize_support_points (order);
87  this->inverse_node_matrix.reinit (n_dofs, n_dofs);
88  this->inverse_node_matrix.fill
89  (FullMatrix<double> (IdentityMatrix (n_dofs)));
90  // From now on, the shape functions
91  // will be the correct ones, not
92  // the raw shape functions anymore.
93 
94  // do not initialize embedding and restriction here. these matrices are
95  // initialized on demand in get_restriction_matrix and
96  // get_prolongation_matrix
97 
98 #ifdef DEBUG_NEDELEC
99  deallog << "Face Embedding" << std::endl;
100 #endif
101  FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
102 
103  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
104  face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face);
105 
106  FETools::compute_face_embedding_matrices<dim,double>
107  (*this, face_embeddings, 0, 0,
108  internal::get_embedding_computation_tolerance(order));
109 
110  switch (dim)
111  {
112  case 1:
113  {
114  this->interface_constraints.reinit (0, 0);
115  break;
116  }
117 
118  case 2:
119  {
120  this->interface_constraints.reinit (2 * this->dofs_per_face,
121  this->dofs_per_face);
122 
123  for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
124  ++i)
125  for (unsigned int j = 0; j < this->dofs_per_face; ++j)
126  for (unsigned int k = 0; k < this->dofs_per_face; ++k)
127  this->interface_constraints (i * this->dofs_per_face + j, k)
128  = face_embeddings[i] (j, k);
129 
130  break;
131  }
132 
133  case 3:
134  {
135  this->interface_constraints.reinit
136  (4 * (this->dofs_per_face - this->degree), this->dofs_per_face);
137 
138  unsigned int target_row = 0;
139 
140  for (unsigned int i = 0; i < 2; ++i)
141  for (unsigned int j = this->degree; j < 2 * this->degree;
142  ++j, ++target_row)
143  for (unsigned int k = 0; k < this->dofs_per_face; ++k)
144  this->interface_constraints (target_row, k)
145  = face_embeddings[2 * i] (j, k);
146 
147  for (unsigned int i = 0; i < 2; ++i)
148  for (unsigned int j = 3 * this->degree;
149  j < GeometryInfo<3>::lines_per_face * this->degree;
150  ++j, ++target_row)
151  for (unsigned int k = 0; k < this->dofs_per_face; ++k)
152  this->interface_constraints (target_row, k)
153  = face_embeddings[i] (j, k);
154 
155  for (unsigned int i = 0; i < 2; ++i)
156  for (unsigned int j = 0; j < 2; ++j)
157  for (unsigned int k = i * this->degree;
158  k < (i + 1) * this->degree; ++k, ++target_row)
159  for (unsigned int l = 0; l < this->dofs_per_face; ++l)
160  this->interface_constraints (target_row, l)
161  = face_embeddings[i + 2 * j] (k, l);
162 
163  for (unsigned int i = 0; i < 2; ++i)
164  for (unsigned int j = 0; j < 2; ++j)
165  for (unsigned int k = (i + 2) * this->degree;
166  k < (i + 3) * this->degree; ++k, ++target_row)
167  for (unsigned int l = 0; l < this->dofs_per_face; ++l)
168  this->interface_constraints (target_row, l)
169  = face_embeddings[2 * i + j] (k, l);
170 
171  for (unsigned int i = 0; i < GeometryInfo<3>::max_children_per_face;
172  ++i)
173  for (unsigned int
174  j = GeometryInfo<3>::lines_per_face * this->degree;
175  j < this->dofs_per_face; ++j, ++target_row)
176  for (unsigned int k = 0; k < this->dofs_per_face; ++k)
177  this->interface_constraints (target_row, k)
178  = face_embeddings[i] (j, k);
179 
180  break;
181  }
182 
183  default:
184  Assert (false, ExcNotImplemented ());
185  }
186 
187 }
188 
189 
190 
191 template <int dim>
192 std::string
193 FE_Nedelec<dim>::get_name () const
194 {
195  // note that the
196  // FETools::get_fe_by_name
197  // function depends on the
198  // particular format of the string
199  // this function returns, so they
200  // have to be kept in synch
201 
202  std::ostringstream namebuf;
203  namebuf << "FE_Nedelec<" << dim << ">(" << this->degree-1 << ")";
204 
205  return namebuf.str();
206 }
207 
208 
209 template <int dim>
210 FiniteElement<dim>
211 *FE_Nedelec<dim>::clone () const
212 {
213  return new FE_Nedelec<dim> (*this);
214 }
215 
216 //---------------------------------------------------------------------------
217 // Auxiliary and internal functions
218 //---------------------------------------------------------------------------
219 
220 
221 
222 // Set the generalized support
223 // points and precompute the
224 // parts of the projection-based
225 // interpolation, which does
226 // not depend on the interpolated
227 // function.
228 template <>
229 void
230 FE_Nedelec<1>::initialize_support_points (const unsigned int)
231 {
232  Assert (false, ExcNotImplemented ());
233 }
234 
235 
236 
237 template <>
238 void
239 FE_Nedelec<2>::initialize_support_points (const unsigned int order)
240 {
241  const int dim = 2;
242 
243  // Create polynomial basis.
244  const std::vector<Polynomials::Polynomial<double> > &lobatto_polynomials
245  = Polynomials::Lobatto::generate_complete_basis (order + 1);
246  std::vector<Polynomials::Polynomial<double> >
247  lobatto_polynomials_grad (order + 1);
248 
249  for (unsigned int i = 0; i < lobatto_polynomials_grad.size (); ++i)
250  lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative ();
251 
252  // Initialize quadratures to obtain
253  // quadrature points later on.
254  const QGauss<dim - 1> reference_edge_quadrature (order + 1);
255  const unsigned int n_edge_points = reference_edge_quadrature.size ();
256  const unsigned int n_boundary_points
257  = GeometryInfo<dim>::lines_per_cell * n_edge_points;
258  const Quadrature<dim> edge_quadrature
259  = QProjector<dim>::project_to_all_faces (reference_edge_quadrature);
260 
261  this->generalized_face_support_points.resize (n_edge_points);
262 
263  // Create face support points.
264  for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
265  this->generalized_face_support_points[q_point]
266  = reference_edge_quadrature.point (q_point);
267 
268  if (order > 0)
269  {
270  // If the polynomial degree is positive
271  // we have support points on the faces
272  // and in the interior of a cell.
273  const QGauss<dim> quadrature (order + 1);
274  const unsigned int &n_interior_points = quadrature.size ();
275 
276  this->generalized_support_points.resize
277  (n_boundary_points + n_interior_points);
278  boundary_weights.reinit (n_edge_points, order);
279 
280  for (unsigned int q_point = 0; q_point < n_edge_points;
281  ++q_point)
282  {
283  for (unsigned int line = 0;
284  line < GeometryInfo<dim>::lines_per_cell; ++line)
285  this->generalized_support_points[line * n_edge_points
286  + q_point]
287  = edge_quadrature.point
288  (QProjector<dim>::DataSetDescriptor::face
289  (line, true, false, false, n_edge_points) + q_point);
290 
291  for (unsigned int i = 0; i < order; ++i)
292  boundary_weights (q_point, i)
293  = reference_edge_quadrature.weight (q_point)
294  * lobatto_polynomials_grad[i + 1].value
295  (this->generalized_face_support_points[q_point] (0));
296  }
297 
298  for (unsigned int q_point = 0; q_point < n_interior_points;
299  ++q_point)
300  this->generalized_support_points[q_point + n_boundary_points]
301  = quadrature.point (q_point);
302  }
303 
304  else
305  {
306  // In this case we only need support points
307  // on the faces of a cell.
308  this->generalized_support_points.resize (n_boundary_points);
309 
310  for (unsigned int line = 0;
311  line < GeometryInfo<dim>::lines_per_cell; ++line)
312  for (unsigned int q_point = 0; q_point < n_edge_points;
313  ++q_point)
314  this->generalized_support_points[line * n_edge_points
315  + q_point]
316  = edge_quadrature.point
317  (QProjector<dim>::DataSetDescriptor::face
318  (line, true, false, false, n_edge_points) + q_point);
319  }
320 }
321 
322 
323 
324 template <>
325 void
326 FE_Nedelec<3>::initialize_support_points (const unsigned int order)
327 {
328  const int dim = 3;
329 
330  // Create polynomial basis.
331  const std::vector<Polynomials::Polynomial<double> > &lobatto_polynomials
332  = Polynomials::Lobatto::generate_complete_basis (order + 1);
333  std::vector<Polynomials::Polynomial<double> >
334  lobatto_polynomials_grad (order + 1);
335 
336  for (unsigned int i = 0; i < lobatto_polynomials_grad.size (); ++i)
337  lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative ();
338 
339  // Initialize quadratures to obtain
340  // quadrature points later on.
341  const QGauss<1> reference_edge_quadrature (order + 1);
342  const unsigned int &n_edge_points = reference_edge_quadrature.size ();
343  const Quadrature<dim - 1>& edge_quadrature
344  = QProjector<dim - 1>::project_to_all_faces
345  (reference_edge_quadrature);
346 
347  if (order > 0)
348  {
349  // If the polynomial order is positive
350  // we have support points on the edges,
351  // faces and in the interior of a cell.
352  const QGauss<dim - 1> reference_face_quadrature (order + 1);
353  const unsigned int &n_face_points
354  = reference_face_quadrature.size ();
355  const unsigned int n_boundary_points
356  = GeometryInfo<dim>::lines_per_cell * n_edge_points
357  + GeometryInfo<dim>::faces_per_cell * n_face_points;
358  const QGauss<dim> quadrature (order + 1);
359  const unsigned int &n_interior_points = quadrature.size ();
360 
361  boundary_weights.reinit (n_edge_points + n_face_points,
362  2 * (order + 1) * order);
363  this->generalized_face_support_points.resize
364  (4 * n_edge_points + n_face_points);
365  this->generalized_support_points.resize
366  (n_boundary_points + n_interior_points);
367 
368  // Create support points on edges.
369  for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
370  {
371  for (unsigned int line = 0;
372  line < GeometryInfo<dim - 1>::lines_per_cell; ++line)
373  this->generalized_face_support_points[line * n_edge_points
374  + q_point]
375  = edge_quadrature.point
376  (QProjector<dim - 1>::DataSetDescriptor::face
377  (line, true, false, false, n_edge_points) + q_point);
378 
379  for (unsigned int i = 0; i < 2; ++i)
380  for (unsigned int j = 0; j < 2; ++j)
381  {
382  this->generalized_support_points
383  [q_point + (i + 4 * j) * n_edge_points]
384  = Point<dim>
385  (i, reference_edge_quadrature.point (q_point) (0),
386  j);
387  this->generalized_support_points
388  [q_point + (i + 4 * j + 2) * n_edge_points]
389  = Point<dim>
390  (reference_edge_quadrature.point (q_point) (0),
391  i, j);
392  this->generalized_support_points
393  [q_point + (i + 2 * (j + 4)) * n_edge_points]
394  = Point<dim>
395  (i, j,
396  reference_edge_quadrature.point (q_point) (0));
397  }
398 
399  for (unsigned int i = 0; i < order; ++i)
400  boundary_weights (q_point, i)
401  = reference_edge_quadrature.weight (q_point)
402  * lobatto_polynomials_grad[i + 1].value
403  (this->generalized_face_support_points[q_point] (1));
404  }
405 
406  // Create support points on faces.
407  for (unsigned int q_point = 0; q_point < n_face_points;
408  ++q_point)
409  {
410  this->generalized_face_support_points[q_point
411  + 4 * n_edge_points]
412  = reference_face_quadrature.point (q_point);
413 
414  for (unsigned int i = 0; i <= order; ++i)
415  for (unsigned int j = 0; j < order; ++j)
416  {
417  boundary_weights (q_point + n_edge_points,
418  2 * (i * order + j))
419  = reference_face_quadrature.weight (q_point)
420  * lobatto_polynomials_grad[i].value
421  (this->generalized_face_support_points
422  [q_point + 4 * n_edge_points] (0))
423  * lobatto_polynomials[j + 2].value
424  (this->generalized_face_support_points
425  [q_point + 4 * n_edge_points] (1));
426  boundary_weights (q_point + n_edge_points,
427  2 * (i * order + j) + 1)
428  = reference_face_quadrature.weight (q_point)
429  * lobatto_polynomials_grad[i].value
430  (this->generalized_face_support_points
431  [q_point + 4 * n_edge_points] (1))
432  * lobatto_polynomials[j + 2].value
433  (this->generalized_face_support_points
434  [q_point + 4 * n_edge_points] (0));
435  }
436  }
437 
438  const Quadrature<dim> &face_quadrature
439  = QProjector<dim>::project_to_all_faces
440  (reference_face_quadrature);
441 
442  for (unsigned int face = 0;
443  face < GeometryInfo<dim>::faces_per_cell; ++face)
444  for (unsigned int q_point = 0; q_point < n_face_points;
445  ++q_point)
446  {
447  this->generalized_support_points
448  [face * n_face_points + q_point
449  + GeometryInfo<dim>::lines_per_cell * n_edge_points]
450  = face_quadrature.point
451  (QProjector<dim>::DataSetDescriptor::face
452  (face, true, false, false, n_face_points) + q_point);
453  }
454 
455  // Create support points in the interior.
456  for (unsigned int q_point = 0; q_point < n_interior_points;
457  ++q_point)
458  this->generalized_support_points[q_point + n_boundary_points]
459  = quadrature.point (q_point);
460  }
461 
462  else
463  {
464  this->generalized_face_support_points.resize (4 * n_edge_points);
465  this->generalized_support_points.resize
466  (GeometryInfo<dim>::lines_per_cell * n_edge_points);
467 
468  for (unsigned int q_point = 0; q_point < n_edge_points;
469  ++q_point)
470  {
471  for (unsigned int line = 0;
472  line < GeometryInfo<dim - 1>::lines_per_cell; ++line)
473  this->generalized_face_support_points[line * n_edge_points
474  + q_point]
475  = edge_quadrature.point
476  (QProjector<dim - 1>::DataSetDescriptor::face
477  (line, true, false, false, n_edge_points) + q_point);
478 
479  for (unsigned int i = 0; i < 2; ++i)
480  for (unsigned int j = 0; j < 2; ++j)
481  {
482  this->generalized_support_points
483  [q_point + (i + 4 * j) * n_edge_points]
484  = Point<dim>
485  (i, reference_edge_quadrature.point (q_point) (0),
486  j);
487  this->generalized_support_points
488  [q_point + (i + 4 * j + 2) * n_edge_points]
489  = Point<dim>
490  (reference_edge_quadrature.point (q_point) (0),
491  i, j);
492  this->generalized_support_points
493  [q_point + (i + 2 * (j + 4)) * n_edge_points]
494  = Point<dim>
495  (i, j,
496  reference_edge_quadrature.point (q_point) (0));
497  }
498  }
499  }
500 }
501 
502 
503 
504 // Set the restriction matrices.
505 template <>
506 void
507 FE_Nedelec<1>::initialize_restriction ()
508 {
509  // there is only one refinement case in 1d,
510  // which is the isotropic one
511  for (unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
512  this->restriction[0][i].reinit(0, 0);
513 }
514 
515 
516 
517 // Restriction operator
518 template <int dim>
519 void
520 FE_Nedelec<dim>::initialize_restriction ()
521 {
522  // This function does the same as the
523  // function interpolate further below.
524  // But since the functions, which we
525  // interpolate here, are discontinuous
526  // we have to use more quadrature
527  // points as in interpolate.
528  const QGauss<1> edge_quadrature (2 * this->degree);
529  const std::vector<Point<1> > &edge_quadrature_points
530  = edge_quadrature.get_points ();
531  const unsigned int &
532  n_edge_quadrature_points = edge_quadrature.size ();
533  const unsigned int
534  index = RefinementCase<dim>::isotropic_refinement - 1;
535 
536  switch (dim)
537  {
538  case 2:
539  {
540  // First interpolate the shape
541  // functions of the child cells
542  // to the lowest order shape
543  // functions of the parent cell.
544  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
545  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
546  ++q_point)
547  {
548  const double weight = 2.0 * edge_quadrature.weight (q_point);
549 
550  if (edge_quadrature_points[q_point] (0) < 0.5)
551  {
552  Point<dim> quadrature_point (0.0,
553  2.0 * edge_quadrature_points[q_point] (0));
554 
555  this->restriction[index][0] (0, dof) += weight
556  * this->shape_value_component
557  (dof,
558  quadrature_point,
559  1);
560  quadrature_point (0) = 1.0;
561  this->restriction[index][1] (this->degree, dof)
562  += weight * this->shape_value_component (dof,
563  quadrature_point,
564  1);
565  quadrature_point (0) = quadrature_point (1);
566  quadrature_point (1) = 0.0;
567  this->restriction[index][0] (2 * this->degree, dof)
568  += weight * this->shape_value_component (dof,
569  quadrature_point,
570  0);
571  quadrature_point (1) = 1.0;
572  this->restriction[index][2] (3 * this->degree, dof)
573  += weight * this->shape_value_component (dof,
574  quadrature_point,
575  0);
576  }
577 
578  else
579  {
580  Point<dim> quadrature_point (0.0,
581  2.0 * edge_quadrature_points[q_point] (0)
582  - 1.0);
583 
584  this->restriction[index][2] (0, dof) += weight
585  * this->shape_value_component
586  (dof,
587  quadrature_point,
588  1);
589  quadrature_point (0) = 1.0;
590  this->restriction[index][3] (this->degree, dof)
591  += weight * this->shape_value_component (dof,
592  quadrature_point,
593  1);
594  quadrature_point (0) = quadrature_point (1);
595  quadrature_point (1) = 0.0;
596  this->restriction[index][1] (2 * this->degree, dof)
597  += weight * this->shape_value_component (dof,
598  quadrature_point,
599  0);
600  quadrature_point (1) = 1.0;
601  this->restriction[index][3] (3 * this->degree, dof)
602  += weight * this->shape_value_component (dof,
603  quadrature_point,
604  0);
605  }
606  }
607 
608  // Then project the shape functions
609  // of the child cells to the higher
610  // order shape functions of the
611  // parent cell.
612  if (this->degree > 1)
613  {
614  const unsigned int deg = this->degree-1;
615  const std::vector<Polynomials::Polynomial<double> > &
616  legendre_polynomials
617  = Polynomials::Legendre::generate_complete_basis (deg);
618  FullMatrix<double> system_matrix_inv (deg, deg);
619 
620  {
621  FullMatrix<double> assembling_matrix (deg,
622  n_edge_quadrature_points);
623 
624  for (unsigned int q_point = 0;
625  q_point < n_edge_quadrature_points; ++q_point)
626  {
627  const double weight
628  = std::sqrt (edge_quadrature.weight (q_point));
629 
630  for (unsigned int i = 0; i < deg; ++i)
631  assembling_matrix (i, q_point) = weight
632  * legendre_polynomials[i + 1].value
633  (edge_quadrature_points[q_point] (0));
634  }
635 
636  FullMatrix<double> system_matrix (deg, deg);
637 
638  assembling_matrix.mTmult (system_matrix, assembling_matrix);
639  system_matrix_inv.invert (system_matrix);
640  }
641 
642  FullMatrix<double> solution (this->degree-1, 4);
643  FullMatrix<double> system_rhs (this->degree-1, 4);
644  Vector<double> tmp (4);
645 
646  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
647  for (unsigned int i = 0; i < 2; ++i)
648  {
649  system_rhs = 0.0;
650 
651  for (unsigned int q_point = 0;
652  q_point < n_edge_quadrature_points; ++q_point)
653  {
654  const double weight
655  = edge_quadrature.weight (q_point);
656  const Point<dim> quadrature_point_0 (i,
657  edge_quadrature_points[q_point] (0));
658  const Point<dim> quadrature_point_1
659  (edge_quadrature_points[q_point] (0),
660  i);
661 
662  if (edge_quadrature_points[q_point] (0) < 0.5)
663  {
664  Point<dim> quadrature_point_2 (i,
665  2.0 * edge_quadrature_points[q_point] (0));
666 
667  tmp (0) = weight
668  * (2.0 * this->shape_value_component
669  (dof, quadrature_point_2, 1)
670  - this->restriction[index][i]
671  (i * this->degree, dof)
672  * this->shape_value_component
673  (i * this->degree,
674  quadrature_point_0, 1));
675  tmp (1) = -1.0 * weight
676  * this->restriction[index][i + 2]
677  (i * this->degree, dof)
678  * this->shape_value_component
679  (i * this->degree,
680  quadrature_point_0, 1);
681  quadrature_point_2
682  = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
683  i);
684  tmp (2) = weight
685  * (2.0 * this->shape_value_component
686  (dof, quadrature_point_2, 0)
687  - this->restriction[index][2 * i]
688  ((i + 2) * this->degree, dof)
689  * this->shape_value_component
690  ((i + 2) * this->degree,
691  quadrature_point_1, 0));
692  tmp (3) = -1.0 * weight
693  * this->restriction[index][2 * i + 1]
694  ((i + 2) * this->degree, dof)
695  * this->shape_value_component
696  ((i + 2) * this->degree,
697  quadrature_point_1, 0);
698  }
699 
700  else
701  {
702  tmp (0) = -1.0 * weight
703  * this->restriction[index][i]
704  (i * this->degree, dof)
705  * this->shape_value_component
706  (i * this->degree,
707  quadrature_point_0, 1);
708 
709  Point<dim> quadrature_point_2 (i,
710  2.0 * edge_quadrature_points[q_point] (0)
711  - 1.0);
712 
713  tmp (1) = weight
714  * (2.0 * this->shape_value_component
715  (dof, quadrature_point_2, 1)
716  - this->restriction[index][i + 2]
717  (i * this->degree, dof)
718  * this->shape_value_component
719  (i * this->degree,
720  quadrature_point_0, 1));
721  tmp (2) = -1.0 * weight
722  * this->restriction[index][2 * i]
723  ((i + 2) * this->degree, dof)
724  * this->shape_value_component
725  ((i + 2) * this->degree,
726  quadrature_point_1, 0);
727  quadrature_point_2
728  = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
729  - 1.0, i);
730  tmp (3) = weight
731  * (2.0 * this->shape_value_component
732  (dof, quadrature_point_2, 0)
733  - this->restriction[index][2 * i + 1]
734  ((i + 2) * this->degree, dof)
735  * this->shape_value_component
736  ((i + 2) * this->degree,
737  quadrature_point_1, 0));
738  }
739 
740  for (unsigned int j = 0; j < this->degree-1; ++j)
741  {
742  const double L_j
743  = legendre_polynomials[j + 1].value
744  (edge_quadrature_points[q_point] (0));
745 
746  for (unsigned int k = 0; k < tmp.size (); ++k)
747  system_rhs (j, k) += tmp (k) * L_j;
748  }
749  }
750 
751  system_matrix_inv.mmult (solution, system_rhs);
752 
753  for (unsigned int j = 0; j < this->degree-1; ++j)
754  for (unsigned int k = 0; k < 2; ++k)
755  {
756  if (std::abs (solution (j, k)) > 1e-14)
757  this->restriction[index][i + 2 * k]
758  (i * this->degree + j + 1, dof)
759  = solution (j, k);
760 
761  if (std::abs (solution (j, k + 2)) > 1e-14)
762  this->restriction[index][2 * i + k]
763  ((i + 2) * this->degree + j + 1, dof)
764  = solution (j, k + 2);
765  }
766  }
767 
768  const QGauss<dim> quadrature (2 * this->degree);
769  const std::vector<Point<dim> > &
770  quadrature_points = quadrature.get_points ();
771  const std::vector<Polynomials::Polynomial<double> > &
772  lobatto_polynomials
773  = Polynomials::Lobatto::generate_complete_basis
774  (this->degree);
775  const unsigned int n_boundary_dofs
776  = GeometryInfo<dim>::faces_per_cell * this->degree;
777  const unsigned int &n_quadrature_points = quadrature.size ();
778 
779  {
780  FullMatrix<double> assembling_matrix ((this->degree-1) * this->degree,
781  n_quadrature_points);
782 
783  for (unsigned int q_point = 0; q_point < n_quadrature_points;
784  ++q_point)
785  {
786  const double weight
787  = std::sqrt (quadrature.weight (q_point));
788 
789  for (unsigned int i = 0; i < this->degree; ++i)
790  {
791  const double L_i = weight
792  * legendre_polynomials[i].value
793  (quadrature_points[q_point] (0));
794 
795  for (unsigned int j = 0; j < this->degree-1; ++j)
796  assembling_matrix (i * (this->degree-1) + j, q_point)
797  = L_i * lobatto_polynomials[j + 2].value
798  (quadrature_points[q_point] (1));
799  }
800  }
801 
802  FullMatrix<double> system_matrix (assembling_matrix.m (),
803  assembling_matrix.m ());
804 
805  assembling_matrix.mTmult (system_matrix, assembling_matrix);
806  system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ());
807  system_matrix_inv.invert (system_matrix);
808  }
809 
810  solution.reinit (system_matrix_inv.m (), 8);
811  system_rhs.reinit (system_matrix_inv.m (), 8);
812  tmp.reinit (8);
813 
814  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
815  {
816  system_rhs = 0.0;
817 
818  for (unsigned int q_point = 0;
819  q_point < n_quadrature_points; ++q_point)
820  {
821  tmp = 0.0;
822 
823  if (quadrature_points[q_point] (0) < 0.5)
824  {
825  if (quadrature_points[q_point] (1) < 0.5)
826  {
827  const Point<dim> quadrature_point
828  (2.0 * quadrature_points[q_point] (0),
829  2.0 * quadrature_points[q_point] (1));
830 
831  tmp (0) += 2.0 * this->shape_value_component
832  (dof, quadrature_point, 0);
833  tmp (1) += 2.0 * this->shape_value_component
834  (dof, quadrature_point, 1);
835  }
836 
837  else
838  {
839  const Point<dim> quadrature_point
840  (2.0 * quadrature_points[q_point] (0),
841  2.0 * quadrature_points[q_point] (1)
842  - 1.0);
843 
844  tmp (4) += 2.0 * this->shape_value_component
845  (dof, quadrature_point, 0);
846  tmp (5) += 2.0 * this->shape_value_component
847  (dof, quadrature_point, 1);
848  }
849  }
850 
851  else if (quadrature_points[q_point] (1) < 0.5)
852  {
853  const Point<dim> quadrature_point
854  (2.0 * quadrature_points[q_point] (0)
855  - 1.0,
856  2.0 * quadrature_points[q_point] (1));
857 
858  tmp (2) += 2.0 * this->shape_value_component
859  (dof, quadrature_point, 0);
860  tmp (3) += 2.0 * this->shape_value_component
861  (dof, quadrature_point, 1);
862  }
863 
864  else
865  {
866  const Point<dim> quadrature_point
867  (2.0 * quadrature_points[q_point] (0)
868  - 1.0,
869  2.0 * quadrature_points[q_point] (1)
870  - 1.0);
871 
872  tmp (6) += 2.0 * this->shape_value_component
873  (dof, quadrature_point, 0);
874  tmp (7) += 2.0 * this->shape_value_component
875  (dof, quadrature_point, 1);
876  }
877 
878  for (unsigned int i = 0; i < 2; ++i)
879  for (unsigned int j = 0; j < this->degree; ++j)
880  {
881  tmp (2 * i) -= this->restriction[index][i]
882  (j + 2 * this->degree, dof)
883  * this->shape_value_component
884  (j + 2 * this->degree,
885  quadrature_points[q_point], 0);
886  tmp (2 * i + 1) -= this->restriction[index][i]
887  (i * this->degree + j, dof)
888  * this->shape_value_component
889  (i * this->degree + j,
890  quadrature_points[q_point], 1);
891  tmp (2 * (i + 2)) -= this->restriction[index][i + 2]
892  (j + 3 * this->degree, dof)
893  * this->shape_value_component
894  (j + 3 * this->degree,
895  quadrature_points[q_point],
896  0);
897  tmp (2 * i + 5) -= this->restriction[index][i + 2]
898  (i * this->degree + j, dof)
899  * this->shape_value_component
900  (i * this->degree + j,
901  quadrature_points[q_point], 1);
902  }
903 
904  tmp *= quadrature.weight (q_point);
905 
906  for (unsigned int i = 0; i < this->degree; ++i)
907  {
908  const double L_i_0
909  = legendre_polynomials[i].value
910  (quadrature_points[q_point] (0));
911  const double L_i_1
912  = legendre_polynomials[i].value
913  (quadrature_points[q_point] (1));
914 
915  for (unsigned int j = 0; j < this->degree-1; ++j)
916  {
917  const double l_j_0
918  = L_i_0 * lobatto_polynomials[j + 2].value
919  (quadrature_points[q_point] (1));
920  const double l_j_1
921  = L_i_1 * lobatto_polynomials[j + 2].value
922  (quadrature_points[q_point] (0));
923 
924  for (unsigned int k = 0; k < 4; ++k)
925  {
926  system_rhs (i * (this->degree-1) + j, 2 * k)
927  += tmp (2 * k) * l_j_0;
928  system_rhs (i * (this->degree-1) + j, 2 * k + 1)
929  += tmp (2 * k + 1) * l_j_1;
930  }
931  }
932  }
933  }
934 
935  system_matrix_inv.mmult (solution, system_rhs);
936 
937  for (unsigned int i = 0; i < this->degree; ++i)
938  for (unsigned int j = 0; j < this->degree-1; ++j)
939  for (unsigned int k = 0; k < 4; ++k)
940  {
941  if (std::abs (solution (i * (this->degree-1) + j, 2 * k))
942  > 1e-14)
943  this->restriction[index][k]
944  (i * (this->degree-1) + j + n_boundary_dofs, dof)
945  = solution (i * (this->degree-1) + j, 2 * k);
946 
947  if (std::abs (solution (i * (this->degree-1) + j, 2 * k + 1))
948  > 1e-14)
949  this->restriction[index][k]
950  (i + (this->degree-1 + j) * this->degree + n_boundary_dofs,
951  dof)
952  = solution (i * (this->degree-1) + j, 2 * k + 1);
953  }
954  }
955  }
956 
957  break;
958  }
959 
960  case 3:
961  {
962  // First interpolate the shape
963  // functions of the child cells
964  // to the lowest order shape
965  // functions of the parent cell.
966  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
967  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
968  ++q_point)
969  {
970  const double weight = 2.0 * edge_quadrature.weight (q_point);
971 
972  if (edge_quadrature_points[q_point] (0) < 0.5)
973  for (unsigned int i = 0; i < 2; ++i)
974  for (unsigned int j = 0; j < 2; ++j)
975  {
976  Point<dim> quadrature_point (i,
977  2.0 * edge_quadrature_points[q_point] (0),
978  j);
979 
980  this->restriction[index][i + 4 * j]
981  ((i + 4 * j) * this->degree, dof)
982  += weight * this->shape_value_component (dof,
983  quadrature_point,
984  1);
985  quadrature_point
986  = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
987  i, j);
988  this->restriction[index][2 * (i + 2 * j)]
989  ((i + 4 * j + 2) * this->degree, dof)
990  += weight * this->shape_value_component (dof,
991  quadrature_point,
992  0);
993  quadrature_point = Point<dim> (i, j,
994  2.0 * edge_quadrature_points[q_point] (0));
995  this->restriction[index][i + 2 * j]
996  ((i + 2 * (j + 4)) * this->degree, dof)
997  += weight * this->shape_value_component (dof,
998  quadrature_point,
999  2);
1000  }
1001 
1002  else
1003  for (unsigned int i = 0; i < 2; ++i)
1004  for (unsigned int j = 0; j < 2; ++j)
1005  {
1006  Point<dim> quadrature_point (i,
1007  2.0 * edge_quadrature_points[q_point] (0)
1008  - 1.0, j);
1009 
1010  this->restriction[index][i + 4 * j + 2]
1011  ((i + 4 * j) * this->degree, dof)
1012  += weight * this->shape_value_component (dof,
1013  quadrature_point,
1014  1);
1015  quadrature_point
1016  = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
1017  - 1.0, i, j);
1018  this->restriction[index][2 * (i + 2 * j) + 1]
1019  ((i + 4 * j + 2) * this->degree, dof)
1020  += weight * this->shape_value_component (dof,
1021  quadrature_point,
1022  0);
1023  quadrature_point = Point<dim> (i, j,
1024  2.0 * edge_quadrature_points[q_point] (0)
1025  - 1.0);
1026  this->restriction[index][i + 2 * (j + 2)]
1027  ((i + 2 * (j + 4)) * this->degree, dof)
1028  += weight * this->shape_value_component (dof,
1029  quadrature_point,
1030  2);
1031  }
1032  }
1033 
1034  // Then project the shape functions
1035  // of the child cells to the higher
1036  // order shape functions of the
1037  // parent cell.
1038  if (this->degree > 1)
1039  {
1040  const unsigned int deg = this->degree-1;
1041  const std::vector<Polynomials::Polynomial<double> > &
1042  legendre_polynomials
1043  = Polynomials::Legendre::generate_complete_basis (deg);
1044  FullMatrix<double> system_matrix_inv (deg, deg);
1045 
1046  {
1047  FullMatrix<double> assembling_matrix (deg,
1048  n_edge_quadrature_points);
1049 
1050  for (unsigned int q_point = 0;
1051  q_point < n_edge_quadrature_points; ++q_point)
1052  {
1053  const double weight = std::sqrt (edge_quadrature.weight
1054  (q_point));
1055 
1056  for (unsigned int i = 0; i < deg; ++i)
1057  assembling_matrix (i, q_point) = weight
1058  * legendre_polynomials[i + 1].value
1059  (edge_quadrature_points[q_point] (0));
1060  }
1061 
1062  FullMatrix<double> system_matrix (deg, deg);
1063 
1064  assembling_matrix.mTmult (system_matrix, assembling_matrix);
1065  system_matrix_inv.invert (system_matrix);
1066  }
1067 
1068  FullMatrix<double> solution (deg, 6);
1069  FullMatrix<double> system_rhs (deg, 6);
1070  Vector<double> tmp (6);
1071 
1072  for (unsigned int i = 0; i < 2; ++i)
1073  for (unsigned int j = 0; j < 2; ++j)
1074  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
1075  {
1076  system_rhs = 0.0;
1077 
1078  for (unsigned int q_point = 0;
1079  q_point < n_edge_quadrature_points; ++q_point)
1080  {
1081  const double weight = edge_quadrature.weight
1082  (q_point);
1083  const Point<dim> quadrature_point_0 (i,
1084  edge_quadrature_points[q_point] (0),
1085  j);
1086  const Point<dim>
1087  quadrature_point_1
1088  (edge_quadrature_points[q_point] (0), i, j);
1089  const Point<dim> quadrature_point_2 (i, j,
1090  edge_quadrature_points[q_point] (0));
1091 
1092  if (edge_quadrature_points[q_point] (0) < 0.5)
1093  {
1094  Point<dim> quadrature_point_3 (i,
1095  2.0 * edge_quadrature_points[q_point] (0),
1096  j);
1097 
1098  tmp (0) = weight
1099  * (2.0 * this->shape_value_component
1100  (dof, quadrature_point_3, 1)
1101  - this->restriction[index][i + 4 * j]
1102  ((i + 4 * j) * this->degree,
1103  dof)
1104  * this->shape_value_component
1105  ((i + 4 * j) * this->degree,
1106  quadrature_point_0, 1));
1107  tmp (1) = -1.0 * weight
1108  * this->restriction[index][i + 4 * j + 2]
1109  ((i + 4 * j) * this->degree,
1110  dof)
1111  * this->shape_value_component
1112  ((i + 4 * j) * this->degree,
1113  quadrature_point_0, 1);
1114  quadrature_point_3
1115  = Point<dim> (2.0 * edge_quadrature_points[q_point] (0),
1116  i, j);
1117  tmp (2) = weight
1118  * (2.0 * this->shape_value_component
1119  (dof, quadrature_point_3, 0)
1120  - this->restriction[index][2 * (i + 2 * j)]
1121  ((i + 4 * j + 2) * this->degree,
1122  dof)
1123  * this->shape_value_component
1124  ((i + 4 * j + 2) * this->degree,
1125  quadrature_point_1, 0));
1126  tmp (3) = -1.0 * weight
1127  * this->restriction[index][2 * (i + 2 * j) + 1]
1128  ((i + 4 * j + 2) * this->degree,
1129  dof)
1130  * this->shape_value_component
1131  ((i + 4 * j + 2) * this->degree,
1132  quadrature_point_1, 0);
1133  quadrature_point_3 = Point<dim> (i, j,
1134  2.0 * edge_quadrature_points[q_point] (0));
1135  tmp (4) = weight
1136  * (2.0 * this->shape_value_component
1137  (dof, quadrature_point_3, 2)
1138  - this->restriction[index][i + 2 * j]
1139  ((i + 2 * (j + 4)) * this->degree,
1140  dof)
1141  * this->shape_value_component
1142  ((i + 2 * (j + 4)) * this->degree,
1143  quadrature_point_2, 2));
1144  tmp (5) = -1.0 * weight
1145  * this->restriction[index][i + 2 * (j + 2)]
1146  ((i + 2 * (j + 4)) * this->degree,
1147  dof)
1148  * this->shape_value_component
1149  ((i + 2 * (j + 4)) * this->degree,
1150  quadrature_point_2, 2);
1151  }
1152 
1153  else
1154  {
1155  tmp (0) = -1.0 * weight
1156  * this->restriction[index][i + 4 * j]
1157  ((i + 4 * j) * this->degree,
1158  dof)
1159  * this->shape_value_component
1160  ((i + 4 * j) * this->degree,
1161  quadrature_point_0, 1);
1162 
1163  Point<dim> quadrature_point_3 (i,
1164  2.0 * edge_quadrature_points[q_point] (0)
1165  - 1.0, j);
1166 
1167  tmp (1) = weight
1168  * (2.0 * this->shape_value_component
1169  (dof, quadrature_point_3, 1)
1170  - this->restriction[index][i + 4 * j + 2]
1171  ((i + 4 * j) * this->degree,
1172  dof)
1173  * this->shape_value_component
1174  ((i + 4 * j) * this->degree,
1175  quadrature_point_0, 1));
1176  tmp (2) = -1.0 * weight
1177  * this->restriction[index][2 * (i + 2 * j)]
1178  ((i + 4 * j + 2) * this->degree,
1179  dof)
1180  * this->shape_value_component
1181  ((i + 4 * j + 2) * this->degree,
1182  quadrature_point_1, 0);
1183  quadrature_point_3
1184  = Point<dim> (2.0 * edge_quadrature_points[q_point] (0)
1185  - 1.0, i, j);
1186  tmp (3) = weight
1187  * (2.0 * this->shape_value_component
1188  (dof, quadrature_point_3, 0)
1189  - this->restriction[index][2 * (i + 2 * j) + 1]
1190  ((i + 4 * j + 2) * this->degree,
1191  dof)
1192  * this->shape_value_component
1193  ((i + 4 * j + 2) * this->degree,
1194  quadrature_point_1, 0));
1195  tmp (4) = -1.0 * weight
1196  * this->restriction[index][i + 2 * j]
1197  ((i + 2 * (j + 4)) * this->degree,
1198  dof)
1199  * this->shape_value_component
1200  ((i + 2 * (j + 4)) * this->degree,
1201  quadrature_point_2, 2);
1202  quadrature_point_3 = Point<dim> (i, j,
1203  2.0 * edge_quadrature_points[q_point] (0)
1204  - 1.0);
1205  tmp (5) = weight
1206  * (2.0 * this->shape_value_component
1207  (dof, quadrature_point_3, 2)
1208  - this->restriction[index][i + 2 * (j + 2)]
1209  ((i + 2 * (j + 4)) * this->degree,
1210  dof)
1211  * this->shape_value_component
1212  ((i + 2 * (j + 4)) * this->degree,
1213  quadrature_point_2, 2));
1214  }
1215 
1216  for (unsigned int k = 0; k < deg; ++k)
1217  {
1218  const double L_k
1219  = legendre_polynomials[k + 1].value
1220  (edge_quadrature_points[q_point] (0));
1221 
1222  for (unsigned int l = 0; l < tmp.size (); ++l)
1223  system_rhs (k, l) += tmp (l) * L_k;
1224  }
1225  }
1226 
1227  system_matrix_inv.mmult (solution, system_rhs);
1228 
1229  for (unsigned int k = 0; k < 2; ++k)
1230  for (unsigned int l = 0; l < deg; ++l)
1231  {
1232  if (std::abs (solution (l, k)) > 1e-14)
1233  this->restriction[index][i + 2 * (2 * j + k)]
1234  ((i + 4 * j) * this->degree + l + 1, dof)
1235  = solution (l, k);
1236 
1237  if (std::abs (solution (l, k + 2)) > 1e-14)
1238  this->restriction[index][2 * (i + 2 * j) + k]
1239  ((i + 4 * j + 2) * this->degree + l + 1, dof)
1240  = solution (l, k + 2);
1241 
1242  if (std::abs (solution (l, k + 4)) > 1e-14)
1243  this->restriction[index][i + 2 * (j + 2 * k)]
1244  ((i + 2 * (j + 4)) * this->degree + l + 1,
1245  dof)
1246  = solution (l, k + 4);
1247  }
1248  }
1249 
1250  const QGauss<2> face_quadrature (2 * this->degree);
1251  const std::vector<Point<2> > &face_quadrature_points
1252  = face_quadrature.get_points ();
1253  const std::vector<Polynomials::Polynomial<double> > &
1254  lobatto_polynomials
1255  = Polynomials::Lobatto::generate_complete_basis
1256  (this->degree);
1257  const unsigned int n_edge_dofs
1258  = GeometryInfo<dim>::lines_per_cell * this->degree;
1259  const unsigned int &n_face_quadrature_points
1260  = face_quadrature.size ();
1261 
1262  {
1263  FullMatrix<double> assembling_matrix
1264  (deg * this->degree,
1265  n_face_quadrature_points);
1266 
1267  for (unsigned int q_point = 0;
1268  q_point < n_face_quadrature_points; ++q_point)
1269  {
1270  const double weight
1271  = std::sqrt (face_quadrature.weight (q_point));
1272 
1273  for (unsigned int i = 0; i <= deg; ++i)
1274  {
1275  const double L_i = weight
1276  * legendre_polynomials[i].value
1277  (face_quadrature_points[q_point] (0));
1278 
1279  for (unsigned int j = 0; j < deg; ++j)
1280  assembling_matrix (i * deg + j, q_point)
1281  = L_i * lobatto_polynomials[j + 2].value
1282  (face_quadrature_points[q_point] (1));
1283  }
1284  }
1285 
1286  FullMatrix<double> system_matrix (assembling_matrix.m (),
1287  assembling_matrix.m ());
1288 
1289  assembling_matrix.mTmult (system_matrix,
1290  assembling_matrix);
1291  system_matrix_inv.reinit (system_matrix.m (),
1292  system_matrix.m ());
1293  system_matrix_inv.invert (system_matrix);
1294  }
1295 
1296  solution.reinit (system_matrix_inv.m (), 24);
1297  system_rhs.reinit (system_matrix_inv.m (), 24);
1298  tmp.reinit (24);
1299 
1300  for (unsigned int i = 0; i < 2; ++i)
1301  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
1302  {
1303  system_rhs = 0.0;
1304 
1305  for (unsigned int q_point = 0;
1306  q_point < n_face_quadrature_points; ++q_point)
1307  {
1308  tmp = 0.0;
1309 
1310  if (face_quadrature_points[q_point] (0) < 0.5)
1311  {
1312  if (face_quadrature_points[q_point] (1) < 0.5)
1313  {
1314  Point<dim> quadrature_point_0 (i,
1315  2.0 * face_quadrature_points[q_point] (0),
1316  2.0 * face_quadrature_points[q_point] (1));
1317 
1318  tmp (0) += 2.0 * this->shape_value_component
1319  (dof, quadrature_point_0, 1);
1320  tmp (1) += 2.0 * this->shape_value_component
1321  (dof, quadrature_point_0, 2);
1322  quadrature_point_0
1323  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
1324  i,
1325  2.0 * face_quadrature_points[q_point] (1));
1326  tmp (8) += 2.0 * this->shape_value_component
1327  (dof, quadrature_point_0, 2);
1328  tmp (9) += 2.0 * this->shape_value_component
1329  (dof, quadrature_point_0, 0);
1330  quadrature_point_0
1331  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
1332  2.0 * face_quadrature_points[q_point] (1),
1333  i);
1334  tmp (16) += 2.0 * this->shape_value_component
1335  (dof, quadrature_point_0, 0);
1336  tmp (17) += 2.0 * this->shape_value_component
1337  (dof, quadrature_point_0, 1);
1338  }
1339 
1340  else
1341  {
1342  Point<dim> quadrature_point_0 (i,
1343  2.0 * face_quadrature_points[q_point] (0),
1344  2.0 * face_quadrature_points[q_point] (1)
1345  - 1.0);
1346 
1347  tmp (2) += 2.0 * this->shape_value_component
1348  (dof, quadrature_point_0, 1);
1349  tmp (3) += 2.0 * this->shape_value_component
1350  (dof, quadrature_point_0, 2);
1351  quadrature_point_0
1352  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
1353  i,
1354  2.0 * face_quadrature_points[q_point] (1)
1355  - 1.0);
1356  tmp (10) += 2.0 * this->shape_value_component
1357  (dof, quadrature_point_0, 2);
1358  tmp (11) += 2.0 * this->shape_value_component
1359  (dof, quadrature_point_0, 0);
1360  quadrature_point_0
1361  = Point<dim> (2.0 * face_quadrature_points[q_point] (0),
1362  2.0 * face_quadrature_points[q_point] (1)
1363  - 1.0, i);
1364  tmp (18) += 2.0 * this->shape_value_component
1365  (dof, quadrature_point_0, 0);
1366  tmp (19) += 2.0 * this->shape_value_component
1367  (dof, quadrature_point_0, 1);
1368  }
1369  }
1370 
1371  else if (face_quadrature_points[q_point] (1) < 0.5)
1372  {
1373  Point<dim> quadrature_point_0 (i,
1374  2.0 * face_quadrature_points[q_point] (0)
1375  - 1.0,
1376  2.0 * face_quadrature_points[q_point] (1));
1377 
1378  tmp (4) += 2.0 * this->shape_value_component
1379  (dof, quadrature_point_0, 1);
1380  tmp (5) += 2.0 * this->shape_value_component
1381  (dof, quadrature_point_0, 2);
1382  quadrature_point_0
1383  = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
1384  - 1.0, i,
1385  2.0 * face_quadrature_points[q_point] (1));
1386  tmp (12) += 2.0 * this->shape_value_component
1387  (dof, quadrature_point_0, 2);
1388  tmp (13) += 2.0 * this->shape_value_component
1389  (dof, quadrature_point_0, 0);
1390  quadrature_point_0
1391  = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
1392  - 1.0,
1393  2.0 * face_quadrature_points[q_point] (1),
1394  i);
1395  tmp (20) += 2.0 * this->shape_value_component
1396  (dof, quadrature_point_0, 0);
1397  tmp (21) += 2.0 * this->shape_value_component
1398  (dof, quadrature_point_0, 1);
1399  }
1400 
1401  else
1402  {
1403  Point<dim> quadrature_point_0 (i,
1404  2.0 * face_quadrature_points[q_point] (0)
1405  - 1.0,
1406  2.0 * face_quadrature_points[q_point] (1)
1407  - 1.0);
1408 
1409  tmp (6) += 2.0 * this->shape_value_component
1410  (dof, quadrature_point_0, 1);
1411  tmp (7) += 2.0 * this->shape_value_component
1412  (dof, quadrature_point_0, 2);
1413  quadrature_point_0
1414  = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
1415  - 1.0, i,
1416  2.0 * face_quadrature_points[q_point] (1)
1417  - 1.0);
1418  tmp (14) += 2.0 * this->shape_value_component
1419  (dof, quadrature_point_0, 2);
1420  tmp (15) += 2.0 * this->shape_value_component
1421  (dof, quadrature_point_0, 0);
1422  quadrature_point_0
1423  = Point<dim> (2.0 * face_quadrature_points[q_point] (0)
1424  - 1.0,
1425  2.0 * face_quadrature_points[q_point] (1)
1426  - 1.0, i);
1427  tmp (22) += 2.0 * this->shape_value_component
1428  (dof, quadrature_point_0, 0);
1429  tmp (23) += 2.0 * this->shape_value_component
1430  (dof, quadrature_point_0, 1);
1431  }
1432 
1433  const Point<dim> quadrature_point_0 (i,
1434  face_quadrature_points[q_point] (0),
1435  face_quadrature_points[q_point] (1));
1436  const Point<dim> quadrature_point_1
1437  (face_quadrature_points[q_point] (0),
1438  i,
1439  face_quadrature_points[q_point] (1));
1440  const Point<dim> quadrature_point_2
1441  (face_quadrature_points[q_point] (0),
1442  face_quadrature_points[q_point] (1),
1443  i);
1444 
1445  for (unsigned int j = 0; j < 2; ++j)
1446  for (unsigned int k = 0; k < 2; ++k)
1447  for (unsigned int l = 0; l <= deg; ++l)
1448  {
1449  tmp (2 * (j + 2 * k))
1450  -= this->restriction[index][i + 2 * (2 * j + k)]
1451  ((i + 4 * j) * this->degree + l, dof)
1452  * this->shape_value_component
1453  ((i + 4 * j) * this->degree + l,
1454  quadrature_point_0, 1);
1455  tmp (2 * (j + 2 * k) + 1)
1456  -= this->restriction[index][i + 2 * (2 * j + k)]
1457  ((i + 2 * (k + 4)) * this->degree + l,
1458  dof)
1459  * this->shape_value_component
1460  ((i + 2 * (k + 4)) * this->degree + l,
1461  quadrature_point_0, 2);
1462  tmp (2 * (j + 2 * (k + 2)))
1463  -= this->restriction[index][2 * (i + 2 * j) + k]
1464  ((2 * (i + 4) + k) * this->degree + l,
1465  dof)
1466  * this->shape_value_component
1467  ((2 * (i + 4) + k) * this->degree + l,
1468  quadrature_point_1, 2);
1469  tmp (2 * (j + 2 * k) + 9)
1470  -= this->restriction[index][2 * (i + 2 * j) + k]
1471  ((i + 4 * j + 2) * this->degree + l,
1472  dof)
1473  * this->shape_value_component
1474  ((i + 4 * j + 2) * this->degree + l,
1475  quadrature_point_1, 0);
1476  tmp (2 * (j + 2 * (k + 4)))
1477  -= this->restriction[index][2 * (2 * i + j) + k]
1478  ((4 * i + j + 2) * this->degree + l,
1479  dof)
1480  * this->shape_value_component
1481  ((4 * i + j + 2) * this->degree + l,
1482  quadrature_point_2, 0);
1483  tmp (2 * (j + 2 * k) + 17)
1484  -= this->restriction[index][2 * (2 * i + j) + k]
1485  ((4 * i + k) * this->degree + l, dof)
1486  * this->shape_value_component
1487  ((4 * i + k) * this->degree + l,
1488  quadrature_point_2, 1);
1489  }
1490 
1491  tmp *= face_quadrature.weight (q_point);
1492 
1493  for (unsigned int j = 0; j <= deg; ++j)
1494  {
1495  const double L_j_0
1496  = legendre_polynomials[j].value
1497  (face_quadrature_points[q_point] (0));
1498  const double L_j_1
1499  = legendre_polynomials[j].value
1500  (face_quadrature_points[q_point] (1));
1501 
1502  for (unsigned int k = 0; k < deg; ++k)
1503  {
1504  const double l_k_0
1505  = L_j_0 * lobatto_polynomials[k + 2].value
1506  (face_quadrature_points[q_point] (1));
1507  const double l_k_1
1508  = L_j_1 * lobatto_polynomials[k + 2].value
1509  (face_quadrature_points[q_point] (0));
1510 
1511  for (unsigned int l = 0; l < 4; ++l)
1512  {
1513  system_rhs (j * deg + k, 2 * l)
1514  += tmp (2 * l) * l_k_0;
1515  system_rhs (j * deg + k, 2 * l + 1)
1516  += tmp (2 * l + 1) * l_k_1;
1517  system_rhs (j * deg + k, 2 * (l + 4))
1518  += tmp (2 * (l + 4)) * l_k_1;
1519  system_rhs (j * deg + k, 2 * l + 9)
1520  += tmp (2 * l + 9) * l_k_0;
1521  system_rhs (j * deg + k, 2 * (l + 8))
1522  += tmp (2 * (l + 8)) * l_k_0;
1523  system_rhs (j * deg + k, 2 * l + 17)
1524  += tmp (2 * l + 17) * l_k_1;
1525  }
1526  }
1527  }
1528  }
1529 
1530  system_matrix_inv.mmult (solution, system_rhs);
1531 
1532  for (unsigned int j = 0; j < 2; ++j)
1533  for (unsigned int k = 0; k < 2; ++k)
1534  for (unsigned int l = 0; l <= deg; ++l)
1535  for (unsigned int m = 0; m < deg; ++m)
1536  {
1537  if (std::abs (solution (l * deg + m,
1538  2 * (j + 2 * k)))
1539  > 1e-14)
1540  this->restriction[index][i + 2 * (2 * j + k)]
1541  ((2 * i * this->degree + l) * deg + m
1542  + n_edge_dofs,
1543  dof) = solution (l * deg + m,
1544  2 * (j + 2 * k));
1545 
1546  if (std::abs (solution (l * deg + m,
1547  2 * (j + 2 * k) + 1))
1548  > 1e-14)
1549  this->restriction[index][i + 2 * (2 * j + k)]
1550  (((2 * i + 1) * deg + m) * this->degree + l
1551  + n_edge_dofs, dof)
1552  = solution (l * deg + m,
1553  2 * (j + 2 * k) + 1);
1554 
1555  if (std::abs (solution (l * deg + m,
1556  2 * (j + 2 * (k + 2))))
1557  > 1e-14)
1558  this->restriction[index][2 * (i + 2 * j) + k]
1559  ((2 * (i + 2) * this->degree + l) * deg + m
1560  + n_edge_dofs,
1561  dof) = solution (l * deg + m,
1562  2 * (j + 2 * (k + 2)));
1563 
1564  if (std::abs (solution (l * deg + m,
1565  2 * (j + 2 * k) + 9))
1566  > 1e-14)
1567  this->restriction[index][2 * (i + 2 * j) + k]
1568  (((2 * i + 5) * deg + m) * this->degree + l
1569  + n_edge_dofs, dof)
1570  = solution (l * deg + m,
1571  2 * (j + 2 * k) + 9);
1572 
1573  if (std::abs (solution (l * deg + m,
1574  2 * (j + 2 * (k + 4))))
1575  > 1e-14)
1576  this->restriction[index][2 * (2 * i + j) + k]
1577  ((2 * (i + 4) * this->degree + l) * deg + m
1578  + n_edge_dofs,
1579  dof) = solution (l * deg + m,
1580  2 * (j + 2 * (k + 4)));
1581 
1582  if (std::abs (solution (l * deg + m,
1583  2 * (j + 2 * k) + 17))
1584  > 1e-14)
1585  this->restriction[index][2 * (2 * i + j) + k]
1586  (((2 * i + 9) * deg + m) * this->degree + l
1587  + n_edge_dofs, dof)
1588  = solution (l * deg + m,
1589  2 * (j + 2 * k) + 17);
1590  }
1591  }
1592 
1593  const QGauss<dim> quadrature (2 * this->degree);
1594  const std::vector<Point<dim> > &
1595  quadrature_points = quadrature.get_points ();
1596  const unsigned int n_boundary_dofs
1597  = 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree
1598  + n_edge_dofs;
1599  const unsigned int &n_quadrature_points = quadrature.size ();
1600 
1601  {
1602  FullMatrix<double>
1603  assembling_matrix (deg * deg * this->degree,
1604  n_quadrature_points);
1605 
1606  for (unsigned int q_point = 0; q_point < n_quadrature_points;
1607  ++q_point)
1608  {
1609  const double weight = std::sqrt (quadrature.weight
1610  (q_point));
1611 
1612  for (unsigned int i = 0; i <= deg; ++i)
1613  {
1614  const double L_i = weight
1615  * legendre_polynomials[i].value
1616  (quadrature_points[q_point] (0));
1617 
1618  for (unsigned int j = 0; j < deg; ++j)
1619  {
1620  const double l_j
1621  = L_i * lobatto_polynomials[j + 2].value
1622  (quadrature_points[q_point] (1));
1623 
1624  for (unsigned int k = 0; k < deg; ++k)
1625  assembling_matrix ((i * deg + j) * deg + k,
1626  q_point)
1627  = l_j * lobatto_polynomials[k + 2].value
1628  (quadrature_points[q_point] (2));
1629  }
1630  }
1631  }
1632 
1633  FullMatrix<double> system_matrix (assembling_matrix.m (),
1634  assembling_matrix.m ());
1635 
1636  assembling_matrix.mTmult (system_matrix, assembling_matrix);
1637  system_matrix_inv.reinit (system_matrix.m (),
1638  system_matrix.m ());
1639  system_matrix_inv.invert (system_matrix);
1640  }
1641 
1642  solution.reinit (system_matrix_inv.m (), 24);
1643  system_rhs.reinit (system_matrix_inv.m (), 24);
1644  tmp.reinit (24);
1645 
1646  for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
1647  {
1648  system_rhs = 0.0;
1649 
1650  for (unsigned int q_point = 0;
1651  q_point < n_quadrature_points; ++q_point)
1652  {
1653  tmp = 0.0;
1654 
1655  if (quadrature_points[q_point] (0) < 0.5)
1656  {
1657  if (quadrature_points[q_point] (1) < 0.5)
1658  {
1659  if (quadrature_points[q_point] (2) < 0.5)
1660  {
1661  const Point<dim> quadrature_point
1662  (2.0 * quadrature_points[q_point] (0),
1663  2.0 * quadrature_points[q_point] (1),
1664  2.0 * quadrature_points[q_point] (2));
1665 
1666  tmp (0) += 2.0 * this->shape_value_component
1667  (dof, quadrature_point, 0);
1668  tmp (1) += 2.0 * this->shape_value_component
1669  (dof, quadrature_point, 1);
1670  tmp (2) += 2.0 * this->shape_value_component
1671  (dof, quadrature_point, 2);
1672  }
1673 
1674  else
1675  {
1676  const Point<dim> quadrature_point
1677  (2.0 * quadrature_points[q_point] (0),
1678  2.0 * quadrature_points[q_point] (1),
1679  2.0 * quadrature_points[q_point] (2)
1680  - 1.0);
1681 
1682  tmp (3) += 2.0 * this->shape_value_component
1683  (dof, quadrature_point, 0);
1684  tmp (4) += 2.0 * this->shape_value_component
1685  (dof, quadrature_point, 1);
1686  tmp (5) += 2.0 * this->shape_value_component
1687  (dof, quadrature_point, 2);
1688  }
1689  }
1690 
1691  else if (quadrature_points[q_point] (2) < 0.5)
1692  {
1693  const Point<dim> quadrature_point
1694  (2.0 * quadrature_points[q_point] (0),
1695  2.0 * quadrature_points[q_point] (1)
1696  - 1.0,
1697  2.0 * quadrature_points[q_point] (2));
1698 
1699  tmp (6) += 2.0 * this->shape_value_component
1700  (dof, quadrature_point, 0);
1701  tmp (7) += 2.0 * this->shape_value_component
1702  (dof, quadrature_point, 1);
1703  tmp (8) += 2.0 * this->shape_value_component
1704  (dof, quadrature_point, 2);
1705  }
1706 
1707  else
1708  {
1709  const Point<dim> quadrature_point
1710  (2.0 * quadrature_points[q_point] (0),
1711  2.0 * quadrature_points[q_point] (1)
1712  - 1.0,
1713  2.0 * quadrature_points[q_point] (2)
1714  - 1.0);
1715 
1716  tmp (9) += 2.0 * this->shape_value_component
1717  (dof, quadrature_point, 0);
1718  tmp (10) += 2.0 * this->shape_value_component
1719  (dof, quadrature_point, 1);
1720  tmp (11) += 2.0 * this->shape_value_component
1721  (dof, quadrature_point, 2);
1722  }
1723  }
1724 
1725  else if (quadrature_points[q_point] (1) < 0.5)
1726  {
1727  if (quadrature_points[q_point] (2) < 0.5)
1728  {
1729  const Point<dim> quadrature_point
1730  (2.0 * quadrature_points[q_point] (0)
1731  - 1.0,
1732  2.0 * quadrature_points[q_point] (1),
1733  2.0 * quadrature_points[q_point] (2));
1734 
1735  tmp (12) += 2.0 * this->shape_value_component
1736  (dof, quadrature_point, 0);
1737  tmp (13) += 2.0 * this->shape_value_component
1738  (dof, quadrature_point, 1);
1739  tmp (14) += 2.0 * this->shape_value_component
1740  (dof, quadrature_point, 2);
1741  }
1742 
1743  else
1744  {
1745  const Point<dim> quadrature_point
1746  (2.0 * quadrature_points[q_point] (0)
1747  - 1.0,
1748  2.0 * quadrature_points[q_point] (1),
1749  2.0 * quadrature_points[q_point] (2)
1750  - 1.0);
1751 
1752  tmp (15) += 2.0 * this->shape_value_component
1753  (dof, quadrature_point, 0);
1754  tmp (16) += 2.0 * this->shape_value_component
1755  (dof, quadrature_point, 1);
1756  tmp (17) += 2.0 * this->shape_value_component
1757  (dof, quadrature_point, 2);
1758  }
1759  }
1760 
1761  else if (quadrature_points[q_point] (2) < 0.5)
1762  {
1763  const Point<dim> quadrature_point
1764  (2.0 * quadrature_points[q_point] (0)
1765  - 1.0,
1766  2.0 * quadrature_points[q_point] (1)
1767  - 1.0,
1768  2.0 * quadrature_points[q_point] (2));
1769 
1770  tmp (18) += 2.0 * this->shape_value_component
1771  (dof, quadrature_point, 0);
1772  tmp (19) += 2.0 * this->shape_value_component
1773  (dof, quadrature_point, 1);
1774  tmp (20) += 2.0 * this->shape_value_component
1775  (dof, quadrature_point, 2);
1776  }
1777 
1778  else
1779  {
1780  const Point<dim> quadrature_point
1781  (2.0 * quadrature_points[q_point] (0)
1782  - 1.0,
1783  2.0 * quadrature_points[q_point] (1)
1784  - 1.0,
1785  2.0 * quadrature_points[q_point] (2)
1786  - 1.0);
1787 
1788  tmp (21) += 2.0 * this->shape_value_component
1789  (dof, quadrature_point, 0);
1790  tmp (22) += 2.0 * this->shape_value_component
1791  (dof, quadrature_point, 1);
1792  tmp (23) += 2.0 * this->shape_value_component
1793  (dof, quadrature_point, 2);
1794  }
1795 
1796  for (unsigned int i = 0; i < 2; ++i)
1797  for (unsigned int j = 0; j < 2; ++j)
1798  for (unsigned int k = 0; k < 2; ++k)
1799  for (unsigned int l = 0; l <= deg; ++l)
1800  {
1801  tmp (3 * (i + 2 * (j + 2 * k)))
1802  -= this->restriction[index][2 * (2 * i + j) + k]
1803  ((4 * i + j + 2) * this->degree + l, dof)
1804  * this->shape_value_component
1805  ((4 * i + j + 2) * this->degree + l,
1806  quadrature_points[q_point], 0);
1807  tmp (3 * (i + 2 * (j + 2 * k)) + 1)
1808  -= this->restriction[index][2 * (2 * i + j) + k]
1809  ((4 * i + k) * this->degree + l, dof)
1810  * this->shape_value_component
1811  ((4 * i + k) * this->degree + l,
1812  quadrature_points[q_point], 1);
1813  tmp (3 * (i + 2 * (j + 2 * k)) + 2)
1814  -= this->restriction[index][2 * (2 * i + j) + k]
1815  ((2 * (j + 4) + k) * this->degree + l,
1816  dof)
1817  * this->shape_value_component
1818  ((2 * (j + 4) + k) * this->degree + l,
1819  quadrature_points[q_point], 2);
1820 
1821  for (unsigned int m = 0; m < deg; ++m)
1822  {
1823  tmp (3 * (i + 2 * (j + 2 * k)))
1824  -= this->restriction[index][2 * (2 * i + j) + k]
1825  (((2 * j + 5) * deg + m)
1826  * this->degree + l + n_edge_dofs,
1827  dof)
1828  * this->shape_value_component
1829  (((2 * j + 5) * deg + m)
1830  * this->degree + l + n_edge_dofs,
1831  quadrature_points[q_point], 0);
1832  tmp (3 * (i + 2 * (j + 2 * k)))
1833  -= this->restriction[index][2 * (2 * i + j) + k]
1834  ((2 * (i + 4) * this->degree + l)
1835  * deg + m + n_edge_dofs, dof)
1836  * this->shape_value_component
1837  ((2 * (i + 4) * this->degree + l)
1838  * deg + m + n_edge_dofs,
1839  quadrature_points[q_point], 0);
1840  tmp (3 * (i + 2 * (j + 2 * k)) + 1)
1841  -= this->restriction[index][2 * (2 * i + j) + k]
1842  ((2 * k * this->degree + l) * deg + m
1843  + n_edge_dofs,
1844  dof)
1845  * this->shape_value_component
1846  ((2 * k * this->degree + l) * deg + m
1847  + n_edge_dofs,
1848  quadrature_points[q_point], 1);
1849  tmp (3 * (i + 2 * (j + 2 * k)) + 1)
1850  -= this->restriction[index][2 * (2 * i + j) + k]
1851  (((2 * i + 9) * deg + m)
1852  * this->degree + l + n_edge_dofs,
1853  dof)
1854  * this->shape_value_component
1855  (((2 * i + 9) * deg + m)
1856  * this->degree + l + n_edge_dofs,
1857  quadrature_points[q_point], 1);
1858  tmp (3 * (i + 2 * (j + 2 * k)) + 2)
1859  -= this->restriction[index][2 * (2 * i + j) + k]
1860  (((2 * k + 1) * deg + m)
1861  * this->degree + l + n_edge_dofs,
1862  dof)
1863  * this->shape_value_component
1864  (((2 * k + 1) * deg + m)
1865  * this->degree + l + n_edge_dofs,
1866  quadrature_points[q_point], 2);
1867  tmp (3 * (i + 2 * (j + 2 * k)) + 2)
1868  -= this->restriction[index][2 * (2 * i + j) + k]
1869  ((2 * (j + 2) * this->degree + l)
1870  * deg + m + n_edge_dofs, dof)
1871  * this->shape_value_component
1872  ((2 * (j + 2) * this->degree + l)
1873  * deg + m + n_edge_dofs,
1874  quadrature_points[q_point], 2);
1875  }
1876  }
1877 
1878  tmp *= quadrature.weight (q_point);
1879 
1880  for (unsigned int i = 0; i <= deg; ++i)
1881  {
1882  const double L_i_0
1883  = legendre_polynomials[i].value
1884  (quadrature_points[q_point] (0));
1885  const double L_i_1
1886  = legendre_polynomials[i].value
1887  (quadrature_points[q_point] (1));
1888  const double L_i_2
1889  = legendre_polynomials[i].value
1890  (quadrature_points[q_point] (2));
1891 
1892  for (unsigned int j = 0; j < deg; ++j)
1893  {
1894  const double l_j_0
1895  = L_i_0 * lobatto_polynomials[j + 2].value
1896  (quadrature_points[q_point] (1));
1897  const double l_j_1
1898  = L_i_1 * lobatto_polynomials[j + 2].value
1899  (quadrature_points[q_point] (0));
1900  const double l_j_2
1901  = L_i_2 * lobatto_polynomials[j + 2].value
1902  (quadrature_points[q_point] (0));
1903 
1904  for (unsigned int k = 0; k < deg; ++k)
1905  {
1906  const double l_k_0
1907  = l_j_0 * lobatto_polynomials[k + 2].value
1908  (quadrature_points[q_point] (2));
1909  const double l_k_1
1910  = l_j_1 * lobatto_polynomials[k + 2].value
1911  (quadrature_points[q_point] (2));
1912  const double l_k_2
1913  = l_j_2 * lobatto_polynomials[k + 2].value
1914  (quadrature_points[q_point] (1));
1915 
1916  for (unsigned int l = 0; l < 8; ++l)
1917  {
1918  system_rhs ((i * deg + j) * deg + k,
1919  3 * l)
1920  += tmp (3 * l) * l_k_0;
1921  system_rhs ((i * deg + j) * deg + k,
1922  3 * l + 1)
1923  += tmp (3 * l + 1) * l_k_1;
1924  system_rhs ((i * deg + j) * deg + k,
1925  3 * l + 2)
1926  += tmp (3 * l + 2) * l_k_2;
1927  }
1928  }
1929  }
1930  }
1931  }
1932 
1933  system_matrix_inv.mmult (solution, system_rhs);
1934 
1935  for (unsigned int i = 0; i < 2; ++i)
1936  for (unsigned int j = 0; j < 2; ++j)
1937  for (unsigned int k = 0; k < 2; ++k)
1938  for (unsigned int l = 0; l <= deg; ++l)
1939  for (unsigned int m = 0; m < deg; ++m)
1940  for (unsigned int n = 0; n < deg; ++n)
1941  {
1942  if (std::abs (solution
1943  ((l * deg + m) * deg + n,
1944  3 * (i + 2 * (j + 2 * k))))
1945  > 1e-14)
1946  this->restriction[index][2 * (2 * i + j) + k]
1947  ((l * deg + m) * deg + n + n_boundary_dofs,
1948  dof) = solution ((l * deg + m) * deg + n,
1949  3 * (i + 2 * (j + 2 * k)));
1950 
1951  if (std::abs (solution
1952  ((l * deg + m) * deg + n,
1953  3 * (i + 2 * (j + 2 * k)) + 1))
1954  > 1e-14)
1955  this->restriction[index][2 * (2 * i + j) + k]
1956  ((l + (m + deg) * this->degree) * deg + n
1957  + n_boundary_dofs,
1958  dof) = solution ((l * deg + m) * deg + n,
1959  3 * (i + 2 * (j + 2 * k)) + 1);
1960 
1961  if (std::abs (solution
1962  ((l * deg + m) * deg + n,
1963  3 * (i + 2 * (j + 2 * k)) + 2))
1964  > 1e-14)
1965  this->restriction[index][2 * (2 * i + j) + k]
1966  (l + ((m + 2 * deg) * deg + n) * this->degree
1967  + n_boundary_dofs, dof)
1968  = solution ((l * deg + m) * deg + n,
1969  3 * (i + 2 * (j + 2 * k)) + 2);
1970  }
1971  }
1972  }
1973 
1974  break;
1975  }
1976 
1977  default:
1978  Assert (false, ExcNotImplemented ());
1979  }
1980 }
1981 
1982 
1983 
1984 template <int dim>
1985 std::vector<unsigned int>
1986 FE_Nedelec<dim>::get_dpo_vector (const unsigned int degree, bool dg)
1987 {
1988  std::vector<unsigned int> dpo (dim + 1);
1989 
1990  if (dg)
1991  {
1992  dpo[dim] = PolynomialsNedelec<dim>::compute_n_pols(degree);
1993  }
1994  else
1995  {
1996  dpo[0] = 0;
1997  dpo[1] = degree + 1;
1998  dpo[2] = 2 * degree * (degree + 1);
1999 
2000  if (dim == 3)
2001  dpo[3] = 3 * degree * degree * (degree + 1);
2002  }
2003 
2004  return dpo;
2005 }
2006 
2007 //---------------------------------------------------------------------------
2008 // Data field initialization
2009 //---------------------------------------------------------------------------
2010 
2011 // Chech wheter a given shape
2012 // function has support on a
2013 // given face.
2014 
2015 // We just switch through the
2016 // faces of the cell and return
2017 // true, if the shape function
2018 // has support on the face
2019 // and false otherwise.
2020 template <int dim>
2021 bool
2022 FE_Nedelec<dim>::has_support_on_face (const unsigned int shape_index,
2023  const unsigned int face_index) const
2024 {
2025  Assert (shape_index < this->dofs_per_cell,
2026  ExcIndexRange (shape_index, 0, this->dofs_per_cell));
2027  Assert (face_index < GeometryInfo<dim>::faces_per_cell,
2028  ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
2029 
2030  const unsigned int deg = this->degree-1;
2031  switch (dim)
2032  {
2033  case 2:
2034  switch (face_index)
2035  {
2036  case 0:
2037  if (!((shape_index > deg) && (shape_index < 2 * this->degree)))
2038  return true;
2039 
2040  else
2041  return false;
2042 
2043  case 1:
2044  if ((shape_index > deg) &&
2045  (shape_index
2046  < GeometryInfo<2>::lines_per_cell * this->degree))
2047  return true;
2048 
2049  else
2050  return false;
2051 
2052  case 2:
2053  if (shape_index < 3 * this->degree)
2054  return true;
2055 
2056  else
2057  return false;
2058 
2059  case 3:
2060  if (!((shape_index >= 2 * this->degree) &&
2061  (shape_index < 3 * this->degree)))
2062  return true;
2063 
2064  else
2065  return false;
2066 
2067  default:
2068  {
2069  Assert (false, ExcNotImplemented ());
2070  return false;
2071  }
2072  }
2073 
2074  case 3:
2075  switch (face_index)
2076  {
2077  case 0:
2078  if (((shape_index > deg) && (shape_index < 2 * this->degree)) ||
2079  ((shape_index >= 5 * this->degree) &&
2080  (shape_index < 6 * this->degree)) ||
2081  ((shape_index >= 9 * this->degree) &&
2082  (shape_index < 10 * this->degree)) ||
2083  ((shape_index >= 11 * this->degree) &&
2084  (shape_index
2085  < GeometryInfo<3>::lines_per_cell * this->degree)) ||
2086  ((shape_index
2087  >= (GeometryInfo<3>::lines_per_cell + 2 * deg)
2088  * this->degree) &&
2089  (shape_index
2090  < (GeometryInfo<3>::lines_per_cell + 4 * deg)
2091  * this->degree)) ||
2092  ((shape_index
2093  >= (GeometryInfo<3>::lines_per_cell + 5 * deg)
2094  * this->degree) &&
2095  (shape_index
2096  < (GeometryInfo<3>::lines_per_cell + 6 * deg)
2097  * this->degree)) ||
2098  ((shape_index
2099  >= (GeometryInfo<3>::lines_per_cell + 7 * deg)
2100  * this->degree) &&
2101  (shape_index
2102  < (GeometryInfo<3>::lines_per_cell + 9 * deg)
2103  * this->degree)) ||
2104  ((shape_index
2105  >= (GeometryInfo<3>::lines_per_cell + 10 * deg)
2106  * this->degree) &&
2107  (shape_index
2108  < (GeometryInfo<3>::lines_per_cell + 11 * deg)
2109  * this->degree)))
2110  return false;
2111 
2112  else
2113  return true;
2114 
2115  case 1:
2116  if (((shape_index > deg) && (shape_index < 4 * this->degree)) ||
2117  ((shape_index >= 5 * this->degree) &&
2118  (shape_index < 8 * this->degree)) ||
2119  ((shape_index >= 9 * this->degree) &&
2120  (shape_index < 10 * this->degree)) ||
2121  ((shape_index >= 11 * this->degree) &&
2122  (shape_index
2123  < GeometryInfo<3>::lines_per_cell * this->degree)) ||
2124  ((shape_index
2125  >= (GeometryInfo<3>::lines_per_cell + 2 * deg)
2126  * this->degree) &&
2127  (shape_index
2128  < (GeometryInfo<3>::lines_per_cell + 5 * deg)
2129  * this->degree)) ||
2130  ((shape_index
2131  >= (GeometryInfo<3>::lines_per_cell + 6 * deg)
2132  * this->degree) &&
2133  (shape_index
2134  < (GeometryInfo<3>::lines_per_cell + 7 * deg)
2135  * this->degree)) ||
2136  ((shape_index
2137  >= (GeometryInfo<3>::lines_per_cell + 9 * deg)
2138  * this->degree) &&
2139  (shape_index
2140  < (GeometryInfo<3>::lines_per_cell + 10 * deg)
2141  * this->degree)) ||
2142  ((shape_index
2143  >= (GeometryInfo<3>::lines_per_cell + 11 * deg)
2144  * this->degree) &&
2145  (shape_index
2146  < (GeometryInfo<3>::lines_per_cell + 12 * deg)
2147  * this->degree)))
2148  return true;
2149 
2150  else
2151  return false;
2152 
2153  case 2:
2154  if ((shape_index < 3 * this->degree) ||
2155  ((shape_index >= 4 * this->degree) &&
2156  (shape_index < 7 * this->degree)) ||
2157  ((shape_index >= 8 * this->degree) &&
2158  (shape_index < 10 * this->degree)) ||
2159  ((shape_index
2160  >= (GeometryInfo<3>::lines_per_cell + deg)
2161  * this->degree) &&
2162  (shape_index
2163  < (GeometryInfo<3>::lines_per_cell + 2 * deg)
2164  * this->degree)) ||
2165  ((shape_index
2166  >= (GeometryInfo<3>::lines_per_cell + 3 * deg)
2167  * this->degree) &&
2168  (shape_index
2169  < (GeometryInfo<3>::lines_per_cell + 6 * deg)
2170  * this->degree)) ||
2171  ((shape_index
2172  >= (GeometryInfo<3>::lines_per_cell + 8 * deg)
2173  * this->degree) &&
2174  (shape_index
2175  < (GeometryInfo<3>::lines_per_cell + 9 * deg)
2176  * this->degree)) ||
2177  ((shape_index
2178  >= (GeometryInfo<3>::lines_per_cell + 10 * deg)
2179  * this->degree) &&
2180  (shape_index
2181  < (GeometryInfo<3>::lines_per_cell + 11 * deg)
2182  * this->degree)))
2183  return true;
2184 
2185  else
2186  return false;
2187 
2188  case 3:
2189  if ((shape_index < 2 * this->degree) ||
2190  ((shape_index >= 3 * this->degree) &&
2191  (shape_index < 6 * this->degree)) ||
2192  ((shape_index >= 7 * this->degree) &&
2193  (shape_index < 8 * this->degree)) ||
2194  ((shape_index >= 10 * this->degree) &&
2195  (shape_index
2196  < GeometryInfo<3>::lines_per_cell * this->degree)) ||
2197  ((shape_index
2198  >= (GeometryInfo<3>::lines_per_cell + deg)
2199  * this->degree) &&
2200  (shape_index
2201  < (GeometryInfo<3>::lines_per_cell + 2 * deg)
2202  * this->degree)) ||
2203  ((shape_index
2204  >= (GeometryInfo<3>::lines_per_cell + 3 * deg)
2205  * this->degree) &&
2206  (shape_index
2207  < (GeometryInfo<3>::lines_per_cell + 4 * deg)
2208  * this->degree)) ||
2209  ((shape_index
2210  >= (GeometryInfo<3>::lines_per_cell + 6 * deg)
2211  * this->degree) &&
2212  (shape_index
2213  < (GeometryInfo<3>::lines_per_cell + 9 * deg)
2214  * this->degree)) ||
2215  ((shape_index
2216  >= (GeometryInfo<3>::lines_per_cell + 10 * deg)
2217  * this->degree) &&
2218  (shape_index
2219  < (GeometryInfo<3>::lines_per_cell + 11 * deg)
2220  * this->degree)))
2221  return true;
2222 
2223  else
2224  return false;
2225 
2226  case 4:
2227  if ((shape_index < 4 * this->degree) ||
2228  ((shape_index >= 8 * this->degree) &&
2229  (shape_index
2230  < (GeometryInfo<3>::lines_per_cell + deg)
2231  * this->degree)) ||
2232  ((shape_index
2233  >= (GeometryInfo<3>::lines_per_cell + 2 * deg)
2234  * this->degree) &&
2235  (shape_index
2236  < (GeometryInfo<3>::lines_per_cell + 3 * deg)
2237  * this->degree)) ||
2238  ((shape_index
2239  >= (GeometryInfo<3>::lines_per_cell + 5 * deg)
2240  * this->degree) &&
2241  (shape_index
2242  < (GeometryInfo<3>::lines_per_cell + 6 * deg)
2243  * this->degree)) ||
2244  ((shape_index
2245  >= (GeometryInfo<3>::lines_per_cell + 7 * deg)
2246  * this->degree) &&
2247  (shape_index
2248  < (GeometryInfo<3>::lines_per_cell + 10 * deg)
2249  * this->degree)))
2250  return true;
2251 
2252  else
2253  return false;
2254 
2255  case 5:
2256  if (((shape_index >= 4 * this->degree) &&
2257  (shape_index
2258  < (GeometryInfo<3>::lines_per_cell + deg)
2259  * this->degree)) ||
2260  ((shape_index
2261  >= (GeometryInfo<3>::lines_per_cell + 2 * deg)
2262  * this->degree) &&
2263  (shape_index
2264  < (GeometryInfo<3>::lines_per_cell + 3 * deg)
2265  * this->degree)) ||
2266  ((shape_index
2267  >= (GeometryInfo<3>::lines_per_cell + 5 * deg)
2268  * this->degree) &&
2269  (shape_index
2270  < (GeometryInfo<3>::lines_per_cell + 6 * deg)
2271  * this->degree)) ||
2272  ((shape_index
2273  >= (GeometryInfo<3>::lines_per_cell + 7 * deg)
2274  * this->degree) &&
2275  (shape_index
2276  < (GeometryInfo<3>::lines_per_cell + 8 * deg)
2277  * this->degree)) ||
2278  ((shape_index
2279  >= (GeometryInfo<3>::lines_per_cell + 10 * deg)
2280  * this->degree) &&
2281  (shape_index
2282  < (GeometryInfo<3>::lines_per_cell + 12 * deg)
2283  * this->degree)))
2284  return true;
2285 
2286  else
2287  return false;
2288 
2289  default:
2290  {
2291  Assert (false, ExcNotImplemented ());
2292  return false;
2293  }
2294  }
2295 
2296  default:
2297  {
2298  Assert (false, ExcNotImplemented ());
2299  return false;
2300  }
2301  }
2302 }
2303 
2304 template <int dim>
2305 FiniteElementDomination::Domination
2306 FE_Nedelec<dim>::compare_for_face_domination (const FiniteElement<dim> &fe_other) const
2307 {
2308  if (const FE_Nedelec<dim> *fe_nedelec_other
2309  = dynamic_cast<const FE_Nedelec<dim>*>(&fe_other))
2310  {
2311  if (this->degree < fe_nedelec_other->degree)
2312  return FiniteElementDomination::this_element_dominates;
2313  else if (this->degree == fe_nedelec_other->degree)
2314  return FiniteElementDomination::either_element_can_dominate;
2315  else
2316  return FiniteElementDomination::other_element_dominates;
2317  }
2318  else if (const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(&fe_other))
2319  {
2320  // TODO: ???
2321  // the FE_Nothing has no
2322  // degrees of
2323  // freedom. nevertheless, we
2324  // say that the FE_Q element
2325  // dominates so that we don't
2326  // have to force the FE_Q side
2327  // to become a zero function
2328  // and rather allow the
2329  // function to be discontinuous
2330  // along the interface
2331 // return FiniteElementDomination::other_element_dominates;
2332  if (fe_nothing->is_dominating())
2333  {
2334  return FiniteElementDomination::other_element_dominates;
2335  }
2336  else
2337  {
2338  // the FE_Nothing has no degrees of freedom and it is typically used in
2339  // a context where we don't require any continuity along the interface
2340  return FiniteElementDomination::no_requirements;
2341  }
2342  }
2343 
2344  Assert (false, ExcNotImplemented());
2345  return FiniteElementDomination::neither_element_dominates;
2346 }
2347 
2348 template <int dim>
2349 bool
2350 FE_Nedelec<dim>::hp_constraints_are_implemented () const
2351 {
2352  return true;
2353 }
2354 
2355 template <int dim>
2356 std::vector<std::pair<unsigned int, unsigned int> >
2357 FE_Nedelec<dim>::hp_vertex_dof_identities (const FiniteElement<dim> &)
2358 const
2359 {
2360  // Nedelec elements do not have any dofs
2361  // on vertices, hence return an empty vector.
2362  return std::vector<std::pair<unsigned int, unsigned int> > ();
2363 }
2364 
2365 template <int dim>
2366 std::vector<std::pair<unsigned int, unsigned int> >
2367 FE_Nedelec<dim>::hp_line_dof_identities (const FiniteElement<dim> &fe_other)
2368 const
2369 {
2370  // we can presently only compute these
2371  // identities if both FEs are
2372  // FE_Nedelec or if the other one is an
2373  // FE_Nothing
2374  if (const FE_Nedelec<dim> *fe_nedelec_other
2375  = dynamic_cast<const FE_Nedelec<dim>*> (&fe_other))
2376  {
2377  // dofs are located on lines, so
2378  // two dofs are identical, if their
2379  // edge shape functions have the
2380  // same polynomial degree.
2381  std::vector<std::pair<unsigned int, unsigned int> > identities;
2382 
2383  for (unsigned int i = 0;
2384  i < std::min (fe_nedelec_other->degree, this->degree); ++i)
2385  identities.push_back (std::make_pair (i, i));
2386 
2387  return identities;
2388  }
2389 
2390  else if (dynamic_cast<const FE_Nothing<dim>*> (&fe_other) != 0)
2391  {
2392  // the FE_Nothing has no
2393  // degrees of freedom, so there
2394  // are no equivalencies to be
2395  // recorded
2396  return std::vector<std::pair<unsigned int, unsigned int> > ();
2397  }
2398 
2399  else
2400  {
2401  Assert (false, ExcNotImplemented ());
2402  return std::vector<std::pair<unsigned int, unsigned int> > ();
2403  }
2404 }
2405 
2406 template <int dim>
2407 std::vector<std::pair<unsigned int, unsigned int> >
2408 FE_Nedelec<dim>::hp_quad_dof_identities (const FiniteElement<dim> &fe_other)
2409 const
2410 {
2411  // we can presently only compute
2412  // these identities if both FEs are
2413  // FE_Nedelec or if the other one is an
2414  // FE_Nothing
2415  if (const FE_Nedelec<dim> *fe_nedelec_other
2416  = dynamic_cast<const FE_Nedelec<dim>*> (&fe_other))
2417  {
2418  // dofs are located on the interior
2419  // of faces, so two dofs are identical,
2420  // if their face shape functions have
2421  // the same polynomial degree.
2422  const unsigned int p = fe_nedelec_other->degree;
2423  const unsigned int q = this->degree;
2424  const unsigned int p_min = std::min (p, q);
2425  std::vector<std::pair<unsigned int, unsigned int> > identities;
2426 
2427  for (unsigned int i = 0; i < p_min; ++i)
2428  for (unsigned int j = 0; j < p_min - 1; ++j)
2429  {
2430  identities.push_back (std::make_pair (i * (q - 1) + j,
2431  i * (p - 1) + j));
2432  identities.push_back (std::make_pair (i + (j + q - 1) * q,
2433  i + (j + p - 1) * p));
2434  }
2435 
2436  return identities;
2437  }
2438 
2439  else if (dynamic_cast<const FE_Nothing<dim>*> (&fe_other) != 0)
2440  {
2441  // the FE_Nothing has no
2442  // degrees of freedom, so there
2443  // are no equivalencies to be
2444  // recorded
2445  return std::vector<std::pair<unsigned int, unsigned int> > ();
2446  }
2447 
2448  else
2449  {
2450  Assert (false, ExcNotImplemented ());
2451  return std::vector<std::pair<unsigned int, unsigned int> > ();
2452  }
2453 }
2454 
2455 // In this function we compute the face
2456 // interpolation matrix. This is usually
2457 // done by projection-based interpolation,
2458 // but, since one can compute the entries
2459 // easy per hand, we save some computation
2460 // time at this point and just fill in the
2461 // correct values.
2462 template <int dim>
2463 void
2464 FE_Nedelec<dim>::get_face_interpolation_matrix
2465 (const FiniteElement<dim> &source, FullMatrix<double> &interpolation_matrix)
2466 const
2467 {
2468  // this is only implemented, if the
2469  // source FE is also a
2470  // Nedelec element
2471  typedef FE_Nedelec<dim> FEN;
2472  typedef FiniteElement<dim> FEL;
2473 
2474  AssertThrow ((source.get_name ().find ("FE_Nedelec<") == 0) ||
2475  (dynamic_cast<const FEN *> (&source) != 0),
2476  typename FEL::ExcInterpolationNotImplemented());
2477  Assert (interpolation_matrix.m () == source.dofs_per_face,
2478  ExcDimensionMismatch (interpolation_matrix.m (),
2479  source.dofs_per_face));
2480  Assert (interpolation_matrix.n () == this->dofs_per_face,
2481  ExcDimensionMismatch (interpolation_matrix.n (),
2482  this->dofs_per_face));
2483 
2484  // ok, source is a Nedelec element, so
2485  // we will be able to do the work
2486  const FE_Nedelec<dim> &source_fe
2487  = dynamic_cast<const FE_Nedelec<dim>&> (source);
2488 
2489  // Make sure, that the element,
2490  // for which the DoFs should be
2491  // constrained is the one with
2492  // the higher polynomial degree.
2493  // Actually the procedure will work
2494  // also if this assertion is not
2495  // satisfied. But the matrices
2496  // produced in that case might
2497  // lead to problems in the
2498  // hp procedures, which use this
2499  // method.
2500  Assert (this->dofs_per_face <= source_fe.dofs_per_face,
2501  typename FEL::ExcInterpolationNotImplemented ());
2502  interpolation_matrix = 0;
2503 
2504  // On lines we can just identify
2505  // all degrees of freedom.
2506  for (unsigned int i = 0; i <this->degree; ++i)
2507  interpolation_matrix (i, i) = 1.0;
2508 
2509  // In 3d we have some lines more
2510  // and a face. The procedure stays
2511  // the same as above, but we have
2512  // to take a bit more care of the
2513  // indices of the degrees of
2514  // freedom.
2515  if (dim == 3)
2516  {
2517  const unsigned int p = source_fe.degree;
2518  const unsigned int q = this->degree;
2519 
2520  for (unsigned int i = 0; i <q; ++i)
2521  {
2522  for (int j = 1; j < (int) GeometryInfo<dim>::lines_per_face; ++j)
2523  interpolation_matrix (j * p + i,
2524  j * q + i) = 1.0;
2525 
2526  for (unsigned int j = 0; j < q-1; ++j)
2527  {
2528  interpolation_matrix (GeometryInfo<dim>::lines_per_face * p + i * (p - 1) + j,
2529  GeometryInfo<dim>::lines_per_face * q + i * (q - 1) + j)
2530  = 1.0;
2531  interpolation_matrix (GeometryInfo<dim>::lines_per_face * p + i + (j + p - 1) * p,
2532  GeometryInfo<dim>::lines_per_face * q + i + (j + q - 1) * q)
2533  = 1.0;
2534  }
2535  }
2536  }
2537 }
2538 
2539 
2540 
2541 template <>
2542 void
2543 FE_Nedelec<1>::get_subface_interpolation_matrix(
2544  const FiniteElement<1,1> &,
2545  const unsigned int,
2546  FullMatrix<double> &) const
2547 {
2548  Assert (false, ExcNotImplemented ());
2549 }
2550 
2551 
2552 
2553 // In this function we compute the
2554 // subface interpolation matrix.
2555 // This is done by a projection-
2556 // based interpolation. Therefore
2557 // we first interpolate the
2558 // shape functions of the higher
2559 // order element on the lowest
2560 // order edge shape functions.
2561 // Then the remaining part of
2562 // the interpolated shape
2563 // functions is projected on the
2564 // higher order edge shape
2565 // functions, the face shape
2566 // functions and the interior
2567 // shape functions (if they all
2568 // exist).
2569 template <int dim>
2570 void
2571 FE_Nedelec<dim>::get_subface_interpolation_matrix(
2572  const FiniteElement<dim> &source,
2573  const unsigned int subface,
2574  FullMatrix<double> &interpolation_matrix) const
2575 {
2576  // this is only implemented, if the
2577  // source FE is also a
2578  // Nedelec element
2579  typedef FE_Nedelec<dim> FEN;
2580  typedef FiniteElement<dim> FEL;
2581 
2582  AssertThrow ((source.get_name ().find ("FE_Nedelec<") == 0) ||
2583  (dynamic_cast<const FEN *> (&source) != 0),
2584  typename FEL::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 FEL::ExcInterpolationNotImplemented ());
2610  interpolation_matrix = 0.0;
2611  // Perform projection-based interpolation
2612  // as usual.
2613  const QGauss<1> edge_quadrature (source_fe.degree);
2614  const std::vector<Point<1> > &
2615  edge_quadrature_points = edge_quadrature.get_points ();
2616  const unsigned int &n_edge_quadrature_points = edge_quadrature.size ();
2617 
2618  switch (dim)
2619  {
2620  case 2:
2621  {
2622  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2623  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2624  ++q_point)
2625  {
2626  const Point<dim> quadrature_point (0.0,
2627  0.5 * (edge_quadrature_points[q_point] (0)
2628  + subface));
2629 
2630  interpolation_matrix (0, dof) += 0.5
2631  * edge_quadrature.weight (q_point)
2632  * this->shape_value_component
2633  (dof, quadrature_point, 1);
2634  }
2635 
2636  if (source_fe.degree > 1)
2637  {
2638  const std::vector<Polynomials::Polynomial<double> > &
2639  legendre_polynomials
2640  = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
2641  FullMatrix<double> system_matrix_inv (source_fe.degree - 1,
2642  source_fe.degree - 1);
2643 
2644  {
2645  FullMatrix<double> assembling_matrix (source_fe.degree - 1,
2646  n_edge_quadrature_points);
2647 
2648  for (unsigned int q_point = 0;
2649  q_point < n_edge_quadrature_points; ++q_point)
2650  {
2651  const double weight
2652  = std::sqrt (edge_quadrature.weight (q_point));
2653 
2654  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2655  assembling_matrix (i, q_point) = weight
2656  * legendre_polynomials[i + 1].value
2657  (edge_quadrature_points[q_point] (0));
2658  }
2659 
2660  FullMatrix<double> system_matrix (source_fe.degree - 1, source_fe.degree - 1);
2661 
2662  assembling_matrix.mTmult (system_matrix, assembling_matrix);
2663  system_matrix_inv.invert (system_matrix);
2664  }
2665 
2666  Vector<double> solution (source_fe.degree - 1);
2667  Vector<double> system_rhs (source_fe.degree - 1);
2668 
2669  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2670  {
2671  system_rhs = 0.0;
2672 
2673  for (unsigned int q_point = 0;
2674  q_point < n_edge_quadrature_points; ++q_point)
2675  {
2676  const Point<dim> quadrature_point_0 (0.0,
2677  0.5 * (edge_quadrature_points[q_point] (0)
2678  + subface));
2679  const Point<dim> quadrature_point_1 (0.0,
2680  edge_quadrature_points[q_point] (0));
2681  const double tmp = edge_quadrature.weight (q_point)
2682  * (0.5 * this->shape_value_component
2683  (dof, quadrature_point_0, 1)
2684  - interpolation_matrix (0,
2685  dof)
2686  * source_fe.shape_value_component
2687  (0, quadrature_point_1, 1));
2688 
2689  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2690  system_rhs (i) += tmp
2691  * legendre_polynomials[i + 1].value
2692  (edge_quadrature_points[q_point] (0));
2693  }
2694 
2695  system_matrix_inv.vmult (solution, system_rhs);
2696 
2697  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2698  if (std::abs (solution (i)) > 1e-14)
2699  interpolation_matrix (i + 1, dof) = solution (i);
2700  }
2701  }
2702 
2703  break;
2704  }
2705 
2706  case 3:
2707  {
2708  const double shifts[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 },
2709  { 0.0, 1.0 }, { 1.0, 1.0 }
2710  };
2711 
2712  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2713  for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2714  ++q_point)
2715  {
2716  const double weight = 0.5 * edge_quadrature.weight (q_point);
2717 
2718  for (unsigned int i = 0; i < 2; ++i)
2719  {
2720  Point<dim>
2721  quadrature_point (0.5 * (i + shifts[subface][0]),
2722  0.5 * (edge_quadrature_points[q_point] (0)
2723  + shifts[subface][1]),
2724  0.0);
2725 
2726  interpolation_matrix (i * source_fe.degree, dof) += weight
2727  * this->shape_value_component
2728  (this->face_to_cell_index (dof, 4),
2729  quadrature_point,
2730  1);
2731  quadrature_point
2732  = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
2733  + shifts[subface][0]),
2734  0.5 * (i + shifts[subface][1]), 0.0);
2735  interpolation_matrix ((i + 2) * source_fe.degree, dof)
2736  += weight * this->shape_value_component
2737  (this->face_to_cell_index (dof, 4),
2738  quadrature_point, 0);
2739  }
2740  }
2741 
2742  if (source_fe.degree > 1)
2743  {
2744  const std::vector<Polynomials::Polynomial<double> > &
2745  legendre_polynomials
2746  = Polynomials::Legendre::generate_complete_basis (source_fe.degree - 1);
2747  FullMatrix<double> system_matrix_inv (source_fe.degree - 1,
2748  source_fe.degree - 1);
2749 
2750  {
2751  FullMatrix<double> assembling_matrix (source_fe.degree - 1,
2752  n_edge_quadrature_points);
2753 
2754  for (unsigned int q_point = 0;
2755  q_point < n_edge_quadrature_points; ++q_point)
2756  {
2757  const double weight
2758  = std::sqrt (edge_quadrature.weight (q_point));
2759 
2760  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2761  assembling_matrix (i, q_point) = weight
2762  * legendre_polynomials[i + 1].value
2763  (edge_quadrature_points[q_point] (0));
2764  }
2765 
2766  FullMatrix<double> system_matrix (source_fe.degree - 1, source_fe.degree - 1);
2767 
2768  assembling_matrix.mTmult (system_matrix, assembling_matrix);
2769  system_matrix_inv.invert (system_matrix);
2770  }
2771 
2772  FullMatrix<double> solution (source_fe.degree - 1,
2773  GeometryInfo<dim>::lines_per_face);
2774  FullMatrix<double> system_rhs (source_fe.degree - 1,
2775  GeometryInfo<dim>::lines_per_face);
2776  Vector<double> tmp (GeometryInfo<dim>::lines_per_face);
2777 
2778  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2779  {
2780  system_rhs = 0.0;
2781 
2782  for (unsigned int q_point = 0;
2783  q_point < n_edge_quadrature_points; ++q_point)
2784  {
2785  const double weight = edge_quadrature.weight (q_point);
2786 
2787  for (unsigned int i = 0; i < 2; ++i)
2788  {
2789  Point<dim>
2790  quadrature_point_0
2791  (0.5 * (i + shifts[subface][0]),
2792  0.5 * (edge_quadrature_points[q_point] (0)
2793  + shifts[subface][1]), 0.0);
2794  Point<dim> quadrature_point_1 (i,
2795  edge_quadrature_points[q_point] (0),
2796  0.0);
2797 
2798  tmp (i) = weight
2799  * (0.5 * this->shape_value_component
2800  (this->face_to_cell_index (dof, 4),
2801  quadrature_point_0, 1)
2802  - interpolation_matrix
2803  (i * source_fe.degree, dof)
2804  * source_fe.shape_value_component
2805  (i * source_fe.degree,
2806  quadrature_point_1, 1));
2807  quadrature_point_0
2808  = Point<dim> (0.5 * (edge_quadrature_points[q_point] (0)
2809  + shifts[subface][0]),
2810  0.5 * (i + shifts[subface][1]),
2811  0.0);
2812  quadrature_point_1
2813  = Point<dim> (edge_quadrature_points[q_point] (0),
2814  i, 0.0);
2815  tmp (i + 2) = weight
2816  * (0.5 * this->shape_value_component
2817  (this->face_to_cell_index (dof, 4),
2818  quadrature_point_0, 0)
2819  - interpolation_matrix
2820  ((i + 2) * source_fe.degree,
2821  dof)
2822  * source_fe.shape_value_component
2823  ((i + 2) * source_fe.degree,
2824  quadrature_point_1, 0));
2825  }
2826 
2827  for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
2828  {
2829  const double L_i
2830  = legendre_polynomials[i + 1].value
2831  (edge_quadrature_points[q_point] (0));
2832 
2833  for (unsigned int j = 0;
2834  j < GeometryInfo<dim>::lines_per_face; ++j)
2835  system_rhs (i, j) += tmp (j) * L_i;
2836  }
2837  }
2838 
2839  system_matrix_inv.mmult (solution, system_rhs);
2840 
2841  for (unsigned int i = 0;
2842  i < GeometryInfo<dim>::lines_per_face; ++i)
2843  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2844  if (std::abs (solution (j, i)) > 1e-14)
2845  interpolation_matrix (i * source_fe.degree + j + 1,
2846  dof) = solution (j, i);
2847  }
2848 
2849  const QGauss<2> quadrature (source_fe.degree);
2850  const std::vector<Point<2> > &
2851  quadrature_points = quadrature.get_points ();
2852  const std::vector<Polynomials::Polynomial<double> > &
2853  lobatto_polynomials
2854  = Polynomials::Lobatto::generate_complete_basis
2855  (source_fe.degree);
2856  const unsigned int n_boundary_dofs
2857  = GeometryInfo<dim>::lines_per_face * source_fe.degree;
2858  const unsigned int &n_quadrature_points = quadrature.size ();
2859 
2860  {
2861  FullMatrix<double>
2862  assembling_matrix (source_fe.degree * (source_fe.degree - 1),
2863  n_quadrature_points);
2864 
2865  for (unsigned int q_point = 0; q_point < n_quadrature_points;
2866  ++q_point)
2867  {
2868  const double weight = std::sqrt (quadrature.weight (q_point));
2869 
2870  for (unsigned int i = 0; i < source_fe.degree; ++i)
2871  {
2872  const double L_i = weight
2873  * legendre_polynomials[i].value
2874  (quadrature_points[q_point] (0));
2875 
2876  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2877  assembling_matrix (i * (source_fe.degree - 1) + j,
2878  q_point)
2879  = L_i * lobatto_polynomials[j + 2].value
2880  (quadrature_points[q_point] (1));
2881  }
2882  }
2883 
2884  FullMatrix<double> system_matrix (assembling_matrix.m (),
2885  assembling_matrix.m ());
2886 
2887  assembling_matrix.mTmult (system_matrix, assembling_matrix);
2888  system_matrix_inv.reinit (system_matrix.m (),
2889  system_matrix.m ());
2890  system_matrix_inv.invert (system_matrix);
2891  }
2892 
2893  solution.reinit (system_matrix_inv.m (), 2);
2894  system_rhs.reinit (system_matrix_inv.m (), 2);
2895  tmp.reinit (2);
2896 
2897  for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof)
2898  {
2899  system_rhs = 0.0;
2900 
2901  for (unsigned int q_point = 0;
2902  q_point < n_quadrature_points; ++q_point)
2903  {
2904  Point<dim>
2905  quadrature_point
2906  (0.5 * (quadrature_points[q_point] (0)
2907  + shifts[subface][0]),
2908  0.5 * (quadrature_points[q_point] (1)
2909  + shifts[subface][1]), 0.0);
2910  tmp (0) = 0.5 * this->shape_value_component
2911  (this->face_to_cell_index (dof, 4),
2912  quadrature_point, 0);
2913  tmp (1) = 0.5 * this->shape_value_component
2914  (this->face_to_cell_index (dof, 4),
2915  quadrature_point, 1);
2916  quadrature_point
2917  = Point<dim> (quadrature_points[q_point] (0),
2918  quadrature_points[q_point] (1), 0.0);
2919 
2920  for (unsigned int i = 0; i < 2; ++i)
2921  for (unsigned int j = 0; j < source_fe.degree; ++j)
2922  {
2923  tmp (0) -= interpolation_matrix
2924  ((i + 2) * source_fe.degree + j, dof)
2925  * source_fe.shape_value_component
2926  ((i + 2) * source_fe.degree + j,
2927  quadrature_point, 0);
2928  tmp (1) -= interpolation_matrix
2929  (i * source_fe.degree + j, dof)
2930  * source_fe.shape_value_component
2931  (i * source_fe.degree + j,
2932  quadrature_point, 1);
2933  }
2934 
2935  tmp *= quadrature.weight (q_point);
2936 
2937  for (unsigned int i = 0; i < source_fe.degree; ++i)
2938  {
2939  const double L_i_0 = legendre_polynomials[i].value
2940  (quadrature_points[q_point] (0));
2941  const double L_i_1 = legendre_polynomials[i].value
2942  (quadrature_points[q_point] (1));
2943 
2944  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2945  {
2946  system_rhs (i * (source_fe.degree - 1) + j, 0)
2947  += tmp (0) * L_i_0
2948  * lobatto_polynomials[j + 2].value
2949  (quadrature_points[q_point] (1));
2950  system_rhs (i * (source_fe.degree - 1) + j, 1)
2951  += tmp (1) * L_i_1
2952  * lobatto_polynomials[j + 2].value
2953  (quadrature_points[q_point] (0));
2954  }
2955  }
2956  }
2957 
2958  system_matrix_inv.mmult (solution, system_rhs);
2959 
2960  for (unsigned int i = 0; i < source_fe.degree; ++i)
2961  for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
2962  {
2963  if (std::abs (solution (i * (source_fe.degree - 1) + j, 0))
2964  > 1e-14)
2965  interpolation_matrix (i * (source_fe.degree - 1)
2966  + j + n_boundary_dofs, dof)
2967  = solution (i * (source_fe.degree - 1) + j, 0);
2968 
2969  if (std::abs (solution (i * (source_fe.degree - 1) + j, 1))
2970  > 1e-14)
2971  interpolation_matrix (i + (j + source_fe.degree - 1)
2972  * source_fe.degree
2973  + n_boundary_dofs, dof)
2974  = solution (i * (source_fe.degree - 1) + j, 1);
2975  }
2976  }
2977  }
2978 
2979  break;
2980  }
2981 
2982  default:
2983  Assert (false, ExcNotImplemented ());
2984  }
2985 }
2986 
2987 template <int dim>
2988 const FullMatrix<double> &
2989 FE_Nedelec<dim>
2990 ::get_prolongation_matrix (const unsigned int child,
2991  const RefinementCase<dim> &refinement_case) const
2992 {
2993  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
2994  ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
2995  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
2996  ExcMessage("Prolongation matrices are only available for refined cells!"));
2997  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
2998  ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
2999 
3000  // initialization upon first request
3001  if (this->prolongation[refinement_case-1][child].n() == 0)
3002  {
3003  Threads::Mutex::ScopedLock lock(this->mutex);
3004 
3005  // if matrix got updated while waiting for the lock
3006  if (this->prolongation[refinement_case-1][child].n() ==
3007  this->dofs_per_cell)
3008  return this->prolongation[refinement_case-1][child];
3009 
3010  // now do the work. need to get a non-const version of data in order to
3011  // be able to modify them inside a const function
3012  FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim>& >(*this);
3013 
3014  // Reinit the vectors of
3015  // restriction and prolongation
3016  // matrices to the right sizes.
3017  // Restriction only for isotropic
3018  // refinement
3019 #ifdef DEBUG_NEDELEC
3020  deallog << "Embedding" << std::endl;
3021 #endif
3022  this_nonconst.reinit_restriction_and_prolongation_matrices ();
3023  // Fill prolongation matrices with embedding operators
3024  FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
3025  internal::get_embedding_computation_tolerance(this->degree));
3026 #ifdef DEBUG_NEDELEC
3027  deallog << "Restriction" << std::endl;
3028 #endif
3029  this_nonconst.initialize_restriction ();
3030  }
3031 
3032  // we use refinement_case-1 here. the -1 takes care of the origin of the
3033  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3034  // available and so the vector indices are shifted
3035  return this->prolongation[refinement_case-1][child];
3036 }
3037 
3038 template <int dim>
3039 const FullMatrix<double> &
3040 FE_Nedelec<dim>
3041 ::get_restriction_matrix (const unsigned int child,
3042  const RefinementCase<dim> &refinement_case) const
3043 {
3044  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
3045  ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
3046  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
3047  ExcMessage("Restriction matrices are only available for refined cells!"));
3048  Assert (child<GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
3049  ExcIndexRange(child,0,GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case))));
3050 
3051  // initialization upon first request
3052  if (this->restriction[refinement_case-1][child].n() == 0)
3053  {
3054  Threads::Mutex::ScopedLock lock(this->mutex);
3055 
3056  // if matrix got updated while waiting for the lock...
3057  if (this->restriction[refinement_case-1][child].n() ==
3058  this->dofs_per_cell)
3059  return this->restriction[refinement_case-1][child];
3060 
3061  // now do the work. need to get a non-const version of data in order to
3062  // be able to modify them inside a const function
3063  FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim>& >(*this);
3064 
3065  // Reinit the vectors of
3066  // restriction and prolongation
3067  // matrices to the right sizes.
3068  // Restriction only for isotropic
3069  // refinement
3070 #ifdef DEBUG_NEDELEC
3071  deallog << "Embedding" << std::endl;
3072 #endif
3073  this_nonconst.reinit_restriction_and_prolongation_matrices ();
3074  // Fill prolongation matrices with embedding operators
3075  FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true,
3076  internal::get_embedding_computation_tolerance(this->degree));
3077 #ifdef DEBUG_NEDELEC
3078  deallog << "Restriction" << std::endl;
3079 #endif
3080  this_nonconst.initialize_restriction ();
3081  }
3082 
3083  // we use refinement_case-1 here. the -1 takes care of the origin of the
3084  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3085  // available and so the vector indices are shifted
3086  return this->restriction[refinement_case-1][child];
3087 }
3088 
3089 
3090 // Interpolate a function, which is given by
3091 // its values at the generalized support
3092 // points in the finite element space on the
3093 // reference cell.
3094 // This is done as usual by projection-based
3095 // interpolation.
3096 template <int dim>
3097 void
3098 FE_Nedelec<dim>::
3099 convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
3100  std::vector<double> &nodal_values) const
3101 {
3102  const unsigned int deg = this->degree-1;
3103  Assert (support_point_values.size () == this->generalized_support_points.size (),
3104  ExcDimensionMismatch (support_point_values.size (),
3105  this->generalized_support_points.size ()));
3106  Assert (support_point_values[0].size () == this->n_components (),
3107  ExcDimensionMismatch (support_point_values[0].size (), this->n_components ()));
3108  Assert (nodal_values.size () == this->dofs_per_cell,
3109  ExcDimensionMismatch (nodal_values.size (), this->dofs_per_cell));
3110  std::fill (nodal_values.begin (), nodal_values.end (), 0.0);
3111 
3112  switch (dim)
3113  {
3114  case 2:
3115  {
3116  // Let us begin with the
3117  // interpolation part.
3118  const QGauss<dim - 1> reference_edge_quadrature (this->degree);
3119  const unsigned int &
3120  n_edge_points = reference_edge_quadrature.size ();
3121 
3122  for (unsigned int i = 0; i < 2; ++i)
3123  for (unsigned int j = 0; j < 2; ++j)
3124  {
3125  for (unsigned int q_point = 0; q_point < n_edge_points;
3126  ++q_point)
3127  nodal_values[(i + 2 * j) * this->degree]
3128  += reference_edge_quadrature.weight (q_point)
3129  * support_point_values[q_point + (i + 2 * j) * n_edge_points][1 - j];
3130 
3131  // Add the computed support_point_values to the resulting vector
3132  // only, if they are not
3133  // too small.
3134  if (std::abs (nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3135  nodal_values[(i + 2 * j) * this->degree] = 0.0;
3136  }
3137 
3138  // If the degree is greater
3139  // than 0, then we have still
3140  // some higher order edge
3141  // shape functions to
3142  // consider.
3143  // Here the projection part
3144  // starts. The dof support_point_values
3145  // are obtained by solving
3146  // a linear system of
3147  // equations.
3148  if (this->degree-1 > 1)
3149  {
3150  // We start with projection
3151  // on the higher order edge
3152  // shape function.
3153  const std::vector<Polynomials::Polynomial<double> > &
3154  lobatto_polynomials
3155  = Polynomials::Lobatto::generate_complete_basis
3156  (this->degree);
3157  FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
3158  std::vector<Polynomials::Polynomial<double> >
3159  lobatto_polynomials_grad (this->degree);
3160 
3161  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
3162  ++i)
3163  lobatto_polynomials_grad[i]
3164  = lobatto_polynomials[i + 1].derivative ();
3165 
3166  // Set up the system matrix.
3167  // This can be used for all
3168  // edges.
3169  for (unsigned int i = 0; i < system_matrix.m (); ++i)
3170  for (unsigned int j = 0; j < system_matrix.n (); ++j)
3171  for (unsigned int q_point = 0; q_point < n_edge_points;
3172  ++q_point)
3173  system_matrix (i, j)
3174  += boundary_weights (q_point, j)
3175  * lobatto_polynomials_grad[i + 1].value
3176  (this->generalized_face_support_points[q_point]
3177  (1));
3178 
3179  FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
3180 
3181  system_matrix_inv.invert (system_matrix);
3182 
3183  const unsigned int
3184  line_coordinate[GeometryInfo<2>::lines_per_cell]
3185  = {1, 1, 0, 0};
3186  Vector<double> system_rhs (system_matrix.m ());
3187  Vector<double> solution (system_rhs.size ());
3188 
3189  for (unsigned int line = 0;
3190  line < GeometryInfo<dim>::lines_per_cell; ++line)
3191  {
3192  // Set up the right hand side.
3193  system_rhs = 0;
3194 
3195  for (unsigned int q_point = 0; q_point < n_edge_points;
3196  ++q_point)
3197  {
3198  const double tmp
3199  = support_point_values[line * n_edge_points + q_point][line_coordinate[line]]
3200  - nodal_values[line * this->degree]
3201  * this->shape_value_component
3202  (line * this->degree,
3203  this->generalized_support_points[line
3204  * n_edge_points
3205  + q_point],
3206  line_coordinate[line]);
3207 
3208  for (unsigned int i = 0; i < system_rhs.size (); ++i)
3209  system_rhs (i) += boundary_weights (q_point, i) * tmp;
3210  }
3211 
3212  system_matrix_inv.vmult (solution, system_rhs);
3213 
3214  // Add the computed support_point_values
3215  // to the resulting vector
3216  // only, if they are not
3217  // too small.
3218  for (unsigned int i = 0; i < solution.size (); ++i)
3219  if (std::abs (solution (i)) > 1e-14)
3220  nodal_values[line * this->degree + i + 1] = solution (i);
3221  }
3222 
3223  // Then we go on to the
3224  // interior shape
3225  // functions. Again we
3226  // set up the system
3227  // matrix and use it
3228  // for both, the
3229  // horizontal and the
3230  // vertical, interior
3231  // shape functions.
3232  const QGauss<dim> reference_quadrature (this->degree);
3233  const unsigned int &
3234  n_interior_points = reference_quadrature.size ();
3235  const std::vector<Polynomials::Polynomial<double> > &
3236  legendre_polynomials
3237  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
3238 
3239  system_matrix.reinit ((this->degree-1) * this->degree,
3240  (this->degree-1) * this->degree);
3241  system_matrix = 0;
3242 
3243  for (unsigned int i = 0; i < this->degree; ++i)
3244  for (unsigned int j = 0; j < this->degree-1; ++j)
3245  for (unsigned int k = 0; k < this->degree; ++k)
3246  for (unsigned int l = 0; l < this->degree-1; ++l)
3247  for (unsigned int q_point = 0;
3248  q_point < n_interior_points; ++q_point)
3249  system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
3250  += reference_quadrature.weight (q_point)
3251  * legendre_polynomials[i].value
3252  (this->generalized_support_points[q_point
3253  + GeometryInfo<dim>::lines_per_cell
3254  * n_edge_points]
3255  (0))
3256  * lobatto_polynomials[j + 2].value
3257  (this->generalized_support_points[q_point
3258  + GeometryInfo<dim>::lines_per_cell
3259  * n_edge_points]
3260  (1))
3261  * lobatto_polynomials_grad[k].value
3262  (this->generalized_support_points[q_point
3263  + GeometryInfo<dim>::lines_per_cell
3264  * n_edge_points]
3265  (0))
3266  * lobatto_polynomials[l + 2].value
3267  (this->generalized_support_points[q_point
3268  + GeometryInfo<dim>::lines_per_cell
3269  * n_edge_points]
3270  (1));
3271 
3272  system_matrix_inv.reinit (system_matrix.m (),
3273  system_matrix.m ());
3274  system_matrix_inv.invert (system_matrix);
3275  // Set up the right hand side
3276  // for the horizontal shape
3277  // functions.
3278  system_rhs.reinit (system_matrix_inv.m ());
3279  system_rhs = 0;
3280 
3281  for (unsigned int q_point = 0; q_point < n_interior_points;
3282  ++q_point)
3283  {
3284  double tmp
3285  = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
3286  * n_edge_points][0];
3287 
3288  for (unsigned int i = 0; i < 2; ++i)
3289  for (unsigned int j = 0; j <= deg; ++j)
3290  tmp -= nodal_values[(i + 2) * this->degree + j]
3291  * this->shape_value_component
3292  ((i + 2) * this->degree + j,
3293  this->generalized_support_points[q_point
3294  + GeometryInfo<dim>::lines_per_cell
3295  * n_edge_points],
3296  0);
3297 
3298  for (unsigned int i = 0; i <= deg; ++i)
3299  for (unsigned int j = 0; j < deg; ++j)
3300  system_rhs (i * deg + j)
3301  += reference_quadrature.weight (q_point) * tmp
3302  * lobatto_polynomials_grad[i].value
3303  (this->generalized_support_points[q_point
3304  + GeometryInfo<dim>::lines_per_cell
3305  * n_edge_points]
3306  (0))
3307  * lobatto_polynomials[j + 2].value
3308  (this->generalized_support_points[q_point
3309  + GeometryInfo<dim>::lines_per_cell
3310  * n_edge_points]
3311  (1));
3312  }
3313 
3314  solution.reinit (system_matrix.m ());
3315  system_matrix_inv.vmult (solution, system_rhs);
3316 
3317  // Add the computed support_point_values
3318  // to the resulting vector
3319  // only, if they are not
3320  // too small.
3321  for (unsigned int i = 0; i <= deg; ++i)
3322  for (unsigned int j = 0; j < deg; ++j)
3323  if (std::abs (solution (i * deg + j)) > 1e-14)
3324  nodal_values[(i + GeometryInfo<dim>::lines_per_cell) * deg
3325  + j + GeometryInfo<dim>::lines_per_cell]
3326  = solution (i * deg + j);
3327 
3328  system_rhs = 0;
3329  // Set up the right hand side
3330  // for the vertical shape
3331  // functions.
3332 
3333  for (unsigned int q_point = 0; q_point < n_interior_points;
3334  ++q_point)
3335  {
3336  double tmp
3337  = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
3338  * n_edge_points][1];
3339 
3340  for (unsigned int i = 0; i < 2; ++i)
3341  for (unsigned int j = 0; j <= deg; ++j)
3342  tmp -= nodal_values[i * this->degree + j]
3343  * this->shape_value_component
3344  (i * this->degree + j,
3345  this->generalized_support_points[q_point
3346  + GeometryInfo<dim>::lines_per_cell
3347  * n_edge_points],
3348  1);
3349 
3350  for (unsigned int i = 0; i <= deg; ++i)
3351  for (unsigned int j = 0; j < deg; ++j)
3352  system_rhs (i * deg + j)
3353  += reference_quadrature.weight (q_point) * tmp
3354  * lobatto_polynomials_grad[i].value
3355  (this->generalized_support_points[q_point
3356  + GeometryInfo<dim>::lines_per_cell
3357  * n_edge_points]
3358  (1))
3359  * lobatto_polynomials[j + 2].value
3360  (this->generalized_support_points[q_point
3361  + GeometryInfo<dim>::lines_per_cell
3362  * n_edge_points]
3363  (0));
3364  }
3365 
3366  system_matrix_inv.vmult (solution, system_rhs);
3367 
3368  // Add the computed support_point_values
3369  // to the resulting vector
3370  // only, if they are not
3371  // too small.
3372  for (unsigned int i = 0; i <= deg; ++i)
3373  for (unsigned int j = 0; j < deg; ++j)
3374  if (std::abs (solution (i * deg + j)) > 1e-14)
3375  nodal_values[i + (j + GeometryInfo<dim>::lines_per_cell
3376  + deg) * this->degree]
3377  = solution (i * deg + j);
3378  }
3379 
3380  break;
3381  }
3382 
3383  case 3:
3384  {
3385  // Let us begin with the
3386  // interpolation part.
3387  const QGauss<1> reference_edge_quadrature (this->degree);
3388  const unsigned int &
3389  n_edge_points = reference_edge_quadrature.size ();
3390 
3391  for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3392  {
3393  for (unsigned int i = 0; i < 4; ++i)
3394  nodal_values[(i + 8) * this->degree]
3395  += reference_edge_quadrature.weight (q_point)
3396  * support_point_values[q_point + (i + 8) * n_edge_points][2];
3397 
3398  for (unsigned int i = 0; i < 2; ++i)
3399  for (unsigned int j = 0; j < 2; ++j)
3400  for (unsigned int k = 0; k < 2; ++k)
3401  nodal_values[(i + 2 * (2 * j + k)) * this->degree]
3402  += reference_edge_quadrature.weight (q_point)
3403  * support_point_values[q_point + (i + 2 * (2 * j + k))
3404  * n_edge_points][1 - k];
3405  }
3406 
3407  // Add the computed support_point_values
3408  // to the resulting vector
3409  // only, if they are not
3410  // too small.
3411  for (unsigned int i = 0; i < 4; ++i)
3412  if (std::abs (nodal_values[(i + 8) * this->degree]) < 1e-14)
3413  nodal_values[(i + 8) * this->degree] = 0.0;
3414 
3415  for (unsigned int i = 0; i < 2; ++i)
3416  for (unsigned int j = 0; j < 2; ++j)
3417  for (unsigned int k = 0; k < 2; ++k)
3418  if (std::abs (nodal_values[(i + 2 * (2 * j + k)) * this->degree])
3419  < 1e-14)
3420  nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3421 
3422  // If the degree is greater
3423  // than 0, then we have still
3424  // some higher order shape
3425  // functions to consider.
3426  // Here the projection part
3427  // starts. The dof support_point_values
3428  // are obtained by solving
3429  // a linear system of
3430  // equations.
3431  if (this->degree > 1)
3432  {
3433  // We start with projection
3434  // on the higher order edge
3435  // shape function.
3436  const std::vector<Polynomials::Polynomial<double> > &
3437  lobatto_polynomials
3438  = Polynomials::Lobatto::generate_complete_basis
3439  (this->degree);
3440  FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
3441  std::vector<Polynomials::Polynomial<double> >
3442  lobatto_polynomials_grad (this->degree);
3443 
3444  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
3445  ++i)
3446  lobatto_polynomials_grad[i]
3447  = lobatto_polynomials[i + 1].derivative ();
3448 
3449  // Set up the system matrix.
3450  // This can be used for all
3451  // edges.
3452  for (unsigned int i = 0; i < system_matrix.m (); ++i)
3453  for (unsigned int j = 0; j < system_matrix.n (); ++j)
3454  for (unsigned int q_point = 0; q_point < n_edge_points;
3455  ++q_point)
3456  system_matrix (i, j)
3457  += boundary_weights (q_point, j)
3458  * lobatto_polynomials_grad[i + 1].value
3459  (this->generalized_face_support_points[q_point]
3460  (1));
3461 
3462  FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
3463 
3464  system_matrix_inv.invert (system_matrix);
3465 
3466  const unsigned int
3467  line_coordinate[GeometryInfo<3>::lines_per_cell]
3468  = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3469  Vector<double> system_rhs (system_matrix.m ());
3470  Vector<double> solution (system_rhs.size ());
3471 
3472  for (unsigned int line = 0;
3473  line < GeometryInfo<dim>::lines_per_cell; ++line)
3474  {
3475  // Set up the right hand side.
3476  system_rhs = 0;
3477 
3478  for (unsigned int q_point = 0; q_point < this->degree; ++q_point)
3479  {
3480  const double tmp
3481  = support_point_values[line * this->degree
3482  + q_point][line_coordinate[line]]
3483  - nodal_values[line * this->degree]
3484  * this->shape_value_component
3485  (line * this->degree,
3486  this->generalized_support_points[line
3487  * this->degree
3488  + q_point],
3489  line_coordinate[line]);
3490 
3491  for (unsigned int i = 0; i < system_rhs.size (); ++i)
3492  system_rhs (i) += boundary_weights (q_point, i)
3493  * tmp;
3494  }
3495 
3496  system_matrix_inv.vmult (solution, system_rhs);
3497 
3498  // Add the computed values
3499  // to the resulting vector
3500  // only, if they are not
3501  // too small.
3502  for (unsigned int i = 0; i < solution.size (); ++i)
3503  if (std::abs (solution (i)) > 1e-14)
3504  nodal_values[line * this->degree + i + 1] = solution (i);
3505  }
3506 
3507  // Then we go on to the
3508  // face shape functions.
3509  // Again we set up the
3510  // system matrix and
3511  // use it for both, the
3512  // horizontal and the
3513  // vertical, shape
3514  // functions.
3515  const std::vector<Polynomials::Polynomial<double> > &
3516  legendre_polynomials
3517  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
3518  const unsigned int n_face_points = n_edge_points * n_edge_points;
3519 
3520  system_matrix.reinit ((this->degree-1) * this->degree,
3521  (this->degree-1) * this->degree);
3522  system_matrix = 0;
3523 
3524  for (unsigned int i = 0; i < this->degree; ++i)
3525  for (unsigned int j = 0; j < this->degree-1; ++j)
3526  for (unsigned int k = 0; k < this->degree; ++k)
3527  for (unsigned int l = 0; l < this->degree-1; ++l)
3528  for (unsigned int q_point = 0; q_point < n_face_points;
3529  ++q_point)
3530  system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
3531  += boundary_weights (q_point + n_edge_points,
3532  2 * (k * (this->degree-1) + l))
3533  * legendre_polynomials[i].value
3534  (this->generalized_face_support_points[q_point
3535  + 4
3536  * n_edge_points]
3537  (0))
3538  * lobatto_polynomials[j + 2].value
3539  (this->generalized_face_support_points[q_point
3540  + 4
3541  * n_edge_points]
3542  (1));
3543 
3544  system_matrix_inv.reinit (system_matrix.m (),
3545  system_matrix.m ());
3546  system_matrix_inv.invert (system_matrix);
3547  solution.reinit (system_matrix.m ());
3548  system_rhs.reinit (system_matrix.m ());
3549 
3550  const unsigned int
3551  face_coordinates[GeometryInfo<3>::faces_per_cell][2]
3552  = {{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3553  const unsigned int
3554  edge_indices[GeometryInfo<3>::faces_per_cell][GeometryInfo<3>::lines_per_face]
3555  = {{0, 4, 8, 10}, {1, 5, 9, 11}, {8, 9, 2, 6},
3556  {10, 11, 3, 7}, {2, 3, 0, 1}, {6, 7, 4, 5}
3557  };
3558 
3559  for (unsigned int face = 0;
3560  face < GeometryInfo<dim>::faces_per_cell; ++face)
3561  {
3562  // Set up the right hand side
3563  // for the horizontal shape
3564  // functions.
3565  system_rhs = 0;
3566 
3567  for (unsigned int q_point = 0; q_point < n_face_points;
3568  ++q_point)
3569  {
3570  double tmp
3571  = support_point_values[q_point
3572  + GeometryInfo<dim>::lines_per_cell
3573  * n_edge_points][face_coordinates[face][0]];
3574 
3575  for (unsigned int i = 0; i < 2; ++i)
3576  for (unsigned int j = 0; j <= deg; ++j)
3577  tmp -= nodal_values[edge_indices[face][i]
3578  * this->degree + j]
3579  * this->shape_value_component
3580  (edge_indices[face][i] * this->degree + j,
3581  this->generalized_support_points[q_point
3582  + GeometryInfo<dim>::lines_per_cell
3583  * n_edge_points],
3584  face_coordinates[face][0]);
3585 
3586  for (unsigned int i = 0; i <= deg; ++i)
3587  for (unsigned int j = 0; j < deg; ++j)
3588  system_rhs (i * deg + j)
3589  += boundary_weights (q_point + n_edge_points,
3590  2 * (i * deg + j)) * tmp;
3591  }
3592 
3593  system_matrix_inv.vmult (solution, system_rhs);
3594 
3595  // Add the computed support_point_values
3596  // to the resulting vector
3597  // only, if they are not
3598  // too small.
3599  for (unsigned int i = 0; i <= deg; ++i)
3600  for (unsigned int j = 0; j < deg; ++j)
3601  if (std::abs (solution (i * deg + j)) > 1e-14)
3602  nodal_values[(2 * face * this->degree + i
3603  + GeometryInfo<dim>::lines_per_cell) * deg
3604  + j + GeometryInfo<dim>::lines_per_cell]
3605  = solution (i * deg + j);
3606 
3607  // Set up the right hand side
3608  // for the vertical shape
3609  // functions.
3610  system_rhs = 0;
3611 
3612  for (unsigned int q_point = 0; q_point < n_face_points;
3613  ++q_point)
3614  {
3615  double tmp
3616  = support_point_values[q_point
3617  + GeometryInfo<dim>::lines_per_cell
3618  * n_edge_points][face_coordinates[face][1]];
3619 
3620  for (int i = 2; i < (int) GeometryInfo<dim>::lines_per_face; ++i)
3621  for (unsigned int j = 0; j <= deg; ++j)
3622  tmp -= nodal_values[edge_indices[face][i]
3623  * this->degree + j]
3624  * this->shape_value_component
3625  (edge_indices[face][i] * this->degree + j,
3626  this->generalized_support_points[q_point
3627  + GeometryInfo<dim>::lines_per_cell
3628  * n_edge_points],
3629  face_coordinates[face][1]);
3630 
3631  for (unsigned int i = 0; i <= deg; ++i)
3632  for (unsigned int j = 0; j < deg; ++j)
3633  system_rhs (i * deg + j)
3634  += boundary_weights (q_point + n_edge_points,
3635  2 * (i * deg + j) + 1)
3636  * tmp;
3637  }
3638 
3639  system_matrix_inv.vmult (solution, system_rhs);
3640 
3641  // Add the computed support_point_values
3642  // to the resulting vector
3643  // only, if they are not
3644  // too small.
3645  for (unsigned int i = 0; i <= deg; ++i)
3646  for (unsigned int j = 0; j < deg; ++j)
3647  if (std::abs (solution (i * deg + j)) > 1e-14)
3648  nodal_values[((2 * face + 1) * deg + j + GeometryInfo<dim>::lines_per_cell)
3649  * this->degree + i]
3650  = solution (i * deg + j);
3651  }
3652 
3653  // Finally we project
3654  // the remaining parts
3655  // of the function on
3656  // the interior shape
3657  // functions.
3658  const QGauss<dim> reference_quadrature (this->degree);
3659  const unsigned int
3660  n_interior_points = reference_quadrature.size ();
3661 
3662  // We create the
3663  // system matrix.
3664  system_matrix.reinit (this->degree * deg * deg,
3665  this->degree * deg * deg);
3666  system_matrix = 0;
3667 
3668  for (unsigned int i = 0; i <= deg; ++i)
3669  for (unsigned int j = 0; j < deg; ++j)
3670  for (unsigned int k = 0; k < deg; ++k)
3671  for (unsigned int l = 0; l <= deg; ++l)
3672  for (unsigned int m = 0; m < deg; ++m)
3673  for (unsigned int n = 0; n < deg; ++n)
3674  for (unsigned int q_point = 0;
3675  q_point < n_interior_points; ++q_point)
3676  system_matrix ((i * deg + j) * deg + k,
3677  (l * deg + m) * deg + n)
3678  += reference_quadrature.weight (q_point)
3679  * legendre_polynomials[i].value
3680  (this->generalized_support_points[q_point
3681  + GeometryInfo<dim>::lines_per_cell
3682  * n_edge_points
3683  + GeometryInfo<dim>::faces_per_cell
3684  * n_face_points]
3685  (0))
3686  * lobatto_polynomials[j + 2].value
3687  (this->generalized_support_points[q_point
3688  + GeometryInfo<dim>::lines_per_cell
3689  * n_edge_points
3690  + GeometryInfo<dim>::faces_per_cell
3691  * n_face_points]
3692  (1))
3693  * lobatto_polynomials[k + 2].value
3694  (this->generalized_support_points[q_point
3695  + GeometryInfo<dim>::lines_per_cell
3696  * n_edge_points
3697  + GeometryInfo<dim>::faces_per_cell
3698  * n_face_points]
3699  (2))
3700  * lobatto_polynomials_grad[l].value
3701  (this->generalized_support_points[q_point
3702  + GeometryInfo<dim>::lines_per_cell
3703  * n_edge_points
3704  + GeometryInfo<dim>::faces_per_cell
3705  * n_face_points]
3706  (0))
3707  * lobatto_polynomials[m + 2].value
3708  (this->generalized_support_points[q_point
3709  + GeometryInfo<dim>::lines_per_cell
3710  * n_edge_points
3711  + GeometryInfo<dim>::faces_per_cell
3712  * n_face_points]
3713  (1))
3714  * lobatto_polynomials[n + 2].value
3715  (this->generalized_support_points[q_point
3716  + GeometryInfo<dim>::lines_per_cell
3717  * n_edge_points
3718  + GeometryInfo<dim>::faces_per_cell
3719  * n_face_points]
3720  (2));
3721 
3722  system_matrix_inv.reinit (system_matrix.m (),
3723  system_matrix.m ());
3724  system_matrix_inv.invert (system_matrix);
3725  // Set up the right hand side.
3726  system_rhs.reinit (system_matrix.m ());
3727  system_rhs = 0;
3728 
3729  for (unsigned int q_point = 0; q_point < n_interior_points;
3730  ++q_point)
3731  {
3732  double tmp
3733  = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
3734  * n_edge_points
3735  + GeometryInfo<dim>::faces_per_cell
3736  * n_face_points][0];
3737 
3738  for (unsigned int i = 0; i <= deg; ++i)
3739  {
3740  for (unsigned int j = 0; j < 2; ++j)
3741  for (unsigned int k = 0; k < 2; ++k)
3742  tmp -= nodal_values[i + (j + 4 * k + 2) * this->degree]
3743  * this->shape_value_component
3744  (i + (j + 4 * k + 2) * this->degree,
3745  this->generalized_support_points[q_point
3746  + GeometryInfo<dim>::lines_per_cell
3747  * n_edge_points
3748  + GeometryInfo<dim>::faces_per_cell
3749  * n_face_points],
3750  0);
3751 
3752  for (unsigned int j = 0; j < deg; ++j)
3753  for (unsigned int k = 0; k < 4; ++k)
3754  tmp -= nodal_values[(i + 2 * (k + 2) * this->degree
3755  + GeometryInfo<dim>::lines_per_cell)
3756  * deg + j
3757  + GeometryInfo<dim>::lines_per_cell]
3758  * this->shape_value_component
3759  ((i + 2 * (k + 2) * this->degree
3760  + GeometryInfo<dim>::lines_per_cell)
3761  * deg + j
3762  + GeometryInfo<dim>::lines_per_cell,
3763  this->generalized_support_points[q_point
3764  + GeometryInfo<dim>::lines_per_cell
3765  * n_edge_points
3766  + GeometryInfo<dim>::faces_per_cell
3767  * n_face_points],
3768  0);
3769  }
3770 
3771  for (unsigned int i = 0; i <= deg; ++i)
3772  for (unsigned int j = 0; j < deg; ++j)
3773  for (unsigned int k = 0; k < deg; ++k)
3774  system_rhs ((i * deg + j) * deg + k)
3775  += reference_quadrature.weight (q_point) * tmp
3776  * lobatto_polynomials_grad[i].value
3777  (this->generalized_support_points[q_point
3778  + GeometryInfo<dim>::lines_per_cell
3779  * n_edge_points
3780  + GeometryInfo<dim>::faces_per_cell
3781  * n_face_points]
3782  (0))
3783  * lobatto_polynomials[j + 2].value
3784  (this->generalized_support_points[q_point
3785  + GeometryInfo<dim>::lines_per_cell
3786  * n_edge_points
3787  + GeometryInfo<dim>::faces_per_cell
3788  * n_face_points]
3789  (1))
3790  * lobatto_polynomials[k + 2].value
3791  (this->generalized_support_points[q_point
3792  + GeometryInfo<dim>::lines_per_cell
3793  * n_edge_points
3794  + GeometryInfo<dim>::faces_per_cell
3795  * n_face_points]
3796  (2));
3797  }
3798 
3799  solution.reinit (system_rhs.size ());
3800  system_matrix_inv.vmult (solution, system_rhs);
3801 
3802  // Add the computed values
3803  // to the resulting vector
3804  // only, if they are not
3805  // too small.
3806  for (unsigned int i = 0; i <= deg; ++i)
3807  for (unsigned int j = 0; j < deg; ++j)
3808  for (unsigned int k = 0; k < deg; ++k)
3809  if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
3810  nodal_values[((i + 2 * GeometryInfo<dim>::faces_per_cell)
3811  * deg + j + GeometryInfo<dim>::lines_per_cell
3812  + 2 * GeometryInfo<dim>::faces_per_cell)
3813  * deg + k + GeometryInfo<dim>::lines_per_cell]
3814  = solution ((i * deg + j) * deg + k);
3815 
3816  // Set up the right hand side.
3817  system_rhs = 0;
3818 
3819  for (unsigned int q_point = 0; q_point < n_interior_points;
3820  ++q_point)
3821  {
3822  double tmp
3823  = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
3824  * n_edge_points
3825  + GeometryInfo<dim>::faces_per_cell
3826  * n_face_points][1];
3827 
3828  for (unsigned int i = 0; i <= deg; ++i)
3829  for (unsigned int j = 0; j < 2; ++j)
3830  {
3831  for (unsigned int k = 0; k < 2; ++k)
3832  tmp -= nodal_values[i + (4 * j + k) * this->degree]
3833  * this->shape_value_component
3834  (i + (4 * j + k) * this->degree,
3835  this->generalized_support_points[q_point
3836  + GeometryInfo<dim>::lines_per_cell
3837  * n_edge_points
3838  + GeometryInfo<dim>::faces_per_cell
3839  * n_face_points],
3840  1);
3841 
3842  for (unsigned int k = 0; k < deg; ++k)
3843  tmp -= nodal_values[(i + 2 * j * this->degree
3844  + GeometryInfo<dim>::lines_per_cell)
3845  * deg + k
3846  + GeometryInfo<dim>::lines_per_cell]
3847  * this->shape_value_component
3848  ((i + 2 * j * this->degree
3849  + GeometryInfo<dim>::lines_per_cell)
3850  * deg + k
3851  + GeometryInfo<dim>::lines_per_cell,
3852  this->generalized_support_points[q_point
3853  + GeometryInfo<dim>::lines_per_cell
3854  * n_edge_points
3855  + GeometryInfo<dim>::faces_per_cell
3856  * n_face_points],
3857  1)
3858  + nodal_values[i + ((2 * j + 9) * deg + k
3859  + GeometryInfo<dim>::lines_per_cell)
3860  * this->degree]
3861  * this->shape_value_component
3862  (i + ((2 * j + 9) * deg + k
3863  + GeometryInfo<dim>::lines_per_cell)
3864  * this->degree,
3865  this->generalized_support_points[q_point
3866  + GeometryInfo<dim>::lines_per_cell
3867  * n_edge_points
3868  + GeometryInfo<dim>::faces_per_cell
3869  * n_face_points],
3870  1);
3871  }
3872 
3873  for (unsigned int i = 0; i <= deg; ++i)
3874  for (unsigned int j = 0; j < deg; ++j)
3875  for (unsigned int k = 0; k < deg; ++k)
3876  system_rhs ((i * deg + j) * deg + k)
3877  += reference_quadrature.weight (q_point) * tmp
3878  * lobatto_polynomials_grad[i].value
3879  (this->generalized_support_points[q_point
3880  + GeometryInfo<dim>::lines_per_cell
3881  * n_edge_points
3882  + GeometryInfo<dim>::faces_per_cell
3883  * n_face_points]
3884  (1))
3885  * lobatto_polynomials[j + 2].value
3886  (this->generalized_support_points[q_point
3887  + GeometryInfo<dim>::lines_per_cell
3888  * n_edge_points
3889  + GeometryInfo<dim>::faces_per_cell
3890  * n_face_points]
3891  (0))
3892  * lobatto_polynomials[k + 2].value
3893  (this->generalized_support_points[q_point
3894  + GeometryInfo<dim>::lines_per_cell
3895  * n_edge_points
3896  + GeometryInfo<dim>::faces_per_cell
3897  * n_face_points]
3898  (2));
3899  }
3900 
3901  system_matrix_inv.vmult (solution, system_rhs);
3902 
3903  // Add the computed support_point_values
3904  // to the resulting vector
3905  // only, if they are not
3906  // too small.
3907  for (unsigned int i = 0; i <= deg; ++i)
3908  for (unsigned int j = 0; j < deg; ++j)
3909  for (unsigned int k = 0; k < deg; ++k)
3910  if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
3911  nodal_values[((i + this->degree + 2
3912  * GeometryInfo<dim>::faces_per_cell) * deg
3913  + j + GeometryInfo<dim>::lines_per_cell + 2
3914  * GeometryInfo<dim>::faces_per_cell) * deg
3915  + k + GeometryInfo<dim>::lines_per_cell]
3916  = solution ((i * deg + j) * deg + k);
3917 
3918  // Set up the right hand side.
3919  system_rhs = 0;
3920 
3921  for (unsigned int q_point = 0; q_point < n_interior_points;
3922  ++q_point)
3923  {
3924  double tmp
3925  = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
3926  * n_edge_points
3927  + GeometryInfo<dim>::faces_per_cell
3928  * n_face_points][2];
3929 
3930  for (unsigned int i = 0; i <= deg; ++i)
3931  for (unsigned int j = 0; j < 4; ++j)
3932  {
3933  tmp -= nodal_values[i + (j + 8) * this->degree]
3934  * this->shape_value_component
3935  (i + (j + 8) * this->degree,
3936  this->generalized_support_points[q_point
3937  + GeometryInfo<dim>::lines_per_cell
3938  * n_edge_points
3939  + GeometryInfo<dim>::faces_per_cell
3940  * n_face_points],
3941  2);
3942 
3943  for (unsigned int k = 0; k < deg; ++k)
3944  tmp -= nodal_values[i + ((2 * j + 1) * deg + k
3945  + GeometryInfo<dim>::lines_per_cell)
3946  * this->degree]
3947  * this->shape_value_component
3948  (i + ((2 * j + 1) * deg + k
3949  + GeometryInfo<dim>::lines_per_cell)
3950  * this->degree,
3951  this->generalized_support_points[q_point
3952  + GeometryInfo<dim>::lines_per_cell
3953  * n_edge_points
3954  + GeometryInfo<dim>::faces_per_cell
3955  * n_face_points],
3956  2);
3957  }
3958 
3959  for (unsigned int i = 0; i <= deg; ++i)
3960  for (unsigned int j = 0; j < deg; ++j)
3961  for (unsigned int k = 0; k < deg; ++k)
3962  system_rhs ((i * deg + j) * deg + k)
3963  += reference_quadrature.weight (q_point) * tmp
3964  * lobatto_polynomials_grad[i].value
3965  (this->generalized_support_points[q_point
3966  + GeometryInfo<dim>::lines_per_cell
3967  * n_edge_points
3968  + GeometryInfo<dim>::faces_per_cell
3969  * n_face_points]
3970  (2))
3971  * lobatto_polynomials[j + 2].value
3972  (this->generalized_support_points[q_point
3973  + GeometryInfo<dim>::lines_per_cell
3974  * n_edge_points
3975  + GeometryInfo<dim>::faces_per_cell
3976  * n_face_points]
3977  (0))
3978  * lobatto_polynomials[k + 2].value
3979  (this->generalized_support_points[q_point
3980  + GeometryInfo<dim>::lines_per_cell
3981  * n_edge_points
3982  + GeometryInfo<dim>::faces_per_cell
3983  * n_face_points]
3984  (1));
3985  }
3986 
3987  system_matrix_inv.vmult (solution, system_rhs);
3988 
3989  // Add the computed support_point_values
3990  // to the resulting vector
3991  // only, if they are not
3992  // too small.
3993  for (unsigned int i = 0; i <= deg; ++i)
3994  for (unsigned int j = 0; j < deg; ++j)
3995  for (unsigned int k = 0; k < deg; ++k)
3996  if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
3997  nodal_values[i + ((j + 2 * (deg
3998  + GeometryInfo<dim>::faces_per_cell))
3999  * deg + k
4000  + GeometryInfo<dim>::lines_per_cell)
4001  * this->degree]
4002  = solution ((i * deg + j) * deg + k);
4003  }
4004 
4005  break;
4006  }
4007 
4008  default:
4009  Assert (false, ExcNotImplemented ());
4010  }
4011 }
4012 
4013 
4014 
4015 
4016 
4017 // Since this is a vector valued element,
4018 // we cannot interpolate a scalar function.
4019 template <int dim>
4020 void FE_Nedelec<dim>::interpolate (std::vector<double> &, const std::vector<double> &) const
4021 {
4022  Assert(false, ExcNotImplemented ());
4023 }
4024 
4025 
4026 // Interpolate a function, which is given by
4027 // its values at the generalized support
4028 // points in the finite element space on the
4029 // reference cell.
4030 // This is done as usual by projection-based
4031 // interpolation.
4032 template <int dim>
4033 void
4034 FE_Nedelec<dim>::interpolate (std::vector<double> &local_dofs,
4035  const std::vector<Vector<double> > &values,
4036  const unsigned int offset) const
4037 {
4038  const unsigned int deg = this->degree-1;
4039 
4040  Assert (values.size () == this->generalized_support_points.size (),
4041  ExcDimensionMismatch (values.size (),
4042  this->generalized_support_points.size ()));
4043  Assert (local_dofs.size () == this->dofs_per_cell,
4044  ExcDimensionMismatch (local_dofs.size (),this->dofs_per_cell));
4045  Assert (values[0].size () >= offset + this->n_components (),
4046  ExcDimensionMismatch (values[0].size (),
4047  offset + this->n_components ()));
4048  std::fill (local_dofs.begin (), local_dofs.end (), 0.);
4049 
4050  if (offset < dim)
4051  switch (dim)
4052  {
4053  case 2:
4054  {
4055  const QGauss<1> reference_edge_quadrature (this->degree);
4056  const unsigned int &n_edge_points
4057  = reference_edge_quadrature.size ();
4058 
4059  // Let us begin with the
4060  // interpolation part.
4061  for (unsigned int i = 0; i < 2; ++i)
4062  {
4063  for (unsigned int q_point = 0; q_point < n_edge_points;
4064  ++q_point)
4065  local_dofs[i * this->degree]
4066  += reference_edge_quadrature.weight (q_point)
4067  * values[q_point + i * n_edge_points] (1);
4068 
4069  // Add the computed values
4070  // to the resulting vector
4071  // only, if they are not
4072  // too small.
4073  if (std::abs (local_dofs[i * this->degree]) < 1e-14)
4074  local_dofs[i * this->degree] = 0.0;
4075  }
4076 
4077  if (offset == 0)
4078  for (unsigned int i = 0; i < 2; ++i)
4079  {
4080  for (unsigned int q_point = 0; q_point < n_edge_points;
4081  ++q_point)
4082  local_dofs[(i + 2) * this->degree]
4083  += reference_edge_quadrature.weight (q_point)
4084  * values[q_point + (i + 2) * n_edge_points] (0);
4085 
4086  if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14)
4087  local_dofs[(i + 2) * this->degree] = 0.0;
4088  }
4089 
4090  // If the degree is greater
4091  // than 0, then we have still
4092  // some higher order edge
4093  // shape functions to
4094  // consider.
4095  // Here the projection part
4096  // starts. The dof values
4097  // are obtained by solving
4098  // a linear system of
4099  // equations.
4100  if (this->degree > 1)
4101  {
4102  // We start with projection
4103  // on the higher order edge
4104  // shape function.
4105  const std::vector<Polynomials::Polynomial<double> > &
4106  lobatto_polynomials
4107  = Polynomials::Lobatto::generate_complete_basis
4108  (this->degree);
4109  const unsigned int
4110  line_coordinate[GeometryInfo<2>::lines_per_cell]
4111  = {1, 1, 0, 0};
4112  std::vector<Polynomials::Polynomial<double> >
4113  lobatto_polynomials_grad (this->degree);
4114 
4115  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
4116  ++i)
4117  lobatto_polynomials_grad[i]
4118  = lobatto_polynomials[i + 1].derivative ();
4119 
4120  // Set up the system matrix.
4121  // This can be used for all
4122  // edges.
4123  FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
4124 
4125  for (unsigned int i = 0; i < system_matrix.m (); ++i)
4126  for (unsigned int j = 0; j < system_matrix.n (); ++j)
4127  for (unsigned int q_point = 0; q_point < n_edge_points;
4128  ++q_point)
4129  system_matrix (i, j)
4130  += boundary_weights (q_point, j)
4131  * lobatto_polynomials_grad[i + 1].value
4132  (this->generalized_face_support_points[q_point]
4133  (1));
4134 
4135  FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
4136 
4137  system_matrix_inv.invert (system_matrix);
4138 
4139  Vector<double> system_rhs (system_matrix.m ());
4140  Vector<double> solution (system_rhs.size ());
4141 
4142  for (unsigned int line = 0;
4143  line < GeometryInfo<dim>::lines_per_cell; ++line)
4144  if ((line < 2) || (offset == 0))
4145  {
4146  // Set up the right hand side.
4147  system_rhs = 0;
4148 
4149  for (unsigned int q_point = 0; q_point < n_edge_points;
4150  ++q_point)
4151  {
4152  const double tmp
4153  = values[line * n_edge_points + q_point]
4154  (line_coordinate[line])
4155  - local_dofs[line * this->degree]
4156  * this->shape_value_component
4157  (line * this->degree,
4158  this->generalized_support_points[line
4159  * n_edge_points
4160  + q_point],
4161  line_coordinate[line]);
4162 
4163  for (unsigned int i = 0; i < system_rhs.size ();
4164  ++i)
4165  system_rhs (i) += boundary_weights (q_point, i)
4166  * tmp;
4167  }
4168 
4169  system_matrix_inv.vmult (solution, system_rhs);
4170 
4171  // Add the computed values
4172  // to the resulting vector
4173  // only, if they are not
4174  // too small.
4175  for (unsigned int i = 0; i < solution.size (); ++i)
4176  if (std::abs (solution (i)) > 1e-14)
4177  local_dofs[line * this->degree + i + 1]
4178  = solution (i);
4179  }
4180 
4181  // Then we go on to the
4182  // interior shape
4183  // functions. Again we
4184  // set up the system
4185  // matrix and use it
4186  // for both, the
4187  // horizontal and the
4188  // vertical, interior
4189  // shape functions.
4190  const QGauss<dim> reference_quadrature (this->degree);
4191  const std::vector<Polynomials::Polynomial<double> > &
4192  legendre_polynomials
4193  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
4194  const unsigned int &n_interior_points
4195  = reference_quadrature.size ();
4196 
4197  system_matrix.reinit ((this->degree-1) * this->degree,
4198  (this->degree-1) * this->degree);
4199  system_matrix = 0;
4200 
4201  for (unsigned int i = 0; i < this->degree; ++i)
4202  for (unsigned int j = 0; j < this->degree-1; ++j)
4203  for (unsigned int k = 0; k < this->degree; ++k)
4204  for (unsigned int l = 0; l < this->degree-1; ++l)
4205  for (unsigned int q_point = 0;
4206  q_point < n_interior_points; ++q_point)
4207  system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
4208  += reference_quadrature.weight (q_point)
4209  * legendre_polynomials[i].value
4210  (this->generalized_support_points[q_point
4211  + GeometryInfo<dim>::lines_per_cell
4212  * n_edge_points]
4213  (0))
4214  * lobatto_polynomials[j + 2].value
4215  (this->generalized_support_points[q_point
4216  + GeometryInfo<dim>::lines_per_cell
4217  * n_edge_points]
4218  (1))
4219  * lobatto_polynomials_grad[k].value
4220  (this->generalized_support_points[q_point
4221  + GeometryInfo<dim>::lines_per_cell
4222  * n_edge_points]
4223  (0))
4224  * lobatto_polynomials[l + 2].value
4225  (this->generalized_support_points[q_point
4226  + GeometryInfo<dim>::lines_per_cell
4227  * n_edge_points]
4228  (1));
4229 
4230  system_matrix_inv.reinit (system_matrix.m (),
4231  system_matrix.m ());
4232  system_matrix_inv.invert (system_matrix);
4233  solution.reinit (system_matrix_inv.m ());
4234  system_rhs.reinit (system_matrix.m ());
4235 
4236  if (offset == 0)
4237  {
4238  // Set up the right hand side
4239  // for the horizontal shape
4240  // functions.
4241  system_rhs = 0;
4242 
4243  for (unsigned int q_point = 0;
4244  q_point < n_interior_points; ++q_point)
4245  {
4246  double tmp
4247  = values[q_point + GeometryInfo<dim>::lines_per_cell
4248  * n_edge_points] (0);
4249 
4250  for (unsigned int i = 0; i < 2; ++i)
4251  for (unsigned int j = 0; j < this->degree; ++j)
4252  tmp -= local_dofs[(i + 2) * this->degree + j]
4253  * this->shape_value_component
4254  ((i + 2) * this->degree + j,
4255  this->generalized_support_points[q_point
4256  + GeometryInfo<dim>::lines_per_cell
4257  * n_edge_points],
4258  0);
4259 
4260  for (unsigned int i = 0; i < this->degree; ++i)
4261  for (unsigned int j = 0; j < this->degree-1; ++j)
4262  system_rhs (i * (this->degree-1) + j)
4263  += reference_quadrature.weight (q_point) * tmp
4264  * lobatto_polynomials_grad[i].value
4265  (this->generalized_support_points[q_point
4266  + GeometryInfo<dim>::lines_per_cell
4267  * n_edge_points]
4268  (0))
4269  * lobatto_polynomials[j + 2].value
4270  (this->generalized_support_points[q_point
4271  + GeometryInfo<dim>::lines_per_cell
4272  * n_edge_points]
4273  (1));
4274  }
4275 
4276  system_matrix_inv.vmult (solution, system_rhs);
4277 
4278  // Add the computed values
4279  // to the resulting vector
4280  // only, if they are not
4281  // too small.
4282  for (unsigned int i = 0; i < this->degree; ++i)
4283  for (unsigned int j = 0; j < this->degree-1; ++j)
4284  if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
4285  local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
4286  * (this->degree-1) + j
4287  + GeometryInfo<dim>::lines_per_cell]
4288  = solution (i * (this->degree-1) + j);
4289  }
4290 
4291  // Set up the right hand side
4292  // for the vertical shape
4293  // functions.
4294  system_rhs = 0;
4295 
4296  for (unsigned int q_point = 0; q_point < n_interior_points;
4297  ++q_point)
4298  {
4299  double tmp
4300  = values[q_point + GeometryInfo<dim>::lines_per_cell
4301  * n_edge_points] (1);
4302 
4303  for (unsigned int i = 0; i < 2; ++i)
4304  for (unsigned int j = 0; j < this->degree; ++j)
4305  tmp -= local_dofs[i * this->degree + j]
4306  * this->shape_value_component
4307  (i * this->degree + j,
4308  this->generalized_support_points[q_point
4309  + GeometryInfo<dim>::lines_per_cell
4310  * n_edge_points],
4311  1);
4312 
4313  for (unsigned int i = 0; i < this->degree; ++i)
4314  for (unsigned int j = 0; j < this->degree-1; ++j)
4315  system_rhs (i * (this->degree-1) + j)
4316  += reference_quadrature.weight (q_point) * tmp
4317  * lobatto_polynomials_grad[i].value
4318  (this->generalized_support_points[q_point
4319  + GeometryInfo<dim>::lines_per_cell
4320  * n_edge_points]
4321  (1))
4322  * lobatto_polynomials[j + 2].value
4323  (this->generalized_support_points[q_point
4324  + GeometryInfo<dim>::lines_per_cell
4325  * n_edge_points]
4326  (0));
4327  }
4328 
4329  system_matrix_inv.vmult (solution, system_rhs);
4330 
4331  // Add the computed values
4332  // to the resulting vector
4333  // only, if they are not
4334  // too small.
4335  for (unsigned int i = 0; i < this->degree; ++i)
4336  for (unsigned int j = 0; j < this->degree-1; ++j)
4337  if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
4338  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
4339  + this->degree-1) * this->degree]
4340  = solution (i * (this->degree-1) + j);
4341  }
4342 
4343  break;
4344  }
4345 
4346  case 3:
4347  {
4348  const QGauss<1>
4349  reference_edge_quadrature (this->degree);
4350  const unsigned int &
4351  n_edge_points = reference_edge_quadrature.size ();
4352 
4353  // Let us begin with the
4354  // interpolation part.
4355  for (unsigned int i = 0; i < 4; ++i)
4356  {
4357  for (unsigned int q_point = 0; q_point < n_edge_points;
4358  ++q_point)
4359  local_dofs[(i + 8) * this->degree]
4360  += reference_edge_quadrature.weight (q_point)
4361  * values[q_point + (i + 8) * n_edge_points] (2);
4362 
4363  // Add the computed values
4364  // to the resulting vector
4365  // only, if they are not
4366  // too small.
4367  if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
4368  local_dofs[(i + 8) * this->degree] = 0.0;
4369  }
4370 
4371  if (offset + 1 < dim)
4372  {
4373  for (unsigned int i = 0; i < 2; ++i)
4374  for (unsigned int j = 0; j < 2; ++j)
4375  {
4376  for (unsigned int q_point = 0; q_point < n_edge_points;
4377  ++q_point)
4378  local_dofs[(i + 4 * j) * this->degree]
4379  += reference_edge_quadrature.weight (q_point)
4380  * values[q_point + (i + 4 * j) * n_edge_points]
4381  (1);
4382 
4383  // Add the computed values
4384  // to the resulting vector
4385  // only, if they are not
4386  // too small.
4387  if (std::abs (local_dofs[(i + 4 * j) * this->degree])
4388  < 1e-14)
4389  local_dofs[(i + 4 * j) * this->degree] = 0.0;
4390  }
4391 
4392  if (offset == 0)
4393  for (unsigned int i = 0; i < 2; ++i)
4394  for (unsigned int j = 0; j < 2; ++j)
4395  {
4396  for (unsigned int q_point = 0;
4397  q_point < n_edge_points; ++q_point)
4398  local_dofs[(i + 4 * j + 2) * this->degree]
4399  += reference_edge_quadrature.weight (q_point)
4400  * values[q_point + (i + 4 * j + 2)
4401  * n_edge_points] (0);
4402 
4403  // Add the computed values
4404  // to the resulting vector
4405  // only, if they are not
4406  // too small.
4407  if (std::abs (local_dofs[(i + 4 * j + 2)
4408  * this->degree]) < 1e-14)
4409  local_dofs[(i + 4 * j + 2) * this->degree] = 0.0;
4410  }
4411  }
4412 
4413  // If the degree is greater
4414  // than 0, then we have still
4415  // some higher order shape
4416  // functions to consider.
4417  // Here the projection part
4418  // starts. The dof values
4419  // are obtained by solving
4420  // a linear system of
4421  // equations.
4422  if (this->degree > 1)
4423  {
4424  // We start with projection
4425  // on the higher order edge
4426  // shape function.
4427  const std::vector<Polynomials::Polynomial<double> > &
4428  lobatto_polynomials
4429  = Polynomials::Lobatto::generate_complete_basis
4430  (this->degree);
4431  const unsigned int
4432  line_coordinate[GeometryInfo<3>::lines_per_cell]
4433  = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
4434  FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
4435  FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
4436  std::vector<Polynomials::Polynomial<double> >
4437  lobatto_polynomials_grad (this->degree);
4438 
4439  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
4440  ++i)
4441  lobatto_polynomials_grad[i]
4442  = lobatto_polynomials[i + 1].derivative ();
4443 
4444  Vector<double> system_rhs (system_matrix.m ());
4445  Vector<double> solution (system_rhs.size ());
4446 
4447  // Set up the system matrix.
4448  // This can be used for all
4449  // edges.
4450  for (unsigned int i = 0; i < system_matrix.m (); ++i)
4451  for (unsigned int j = 0; j < system_matrix.n (); ++j)
4452  for (unsigned int q_point = 0; q_point < n_edge_points;
4453  ++q_point)
4454  system_matrix (i, j)
4455  += boundary_weights (q_point, j)
4456  * lobatto_polynomials_grad[i + 1].value
4457  (this->generalized_face_support_points[q_point]
4458  (1));
4459 
4460  system_matrix_inv.invert (system_matrix);
4461 
4462  for (unsigned int line = 0;
4463  line < GeometryInfo<dim>::lines_per_cell; ++line)
4464  {
4465  // Set up the right hand side.
4466  system_rhs = 0;
4467 
4468  if ((((line == 0) || (line == 1) || (line == 4) ||
4469  (line == 5)) && (offset + 1 < dim)) ||
4470  (((line == 2) || (line == 3) || (line == 6) ||
4471  (line == 7)) && (offset == 0)) || (line > 7))
4472  {
4473  for (unsigned int q_point = 0; q_point < n_edge_points;
4474  ++q_point)
4475  {
4476  double tmp
4477  = values[line * n_edge_points + q_point]
4478  (line_coordinate[line])
4479  - local_dofs[line * this->degree]
4480  * this->shape_value_component
4481  (line * this->degree,
4482  this->generalized_support_points[line
4483  * this->degree
4484  + q_point],
4485  line_coordinate[line]);
4486 
4487  for (unsigned int i = 0; i < system_rhs.size ();
4488  ++i)
4489  system_rhs (i)
4490  += boundary_weights (q_point, i) * tmp;
4491  }
4492 
4493  system_matrix_inv.vmult (solution, system_rhs);
4494 
4495  // Add the computed values
4496  // to the resulting vector
4497  // only, if they are not
4498  // too small.
4499  for (unsigned int i = 0; i < solution.size (); ++i)
4500  if (std::abs (solution (i)) > 1e-14)
4501  local_dofs[line * this->degree + i + 1]
4502  = solution (i);
4503  }
4504  }
4505 
4506  // Then we go on to the
4507  // face shape functions.
4508  // Again we set up the
4509  // system matrix and
4510  // use it for both, the
4511  // horizontal and the
4512  // vertical, shape
4513  // functions.
4514  const std::vector<Polynomials::Polynomial<double> > &
4515  legendre_polynomials
4516  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
4517  const unsigned int
4518  n_face_points = n_edge_points * n_edge_points;
4519 
4520  system_matrix.reinit ((this->degree-1) * this->degree,
4521  (this->degree-1) * this->degree);
4522  system_matrix = 0;
4523 
4524  for (unsigned int i = 0; i < this->degree; ++i)
4525  for (unsigned int j = 0; j < this->degree-1; ++j)
4526  for (unsigned int k = 0; k < this->degree; ++k)
4527  for (unsigned int l = 0; l < this->degree-1; ++l)
4528  for (unsigned int q_point = 0; q_point < n_face_points;
4529  ++q_point)
4530  system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
4531  += boundary_weights (q_point + n_edge_points,
4532  2 * (k * (this->degree-1) + l))
4533  * legendre_polynomials[i].value
4534  (this->generalized_face_support_points[q_point
4535  + 4
4536  * n_edge_points]
4537  (0))
4538  * lobatto_polynomials[j + 2].value
4539  (this->generalized_face_support_points[q_point
4540  + 4
4541  * n_edge_points]
4542  (1));
4543 
4544  system_matrix_inv.reinit (system_matrix.m (),
4545  system_matrix.n ());
4546  system_matrix_inv.invert (system_matrix);
4547  solution.reinit (system_matrix.m ());
4548  system_rhs.reinit (system_matrix.m ());
4549 
4550  for (unsigned int face = 0;
4551  face < GeometryInfo<dim>::faces_per_cell; ++face)
4552  {
4553  switch (face)
4554  {
4555  case 0:
4556  {
4557  if (offset + 1 < dim)
4558  {
4559  // Set up the right hand side
4560  // for the horizontal shape
4561  // functions.
4562  system_rhs = 0;
4563 
4564  for (unsigned int q_point = 0;
4565  q_point < n_face_points; ++q_point)
4566  {
4567  double tmp
4568  = values[q_point
4569  + GeometryInfo<dim>::lines_per_cell
4570  * n_edge_points] (1);
4571 
4572  for (unsigned int i = 0; i < 2; ++i)
4573  for (unsigned int j = 0; j < this->degree; ++j)
4574  tmp
4575  -= local_dofs[4 * i * this->degree
4576  + j]
4577  * this->shape_value_component
4578  (4 * i * this->degree + j,
4579  this->generalized_support_points[q_point
4580  + GeometryInfo<dim>::lines_per_cell
4581  * n_edge_points],
4582  1);
4583 
4584  for (unsigned int i = 0; i < this->degree; ++i)
4585  for (unsigned int j = 0; j < this->degree-1; ++j)
4586  system_rhs (i * (this->degree-1) + j)
4587  += boundary_weights
4588  (q_point + n_edge_points,
4589  2 * (i * (this->degree-1) + j)) * tmp;
4590  }
4591 
4592  system_matrix_inv.vmult (solution, system_rhs);
4593 
4594  // Add the computed values
4595  // to the resulting vector
4596  // only, if they are not
4597  // too small.
4598  for (unsigned int i = 0; i < this->degree; ++i)
4599  for (unsigned int j = 0; j < this->degree-1; ++j)
4600  if (std::abs (solution (i * (this->degree-1) + j))
4601  > 1e-14)
4602  local_dofs[(i
4603  + GeometryInfo<dim>::lines_per_cell)
4604  * (this->degree-1) + j
4605  + GeometryInfo<dim>::lines_per_cell]
4606  = solution (i * (this->degree-1) + j);
4607  }
4608 
4609  // Set up the right hand side
4610  // for the vertical shape
4611  // functions.
4612  system_rhs = 0;
4613 
4614  for (unsigned int q_point = 0;
4615  q_point < n_face_points; ++q_point)
4616  {
4617  double tmp
4618  = values[q_point
4619  + GeometryInfo<dim>::lines_per_cell
4620  * n_edge_points] (2);
4621 
4622  for (unsigned int i = 0; i < 2; ++i)
4623  for (unsigned int j = 0; j < this->degree; ++j)
4624  tmp -= local_dofs[2 * (i + 4)
4625  * this->degree + j]
4626  * this->shape_value_component
4627  (2 * (i + 4) * this->degree + j,
4628  this->generalized_support_points[q_point
4629  + GeometryInfo<dim>::lines_per_cell
4630  * n_edge_points],
4631  2);
4632 
4633  for (unsigned int i = 0; i < this->degree; ++i)
4634  for (unsigned int j = 0; j < this->degree-1; ++j)
4635  system_rhs (i * (this->degree-1) + j)
4636  += boundary_weights
4637  (q_point + n_edge_points,
4638  2 * (i * (this->degree-1) + j) + 1)
4639  * tmp;
4640  }
4641 
4642  system_matrix_inv.vmult (solution, system_rhs);
4643 
4644  // Add the computed values
4645  // to the resulting vector
4646  // only, if they are not
4647  // too small.
4648  for (unsigned int i = 0; i < this->degree; ++i)
4649  for (unsigned int j = 0; j < this->degree-1; ++j)
4650  if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
4651  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
4652  + this->degree-1)
4653  * this->degree]
4654  = solution (i * (this->degree-1) + j);
4655 
4656  break;
4657  }
4658 
4659  case 1:
4660  {
4661  if (offset + 1 < dim)
4662  {
4663  // Set up the right hand side
4664  // for the horizontal shape
4665  // functions.
4666  system_rhs = 0;
4667 
4668  for (unsigned int q_point = 0;
4669  q_point < n_face_points; ++q_point)
4670  {
4671  double tmp
4672  = values[q_point
4673  + GeometryInfo<dim>::lines_per_cell
4674  * n_edge_points
4675  + n_face_points] (1);
4676 
4677  for (unsigned int i = 0; i < 2; ++i)
4678  for (unsigned int j = 0; j <= deg; ++j)
4679  tmp -= local_dofs[(4 * i + 1)
4680  * this->degree + j]
4681  * this->shape_value_component
4682  ((4 * i + 1) * this->degree
4683  + j,
4684  this->generalized_support_points[q_point
4685  + GeometryInfo<dim>::lines_per_cell
4686  * n_edge_points
4687  + n_face_points],
4688  1);
4689 
4690  for (unsigned int i = 0; i <= deg; ++i)
4691  for (unsigned int j = 0; j < deg; ++j)
4692  system_rhs (i * deg + j)
4693  += boundary_weights
4694  (q_point + n_edge_points,
4695  2 * (i * deg + j)) * tmp;
4696  }
4697 
4698  system_matrix_inv.vmult (solution, system_rhs);
4699 
4700  // Add the computed values
4701  // to the resulting vector
4702  // only, if they are not
4703  // too small.
4704  for (unsigned int i = 0; i <= deg; ++i)
4705  for (unsigned int j = 0; j < deg; ++j)
4706  if (std::abs (solution (i * deg + j))
4707  > 1e-14)
4708  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
4709  + 2 * this->degree) * deg + j
4710  + GeometryInfo<dim>::lines_per_cell]
4711  = solution (i * deg + j);
4712  }
4713 
4714  // Set up the right hand side
4715  // for the vertical shape
4716  // functions.
4717  system_rhs = 0;
4718 
4719  for (unsigned int q_point = 0;
4720  q_point < n_face_points; ++q_point)
4721  {
4722  double tmp
4723  = values[q_point
4724  + GeometryInfo<dim>::lines_per_cell
4725  * n_edge_points + n_face_points]
4726  (2);
4727 
4728  for (unsigned int i = 0; i < 2; ++i)
4729  for (unsigned int j = 0; j <= deg; ++j)
4730  tmp -= local_dofs[(2 * (i + 4) + 1)
4731  * this->degree + j]
4732  * this->shape_value_component
4733  ((2 * (i + 4) + 1) * this->degree
4734  + j,
4735  this->generalized_support_points[q_point
4736  + GeometryInfo<dim>::lines_per_cell
4737  * n_edge_points
4738  + n_face_points],
4739  2);
4740 
4741  for (unsigned int i = 0; i <= deg; ++i)
4742  for (unsigned int j = 0; j < deg; ++j)
4743  system_rhs (i * deg + j)
4744  += boundary_weights
4745  (q_point + n_edge_points,
4746  2 * (i * deg + j) + 1) * tmp;
4747  }
4748 
4749  system_matrix_inv.vmult (solution, system_rhs);
4750 
4751  // Add the computed values
4752  // to the resulting vector
4753  // only, if they are not
4754  // too small.
4755  for (unsigned int i = 0; i <= deg; ++i)
4756  for (unsigned int j = 0; j < deg; ++j)
4757  if (std::abs (solution (i * deg + j)) > 1e-14)
4758  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
4759  + 3 * deg)
4760  * this->degree]
4761  = solution (i * deg + j);
4762 
4763  break;
4764  }
4765 
4766  case 2:
4767  {
4768  if (offset == 0)
4769  {
4770  // Set up the right hand side
4771  // for the horizontal shape
4772  // functions.
4773  system_rhs = 0;
4774 
4775  for (unsigned int q_point = 0;
4776  q_point < n_face_points; ++q_point)
4777  {
4778  double tmp
4779  = values[q_point
4780  + GeometryInfo<dim>::lines_per_cell
4781  * n_edge_points + 2 * n_face_points]
4782  (2);
4783 
4784  for (unsigned int i = 0; i < 2; ++i)
4785  for (unsigned int j = 0; j <= deg; ++j)
4786  tmp -= local_dofs[(i + 8) * this->degree
4787  + j]
4788  * this->shape_value_component
4789  ((i + 8) * this->degree + j,
4790  this->generalized_support_points[q_point
4791  + GeometryInfo<dim>::lines_per_cell
4792  * n_edge_points
4793  + 2
4794  * n_face_points],
4795  2);
4796 
4797  for (unsigned int i = 0; i <= deg; ++i)
4798  for (unsigned int j = 0; j < deg; ++j)
4799  system_rhs (i * deg + j)
4800  += boundary_weights
4801  (q_point + n_edge_points,
4802  2 * (i * deg + j)) * tmp;
4803  }
4804 
4805  system_matrix_inv.vmult (solution, system_rhs);
4806 
4807  // Add the computed values
4808  // to the resulting vector
4809  // only, if they are not
4810  // too small.
4811  for (unsigned int i = 0; i <= deg; ++i)
4812  for (unsigned int j = 0; j < deg; ++j)
4813  if (std::abs (solution (i * deg + j))
4814  > 1e-14)
4815  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
4816  + 4 * this->degree) * deg
4817  + j
4818  + GeometryInfo<dim>::lines_per_cell]
4819  = solution (i * deg + j);
4820  }
4821 
4822  // Set up the right hand side
4823  // for the vertical shape
4824  // functions.
4825  system_rhs = 0;
4826 
4827  for (unsigned int q_point = 0;
4828  q_point < n_face_points; ++q_point)
4829  {
4830  double tmp
4831  = values[q_point
4832  + GeometryInfo<dim>::lines_per_cell
4833  * n_edge_points
4834  + 2 * n_face_points] (0);
4835 
4836  for (unsigned int i = 0; i < 2; ++i)
4837  for (unsigned int j = 0; j <= deg; ++j)
4838  tmp -= local_dofs[(4 * i + 2)
4839  * this->degree + j]
4840  * this->shape_value_component
4841  ((4 * i + 2) * this->degree
4842  + j,
4843  this->generalized_support_points[q_point
4844  + GeometryInfo<dim>::lines_per_cell
4845  * n_edge_points
4846  + 2
4847  * n_face_points],
4848  0);
4849 
4850  for (unsigned int i = 0; i <= deg; ++i)
4851  for (unsigned int j = 0; j < deg; ++j)
4852  system_rhs (i * deg + j)
4853  += boundary_weights
4854  (q_point + n_edge_points,
4855  2 * (i * deg + j) + 1) * tmp;
4856  }
4857 
4858  system_matrix_inv.vmult (solution, system_rhs);
4859 
4860  // Add the computed values
4861  // to the resulting vector
4862  // only, if they are not
4863  // too small.
4864  for (unsigned int i = 0; i <= deg; ++i)
4865  for (unsigned int j = 0; j < deg; ++j)
4866  if (std::abs (solution (i * deg + j)) > 1e-14)
4867  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
4868  + 5 * deg) * this->degree]
4869  = solution (i * deg + j);
4870 
4871  break;
4872  }
4873 
4874  case 3:
4875  {
4876  if (offset == 0)
4877  {
4878  // Set up the right hand side
4879  // for the horizontal shape
4880  // functions.
4881  system_rhs = 0;
4882 
4883  for (unsigned int q_point = 0;
4884  q_point < n_face_points; ++q_point)
4885  {
4886  double tmp
4887  = values[q_point
4888  + GeometryInfo<dim>::lines_per_cell
4889  * n_edge_points + 3 * n_face_points]
4890  (2);
4891 
4892  for (unsigned int i = 0; i < 2; ++i)
4893  for (unsigned int j = 0; j <= deg; ++j)
4894  tmp -= local_dofs[(i + 10) * this->degree
4895  + j]
4896  * this->shape_value_component
4897  ((i + 10) * this->degree + j,
4898  this->generalized_support_points[q_point
4899  + GeometryInfo<dim>::lines_per_cell
4900  * n_edge_points
4901  + 3
4902  * n_face_points],
4903  2);
4904 
4905  for (unsigned int i = 0; i <= deg; ++i)
4906  for (unsigned int j = 0; j < deg; ++j)
4907  system_rhs (i * deg + j)
4908  += boundary_weights
4909  (q_point + n_edge_points,
4910  2 * (i * deg + j)) * tmp;
4911  }
4912 
4913  system_matrix_inv.vmult (solution, system_rhs);
4914 
4915  // Add the computed values
4916  // to the resulting vector
4917  // only, if they are not
4918  // too small.
4919  for (unsigned int i = 0; i <= deg; ++i)
4920  for (unsigned int j = 0; j < deg; ++j)
4921  if (std::abs (solution (i * deg + j))
4922  > 1e-14)
4923  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
4924  + 6 * this->degree) * deg + j
4925  + GeometryInfo<dim>::lines_per_cell]
4926  = solution (i * deg + j);
4927  }
4928 
4929  // Set up the right hand side
4930  // for the vertical shape
4931  // functions.
4932  system_rhs = 0;
4933 
4934  for (unsigned int q_point = 0;
4935  q_point < n_face_points; ++q_point)
4936  {
4937  double tmp
4938  = values[q_point
4939  + GeometryInfo<dim>::lines_per_cell
4940  * n_edge_points + 3
4941  * n_face_points] (0);
4942 
4943  for (unsigned int i = 0; i < 2; ++i)
4944  for (unsigned int j = 0; j <= deg; ++j)
4945  tmp -= local_dofs[(4 * i + 3)
4946  * this->degree + j]
4947  * this->shape_value_component
4948  ((4 * i + 3) * this->degree
4949  + j,
4950  this->generalized_support_points[q_point
4951  + GeometryInfo<dim>::lines_per_cell
4952  * n_edge_points
4953  + 3
4954  * n_face_points],
4955  0);
4956 
4957  for (unsigned int i = 0; i <= deg; ++i)
4958  for (unsigned int j = 0; j < deg; ++j)
4959  system_rhs (i * deg + j)
4960  += boundary_weights
4961  (q_point + n_edge_points,
4962  2 * (i * deg + j) + 1) * tmp;
4963  }
4964 
4965  system_matrix_inv.vmult (solution, system_rhs);
4966 
4967  // Add the computed values
4968  // to the resulting vector
4969  // only, if they are not
4970  // too small.
4971  for (unsigned int i = 0; i <= deg; ++i)
4972  for (unsigned int j = 0; j < deg; ++j)
4973  if (std::abs (solution (i * deg + j)) > 1e-14)
4974  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
4975  + 7 * deg) * this->degree]
4976  = solution (i * deg + j);
4977 
4978  break;
4979  }
4980 
4981  case 4:
4982  {
4983  if (offset + 1 < dim)
4984  {
4985  // Set up the right hand side
4986  // for the horizontal shape
4987  // functions.
4988  if (offset == 0)
4989  {
4990  system_rhs = 0;
4991 
4992  for (unsigned int q_point = 0;
4993  q_point < n_face_points; ++q_point)
4994  {
4995  double tmp
4996  = values[q_point
4997  + GeometryInfo<dim>::lines_per_cell
4998  * n_edge_points + 4
4999  * n_face_points] (0);
5000 
5001  for (unsigned int i = 0; i < 2; ++i)
5002  for (unsigned int j = 0; j <= deg; ++j)
5003  tmp -= local_dofs[(i + 2)
5004  * this->degree
5005  + j]
5006  * this->shape_value_component
5007  ((i + 2) * this->degree
5008  + j,
5009  this->generalized_support_points[q_point
5010  + GeometryInfo<dim>::lines_per_cell
5011  * n_edge_points
5012  + 4
5013  * n_face_points],
5014  0);
5015 
5016  for (unsigned int i = 0; i <= deg; ++i)
5017  for (unsigned int j = 0; j < deg; ++j)
5018  system_rhs (i * deg + j)
5019  += boundary_weights
5020  (q_point + n_edge_points,
5021  2 * (i * deg + j)) * tmp;
5022  }
5023 
5024  system_matrix_inv.vmult
5025  (solution, system_rhs);
5026 
5027  // Add the computed values
5028  // to the resulting vector
5029  // only, if they are not
5030  // too small.
5031  for (unsigned int i = 0; i <= deg; ++i)
5032  for (unsigned int j = 0; j < deg; ++j)
5033  if (std::abs (solution (i * deg + j))
5034  > 1e-14)
5035  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
5036  + 8 * this->degree) * deg
5037  + j
5038  + GeometryInfo<dim>::lines_per_cell]
5039  = solution (i * deg + j);
5040  }
5041 
5042  // Set up the right hand side
5043  // for the vertical shape
5044  // functions.
5045  system_rhs = 0;
5046 
5047  for (unsigned int q_point = 0;
5048  q_point < n_face_points; ++q_point)
5049  {
5050  double tmp
5051  = values[q_point
5052  + GeometryInfo<dim>::lines_per_cell
5053  * n_edge_points + 4
5054  * n_face_points] (1);
5055 
5056  for (unsigned int i = 0; i < 2; ++i)
5057  for (unsigned int j = 0; j <= deg; ++j)
5058  tmp -= local_dofs[i * this->degree + j]
5059  * this->shape_value_component
5060  (i * this->degree + j,
5061  this->generalized_support_points[q_point
5062  + GeometryInfo<dim>::lines_per_cell
5063  * n_edge_points
5064  + 4
5065  * n_face_points],
5066  1);
5067 
5068  for (unsigned int i = 0; i <= deg; ++i)
5069  for (unsigned int j = 0; j < deg; ++j)
5070  system_rhs (i * deg + j)
5071  += boundary_weights
5072  (q_point + n_edge_points,
5073  2 * (i * deg + j) + 1) * tmp;
5074  }
5075 
5076  system_matrix_inv.vmult (solution, system_rhs);
5077 
5078  // Add the computed values
5079  // to the resulting vector
5080  // only, if they are not
5081  // too small.
5082  for (unsigned int i = 0; i <= deg; ++i)
5083  for (unsigned int j = 0; j < deg; ++j)
5084  if (std::abs (solution (i * deg + j))
5085  > 1e-14)
5086  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
5087  + 9 * deg)
5088  * this->degree]
5089  = solution (i * deg + j);
5090  }
5091 
5092  break;
5093  }
5094 
5095  default:
5096  if (offset + 1 < dim)
5097  {
5098  // Set up the right hand side
5099  // for the horizontal shape
5100  // functions.
5101  if (offset == 0)
5102  {
5103  system_rhs = 0;
5104 
5105  for (unsigned int q_point = 0;
5106  q_point < n_face_points; ++q_point)
5107  {
5108  double tmp
5109  = values[q_point
5110  + GeometryInfo<dim>::lines_per_cell
5111  * n_edge_points
5112  + 5 * n_face_points] (0);
5113 
5114  for (unsigned int i = 0; i < 2; ++i)
5115  for (unsigned int j = 0; j <= deg; ++j)
5116  tmp -= local_dofs[(i + 6)
5117  * this->degree + j]
5118  * this->shape_value_component
5119  ((i + 6) * this->degree + j,
5120  this->generalized_support_points[q_point
5121  + GeometryInfo<dim>::lines_per_cell
5122  * n_edge_points
5123  + 5
5124  * n_face_points],
5125  0);
5126 
5127  for (unsigned int i = 0; i <= deg; ++i)
5128  for (unsigned int j = 0; j < deg; ++j)
5129  system_rhs (i * deg + j)
5130  += boundary_weights
5131  (q_point + n_edge_points,
5132  2 * (i * deg + j)) * tmp;
5133  }
5134 
5135  system_matrix_inv.vmult
5136  (solution, system_rhs);
5137 
5138  // Add the computed values
5139  // to the resulting vector
5140  // only, if they are not
5141  // too small.
5142  for (unsigned int i = 0; i <= deg; ++i)
5143  for (unsigned int j = 0; j < deg; ++j)
5144  if (std::abs (solution (i * deg + j))
5145  > 1e-14)
5146  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
5147  + 10 * this->degree)
5148  * deg + j
5149  + GeometryInfo<dim>::lines_per_cell]
5150  = solution (i * deg + j);
5151  }
5152 
5153  // Set up the right hand side
5154  // for the vertical shape
5155  // functions.
5156  system_rhs = 0;
5157 
5158  for (unsigned int q_point = 0;
5159  q_point < n_face_points; ++q_point)
5160  {
5161  double tmp
5162  = values[q_point
5163  + GeometryInfo<dim>::lines_per_cell
5164  * n_edge_points + 5
5165  * n_face_points] (1);
5166 
5167  for (unsigned int i = 0; i < 2; ++i)
5168  for (unsigned int j = 0; j <= deg; ++j)
5169  tmp -= local_dofs[(i + 4)
5170  * this->degree + j]
5171  * this->shape_value_component
5172  ((i + 4) * this->degree + j,
5173  this->generalized_support_points[q_point
5174  + GeometryInfo<dim>::lines_per_cell
5175  * n_edge_points
5176  + 5
5177  * n_face_points],
5178  1);
5179 
5180  for (unsigned int i = 0; i <= deg; ++i)
5181  for (unsigned int j = 0; j < deg; ++j)
5182  system_rhs (i * deg + j)
5183  += boundary_weights
5184  (q_point + n_edge_points,
5185  2 * (i * deg + j) + 1) * tmp;
5186  }
5187 
5188  system_matrix_inv.vmult (solution, system_rhs);
5189 
5190  // Add the computed values
5191  // to the resulting vector
5192  // only, if they are not
5193  // too small.
5194  for (unsigned int i = 0; i <= deg; ++i)
5195  for (unsigned int j = 0; j < deg; ++j)
5196  if (std::abs (solution (i * deg + j))
5197  > 1e-14)
5198  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
5199  + 11 * deg) * this->degree]
5200  = solution (i * deg + j);
5201  }
5202  }
5203  }
5204 
5205  // Finally we project
5206  // the remaining parts
5207  // of the function on
5208  // the interior shape
5209  // functions.
5210  const QGauss<dim> reference_quadrature (this->degree);
5211  const unsigned int &
5212  n_interior_points = reference_quadrature.size ();
5213 
5214  // We create the
5215  // system matrix.
5216  system_matrix.reinit (this->degree * deg * deg,
5217  this->degree * deg * deg);
5218  system_matrix = 0;
5219 
5220  for (unsigned int i = 0; i <= deg; ++i)
5221  for (unsigned int j = 0; j < deg; ++j)
5222  for (unsigned int k = 0; k < deg; ++k)
5223  for (unsigned int l = 0; l <= deg; ++l)
5224  for (unsigned int m = 0; m < deg; ++m)
5225  for (unsigned int n = 0; n < deg; ++n)
5226  for (unsigned int q_point = 0;
5227  q_point < n_interior_points; ++q_point)
5228  system_matrix ((i * deg + j) * deg + k,
5229  (l * deg + m) * deg + n)
5230  += reference_quadrature.weight (q_point)
5231  * legendre_polynomials[i].value
5232  (this->generalized_support_points[q_point
5233  + GeometryInfo<dim>::lines_per_cell
5234  * n_edge_points
5235  + GeometryInfo<dim>::faces_per_cell
5236  * n_face_points]
5237  (0)) * lobatto_polynomials[j + 2].value
5238  (this->generalized_support_points[q_point
5239  + GeometryInfo<dim>::lines_per_cell
5240  * n_edge_points
5241  + GeometryInfo<dim>::faces_per_cell
5242  * n_face_points]
5243  (1))
5244  * lobatto_polynomials[k + 2].value
5245  (this->generalized_support_points[q_point
5246  + GeometryInfo<dim>::lines_per_cell
5247  * n_edge_points
5248  + GeometryInfo<dim>::faces_per_cell
5249  * n_face_points]
5250  (2))
5251  * lobatto_polynomials_grad[l].value
5252  (this->generalized_support_points[q_point
5253  + GeometryInfo<dim>::lines_per_cell
5254  * n_edge_points
5255  + GeometryInfo<dim>::faces_per_cell
5256  * n_face_points]
5257  (0))
5258  * lobatto_polynomials[m + 2].value
5259  (this->generalized_support_points[q_point
5260  + GeometryInfo<dim>::lines_per_cell
5261  * n_edge_points
5262  + GeometryInfo<dim>::faces_per_cell
5263  * n_face_points]
5264  (1))
5265  * lobatto_polynomials[n + 2].value
5266  (this->generalized_support_points[q_point
5267  + GeometryInfo<dim>::lines_per_cell
5268  * n_edge_points
5269  + GeometryInfo<dim>::faces_per_cell
5270  * n_face_points]
5271  (2));
5272 
5273  system_matrix_inv.reinit (system_matrix.m (),
5274  system_matrix.m ());
5275  system_matrix_inv.invert (system_matrix);
5276  system_rhs.reinit (system_matrix_inv.m ());
5277  solution.reinit (system_matrix.m ());
5278 
5279  if (offset + 1 < dim)
5280  {
5281  if (offset == 0)
5282  {
5283  // Set up the right hand side.
5284  system_rhs = 0;
5285 
5286  for (unsigned int q_point = 0;
5287  q_point < n_interior_points; ++q_point)
5288  {
5289  double tmp
5290  = values[q_point
5291  + GeometryInfo<dim>::lines_per_cell
5292  * n_edge_points
5293  + GeometryInfo<dim>::faces_per_cell
5294  * n_face_points] (0);
5295 
5296  for (unsigned int i = 0; i <= deg; ++i)
5297  {
5298  for (unsigned int j = 0; j < 2; ++j)
5299  for (unsigned int k = 0; k < 2; ++k)
5300  tmp -= local_dofs[i + (j + 4 * k + 2)
5301  * this->degree]
5302  * this->shape_value_component
5303  (i + (j + 4 * k + 2)
5304  * this->degree,
5305  this->generalized_support_points[q_point
5306  + GeometryInfo<dim>::lines_per_cell
5307  * n_edge_points
5308  + GeometryInfo<dim>::faces_per_cell
5309  * n_face_points],
5310  0);
5311 
5312  for (unsigned int j = 0; j < deg; ++j)
5313  for (unsigned int k = 0; k < 4; ++k)
5314  tmp -= local_dofs[(i + 2 * (k + 2)
5315  * this->degree
5316  + GeometryInfo<dim>::lines_per_cell)
5317  * deg + j
5318  + GeometryInfo<dim>::lines_per_cell]
5319  * this->shape_value_component
5320  ((i + 2 * (k + 2) * this->degree
5321  + GeometryInfo<dim>::lines_per_cell)
5322  * deg + j
5323  + GeometryInfo<dim>::lines_per_cell,
5324  this->generalized_support_points[q_point
5325  + GeometryInfo<dim>::lines_per_cell
5326  * n_edge_points
5327  + GeometryInfo<dim>::faces_per_cell
5328  * n_face_points],
5329  0);
5330  }
5331 
5332  for (unsigned int i = 0; i <= deg; ++i)
5333  for (unsigned int j = 0; j < deg; ++j)
5334  for (unsigned int k = 0; k < deg; ++k)
5335  system_rhs ((i * deg + j) * deg + k)
5336  += reference_quadrature.weight (q_point)
5337  * tmp
5338  * lobatto_polynomials_grad[i].value
5339  (this->generalized_support_points[q_point
5340  + GeometryInfo<dim>::lines_per_cell
5341  * n_edge_points
5342  + GeometryInfo<dim>::faces_per_cell
5343  * n_face_points]
5344  (0))
5345  * lobatto_polynomials[j + 2].value
5346  (this->generalized_support_points[q_point
5347  + GeometryInfo<dim>::lines_per_cell
5348  * n_edge_points
5349  + GeometryInfo<dim>::faces_per_cell
5350  * n_face_points]
5351  (1))
5352  * lobatto_polynomials[k + 2].value
5353  (this->generalized_support_points[q_point
5354  + GeometryInfo<dim>::lines_per_cell
5355  * n_edge_points
5356  + GeometryInfo<dim>::faces_per_cell
5357  * n_face_points]
5358  (2));
5359  }
5360 
5361  system_matrix_inv.vmult (solution, system_rhs);
5362 
5363  // Add the computed values
5364  // to the resulting vector
5365  // only, if they are not
5366  // too small.
5367  for (unsigned int i = 0; i <= deg; ++i)
5368  for (unsigned int j = 0; j < deg; ++j)
5369  for (unsigned int k = 0; k < deg; ++k)
5370  if (std::abs (solution ((i * deg + j) * deg + k))
5371  > 1e-14)
5372  local_dofs[((i + 2
5373  * GeometryInfo<dim>::faces_per_cell)
5374  * deg + j
5375  + GeometryInfo<dim>::lines_per_cell
5376  + 2
5377  * GeometryInfo<dim>::faces_per_cell)
5378  * deg + k
5379  + GeometryInfo<dim>::lines_per_cell]
5380  = solution ((i * deg + j) * deg + k);
5381  }
5382 
5383  // Set up the right hand side.
5384  system_rhs = 0;
5385 
5386  for (unsigned int q_point = 0; q_point < n_interior_points;
5387  ++q_point)
5388  {
5389  double tmp
5390  = values[q_point + GeometryInfo<dim>::lines_per_cell
5391  * n_edge_points
5392  + GeometryInfo<dim>::faces_per_cell
5393  * n_face_points] (1);
5394 
5395  for (unsigned int i = 0; i <= deg; ++i)
5396  for (unsigned int j = 0; j < 2; ++j)
5397  {
5398  for (unsigned int k = 0; k < 2; ++k)
5399  tmp -= local_dofs[i + (4 * j + k)
5400  * this->degree]
5401  * this->shape_value_component
5402  (i + (4 * j + k) * this->degree,
5403  this->generalized_support_points[q_point
5404  + GeometryInfo<dim>::lines_per_cell
5405  * n_edge_points
5406  + GeometryInfo<dim>::faces_per_cell
5407  * n_face_points],
5408  1);
5409 
5410  for (unsigned int k = 0; k < deg; ++k)
5411  tmp -= local_dofs[(i + 2 * j * this->degree
5412  + GeometryInfo<dim>::lines_per_cell)
5413  * deg + k
5414  + GeometryInfo<dim>::lines_per_cell]
5415  * this->shape_value_component
5416  ((i + 2 * j * this->degree
5417  + GeometryInfo<dim>::lines_per_cell)
5418  * deg + k
5419  + GeometryInfo<dim>::lines_per_cell,
5420  this->generalized_support_points[q_point
5421  + GeometryInfo<dim>::lines_per_cell
5422  * n_edge_points
5423  + GeometryInfo<dim>::faces_per_cell
5424  * n_face_points],
5425  1)
5426  + local_dofs[i + ((2 * j + 9) * deg + k
5427  + GeometryInfo<dim>::lines_per_cell)
5428  * this->degree]
5429  * this->shape_value_component
5430  (i + ((2 * j + 9) * deg + k
5431  + GeometryInfo<dim>::lines_per_cell)
5432  * this->degree,
5433  this->generalized_support_points[q_point
5434  + GeometryInfo<dim>::lines_per_cell
5435  * n_edge_points
5436  + GeometryInfo<dim>::faces_per_cell
5437  * n_face_points],
5438  1);
5439  }
5440 
5441  for (unsigned int i = 0; i <= deg; ++i)
5442  for (unsigned int j = 0; j < deg; ++j)
5443  for (unsigned int k = 0; k < deg; ++k)
5444  system_rhs ((i * deg + j) * deg + k)
5445  += reference_quadrature.weight (q_point) * tmp
5446  * lobatto_polynomials_grad[i].value
5447  (this->generalized_support_points[q_point
5448  + GeometryInfo<dim>::lines_per_cell
5449  * n_edge_points
5450  + GeometryInfo<dim>::faces_per_cell
5451  * n_face_points]
5452  (1))
5453  * lobatto_polynomials[j + 2].value
5454  (this->generalized_support_points[q_point
5455  + GeometryInfo<dim>::lines_per_cell
5456  * n_edge_points
5457  + GeometryInfo<dim>::faces_per_cell
5458  * n_face_points]
5459  (0))
5460  * lobatto_polynomials[k + 2].value
5461  (this->generalized_support_points[q_point
5462  + GeometryInfo<dim>::lines_per_cell
5463  * n_edge_points
5464  + GeometryInfo<dim>::faces_per_cell
5465  * n_face_points]
5466  (2));
5467  }
5468 
5469  system_matrix_inv.vmult (solution, system_rhs);
5470 
5471  // Add the computed values
5472  // to the resulting vector
5473  // only, if they are not
5474  // too small.
5475  for (unsigned int i = 0; i <= deg; ++i)
5476  for (unsigned int j = 0; j < deg; ++j)
5477  for (unsigned int k = 0; k < deg; ++k)
5478  if (std::abs (solution ((i * deg + j) * deg + k))
5479  > 1e-14)
5480  local_dofs[((i + this->degree + 2
5481  * GeometryInfo<dim>::faces_per_cell)
5482  * deg + j
5483  + GeometryInfo<dim>::lines_per_cell + 2
5484  * GeometryInfo<dim>::faces_per_cell)
5485  * deg + k
5486  + GeometryInfo<dim>::lines_per_cell]
5487  = solution ((i * deg + j) * deg + k);
5488  }
5489 
5490  // Set up the right hand side.
5491  system_rhs = 0;
5492 
5493  for (unsigned int q_point = 0; q_point < n_interior_points;
5494  ++q_point)
5495  {
5496  double tmp
5497  = values[q_point + GeometryInfo<dim>::lines_per_cell
5498  * n_edge_points
5499  + GeometryInfo<dim>::faces_per_cell
5500  * n_face_points] (2);
5501 
5502  for (unsigned int i = 0; i <= deg; ++i)
5503  for (unsigned int j = 0; j < 4; ++j)
5504  {
5505  tmp -= local_dofs[i + (j + 8) * this->degree]
5506  * this->shape_value_component
5507  (i + (j + 8) * this->degree,
5508  this->generalized_support_points[q_point
5509  + GeometryInfo<dim>::lines_per_cell
5510  * n_edge_points
5511  + GeometryInfo<dim>::faces_per_cell
5512  * n_face_points],
5513  2);
5514 
5515  for (unsigned int k = 0; k < deg; ++k)
5516  tmp -= local_dofs[i + ((2 * j + 1) * deg + k
5517  + GeometryInfo<dim>::lines_per_cell)
5518  * this->degree]
5519  * this->shape_value_component
5520  (i + ((2 * j + 1) * deg + k
5521  + GeometryInfo<dim>::lines_per_cell)
5522  * this->degree,
5523  this->generalized_support_points[q_point
5524  + GeometryInfo<dim>::lines_per_cell
5525  * n_edge_points
5526  + GeometryInfo<dim>::faces_per_cell
5527  * n_face_points],
5528  2);
5529  }
5530 
5531  for (unsigned int i = 0; i <= deg; ++i)
5532  for (unsigned int j = 0; j < deg; ++j)
5533  for (unsigned int k = 0; k < deg; ++k)
5534  system_rhs ((i * deg + j) * deg + k)
5535  += reference_quadrature.weight (q_point) * tmp
5536  * lobatto_polynomials_grad[i].value
5537  (this->generalized_support_points[q_point
5538  + GeometryInfo<dim>::lines_per_cell
5539  * n_edge_points
5540  + GeometryInfo<dim>::faces_per_cell
5541  * n_face_points]
5542  (2))
5543  * lobatto_polynomials[j + 2].value
5544  (this->generalized_support_points[q_point
5545  + GeometryInfo<dim>::lines_per_cell
5546  * n_edge_points
5547  + GeometryInfo<dim>::faces_per_cell
5548  * n_face_points]
5549  (0))
5550  * lobatto_polynomials[k + 2].value
5551  (this->generalized_support_points[q_point
5552  + GeometryInfo<dim>::lines_per_cell
5553  * n_edge_points
5554  + GeometryInfo<dim>::faces_per_cell
5555  * n_face_points]
5556  (1));
5557  }
5558 
5559  system_matrix_inv.vmult (solution, system_rhs);
5560 
5561  // Add the computed values
5562  // to the resulting vector
5563  // only, if they are not
5564  // too small.
5565  for (unsigned int i = 0; i <= deg; ++i)
5566  for (unsigned int j = 0; j < deg; ++j)
5567  for (unsigned int k = 0; k < deg; ++k)
5568  if (std::abs (solution ((i * deg + j) * deg + k))
5569  > 1e-14)
5570  local_dofs[i + ((j + 2
5571  * (deg + GeometryInfo<dim>::faces_per_cell))
5572  * deg + k
5573  + GeometryInfo<dim>::lines_per_cell)
5574  * this->degree]
5575  = solution ((i * deg + j) * deg + k);
5576  }
5577 
5578  break;
5579  }
5580 
5581  default:
5582  Assert (false, ExcNotImplemented ());
5583  }
5584 }
5585 
5586 
5587 // Interpolate a function, which is given by
5588 // its values at the generalized support
5589 // points in the finite element space on the
5590 // reference cell.
5591 // This is done as usual by projection-based
5592 // interpolation.
5593 template <int dim>
5594 void
5595 FE_Nedelec<dim>::interpolate (std::vector<double> &local_dofs,
5596  const VectorSlice<const std::vector<std::vector<double> > > &values)
5597 const
5598 {
5599  const unsigned int deg = this->degree-1;
5600  Assert (values.size () == this->n_components (),
5601  ExcDimensionMismatch (values.size (), this->n_components ()));
5602  Assert (values[0].size () == this->generalized_support_points.size (),
5603  ExcDimensionMismatch (values[0].size (),
5604  this->generalized_support_points.size ()));
5605  Assert (local_dofs.size () == this->dofs_per_cell,
5606  ExcDimensionMismatch (local_dofs.size (), this->dofs_per_cell));
5607  std::fill (local_dofs.begin (), local_dofs.end (), 0.0);
5608 
5609  switch (dim)
5610  {
5611  case 2:
5612  {
5613  // Let us begin with the
5614  // interpolation part.
5615  const QGauss<dim - 1> reference_edge_quadrature (this->degree);
5616  const unsigned int &
5617  n_edge_points = reference_edge_quadrature.size ();
5618 
5619  for (unsigned int i = 0; i < 2; ++i)
5620  for (unsigned int j = 0; j < 2; ++j)
5621  {
5622  for (unsigned int q_point = 0; q_point < n_edge_points;
5623  ++q_point)
5624  local_dofs[(i + 2 * j) * this->degree]
5625  += reference_edge_quadrature.weight (q_point)
5626  * values[1 - j][q_point + (i + 2 * j) * n_edge_points];
5627 
5628  // Add the computed values
5629  // to the resulting vector
5630  // only, if they are not
5631  // too small.
5632  if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14)
5633  local_dofs[(i + 2 * j) * this->degree] = 0.0;
5634  }
5635 
5636  // If the degree is greater
5637  // than 0, then we have still
5638  // some higher order edge
5639  // shape functions to
5640  // consider.
5641  // Here the projection part
5642  // starts. The dof values
5643  // are obtained by solving
5644  // a linear system of
5645  // equations.
5646  if (this->degree-1 > 1)
5647  {
5648  // We start with projection
5649  // on the higher order edge
5650  // shape function.
5651  const std::vector<Polynomials::Polynomial<double> > &
5652  lobatto_polynomials
5653  = Polynomials::Lobatto::generate_complete_basis
5654  (this->degree);
5655  FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
5656  std::vector<Polynomials::Polynomial<double> >
5657  lobatto_polynomials_grad (this->degree);
5658 
5659  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
5660  ++i)
5661  lobatto_polynomials_grad[i]
5662  = lobatto_polynomials[i + 1].derivative ();
5663 
5664  // Set up the system matrix.
5665  // This can be used for all
5666  // edges.
5667  for (unsigned int i = 0; i < system_matrix.m (); ++i)
5668  for (unsigned int j = 0; j < system_matrix.n (); ++j)
5669  for (unsigned int q_point = 0; q_point < n_edge_points;
5670  ++q_point)
5671  system_matrix (i, j)
5672  += boundary_weights (q_point, j)
5673  * lobatto_polynomials_grad[i + 1].value
5674  (this->generalized_face_support_points[q_point]
5675  (1));
5676 
5677  FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
5678 
5679  system_matrix_inv.invert (system_matrix);
5680 
5681  const unsigned int
5682  line_coordinate[GeometryInfo<2>::lines_per_cell]
5683  = {1, 1, 0, 0};
5684  Vector<double> system_rhs (system_matrix.m ());
5685  Vector<double> solution (system_rhs.size ());
5686 
5687  for (unsigned int line = 0;
5688  line < GeometryInfo<dim>::lines_per_cell; ++line)
5689  {
5690  // Set up the right hand side.
5691  system_rhs = 0;
5692 
5693  for (unsigned int q_point = 0; q_point < n_edge_points;
5694  ++q_point)
5695  {
5696  const double tmp
5697  = values[line_coordinate[line]][line * n_edge_points
5698  + q_point]
5699  - local_dofs[line * this->degree]
5700  * this->shape_value_component
5701  (line * this->degree,
5702  this->generalized_support_points[line
5703  * n_edge_points
5704  + q_point],
5705  line_coordinate[line]);
5706 
5707  for (unsigned int i = 0; i < system_rhs.size (); ++i)
5708  system_rhs (i) += boundary_weights (q_point, i) * tmp;
5709  }
5710 
5711  system_matrix_inv.vmult (solution, system_rhs);
5712 
5713  // Add the computed values
5714  // to the resulting vector
5715  // only, if they are not
5716  // too small.
5717  for (unsigned int i = 0; i < solution.size (); ++i)
5718  if (std::abs (solution (i)) > 1e-14)
5719  local_dofs[line * this->degree + i + 1] = solution (i);
5720  }
5721 
5722  // Then we go on to the
5723  // interior shape
5724  // functions. Again we
5725  // set up the system
5726  // matrix and use it
5727  // for both, the
5728  // horizontal and the
5729  // vertical, interior
5730  // shape functions.
5731  const QGauss<dim> reference_quadrature (this->degree);
5732  const unsigned int &
5733  n_interior_points = reference_quadrature.size ();
5734  const std::vector<Polynomials::Polynomial<double> > &
5735  legendre_polynomials
5736  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
5737 
5738  system_matrix.reinit ((this->degree-1) * this->degree,
5739  (this->degree-1) * this->degree);
5740  system_matrix = 0;
5741 
5742  for (unsigned int i = 0; i < this->degree; ++i)
5743  for (unsigned int j = 0; j < this->degree-1; ++j)
5744  for (unsigned int k = 0; k < this->degree; ++k)
5745  for (unsigned int l = 0; l < this->degree-1; ++l)
5746  for (unsigned int q_point = 0;
5747  q_point < n_interior_points; ++q_point)
5748  system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
5749  += reference_quadrature.weight (q_point)
5750  * legendre_polynomials[i].value
5751  (this->generalized_support_points[q_point
5752  + GeometryInfo<dim>::lines_per_cell
5753  * n_edge_points]
5754  (0))
5755  * lobatto_polynomials[j + 2].value
5756  (this->generalized_support_points[q_point
5757  + GeometryInfo<dim>::lines_per_cell
5758  * n_edge_points]
5759  (1))
5760  * lobatto_polynomials_grad[k].value
5761  (this->generalized_support_points[q_point
5762  + GeometryInfo<dim>::lines_per_cell
5763  * n_edge_points]
5764  (0))
5765  * lobatto_polynomials[l + 2].value
5766  (this->generalized_support_points[q_point
5767  + GeometryInfo<dim>::lines_per_cell
5768  * n_edge_points]
5769  (1));
5770 
5771  system_matrix_inv.reinit (system_matrix.m (),
5772  system_matrix.m ());
5773  system_matrix_inv.invert (system_matrix);
5774  // Set up the right hand side
5775  // for the horizontal shape
5776  // functions.
5777  system_rhs.reinit (system_matrix_inv.m ());
5778  system_rhs = 0;
5779 
5780  for (unsigned int q_point = 0; q_point < n_interior_points;
5781  ++q_point)
5782  {
5783  double tmp
5784  = values[0][q_point + GeometryInfo<dim>::lines_per_cell
5785  * n_edge_points];
5786 
5787  for (unsigned int i = 0; i < 2; ++i)
5788  for (unsigned int j = 0; j <= deg; ++j)
5789  tmp -= local_dofs[(i + 2) * this->degree + j]
5790  * this->shape_value_component
5791  ((i + 2) * this->degree + j,
5792  this->generalized_support_points[q_point
5793  + GeometryInfo<dim>::lines_per_cell
5794  * n_edge_points],
5795  0);
5796 
5797  for (unsigned int i = 0; i <= deg; ++i)
5798  for (unsigned int j = 0; j < deg; ++j)
5799  system_rhs (i * deg + j)
5800  += reference_quadrature.weight (q_point) * tmp
5801  * lobatto_polynomials_grad[i].value
5802  (this->generalized_support_points[q_point
5803  + GeometryInfo<dim>::lines_per_cell
5804  * n_edge_points]
5805  (0))
5806  * lobatto_polynomials[j + 2].value
5807  (this->generalized_support_points[q_point
5808  + GeometryInfo<dim>::lines_per_cell
5809  * n_edge_points]
5810  (1));
5811  }
5812 
5813  solution.reinit (system_matrix.m ());
5814  system_matrix_inv.vmult (solution, system_rhs);
5815 
5816  // Add the computed values
5817  // to the resulting vector
5818  // only, if they are not
5819  // too small.
5820  for (unsigned int i = 0; i <= deg; ++i)
5821  for (unsigned int j = 0; j < deg; ++j)
5822  if (std::abs (solution (i * deg + j)) > 1e-14)
5823  local_dofs[(i + GeometryInfo<dim>::lines_per_cell) * deg
5824  + j + GeometryInfo<dim>::lines_per_cell]
5825  = solution (i * deg + j);
5826 
5827  system_rhs = 0;
5828  // Set up the right hand side
5829  // for the vertical shape
5830  // functions.
5831 
5832  for (unsigned int q_point = 0; q_point < n_interior_points;
5833  ++q_point)
5834  {
5835  double tmp
5836  = values[1][q_point + GeometryInfo<dim>::lines_per_cell
5837  * n_edge_points];
5838 
5839  for (unsigned int i = 0; i < 2; ++i)
5840  for (unsigned int j = 0; j <= deg; ++j)
5841  tmp -= local_dofs[i * this->degree + j]
5842  * this->shape_value_component
5843  (i * this->degree + j,
5844  this->generalized_support_points[q_point
5845  + GeometryInfo<dim>::lines_per_cell
5846  * n_edge_points],
5847  1);
5848 
5849  for (unsigned int i = 0; i <= deg; ++i)
5850  for (unsigned int j = 0; j < deg; ++j)
5851  system_rhs (i * deg + j)
5852  += reference_quadrature.weight (q_point) * tmp
5853  * lobatto_polynomials_grad[i].value
5854  (this->generalized_support_points[q_point
5855  + GeometryInfo<dim>::lines_per_cell
5856  * n_edge_points]
5857  (1))
5858  * lobatto_polynomials[j + 2].value
5859  (this->generalized_support_points[q_point
5860  + GeometryInfo<dim>::lines_per_cell
5861  * n_edge_points]
5862  (0));
5863  }
5864 
5865  system_matrix_inv.vmult (solution, system_rhs);
5866 
5867  // Add the computed values
5868  // to the resulting vector
5869  // only, if they are not
5870  // too small.
5871  for (unsigned int i = 0; i <= deg; ++i)
5872  for (unsigned int j = 0; j < deg; ++j)
5873  if (std::abs (solution (i * deg + j)) > 1e-14)
5874  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
5875  + deg) * this->degree]
5876  = solution (i * deg + j);
5877  }
5878 
5879  break;
5880  }
5881 
5882  case 3:
5883  {
5884  // Let us begin with the
5885  // interpolation part.
5886  const QGauss<1> reference_edge_quadrature (this->degree);
5887  const unsigned int &
5888  n_edge_points = reference_edge_quadrature.size ();
5889 
5890  for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
5891  {
5892  for (unsigned int i = 0; i < 4; ++i)
5893  local_dofs[(i + 8) * this->degree]
5894  += reference_edge_quadrature.weight (q_point)
5895  * values[2][q_point + (i + 8) * n_edge_points];
5896 
5897  for (unsigned int i = 0; i < 2; ++i)
5898  for (unsigned int j = 0; j < 2; ++j)
5899  for (unsigned int k = 0; k < 2; ++k)
5900  local_dofs[(i + 2 * (2 * j + k)) * this->degree]
5901  += reference_edge_quadrature.weight (q_point)
5902  * values[1 - k][q_point + (i + 2 * (2 * j + k))
5903  * n_edge_points];
5904  }
5905 
5906  // Add the computed values
5907  // to the resulting vector
5908  // only, if they are not
5909  // too small.
5910  for (unsigned int i = 0; i < 4; ++i)
5911  if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
5912  local_dofs[(i + 8) * this->degree] = 0.0;
5913 
5914  for (unsigned int i = 0; i < 2; ++i)
5915  for (unsigned int j = 0; j < 2; ++j)
5916  for (unsigned int k = 0; k < 2; ++k)
5917  if (std::abs (local_dofs[(i + 2 * (2 * j + k)) * this->degree])
5918  < 1e-14)
5919  local_dofs[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
5920 
5921  // If the degree is greater
5922  // than 0, then we have still
5923  // some higher order shape
5924  // functions to consider.
5925  // Here the projection part
5926  // starts. The dof values
5927  // are obtained by solving
5928  // a linear system of
5929  // equations.
5930  if (this->degree > 1)
5931  {
5932  // We start with projection
5933  // on the higher order edge
5934  // shape function.
5935  const std::vector<Polynomials::Polynomial<double> > &
5936  lobatto_polynomials
5937  = Polynomials::Lobatto::generate_complete_basis
5938  (this->degree);
5939  FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
5940  std::vector<Polynomials::Polynomial<double> >
5941  lobatto_polynomials_grad (this->degree);
5942 
5943  for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
5944  ++i)
5945  lobatto_polynomials_grad[i]
5946  = lobatto_polynomials[i + 1].derivative ();
5947 
5948  // Set up the system matrix.
5949  // This can be used for all
5950  // edges.
5951  for (unsigned int i = 0; i < system_matrix.m (); ++i)
5952  for (unsigned int j = 0; j < system_matrix.n (); ++j)
5953  for (unsigned int q_point = 0; q_point < n_edge_points;
5954  ++q_point)
5955  system_matrix (i, j)
5956  += boundary_weights (q_point, j)
5957  * lobatto_polynomials_grad[i + 1].value
5958  (this->generalized_face_support_points[q_point]
5959  (1));
5960 
5961  FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
5962 
5963  system_matrix_inv.invert (system_matrix);
5964 
5965  const unsigned int
5966  line_coordinate[GeometryInfo<3>::lines_per_cell]
5967  = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
5968  Vector<double> system_rhs (system_matrix.m ());
5969  Vector<double> solution (system_rhs.size ());
5970 
5971  for (unsigned int line = 0;
5972  line < GeometryInfo<dim>::lines_per_cell; ++line)
5973  {
5974  // Set up the right hand side.
5975  system_rhs = 0;
5976 
5977  for (unsigned int q_point = 0; q_point < this->degree; ++q_point)
5978  {
5979  const double tmp
5980  = values[line_coordinate[line]][line * this->degree
5981  + q_point]
5982  - local_dofs[line * this->degree]
5983  * this->shape_value_component
5984  (line * this->degree,
5985  this->generalized_support_points[line
5986  * this->degree
5987  + q_point],
5988  line_coordinate[line]);
5989 
5990  for (unsigned int i = 0; i < system_rhs.size (); ++i)
5991  system_rhs (i) += boundary_weights (q_point, i)
5992  * tmp;
5993  }
5994 
5995  system_matrix_inv.vmult (solution, system_rhs);
5996 
5997  // Add the computed values
5998  // to the resulting vector
5999  // only, if they are not
6000  // too small.
6001  for (unsigned int i = 0; i < solution.size (); ++i)
6002  if (std::abs (solution (i)) > 1e-14)
6003  local_dofs[line * this->degree + i + 1] = solution (i);
6004  }
6005 
6006  // Then we go on to the
6007  // face shape functions.
6008  // Again we set up the
6009  // system matrix and
6010  // use it for both, the
6011  // horizontal and the
6012  // vertical, shape
6013  // functions.
6014  const std::vector<Polynomials::Polynomial<double> > &
6015  legendre_polynomials
6016  = Polynomials::Legendre::generate_complete_basis (this->degree-1);
6017  const unsigned int n_face_points = n_edge_points * n_edge_points;
6018 
6019  system_matrix.reinit ((this->degree-1) * this->degree,
6020  (this->degree-1) * this->degree);
6021  system_matrix = 0;
6022 
6023  for (unsigned int i = 0; i < this->degree; ++i)
6024  for (unsigned int j = 0; j < this->degree-1; ++j)
6025  for (unsigned int k = 0; k < this->degree; ++k)
6026  for (unsigned int l = 0; l < this->degree-1; ++l)
6027  for (unsigned int q_point = 0; q_point < n_face_points;
6028  ++q_point)
6029  system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
6030  += boundary_weights (q_point + n_edge_points,
6031  2 * (k * (this->degree-1) + l))
6032  * legendre_polynomials[i].value
6033  (this->generalized_face_support_points[q_point
6034  + 4
6035  * n_edge_points]
6036  (0))
6037  * lobatto_polynomials[j + 2].value
6038  (this->generalized_face_support_points[q_point
6039  + 4
6040  * n_edge_points]
6041  (1));
6042 
6043  system_matrix_inv.reinit (system_matrix.m (),
6044  system_matrix.m ());
6045  system_matrix_inv.invert (system_matrix);
6046  solution.reinit (system_matrix.m ());
6047  system_rhs.reinit (system_matrix.m ());
6048 
6049  const unsigned int
6050  face_coordinates[GeometryInfo<3>::faces_per_cell][2]
6051  = {{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
6052  const unsigned int
6053  edge_indices[GeometryInfo<3>::faces_per_cell][GeometryInfo<3>::lines_per_face]
6054  = {{0, 4, 8, 10}, {1, 5, 9, 11}, {8, 9, 2, 6},
6055  {10, 11, 3, 7}, {2, 3, 0, 1}, {6, 7, 4, 5}
6056  };
6057 
6058  for (unsigned int face = 0;
6059  face < GeometryInfo<dim>::faces_per_cell; ++face)
6060  {
6061  // Set up the right hand side
6062  // for the horizontal shape
6063  // functions.
6064  system_rhs = 0;
6065 
6066  for (unsigned int q_point = 0; q_point < n_face_points;
6067  ++q_point)
6068  {
6069  double tmp
6070  = values[face_coordinates[face][0]][q_point
6071  + GeometryInfo<dim>::lines_per_cell
6072  * n_edge_points];
6073 
6074  for (unsigned int i = 0; i < 2; ++i)
6075  for (unsigned int j = 0; j <= deg; ++j)
6076  tmp -= local_dofs[edge_indices[face][i]
6077  * this->degree + j]
6078  * this->shape_value_component
6079  (edge_indices[face][i] * this->degree + j,
6080  this->generalized_support_points[q_point
6081  + GeometryInfo<dim>::lines_per_cell
6082  * n_edge_points],
6083  face_coordinates[face][0]);
6084 
6085  for (unsigned int i = 0; i <= deg; ++i)
6086  for (unsigned int j = 0; j < deg; ++j)
6087  system_rhs (i * deg + j)
6088  += boundary_weights (q_point + n_edge_points,
6089  2 * (i * deg + j)) * tmp;
6090  }
6091 
6092  system_matrix_inv.vmult (solution, system_rhs);
6093 
6094  // Add the computed values
6095  // to the resulting vector
6096  // only, if they are not
6097  // too small.
6098  for (unsigned int i = 0; i <= deg; ++i)
6099  for (unsigned int j = 0; j < deg; ++j)
6100  if (std::abs (solution (i * deg + j)) > 1e-14)
6101  local_dofs[(2 * face * this->degree + i
6102  + GeometryInfo<dim>::lines_per_cell) * deg
6103  + j + GeometryInfo<dim>::lines_per_cell]
6104  = solution (i * deg + j);
6105 
6106  // Set up the right hand side
6107  // for the vertical shape
6108  // functions.
6109  system_rhs = 0;
6110 
6111  for (unsigned int q_point = 0; q_point < n_face_points;
6112  ++q_point)
6113  {
6114  double tmp
6115  = values[face_coordinates[face][1]][q_point
6116  + GeometryInfo<dim>::lines_per_cell
6117  * n_edge_points];
6118 
6119  for (int i = 2; i < (int) GeometryInfo<dim>::lines_per_face; ++i)
6120  for (unsigned int j = 0; j <= deg; ++j)
6121  tmp -= local_dofs[edge_indices[face][i]
6122  * this->degree + j]
6123  * this->shape_value_component
6124  (edge_indices[face][i] * this->degree + j,
6125  this->generalized_support_points[q_point
6126  + GeometryInfo<dim>::lines_per_cell
6127  * n_edge_points],
6128  face_coordinates[face][1]);
6129 
6130  for (unsigned int i = 0; i <= deg; ++i)
6131  for (unsigned int j = 0; j < deg; ++j)
6132  system_rhs (i * deg + j)
6133  += boundary_weights (q_point + n_edge_points,
6134  2 * (i * deg + j) + 1)
6135  * tmp;
6136  }
6137 
6138  system_matrix_inv.vmult (solution, system_rhs);
6139 
6140  // Add the computed values
6141  // to the resulting vector
6142  // only, if they are not
6143  // too small.
6144  for (unsigned int i = 0; i <= deg; ++i)
6145  for (unsigned int j = 0; j < deg; ++j)
6146  if (std::abs (solution (i * deg + j)) > 1e-14)
6147  local_dofs[((2 * face + 1) * deg + j + GeometryInfo<dim>::lines_per_cell)
6148  * this->degree + i]
6149  = solution (i * deg + j);
6150  }
6151 
6152  // Finally we project
6153  // the remaining parts
6154  // of the function on
6155  // the interior shape
6156  // functions.
6157  const QGauss<dim> reference_quadrature (this->degree);
6158  const unsigned int
6159  n_interior_points = reference_quadrature.size ();
6160 
6161  // We create the
6162  // system matrix.
6163  system_matrix.reinit (this->degree * deg * deg,
6164  this->degree * deg * deg);
6165  system_matrix = 0;
6166 
6167  for (unsigned int i = 0; i <= deg; ++i)
6168  for (unsigned int j = 0; j < deg; ++j)
6169  for (unsigned int k = 0; k < deg; ++k)
6170  for (unsigned int l = 0; l <= deg; ++l)
6171  for (unsigned int m = 0; m < deg; ++m)
6172  for (unsigned int n = 0; n < deg; ++n)
6173  for (unsigned int q_point = 0;
6174  q_point < n_interior_points; ++q_point)
6175  system_matrix ((i * deg + j) * deg + k,
6176  (l * deg + m) * deg + n)
6177  += reference_quadrature.weight (q_point)
6178  * legendre_polynomials[i].value
6179  (this->generalized_support_points[q_point
6180  + GeometryInfo<dim>::lines_per_cell
6181  * n_edge_points
6182  + GeometryInfo<dim>::faces_per_cell
6183  * n_face_points]
6184  (0))
6185  * lobatto_polynomials[j + 2].value
6186  (this->generalized_support_points[q_point
6187  + GeometryInfo<dim>::lines_per_cell
6188  * n_edge_points
6189  + GeometryInfo<dim>::faces_per_cell
6190  * n_face_points]
6191  (1))
6192  * lobatto_polynomials[k + 2].value
6193  (this->generalized_support_points[q_point
6194  + GeometryInfo<dim>::lines_per_cell
6195  * n_edge_points
6196  + GeometryInfo<dim>::faces_per_cell
6197  * n_face_points]
6198  (2))
6199  * lobatto_polynomials_grad[l].value
6200  (this->generalized_support_points[q_point
6201  + GeometryInfo<dim>::lines_per_cell
6202  * n_edge_points
6203  + GeometryInfo<dim>::faces_per_cell
6204  * n_face_points]
6205  (0))
6206  * lobatto_polynomials[m + 2].value
6207  (this->generalized_support_points[q_point
6208  + GeometryInfo<dim>::lines_per_cell
6209  * n_edge_points
6210  + GeometryInfo<dim>::faces_per_cell
6211  * n_face_points]
6212  (1))
6213  * lobatto_polynomials[n + 2].value
6214  (this->generalized_support_points[q_point
6215  + GeometryInfo<dim>::lines_per_cell
6216  * n_edge_points
6217  + GeometryInfo<dim>::faces_per_cell
6218  * n_face_points]
6219  (2));
6220 
6221  system_matrix_inv.reinit (system_matrix.m (),
6222  system_matrix.m ());
6223  system_matrix_inv.invert (system_matrix);
6224  // Set up the right hand side.
6225  system_rhs.reinit (system_matrix.m ());
6226  system_rhs = 0;
6227 
6228  for (unsigned int q_point = 0; q_point < n_interior_points;
6229  ++q_point)
6230  {
6231  double tmp
6232  = values[0][q_point + GeometryInfo<dim>::lines_per_cell
6233  * n_edge_points
6234  + GeometryInfo<dim>::faces_per_cell
6235  * n_face_points];
6236 
6237  for (unsigned int i = 0; i <= deg; ++i)
6238  {
6239  for (unsigned int j = 0; j < 2; ++j)
6240  for (unsigned int k = 0; k < 2; ++k)
6241  tmp -= local_dofs[i + (j + 4 * k + 2) * this->degree]
6242  * this->shape_value_component
6243  (i + (j + 4 * k + 2) * this->degree,
6244  this->generalized_support_points[q_point
6245  + GeometryInfo<dim>::lines_per_cell
6246  * n_edge_points
6247  + GeometryInfo<dim>::faces_per_cell
6248  * n_face_points],
6249  0);
6250 
6251  for (unsigned int j = 0; j < deg; ++j)
6252  for (unsigned int k = 0; k < 4; ++k)
6253  tmp -= local_dofs[(i + 2 * (k + 2) * this->degree
6254  + GeometryInfo<dim>::lines_per_cell)
6255  * deg + j
6256  + GeometryInfo<dim>::lines_per_cell]
6257  * this->shape_value_component
6258  ((i + 2 * (k + 2) * this->degree
6259  + GeometryInfo<dim>::lines_per_cell)
6260  * deg + j
6261  + GeometryInfo<dim>::lines_per_cell,
6262  this->generalized_support_points[q_point
6263  + GeometryInfo<dim>::lines_per_cell
6264  * n_edge_points
6265  + GeometryInfo<dim>::faces_per_cell
6266  * n_face_points],
6267  0);
6268  }
6269 
6270  for (unsigned int i = 0; i <= deg; ++i)
6271  for (unsigned int j = 0; j < deg; ++j)
6272  for (unsigned int k = 0; k < deg; ++k)
6273  system_rhs ((i * deg + j) * deg + k)
6274  += reference_quadrature.weight (q_point) * tmp
6275  * lobatto_polynomials_grad[i].value
6276  (this->generalized_support_points[q_point
6277  + GeometryInfo<dim>::lines_per_cell
6278  * n_edge_points
6279  + GeometryInfo<dim>::faces_per_cell
6280  * n_face_points]
6281  (0))
6282  * lobatto_polynomials[j + 2].value
6283  (this->generalized_support_points[q_point
6284  + GeometryInfo<dim>::lines_per_cell
6285  * n_edge_points
6286  + GeometryInfo<dim>::faces_per_cell
6287  * n_face_points]
6288  (1))
6289  * lobatto_polynomials[k + 2].value
6290  (this->generalized_support_points[q_point
6291  + GeometryInfo<dim>::lines_per_cell
6292  * n_edge_points
6293  + GeometryInfo<dim>::faces_per_cell
6294  * n_face_points]
6295  (2));
6296  }
6297 
6298  solution.reinit (system_rhs.size ());
6299  system_matrix_inv.vmult (solution, system_rhs);
6300 
6301  // Add the computed values
6302  // to the resulting vector
6303  // only, if they are not
6304  // too small.
6305  for (unsigned int i = 0; i <= deg; ++i)
6306  for (unsigned int j = 0; j < deg; ++j)
6307  for (unsigned int k = 0; k < deg; ++k)
6308  if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
6309  local_dofs[((i + 2 * GeometryInfo<dim>::faces_per_cell)
6310  * deg + j + GeometryInfo<dim>::lines_per_cell
6311  + 2 * GeometryInfo<dim>::faces_per_cell)
6312  * deg + k + GeometryInfo<dim>::lines_per_cell]
6313  = solution ((i * deg + j) * deg + k);
6314 
6315  // Set up the right hand side.
6316  system_rhs = 0;
6317 
6318  for (unsigned int q_point = 0; q_point < n_interior_points;
6319  ++q_point)
6320  {
6321  double tmp
6322  = values[1][q_point + GeometryInfo<dim>::lines_per_cell
6323  * n_edge_points
6324  + GeometryInfo<dim>::faces_per_cell
6325  * n_face_points];
6326 
6327  for (unsigned int i = 0; i <= deg; ++i)
6328  for (unsigned int j = 0; j < 2; ++j)
6329  {
6330  for (unsigned int k = 0; k < 2; ++k)
6331  tmp -= local_dofs[i + (4 * j + k) * this->degree]
6332  * this->shape_value_component
6333  (i + (4 * j + k) * this->degree,
6334  this->generalized_support_points[q_point
6335  + GeometryInfo<dim>::lines_per_cell
6336  * n_edge_points
6337  + GeometryInfo<dim>::faces_per_cell
6338  * n_face_points],
6339  1);
6340 
6341  for (unsigned int k = 0; k < deg; ++k)
6342  tmp -= local_dofs[(i + 2 * j * this->degree
6343  + GeometryInfo<dim>::lines_per_cell)
6344  * deg + k
6345  + GeometryInfo<dim>::lines_per_cell]
6346  * this->shape_value_component
6347  ((i + 2 * j * this->degree
6348  + GeometryInfo<dim>::lines_per_cell)
6349  * deg + k
6350  + GeometryInfo<dim>::lines_per_cell,
6351  this->generalized_support_points[q_point
6352  + GeometryInfo<dim>::lines_per_cell
6353  * n_edge_points
6354  + GeometryInfo<dim>::faces_per_cell
6355  * n_face_points],
6356  1)
6357  + local_dofs[i + ((2 * j + 9) * deg + k
6358  + GeometryInfo<dim>::lines_per_cell)
6359  * this->degree]
6360  * this->shape_value_component
6361  (i + ((2 * j + 9) * deg + k
6362  + GeometryInfo<dim>::lines_per_cell)
6363  * this->degree,
6364  this->generalized_support_points[q_point
6365  + GeometryInfo<dim>::lines_per_cell
6366  * n_edge_points
6367  + GeometryInfo<dim>::faces_per_cell
6368  * n_face_points],
6369  1);
6370  }
6371 
6372  for (unsigned int i = 0; i <= deg; ++i)
6373  for (unsigned int j = 0; j < deg; ++j)
6374  for (unsigned int k = 0; k < deg; ++k)
6375  system_rhs ((i * deg + j) * deg + k)
6376  += reference_quadrature.weight (q_point) * tmp
6377  * lobatto_polynomials_grad[i].value
6378  (this->generalized_support_points[q_point
6379  + GeometryInfo<dim>::lines_per_cell
6380  * n_edge_points
6381  + GeometryInfo<dim>::faces_per_cell
6382  * n_face_points]
6383  (1))
6384  * lobatto_polynomials[j + 2].value
6385  (this->generalized_support_points[q_point
6386  + GeometryInfo<dim>::lines_per_cell
6387  * n_edge_points
6388  + GeometryInfo<dim>::faces_per_cell
6389  * n_face_points]
6390  (0))
6391  * lobatto_polynomials[k + 2].value
6392  (this->generalized_support_points[q_point
6393  + GeometryInfo<dim>::lines_per_cell
6394  * n_edge_points
6395  + GeometryInfo<dim>::faces_per_cell
6396  * n_face_points]
6397  (2));
6398  }
6399 
6400  system_matrix_inv.vmult (solution, system_rhs);
6401 
6402  // Add the computed values
6403  // to the resulting vector
6404  // only, if they are not
6405  // too small.
6406  for (unsigned int i = 0; i <= deg; ++i)
6407  for (unsigned int j = 0; j < deg; ++j)
6408  for (unsigned int k = 0; k < deg; ++k)
6409  if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
6410  local_dofs[((i + this->degree + 2
6411  * GeometryInfo<dim>::faces_per_cell) * deg
6412  + j + GeometryInfo<dim>::lines_per_cell + 2
6413  * GeometryInfo<dim>::faces_per_cell) * deg
6414  + k + GeometryInfo<dim>::lines_per_cell]
6415  = solution ((i * deg + j) * deg + k);
6416 
6417  // Set up the right hand side.
6418  system_rhs = 0;
6419 
6420  for (unsigned int q_point = 0; q_point < n_interior_points;
6421  ++q_point)
6422  {
6423  double tmp
6424  = values[2][q_point + GeometryInfo<dim>::lines_per_cell
6425  * n_edge_points
6426  + GeometryInfo<dim>::faces_per_cell
6427  * n_face_points];
6428 
6429  for (unsigned int i = 0; i <= deg; ++i)
6430  for (unsigned int j = 0; j < 4; ++j)
6431  {
6432  tmp -= local_dofs[i + (j + 8) * this->degree]
6433  * this->shape_value_component
6434  (i + (j + 8) * this->degree,
6435  this->generalized_support_points[q_point
6436  + GeometryInfo<dim>::lines_per_cell
6437  * n_edge_points
6438  + GeometryInfo<dim>::faces_per_cell
6439  * n_face_points],
6440  2);
6441 
6442  for (unsigned int k = 0; k < deg; ++k)
6443  tmp -= local_dofs[i + ((2 * j + 1) * deg + k
6444  + GeometryInfo<dim>::lines_per_cell)
6445  * this->degree]
6446  * this->shape_value_component
6447  (i + ((2 * j + 1) * deg + k
6448  + GeometryInfo<dim>::lines_per_cell)
6449  * this->degree,
6450  this->generalized_support_points[q_point
6451  + GeometryInfo<dim>::lines_per_cell
6452  * n_edge_points
6453  + GeometryInfo<dim>::faces_per_cell
6454  * n_face_points],
6455  2);
6456  }
6457 
6458  for (unsigned int i = 0; i <= deg; ++i)
6459  for (unsigned int j = 0; j < deg; ++j)
6460  for (unsigned int k = 0; k < deg; ++k)
6461  system_rhs ((i * deg + j) * deg + k)
6462  += reference_quadrature.weight (q_point) * tmp
6463  * lobatto_polynomials_grad[i].value
6464  (this->generalized_support_points[q_point
6465  + GeometryInfo<dim>::lines_per_cell
6466  * n_edge_points
6467  + GeometryInfo<dim>::faces_per_cell
6468  * n_face_points]
6469  (2))
6470  * lobatto_polynomials[j + 2].value
6471  (this->generalized_support_points[q_point
6472  + GeometryInfo<dim>::lines_per_cell
6473  * n_edge_points
6474  + GeometryInfo<dim>::faces_per_cell
6475  * n_face_points]
6476  (0))
6477  * lobatto_polynomials[k + 2].value
6478  (this->generalized_support_points[q_point
6479  + GeometryInfo<dim>::lines_per_cell
6480  * n_edge_points
6481  + GeometryInfo<dim>::faces_per_cell
6482  * n_face_points]
6483  (1));
6484  }
6485 
6486  system_matrix_inv.vmult (solution, system_rhs);
6487 
6488  // Add the computed values
6489  // to the resulting vector
6490  // only, if they are not
6491  // too small.
6492  for (unsigned int i = 0; i <= deg; ++i)
6493  for (unsigned int j = 0; j < deg; ++j)
6494  for (unsigned int k = 0; k < deg; ++k)
6495  if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
6496  local_dofs[i + ((j + 2 * (deg
6497  + GeometryInfo<dim>::faces_per_cell))
6498  * deg + k
6499  + GeometryInfo<dim>::lines_per_cell)
6500  * this->degree]
6501  = solution ((i * deg + j) * deg + k);
6502  }
6503 
6504  break;
6505  }
6506 
6507  default:
6508  Assert (false, ExcNotImplemented ());
6509  }
6510 }
6511 
6512 
6513 
6514 template <int dim>
6515 std::pair<Table<2,bool>, std::vector<unsigned int> >
6516 FE_Nedelec<dim>::get_constant_modes() const
6517 {
6518  Table<2,bool> constant_modes(dim, this->dofs_per_cell);
6519  for (unsigned int d=0; d<dim; ++d)
6520  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
6521  constant_modes(d,i) = true;
6522  std::vector<unsigned int> components;
6523  for (unsigned int d=0; d<dim; ++d)
6524  components.push_back(d);
6525  return std::pair<Table<2,bool>, std::vector<unsigned int> >
6526  (constant_modes, components);
6527 }
6528 
6529 
6530 template <int dim>
6531 std::size_t
6532 FE_Nedelec<dim>::memory_consumption () const
6533 {
6534  Assert (false, ExcNotImplemented ());
6535  return 0;
6536 }
6537 
6538 
6539 //----------------------------------------------------------------------//
6540 
6541 
6542 // explicit instantiations
6543 #include "fe_nedelec.inst"
6544 
6545 
6546 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:2022
FE_Nedelec::FE_Nedelec
friend class FE_Nedelec
Definition: fe_nedelec.h:329
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:2211
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:3041
FE_Nedelec::compare_for_face_domination
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim > &fe_other) const
Definition: fe_nedelec.cc:2306
Polynomials::Lobatto::generate_complete_basis
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Definition: polynomial.cc:959
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:873
FE_Nedelec::memory_consumption
virtual std::size_t memory_consumption() const
Definition: fe_nedelec.cc:6532
FE_PolyTensor< PolynomialsNedelec< dim >, dim >::mapping_type
MappingType mapping_type
Definition: fe_poly_tensor.h:205
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:2571
std
STL namespace.
mapping_nedelec
Definition: mapping.h:101
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)
Definition: fe_tools.cc:1897
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
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:2367
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:1986
FE_PolyTensor
Definition: fe_poly_tensor.h:139
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:193
FullMatrix::n
size_type n() const
FiniteElement< dim, dim >::prolongation
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2199
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:313
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Vector::size
std::size_t size() const
FiniteElement::get_name
virtual std::string get_name() const =0
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:2990
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:2408
PolynomialsNedelec::compute_n_pols
static unsigned int compute_n_pols(unsigned int degree)
Definition: polynomials_nedelec.cc:1388
FE_Nedelec::interpolate
virtual void interpolate(std::vector< double > &local_dofs, const std::vector< double > &values) const 1
Definition: fe_nedelec.cc:4020
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:2357
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
FE_Nedelec::clone
virtual FiniteElement< dim > * clone() const
Definition: fe_nedelec.cc:211
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:277
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:176
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
FE_Nedelec::hp_constraints_are_implemented
virtual bool hp_constraints_are_implemented() const
Definition: fe_nedelec.cc:2350
FE_Nedelec::initialize_restriction
void initialize_restriction()
Definition: fe_nedelec.cc:520
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:432
FE_Nedelec::convert_generalized_support_point_values_to_nodal_values
virtual void convert_generalized_support_point_values_to_nodal_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const
Definition: fe_nedelec.cc:3099
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:6516
FE_Nedelec::get_face_interpolation_matrix
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const
Definition: fe_nedelec.cc:2465
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