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