]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added FE_NedelecSZ, a H(curl)-conforming element which overcomes the sign-conflict... 2240/head
authorRoss Kynch <rkynch@gmail.com>
Sat, 28 Apr 2018 20:48:37 +0000 (21:48 +0100)
committerRoss Kynch <rkynch@gmail.com>
Sat, 28 Apr 2018 20:48:37 +0000 (21:48 +0100)
14 files changed:
doc/news/changes/major/20180428RossKynch [new file with mode: 0644]
include/deal.II/fe/fe_nedelec_sz.h [new file with mode: 0644]
include/deal.II/numerics/vector_tools.templates.h
source/fe/CMakeLists.txt
source/fe/fe_nedelec_sz.cc [new file with mode: 0644]
source/fe/fe_nedelec_sz.inst.in [new file with mode: 0644]
tests/fe/fe_nedelec_sz_face_orientation.cc [new file with mode: 0644]
tests/fe/fe_nedelec_sz_face_orientation.output [new file with mode: 0644]
tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc [new file with mode: 0644]
tests/fe/fe_nedelec_sz_gradient_divergence_theorem.output [new file with mode: 0644]
tests/fe/fe_nedelec_sz_non_rect_2d.cc [new file with mode: 0644]
tests/fe/fe_nedelec_sz_non_rect_2d.output [new file with mode: 0644]
tests/fe/fe_nedelec_sz_non_rect_face.cc [new file with mode: 0644]
tests/fe/fe_nedelec_sz_non_rect_face.with_umfpack=on.output [new file with mode: 0644]

diff --git a/doc/news/changes/major/20180428RossKynch b/doc/news/changes/major/20180428RossKynch
new file mode 100644 (file)
index 0000000..0d86c7d
--- /dev/null
@@ -0,0 +1,8 @@
+New: FE_NedelecSZ is a new H(curl)-conforming element which overcomes the sign conflict problem
+which can be encountered for Nedelec elements on hexahedral meshes. This is an implementation of
+the method described in the PhD thesis of Sabine Zaglmayr and should be considered over the
+FE_Nedelec element where meshes with non-standard orientation are used. Note that not all
+functionality (e.g. Hessians), has been implemented for this element.
+<br>
+(Ross Kynch, 2018/04/28)
+
diff --git a/include/deal.II/fe/fe_nedelec_sz.h b/include/deal.II/fe/fe_nedelec_sz.h
new file mode 100644 (file)
index 0000000..5575de3
--- /dev/null
@@ -0,0 +1,427 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__fe_nedelec_sz_h
+#define dealii__fe_nedelec_sz_h
+
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/polynomials_integrated_legendre_sz.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/*!@addtogroup fe */
+/*@{*/
+
+/**
+* This class represents an implementation of the
+* H<sup>curl</sup>-conforming Nédélec element described in the PhD
+* thesis of S. Zaglmayr, <b>High Order Finite Element Methods
+* for Electromagnetic Field Computation</b>, Johannes Kepler
+* Univerität Linz, 2006.
+*
+* This element overcomes the sign conflict issues present in
+* traditional Nédélec elements which arise from the edge and face
+* parameterisations used in the basis functions. Therefore, this
+* element should provide consistent results for general quadrilateral
+* and hexahedral elements.
+*
+* The way this element addresses the sign conflict problem is to assign
+* local edges and faces a globally defined orientation for all cells. The
+* local edge orientation is always chosen such that the first vertex defining
+* the edge is always chosen to be that which has the highest global vertex
+* numbering, with the second being that which has the lowest global vertex
+* numbering.
+*
+* Similarly, the face orientation is always chosen such that the first
+* vertex is chosen to be that with the highest global vertex numbering of the
+* four vertices making up the face. The third vertex is then chosen to be that
+* which is geometrically opposite the first vertex, and the second and fourth
+* vertices are decided such that the second has a higher global vertex numbering
+* than the fourth.
+*
+* Note that this element does not support non-conforming meshes at this time.
+*
+* Further details on this element, including some benchmarking, can be
+* in the paper R. Kynch, P. Ledger: <b>Resolving the sign conflict
+* problem for hp\96hexahedral Nédélec elements with application to eddy
+* current problems</b>, Computers & Structures 181, 41-54, 2017 (see
+* https://doi.org/10.1016/j.compstruc.2016.05.021).
+*
+* @author Ross Kynch
+* @date 2016, 2017, 2018
+*/
+template <int dim>
+class FE_NedelecSZ : public FiniteElement<dim,dim>
+{
+public:
+  /**
+  * Constructor for an element of given @p degree.
+  */
+  FE_NedelecSZ (const unsigned int degree);
+
+  virtual
+  UpdateFlags
+  requires_update_flags (const UpdateFlags update_flags) const;
+
+  virtual std::string get_name () const;
+
+  virtual std::unique_ptr<FiniteElement<dim,dim> > clone() const;
+
+  /**
+  * This element is vector-valued so this function will
+  * throw an exception.
+  */
+  virtual double shape_value (const unsigned int i,
+                              const Point<dim>  &p) const;
+
+  /**
+  * Not implemented.
+  */
+  virtual double shape_value_component (const unsigned int i,
+                                        const Point<dim>  &p,
+                                        const unsigned int component) const;
+
+  /**
+  * This element is vector-valued so this function will
+  * throw an exception.
+  */
+  virtual Tensor<1,dim> shape_grad (const unsigned int i,
+                                    const Point<dim>  &p) const;
+
+  /**
+  * Not implemented.
+  */
+  virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+                                              const Point<dim>  &p,
+                                              const unsigned int component) const;
+
+  /**
+  * This element is vector-valued so this function will
+  * throw an exception.
+  */
+  virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
+                                         const Point<dim>  &p) const;
+
+  /**
+  * Not implemented.
+  */
+  virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+                                                   const Point<dim>  &p,
+                                                   const unsigned int component) const;
+
+  /**
+   * Given <tt>flags</tt>, determines the values which must be computed only
+   * for the reference cell. Make sure, that #mapping_type is set by the
+   * derived class, such that this function can operate correctly.
+   */
+  UpdateFlags update_once (const UpdateFlags flags) const;
+
+  /**
+   * Given <tt>flags</tt>, determines the values which must be computed in
+   * each cell cell. Make sure, that #mapping_type is set by the derived
+   * class, such that this function can operate correctly.
+   */
+  UpdateFlags update_each (const UpdateFlags flags) const;
+
+protected:
+  /**
+   * The mapping type to be used to map shape functions from the reference
+   * cell to the mesh cell.
+   */
+  MappingType mapping_type;
+
+  virtual
+  std::unique_ptr<typename FiniteElement<dim,dim>::InternalDataBase>
+  get_data (const UpdateFlags                                                             update_flags,
+            const Mapping<dim,dim>                                                       &mapping,
+            const Quadrature<dim>                                                        &quadrature,
+            dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &output_data) const;
+
+  /**
+  * Compute information about the shape functions on the cell denoted by the
+  * first argument. Note that this function must recompute the cell-dependent
+  * degrees of freedom, and so is not thread-safe at this time.
+  */
+  virtual void
+  fill_fe_values (const typename Triangulation<dim,dim>::cell_iterator                         &cell,
+                  const CellSimilarity::Similarity                                              cell_similarity,
+                  const Quadrature<dim>                                                        &quadrature,
+                  const Mapping<dim,dim>                                                       &mapping,
+                  const typename Mapping<dim,dim>::InternalDataBase                            &mapping_internal,
+                  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
+                  const typename FiniteElement<dim,dim>::InternalDataBase                      &fedata,
+                  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &data) const;
+
+  /**
+  * Compute information about the shape functions on the cell and face denoted
+  * by the first two arguments. Note that this function must recompute the
+  * cell-dependent degrees of freedom, and so is not thread-safe at this time.
+  */
+  virtual void
+  fill_fe_face_values (const typename Triangulation<dim,dim>::cell_iterator                         &cell,
+                       const unsigned int                                                            face_no,
+                       const Quadrature<dim-1>                                                      &quadrature,
+                       const Mapping<dim,dim>                                                       &mapping,
+                       const typename Mapping<dim,dim>::InternalDataBase                            &mapping_internal,
+                       const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
+                       const typename FiniteElement<dim,dim>::InternalDataBase                      &fedata,
+                       dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &data) const;
+
+  /**
+  * Not implemented.
+  */
+  virtual void
+  fill_fe_subface_values (const typename Triangulation<dim,dim>::cell_iterator                         &cell,
+                          const unsigned int                                                            face_no,
+                          const unsigned int                                                            sub_no,
+                          const Quadrature<dim-1>                                                      &quadrature,
+                          const Mapping<dim,dim>                                                       &mapping,
+                          const typename Mapping<dim,dim>::InternalDataBase                            &mapping_internal,
+                          const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
+                          const typename FiniteElement<dim,dim>::InternalDataBase                      &fedata,
+                          dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &data) const;
+
+  /**
+  * Derived Internal data which is used to store cell-independent data.
+  * Note that due to the nature of this element, a number of useful
+  * pre-computed quantities are stored for the computation of cell-dependent
+  * shape functions.
+  *
+  * The main quantities which are stored are associated with edge and face
+  * parameterisations. These are:
+  * <ul>
+  * <li> \lambda_{i} - trilinear function, equal to one at the i-th vertex
+  * and zero at all other vertices.</li>
+  * <li> \sigma_{i} - linear functional associated with the i-th vertex.</li>
+  * </ul>
+  *
+  * The definitions of these functionals, as well as the edge and face
+  * parametrisations and edge and face extension parameters, can be found on
+  * page 82 of Zaglmayr's thesis. The details of the definition of the
+  * globally-defined edge and face orientations can be found on page 67.
+  */
+  class InternalData : public FiniteElement<dim,dim>::InternalDataBase
+  {
+  public:
+    /**
+    * Storage for shape functions on the reference element. We only pre-compute
+    * cell-based DoFs, as the edge- and face-based DoFs depend on the cell.
+    *
+    * Due to the cell-dependent DoFs, this variable is declared mutable.
+    */
+    mutable std::vector<std::vector<Tensor<1,dim> > > shape_values;
+
+    /**
+    * Storage for shape function gradients on the reference element. We only
+    * pre-compute cell-based DoFs, as the edge- and face-based DoFs depend on
+    * the cell.
+    *
+    * Due to the cell-dependent DoFs, this variable is declared mutable.
+    */
+    mutable std::vector<std::vector<DerivativeForm<1, dim, dim> > > shape_grads;
+
+    /**
+    * Storage for all possible edge parameterisation between vertices. These
+    * are required in the computation of edge- and face-based DoFs, which are
+    * cell-dependent.
+    *
+    * The edge parameterisation of an edge, E, starting at vertex i and ending
+    * at vertex j is given by \sigma_{E} = \sigma_{i} - \sigma{j}.
+    *
+    * sigma_imj_values[q][i][j] stores the value of the edge parametrisation
+    * connected by vertices i and j at the q-th quadrature point.
+    *
+    * Note that not all of the i and j combinations result in valid edges on
+    * the hexahedral cell, but they are computed in this fashion for use with
+    * non-standard edge and face orientations.
+    */
+    std::vector<std::vector<std::vector<double> > > sigma_imj_values;
+
+    /**
+    * Storage for gradients of all possible edge parameterisations between
+    * vertices. These are required in the computation of edge- and face-based
+    * DoFs, which are cell-dependent. Note that the components of the gradient
+    * are constant.
+    *
+    * The edge parameterisation of an edge, E, starting at vertex i and ending
+    * at vertex j is given by \sigma_{E} = \sigma_{i} - \sigma{j}.
+    *
+    * sigma_imj_grads[i][j][d] stores the gradient of the edge parametrisation
+    * connected by vertices i and j in component d.
+    *
+    * Note that the gradient of the edge parametrisation is constant on an
+    * edge, so we do not need to store it at every quadrature point.
+    */
+    std::vector<std::vector<std::vector<double> > > sigma_imj_grads;
+
+    /**
+    * Storage for values of edge parametrisations at quadrature points. These
+    * are stored for the 12 edges such that the global vertex numbering would
+    * follow the order defined by the "standard" deal.II cell.
+    *
+    * edge_sigma_values[m][q] stores the edge parameterisation value at the
+    * q-th quadrature point on edge m.
+    *
+    * These values change with the orientation of the edges of a physical cell
+    * and so must take the "sign" into account when used for computation.
+    */
+    std::vector<std::vector<double> > edge_sigma_values;
+
+    /**
+    * Storage for gradients of edge parametrisation at quadrature points.
+    * These are stored for the 12 edges such that the global vertex numbering
+    * would follow the order defined by the "standard" deal.II cell.
+    *
+    * edge_sigma_grads[m][d] stores the gradient of the edge parameterisation
+    * for component d on edge m.
+    *
+    * These values change with the orientation of the edges of a physical cell
+    * and so must take the "sign" into account when used for computation.
+    *
+    */
+    std::vector<std::vector<double> > edge_sigma_grads;
+
+    /**
+    * Storage for edge extension parameters at quadrature points. These are
+    * stored for the 12 edges such that the global vertex numbering would
+    * follow the order defined by the "standard" deal.II cell.
+    *
+    * The edge extension parameter of an edge, E, starting at vertex i and
+    * ending at vertex j is given by \lambda_{E} = \lambda_{i} + \lambda_{j}.
+    *
+    * Note that under this definition, the values of \lambda_{E} do not change
+    * with the orientation of the edge.
+    *
+    * edge_lambda_values[m][q] stores the edge extension parameter value at
+    * the q-th quadrature point on edge m.
+    */
+    std::vector<std::vector<double> > edge_lambda_values;
+
+    /**
+    * Storage for gradients of edge extension parameters in 2D. In this case
+    * they are constant. These are stored for the 12 edges such that the global
+    * vertex numbering* would follow the order defined by the "standard"
+    * deal.II cell.
+    *
+    * edge_lambda_grads_2d[m][d] stores the gradient of the edge extension
+    * parameter for component d on edge m.
+    */
+    std::vector<std::vector<double> > edge_lambda_grads_2d;
+
+    /**
+    * Storage for gradients of edge extension parameters in 3D. In this case
+    * they are non-constant. These are stored for the 12 edges such that the
+    * global vertex numbering* would follow the order defined by the
+    * "standard" deal.II cell.
+    *
+    * edge_lambda_grads_3d[m][q][d] stores the gradient of the edge extension
+    * parameter for component d at the q-th quadrature point on edge m.
+    */
+    std::vector<std::vector<std::vector<double> > > edge_lambda_grads_3d;
+
+    /**
+    * Storage for 2nd derivatives of edge extension parameters in 3D, which are
+    * constant across the cell. These are stored for the 12 edges such that the
+    * global vertex numbering* would follow the order defined by the
+    * "standard" deal.II cell.
+    *
+    * edge_lambda_gradgrads_3d[m][d1][d2] stores the 2nd derivatives of the
+    * edge extension parameters with respect to components d1 and d2 on edge m.
+    */
+    std::vector<std::vector<std::vector<double> > > edge_lambda_gradgrads_3d;
+
+    /**
+    * Storage for the face extension parameters. These are stored for the 6
+    * faces such that the global vertex numbering would follow the order
+    * defined by the "standard" deal.II cell.
+    *
+    * The face extension parameter of a face, F, defined by the vertices
+    * v1, v2, v3, v4 is given by
+    * \lambda_{F} = \lambda_{v1} + \lambda_{v2} + \lambda_{v3} + \lambda_{v4}.
+    *
+    * Note that under this definition, the values of \lambda_{F} do not change
+    * with the orientation of the face.
+    *
+    * face_lambda_values[m][q] stores the face extension parameter value at
+    * the q-th quadrature point on face m.
+    */
+    std::vector<std::vector<double> > face_lambda_values;
+
+    /**
+    * Storage for gradients of face extension parameters. These are stored for
+    * the 6 faces such that the global vertex numbering would follow the order
+    * defined by the "standard" deal.II cell.
+    *
+    * face_lambda_grads[m][d] stores the gradient of the face extension
+    * parameters for component d on face m.
+    */
+    std::vector<std::vector<double> > face_lambda_grads;
+  };
+
+private:
+  /**
+  * Internal function to return a vector of "dofs per object"
+  * where the components of the returned vector refer to:
+  * 0 = vertex
+  * 1 = edge
+  * 2 = face (which is a cell in 2D)
+  * 3 = cell
+  */
+  static std::vector<unsigned int> get_dpo_vector (const unsigned int degree);
+
+  /**
+   * Internal storage for all required integrated Legendre polynomials.
+   */
+  std::vector<Polynomials::Polynomial<double> > IntegratedLegendrePolynomials;
+
+  /**
+  * Internal function to populate the internal array of integrated Legendre
+  * polynomials.
+  */
+  void create_polynomials (const unsigned int degree);
+
+  /**
+  * Returns the number of DoFs in the basis set.
+  */
+  unsigned int compute_num_dofs (const unsigned int degree) const;
+
+  /**
+  * Populates cell-dependent edge-based shape functions on the given
+  * InternalData object.
+  */
+  void fill_edge_values(const typename Triangulation<dim,dim>::cell_iterator &cell,
+                        const Quadrature<dim>                                &quadrature,
+                        const InternalData                                   &fedata) const;
+
+  /**
+  * Populates the cell-dependent face-based shape functions on the given
+  * InternalData object.
+  */
+  void fill_face_values(const typename Triangulation<dim,dim>::cell_iterator &cell,
+                        const Quadrature<dim>                                &quadrature,
+                        const InternalData                                   &fedata) const;
+};
+
+
+
+/*@}*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 008e03b53ee361c7a917cff3daa02a6d848c9e1f..14a1aa8eaf4af3eaf8556ec794a1553d3b24f600 100644 (file)
@@ -59,6 +59,7 @@
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_q_dg0.h>
 #include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_nedelec_sz.h>
 #include <deal.II/fe/fe_raviart_thomas.h>
 #include <deal.II/fe/fe_system.h>
 #include <deal.II/fe/fe_tools.h>
@@ -4787,7 +4788,8 @@ namespace VectorTools
                && (fe.system_to_base_index (cell_idx).first == base_indices)
                && (lower_bound <= fe.system_to_base_index (cell_idx).second)
                && (fe.system_to_base_index (cell_idx).second <= upper_bound))
-              || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)
+              || (((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)
+                   || (dynamic_cast<const FE_NedelecSZ<dim>*> (&fe) != nullptr))
                   && (line * (degree + 1) <= face_idx)
                   && (face_idx <= (line + 1) * (degree + 1) - 1)))
             {
@@ -5005,7 +5007,8 @@ namespace VectorTools
               const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face);
               if (((dynamic_cast<const FESystem<dim>*> (&fe) != nullptr)
                    && (fe.system_to_base_index (cell_idx).first == base_indices))
-                  || (dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr))
+                  || (dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)
+                  || (dynamic_cast<const FE_NedelecSZ<dim>*> (&fe) != nullptr))
                 {
                   associated_edge_dof_to_face_dof[associated_edge_dof_index] = face_idx;
                   ++associated_edge_dof_index;
@@ -5122,7 +5125,8 @@ namespace VectorTools
                        && (fe.system_to_base_index (cell_idx).first == base_indices)
                        && (lower_bound <= fe.system_to_base_index (cell_idx).second)
                        && (fe.system_to_base_index (cell_idx).second <= upper_bound))
-                      || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)
+                      || (((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)
+                           || (dynamic_cast<const FE_NedelecSZ<dim>*> (&fe) != nullptr))
                           && (line * (degree + 1) <= face_idx)
                           && (face_idx <= (line + 1) * (degree + 1) - 1)))
                     {
@@ -5159,7 +5163,8 @@ namespace VectorTools
               const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face);
               if (((dynamic_cast<const FESystem<dim>*> (&fe) != nullptr)
                    && (fe.system_to_base_index (cell_idx).first == base_indices))
-                  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)))
+                  || (dynamic_cast<const FE_Nedelec<dim>*> (&fe) != nullptr)
+                  || (dynamic_cast<const FE_NedelecSZ<dim>*> (&fe) != nullptr))
                 {
                   associated_face_dof_to_face_dof[associated_face_dof_index] = face_idx;
                   ++associated_face_dof_index;
@@ -5357,7 +5362,8 @@ namespace VectorTools
                           // If the FE is a FESystem we cannot check this.
                           if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == nullptr)
                             {
-                              AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != nullptr,
+                              AssertThrow ((dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != nullptr)
+                                           || (dynamic_cast<const FE_NedelecSZ<dim>*> (&cell->get_fe ()) != nullptr),
                                            typename FiniteElement<dim>::ExcInterpolationNotImplemented ());
 
                             }
@@ -5453,7 +5459,8 @@ namespace VectorTools
                           // If the FE is a FESystem we cannot check this.
                           if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == nullptr)
                             {
-                              AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != nullptr,
+                              AssertThrow ((dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != nullptr)
+                                           || (dynamic_cast<const FE_NedelecSZ<dim>*> (&cell->get_fe ()) != nullptr),
                                            typename FiniteElement<dim>::ExcInterpolationNotImplemented ());
                             }
 
index f8fd6112534924071a5cba7d170df884e9e4674e..8ec617c984367ffa152b3aaa05944bd1b66f6b1f 100644 (file)
@@ -30,6 +30,7 @@ SET(_unity_include_src
   fe_dg_vector.cc
   fe_face.cc
   fe_nedelec.cc
+  fe_nedelec_sz.cc
   fe_nothing.cc
   fe_poly.cc
   fe_poly_tensor.cc
@@ -94,6 +95,7 @@ SET(_inst
   fe_face.inst.in
   fe.inst.in
   fe_nedelec.inst.in
+  fe_nedelec_sz.inst.in
   fe_nothing.inst.in
   fe_poly.inst.in
   fe_poly_tensor.inst.in
diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc
new file mode 100644 (file)
index 0000000..ffe25cf
--- /dev/null
@@ -0,0 +1,1839 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/fe/fe_nedelec_sz.h>
+
+#include <deal.II/base/std_cxx14/memory.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Constructor:
+template <int dim>
+FE_NedelecSZ<dim>::FE_NedelecSZ (const unsigned int degree)
+  :
+  FiniteElement<dim,dim>
+  (FiniteElementData<dim> (get_dpo_vector (degree), dim, degree + 1,
+                           FiniteElementData<dim>::Hcurl),
+   std::vector<bool> (compute_num_dofs(degree), true),
+   std::vector<ComponentMask> (compute_num_dofs(degree),
+                               std::vector<bool> (dim, true) ) )
+{
+  Assert (dim >= 2, ExcImpossibleInDim(dim));
+
+  this->mapping_type = mapping_nedelec;
+  // Set up the table converting components to base components. Since we have only
+  // one base element, everything remains zero except the component in the base,
+  // which is the component itself.
+  for (unsigned int comp=0; comp<this->n_components() ; ++comp)
+    {
+      this->component_to_base_table[comp].first.second = comp;
+    }
+
+  // Generate the 1-D polynomial basis.
+  create_polynomials(degree);
+}
+
+// Shape functions:
+template <int dim>
+double
+FE_NedelecSZ<dim>::shape_value (const unsigned int /*i*/,
+                                const Point<dim>  &/*p*/) const
+{
+  typedef    FiniteElement<dim,dim> FEE;
+  Assert(false, typename FEE::ExcFENotPrimitive());
+  return 0.;
+}
+
+
+
+template <int dim>
+double
+FE_NedelecSZ<dim>::shape_value_component (const unsigned int /*i*/,
+                                          const Point<dim>  &/*p*/,
+                                          const unsigned int /*component*/) const
+{
+  // Not implemented yet:
+  Assert (false, ExcNotImplemented ());
+  return 0.;
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_NedelecSZ<dim>::shape_grad (const unsigned int /*i*/,
+                               const Point<dim>  &/*p*/) const
+{
+  typedef    FiniteElement<dim,dim> FEE;
+  Assert(false, typename FEE::ExcFENotPrimitive());
+  return Tensor<1,dim>();
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_NedelecSZ<dim>::shape_grad_component (const unsigned int /*i*/,
+                                         const Point<dim>  &/*p*/,
+                                         const unsigned int /*component*/) const
+{
+  Assert (false, ExcNotImplemented ());
+  return Tensor<1,dim> ();
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_NedelecSZ<dim>::shape_grad_grad (const unsigned int /*i*/,
+                                    const Point<dim>  &/*p*/) const
+{
+  typedef    FiniteElement<dim,dim> FEE;
+  Assert(false, typename FEE::ExcFENotPrimitive());
+  return Tensor<2,dim>();
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_NedelecSZ<dim>::shape_grad_grad_component (const unsigned int /*i*/,
+                                              const Point<dim>  &/*p*/,
+                                              const unsigned int /*component*/) const
+{
+  Assert (false, ExcNotImplemented ());
+  return Tensor<2,dim>();
+}
+
+
+
+template <int dim>
+std::unique_ptr<typename FiniteElement<dim,dim>::InternalDataBase>
+FE_NedelecSZ<dim>::get_data (
+  const UpdateFlags                                                             update_flags,
+  const Mapping<dim,dim>                                                       &/*mapping*/,
+  const Quadrature<dim>                                                        &quadrature,
+  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &/*output_data*/) const
+{
+  std::unique_ptr<InternalData> data(new InternalData);
+  data->update_each = update_each(update_flags) | update_once(update_flags);
+
+  // Useful quantities:
+  const unsigned int degree (this->degree-1); // Note: FE holds input degree+1
+
+  const unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
+  const unsigned int lines_per_cell = GeometryInfo<dim>::lines_per_cell;
+  const unsigned int faces_per_cell = GeometryInfo<dim>::faces_per_cell;
+
+  const unsigned int n_line_dofs = this->dofs_per_line*lines_per_cell;
+  const unsigned int n_face_dofs = this->dofs_per_quad*faces_per_cell;
+
+  const UpdateFlags flags(data->update_each);
+  const unsigned int n_q_points = quadrature.size();
+
+  // Resize the internal data storage:
+  data->sigma_imj_values.resize(
+    n_q_points,
+    std::vector<std::vector<double>> (vertices_per_cell,
+                                      std::vector<double> (vertices_per_cell)));
+
+  data->sigma_imj_grads.resize(
+    vertices_per_cell,
+    std::vector<std::vector<double>> (vertices_per_cell,
+                                      std::vector<double> (dim)));
+
+  // Resize shape function arrays according to update flags:
+  if (flags & update_values)
+    {
+      data->shape_values.resize(
+        this->dofs_per_cell, std::vector<Tensor<1,dim>> (n_q_points));
+    }
+
+  if (flags & update_gradients)
+    {
+      data->shape_grads.resize(
+        this->dofs_per_cell, std::vector<DerivativeForm<1, dim, dim> > (n_q_points));
+    }
+  // Not implementing second derivatives yet:
+  if (flags & update_hessians)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+  std::vector< Point<dim> > p_list(n_q_points);
+  p_list = quadrature.get_points();
+
+
+  switch (dim)
+    {
+    case 2:
+    {
+      // Compute values of sigma & lambda and the sigma differences and lambda additions.
+      std::vector<std::vector<double>> sigma(n_q_points, std::vector<double> (lines_per_cell));
+      std::vector<std::vector<double>> lambda(n_q_points, std::vector<double> (lines_per_cell));
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
+          sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
+          sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
+          sigma[q][3] = p_list[q][0] + p_list[q][1];
+
+          lambda[q][0] = (1.0 - p_list[q][0])*(1.0 - p_list[q][1]);
+          lambda[q][1] = p_list[q][0]*(1.0 - p_list[q][1]);
+          lambda[q][2] = (1.0 - p_list[q][0])*p_list[q][1];
+          lambda[q][3] = p_list[q][0]*p_list[q][1];
+          for (unsigned int i=0; i<vertices_per_cell; ++i)
+            {
+              for (unsigned int j=0; j<vertices_per_cell; ++j)
+                {
+                  data->sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
+                }
+            }
+        }
+
+      // Calulate the gradient of sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j]
+      //   - this depends on the component and the direction of the corresponding edge.
+      //   - the direction of the edge is determined by sigma_imj_sign[i][j].
+      // Helper arrays:
+      const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2]
+      = { {-1,-1}, {1,-1}, {-1, 1}, {1, 1} };
+      int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
+      unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell];
+
+      for (unsigned int i=0; i<vertices_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<vertices_per_cell; ++j)
+            {
+              // sigma_imj_sign is the sign (+/-) of the coefficient of x/y/z in sigma_imj_values
+              // Due to the numbering of vertices on the reference element it is easy to find
+              // edges in the positive direction are from smaller to higher local vertex numbering.
+              sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
+              sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
+
+              // Now store the component which the sigma_i - sigma_j corresponds to:
+              sigma_imj_component[i][j] = 0;
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
+                  // Only interested in the first non-zero
+                  // as if there is a second, it can not be a valid edge.
+                  if (temp_imj != 0)
+                    {
+                      sigma_imj_component[i][j] = d;
+                      break;
+                    }
+                }
+              // Can now calculate the gradient, only non-zero in the component given:
+              // Note some i,j combinations will be incorrect, but only on invalid edges.
+              data->sigma_imj_grads[i][j][sigma_imj_component[i][j]]
+                = 2.0*(double)sigma_imj_sign[i][j];
+            }
+        }
+
+      // Now compute the edge parameterisations for a single element
+      // with global numbering matching that of the reference element:
+
+      // Resize the edge parameterisations
+      data->edge_sigma_values.resize(lines_per_cell);
+      data->edge_sigma_grads.resize(lines_per_cell);
+      for (unsigned int m=0; m<lines_per_cell; ++m)
+        {
+          data->edge_sigma_values[m].resize(n_q_points);
+
+          // sigma grads are constant in a cell (no need for quad points)
+          data->edge_sigma_grads[m].resize(dim);
+        }
+
+      // Fill the values for edge lambda and edge sigma:
+      const unsigned int
+      edge_sigma_direction[GeometryInfo<2>::lines_per_cell]
+        = {1, 1, 0, 0};
+
+      data->edge_lambda_values.resize(lines_per_cell, std::vector<double> (n_q_points));
+      data->edge_lambda_grads_2d.resize(lines_per_cell, std::vector<double> (dim));
+      for (unsigned int m=0; m<lines_per_cell; ++m)
+        {
+          // e1=max(reference vertex numbering on this edge)
+          // e2=min(reference vertex numbering on this edge)
+          // Which is guaranteed to be:
+          const unsigned int e1 (GeometryInfo<dim>::line_to_cell_vertices(m,1));
+          const unsigned int e2 (GeometryInfo<dim>::line_to_cell_vertices(m,0));
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              data->edge_sigma_values[m][q] = data->sigma_imj_values[q][e2][e1];
+              data->edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
+            }
+
+          data->edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
+        }
+      data->edge_lambda_grads_2d[0] = {-1.0, 0.0};
+      data->edge_lambda_grads_2d[1] = {1.0, 0.0};
+      data->edge_lambda_grads_2d[2] = {0.0, -1.0};
+      data->edge_lambda_grads_2d[3] = {0.0, 1.0};
+
+      // If the polynomial order is 0, then no more work to do:
+      if (degree < 1)
+        {
+          break;
+        }
+
+      // Otherwise, we can compute the non-cell dependent shape functions.
+      //
+      // Note: the local dof numberings follow the usual order of lines -> faces -> cells
+      //       (we have no vertex-based DoFs in this element).
+      // For a given cell we have:
+      //      n_line_dofs = dofs_per_line*lines_per_cell.
+      //      n_face_dofs = dofs_per_face*faces_per_cell.
+      //      n_cell_dofs = dofs_per_quad (2D)
+      //                  = dofs_per_hex (3D)
+      //
+      // i.e. For the local dof numbering:
+      //      the first line dof is 0,
+      //      the first face dof is n_line_dofs,
+      //      the first cell dof is n_line_dofs + n_face_dofs.
+      //
+      // On a line, DoFs are ordered first by line_dof and then line_index:
+      // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
+      //
+      // and similarly for faces:
+      // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
+      //
+      // HOWEVER, we have different types of DoFs on a line/face/cell.
+      // On a line we have two types, lowest order and higher order gradients.
+      //    - The numbering is such the lowest order is first, then higher order.
+      //      This is simple enough as there is only 1 lowest order and degree higher
+      //      orders DoFs per line.
+      //
+      // On a 2D cell, we have 3 types: Type 1/2/3:
+      //    - The ordering done by type:
+      //      - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
+      //        Numbered: ij1 = i1 + j1*(degree).        i.e. cell_dof_index = ij1.
+      //      - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
+      //        Numbered: ij2 = i2 + j2*(degree).        i.e. cell_dof_index = degree^2 + ij2
+      //      - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
+      //        Numbered: ij3 = i3.                      i.e. cell_dof_index =  2*(degree^2) + ij3.
+      //
+      // These then fit into the local dof numbering described above:
+      // - local dof numberings are:
+      //   line_dofs: local_dof = line_dof_index.    0 <= local_dof < dofs_per_line*lines_per_cell
+      //   face_dofs: local_dof = n_line_dofs*lines_per_cell + face_dof_index.
+      //   cell dofs: local_dof = n_lines_dof + n_face_dofs + cell_dof_index.
+      //
+      // The cell-based shape functions are:
+      //
+      // Type 1 (gradients):
+      // \phi^{C_{1}}_{ij) = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
+      //
+      // 0 <= i,j < degree.
+      //
+      // NOTE: The derivative produced by IntegratedLegendrePolynomials does not account for the
+      //       (2*x-1) or (2*y-1) so we must take this into account when taking
+      //       derivatives.
+      const unsigned int cell_type1_offset = n_line_dofs;
+
+      // Type 2:
+      // \phi^{C_{2}}_{ij) = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
+      //                     - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
+      //
+      // 0 <= i,j < degree.
+      const unsigned int cell_type2_offset = cell_type1_offset + degree*degree;
+
+      // Type 3 (two subtypes):
+      // \phi^{C_{3}}_{j)        = L_{j+2}(2y-1) \mathbf{e}_{x}
+      //
+      // \phi^{C_{3}}_{j+degree) = L_{j+2}(2x-1) \mathbf{e}_{y}
+      //
+      // 0 <= j < degree
+      const unsigned int cell_type3_offset1 = cell_type2_offset + degree*degree;
+      const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
+
+      if (flags & (update_values | update_gradients))
+        {
+          // compute all points we must evaluate the 1d polynomials at:
+          std::vector< Point<dim> > cell_points (n_q_points);
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  cell_points[q][d] = 2.0*p_list[q][d] - 1.0;
+                }
+            }
+
+          // Loop through quad points:
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              // pre-compute values & required derivatives at this quad point (x,y):
+              // polyx = L_{i+2}(2x-1),
+              // polyy = L_{j+2}(2y-1),
+              //
+              // for each polyc[d], c=x,y, contains the d-th derivative with respect to the
+              // co-ordinate c.
+
+              // We only need poly values and 1st derivative for update_values, but need the
+              // 2nd derivative too for update_gradients.
+              //
+              // Note that this will need to be updated if we're supporting update_hessians.
+              const unsigned int poly_length( (flags & update_gradients) ? 3 : 2);
+
+              std::vector<std::vector<double> > polyx(degree, std::vector<double> (poly_length));
+              std::vector<std::vector<double> > polyy(degree, std::vector<double> (poly_length));
+              for (unsigned int i=0; i<degree; ++i)
+                {
+                  // Compute all required 1d polynomials and their derivatives,
+                  // starting at degree 2. e.g. to access L'_{3}(2x-1) use polyx[1][1].
+                  IntegratedLegendrePolynomials[i+2].value(cell_points[q][0], polyx[i]);
+                  IntegratedLegendrePolynomials[i+2].value(cell_points[q][1], polyy[i]);
+                }
+              // Now use these to compute the shape functions:
+              if (flags & update_values)
+                {
+                  for (unsigned int j=0; j<degree; ++j)
+                    {
+                      const unsigned int shift_j(j*degree);
+                      for (unsigned int i=0; i<degree; ++i)
+                        {
+                          const unsigned int shift_ij (i+shift_j);
+
+                          // Type 1:
+                          const unsigned int dof_index1 (cell_type1_offset + shift_ij);
+                          data->shape_values[dof_index1][q][0] = 2.0*polyx[i][1]*polyy[j][0];
+                          data->shape_values[dof_index1][q][1] = 2.0*polyx[i][0]*polyy[j][1];
+
+                          // Type 2:
+                          const unsigned int dof_index2 (cell_type2_offset + shift_ij);
+                          data->shape_values[dof_index2][q][0] =      data->shape_values[dof_index1][q][0];
+                          data->shape_values[dof_index2][q][1] = -1.0*data->shape_values[dof_index1][q][1];
+                        }
+                      // Type 3:
+                      const unsigned int dof_index3_1 (cell_type3_offset1 + j);
+                      data->shape_values[dof_index3_1][q][0] = polyy[j][0];
+                      data->shape_values[dof_index3_1][q][1] = 0.0;
+
+                      const unsigned int dof_index3_2 (cell_type3_offset2 + j);
+                      data->shape_values[dof_index3_2][q][0] = 0.0;
+                      data->shape_values[dof_index3_2][q][1] = polyx[j][0];
+                    }
+                }
+              if (flags & update_gradients)
+                {
+                  for (unsigned int j=0; j<degree; ++j)
+                    {
+                      const unsigned int shift_j(j*degree);
+                      for (unsigned int i=0; i<degree; ++i)
+                        {
+                          const unsigned int shift_ij (i+shift_j);
+
+                          // Type 1:
+                          const unsigned int dof_index1 (cell_type1_offset + shift_ij);
+                          data->shape_grads[dof_index1][q][0][0] = 4.0*polyx[i][2]*polyy[j][0];
+                          data->shape_grads[dof_index1][q][0][1] = 4.0*polyx[i][1]*polyy[j][1];
+                          data->shape_grads[dof_index1][q][1][0] = data->shape_grads[dof_index1][q][0][1];
+                          data->shape_grads[dof_index1][q][1][1] = 4.0*polyx[i][0]*polyy[j][2];
+
+                          // Type 2:
+                          const unsigned int dof_index2 (cell_type2_offset + shift_ij);
+                          data->shape_grads[dof_index2][q][0][0] = data->shape_grads[dof_index1][q][0][0];
+                          data->shape_grads[dof_index2][q][0][1] = data->shape_grads[dof_index1][q][0][1];
+                          data->shape_grads[dof_index2][q][1][0] = -1.0*data->shape_grads[dof_index1][q][1][0];
+                          data->shape_grads[dof_index2][q][1][1] = -1.0*data->shape_grads[dof_index1][q][1][1];
+                        }
+                      // Type 3:
+                      const unsigned int dof_index3_1 (cell_type3_offset1 + j);
+                      data->shape_grads[dof_index3_1][q][0][0] = 0.0;
+                      data->shape_grads[dof_index3_1][q][0][1] = 2.0*polyy[j][1];
+                      data->shape_grads[dof_index3_1][q][1][0] = 0.0;
+                      data->shape_grads[dof_index3_1][q][1][1] = 0.0;
+
+                      const unsigned int dof_index3_2 (cell_type3_offset2 + j);
+                      data->shape_grads[dof_index3_2][q][0][0] = 0.0;
+                      data->shape_grads[dof_index3_2][q][0][1] = 0.0;
+                      data->shape_grads[dof_index3_2][q][1][0] = 2.0*polyx[j][1];
+                      data->shape_grads[dof_index3_2][q][1][1] = 0.0;
+                    }
+                }
+            }
+        }
+      break;
+    }
+    case 3:
+    {
+      // Compute values of sigma & lambda and the sigma differences and lambda additions.
+      std::vector<std::vector<double>> sigma(n_q_points, std::vector<double> (lines_per_cell));
+      std::vector<std::vector<double>> lambda(n_q_points, std::vector<double> (lines_per_cell));
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
+          sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
+          sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
+          sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
+          sigma[q][4] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
+          sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
+          sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
+          sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
+
+          lambda[q][0] = (1.0 - p_list[q][0])*(1.0 - p_list[q][1])*(1.0 - p_list[q][2]);
+          lambda[q][1] = p_list[q][0]*(1.0 - p_list[q][1])*(1.0 - p_list[q][2]);
+          lambda[q][2] = (1.0 - p_list[q][0])*p_list[q][1]*(1.0 - p_list[q][2]);
+          lambda[q][3] = p_list[q][0]*p_list[q][1]*(1.0 - p_list[q][2]);
+          lambda[q][4] = (1.0 - p_list[q][0])*(1.0 - p_list[q][1])*p_list[q][2];
+          lambda[q][5] = p_list[q][0]*(1.0 - p_list[q][1])*p_list[q][2];
+          lambda[q][6] = (1.0 - p_list[q][0])*p_list[q][1]*p_list[q][2];
+          lambda[q][7] = p_list[q][0]*p_list[q][1]*p_list[q][2];
+
+          // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
+          // and lambda_ipj = \lambda_{i} + \lambda_{j}.
+          for (unsigned int i=0; i<vertices_per_cell; ++i)
+            {
+              for (unsigned int j=0; j<vertices_per_cell; ++j)
+                {
+                  data->sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
+                }
+            }
+        }
+
+      // We now want some additional information about sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j]
+      // In order to calculate values & derivatives of the shape functions we need to know:
+      // - The component the sigma_imj value corresponds to - this varies with i & j.
+      // - The gradient of the sigma_imj value
+      //   - this depends on the component and the direction of the corresponding edge.
+      //   - the direction of the edge is determined by sigma_imj_sign[i][j].
+      //
+      // Note that not every i,j combination is a valid edge (there are only 12 valid edges in 3D),
+      // but we compute them all as it simplifies things.
+
+      // store the sign of each component x, y, z in the sigma list.
+      // can use this to fill in the sigma_imj_component data.
+      const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3]
+      = { {-1, -1, -1}, {1, -1, -1}, {-1, 1, -1}, {1, 1, -1},
+        {-1, -1, 1}, {1, -1, 1}, {-1, 1, 1}, {1, 1, 1}
+      };
+
+      int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
+      unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell];
+
+      for (unsigned int i=0; i<vertices_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<vertices_per_cell; ++j)
+            {
+              // sigma_imj_sign is the sign (+/-) of the coefficient of x/y/z in sigma_imj.
+              // Due to the numbering of vertices on the reference element
+              // this is easy to work out because edges in the positive direction
+              // go from smaller to higher local vertex numbering.
+              sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
+              sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
+
+              // Now store the component which the sigma_i - sigma_j corresponds to:
+              sigma_imj_component[i][j] = 0;
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
+                  // Only interested in the first non-zero
+                  // as if there is a second, it will not be a valid edge.
+                  if (temp_imj !=0)
+                    {
+                      sigma_imj_component[i][j] = d;
+                      break;
+                    }
+                }
+              // Can now calculate the gradient, only non-zero in the component given:
+              // Note some i,j combinations will be incorrect, but only on invalid edges.
+              data->sigma_imj_grads[i][j][sigma_imj_component[i][j]] = 2.0*(double)sigma_imj_sign[i][j];
+            }
+        }
+      // Now compute the edge parameterisations for a single element
+      // with global numbering matching that of the reference element:
+
+      // resize the edge parameterisations
+      data->edge_sigma_values.resize(lines_per_cell);
+      data->edge_lambda_values.resize(lines_per_cell);
+      data->edge_sigma_grads.resize(lines_per_cell);
+      data->edge_lambda_grads_3d.resize(lines_per_cell);
+      data->edge_lambda_gradgrads_3d.resize(lines_per_cell);
+      for (unsigned int m=0; m<lines_per_cell; ++m)
+        {
+          data->edge_sigma_values[m].resize(n_q_points);
+          data->edge_lambda_values[m].resize(n_q_points);
+
+          // sigma grads are constant in a cell (no need for quad points)
+          data->edge_sigma_grads[m].resize(dim);
+
+          data->edge_lambda_grads_3d[m].resize(n_q_points);
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              data->edge_lambda_grads_3d[m][q].resize(dim);
+            }
+          //lambda_gradgrads are constant in a cell (no need for quad points)
+          data->edge_lambda_gradgrads_3d[m].resize(dim);
+          for (unsigned int d=0; d<dim; ++d)
+            {
+              data->edge_lambda_gradgrads_3d[m][d].resize(dim);
+            }
+        }
+
+      // Fill the values:
+      const unsigned int
+      edge_sigma_direction[GeometryInfo<3>::lines_per_cell]
+        = {1, 1, 0, 0,
+           1, 1, 0, 0,
+           2, 2, 2, 2
+          };
+
+      for (unsigned int m=0; m<lines_per_cell; ++m)
+        {
+          // e1=max(reference vertex numbering on this edge)
+          // e2=min(reference vertex numbering on this edge)
+          // Which is guaranteed to be:
+          const unsigned int e1 (GeometryInfo<dim>::line_to_cell_vertices(m,1));
+          const unsigned int e2 (GeometryInfo<dim>::line_to_cell_vertices(m,0));
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              data->edge_sigma_values[m][q] = data->sigma_imj_values[q][e2][e1];
+              data->edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
+            }
+
+          data->edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
+        }
+      // edge_lambda_grads
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          double x(p_list[q][0]);
+          double y(p_list[q][1]);
+          double z(p_list[q][2]);
+          data->edge_lambda_grads_3d[0][q] = {z-1.0, 0.0, x-1.0};
+          data->edge_lambda_grads_3d[1][q] = {1.0-z, 0.0, -x };
+          data->edge_lambda_grads_3d[2][q] = {0.0, z-1.0, y-1.0};
+          data->edge_lambda_grads_3d[3][q] = {0.0, 1.0-z, -y};
+          data->edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0-x};
+          data->edge_lambda_grads_3d[5][q] = {z, 0.0, x};
+          data->edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0-y};
+          data->edge_lambda_grads_3d[7][q] = {0.0, z, y};
+          data->edge_lambda_grads_3d[8][q] = {y-1.0, x-1.0, 0.0};
+          data->edge_lambda_grads_3d[9][q] = {1.0-y, -x, 0.0};
+          data->edge_lambda_grads_3d[10][q] = {-y, 1.0-x, 0.0};
+          data->edge_lambda_grads_3d[11][q] = {y, x, 0.0};
+        }
+      // edge_lambda gradgrads:
+      const int
+      edge_lambda_sign[GeometryInfo<3>::lines_per_cell]
+        = {1, -1, 1, -1,
+           -1, 1, -1, 1,
+           1, -1, -1, 1
+          }; // sign of the 2nd derivative for each edge.
+      const unsigned int
+      edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2]
+      = { {0,2}, {0,2}, {1,2}, {1,2},
+        {0,2}, {0,2}, {1,2}, {1,2},
+        {0,1}, {0,1}, {0,1}, {0,1}
+      }; // component which edge_lambda[m] depends on.
+      for (unsigned int m=0; m<lines_per_cell; ++m)
+        {
+          data->edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]][edge_lambda_directions[m][1]]
+            = (double) edge_lambda_sign[m]*1.0;
+          data->edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]][edge_lambda_directions[m][0]]
+            = (double) edge_lambda_sign[m]*1.0;
+        }
+      // Precomputation for higher order shape functions,
+      // and the face parameterisation.
+      if (degree > 0)
+        {
+          // resize required data:
+          data->face_lambda_values.resize(faces_per_cell);
+          data->face_lambda_grads.resize(faces_per_cell);
+          // for face-based shape functions:
+          for (unsigned int m=0; m < faces_per_cell; ++m)
+            {
+              data->face_lambda_values[m].resize(n_q_points);
+              data->face_lambda_grads[m].resize(3);
+            }
+          // Fill in the values (these don't change between cells).
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              double x(p_list[q][0]);
+              double y(p_list[q][1]);
+              double z(p_list[q][2]);
+              data->face_lambda_values[0][q] = 1.0-x;
+              data->face_lambda_values[1][q] = x;
+              data->face_lambda_values[2][q] = 1.0-y;
+              data->face_lambda_values[3][q] = y;
+              data->face_lambda_values[4][q] = 1.0-z;
+              data->face_lambda_values[5][q] = z;
+            }
+          // gradients are constant:
+          data->face_lambda_grads[0] = {-1.0, 0.0, 0.0};
+          data->face_lambda_grads[1] = {1.0, 0.0, 0.0};
+          data->face_lambda_grads[2] = {0.0, -1.0, 0.0};
+          data->face_lambda_grads[3] = {0.0, 1.0, 0.0};
+          data->face_lambda_grads[4] = {0.0, 0.0, -1.0};
+          data->face_lambda_grads[5] = {0.0, 0.0, 1.0};
+
+          // for cell-based shape functions:
+          // these don't depend on the cell, so can precompute all here:
+          if (flags & (update_values | update_gradients))
+            {
+              // Cell-based shape functions:
+              //
+              // Type-1 (gradients):
+              // \phi^{C_{1}}_{ijk} = grad( L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
+              //
+              // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
+              const unsigned int cell_type1_offset(n_line_dofs + n_face_dofs);
+              // Type-2:
+              //
+              // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
+              // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1, -1)\phi^{C_{1}}_{ijk}
+              //
+              // 0 <= i,j,k < degree. (subtypes in groups of degree*degree*degree)
+              //
+              // here we order so that all of subtype 1 comes first, then subtype 2.
+              const unsigned int cell_type2_offset1(cell_type1_offset + degree*degree*degree);
+              const unsigned int cell_type2_offset2(cell_type2_offset1 + degree*degree*degree);
+              // Type-3
+              // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
+              // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
+              // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
+              //
+              // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
+              //
+              // again we order so we compute all of subtype 1 first, then subtype 2, etc.
+              const unsigned int cell_type3_offset1(cell_type2_offset2 + degree*degree*degree);
+              const unsigned int cell_type3_offset2(cell_type3_offset1 + degree*degree);
+              const unsigned int cell_type3_offset3(cell_type3_offset2 + degree*degree);
+
+              // compute all points we must evaluate the 1d polynomials at:
+              std::vector< Point<dim> > cell_points (n_q_points);
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  for (unsigned int d=0; d<dim; ++d)
+                    {
+                      cell_points[q][d] = 2.0*p_list[q][d] - 1.0;
+                    }
+                }
+
+
+              // only need poly values and 1st derivative for update_values,
+              // but need 2nd derivative too for update_gradients.
+              const unsigned int poly_length( (flags & update_gradients) ? 3 : 2);
+              // Loop through quad points:
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  // pre-compute values & required derivatives at this quad point, (x,y,z):
+                  // polyx = L_{i+2}(2x-1),
+                  // polyy = L_{j+2}(2y-1),
+                  // polyz = L_{k+2}(2z-1).
+                  //
+                  // for each polyc[d], c=x,y,z, contains the d-th derivative with respect to the
+                  // co-ordinate c.
+                  std::vector<std::vector<double> > polyx(degree, std::vector<double> (poly_length));
+                  std::vector<std::vector<double> > polyy(degree, std::vector<double> (poly_length));
+                  std::vector<std::vector<double> > polyz(degree, std::vector<double> (poly_length));
+                  for (unsigned int i=0; i<degree; ++i)
+                    {
+                      // compute all required 1d polynomials for i
+                      IntegratedLegendrePolynomials[i+2].value(cell_points[q][0], polyx[i]);
+                      IntegratedLegendrePolynomials[i+2].value(cell_points[q][1], polyy[i]);
+                      IntegratedLegendrePolynomials[i+2].value(cell_points[q][2], polyz[i]);
+                    }
+                  // Now use these to compute the shape functions:
+                  if (flags & update_values)
+                    {
+                      for (unsigned int k=0; k<degree; ++k)
+                        {
+                          const unsigned int shift_k(k*degree*degree);
+                          const unsigned int shift_j(k*degree); // Used below when subbing k for j (type 3)
+                          for (unsigned int j=0; j<degree; ++j)
+                            {
+                              const unsigned int shift_jk(j*degree + shift_k);
+                              for (unsigned int i=0; i<degree; ++i)
+                                {
+                                  const unsigned int shift_ijk(shift_jk + i);
+
+                                  // Type 1:
+                                  const unsigned int dof_index1(cell_type1_offset + shift_ijk);
+
+                                  data->shape_values[dof_index1][q][0] = 2.0*polyx[i][1]*polyy[j][0]*polyz[k][0];
+                                  data->shape_values[dof_index1][q][1] = 2.0*polyx[i][0]*polyy[j][1]*polyz[k][0];
+                                  data->shape_values[dof_index1][q][2] = 2.0*polyx[i][0]*polyy[j][0]*polyz[k][1];
+
+                                  // Type 2:
+                                  const unsigned int dof_index2_1(cell_type2_offset1 + shift_ijk);
+                                  const unsigned int dof_index2_2(cell_type2_offset2 + shift_ijk);
+
+                                  data->shape_values[dof_index2_1][q][0] = data->shape_values[dof_index1][q][0];
+                                  data->shape_values[dof_index2_1][q][1] = -1.0*data->shape_values[dof_index1][q][1];
+                                  data->shape_values[dof_index2_1][q][2] = data->shape_values[dof_index1][q][2];
+
+                                  data->shape_values[dof_index2_2][q][0] = data->shape_values[dof_index1][q][0];
+                                  data->shape_values[dof_index2_2][q][1] = -1.0*data->shape_values[dof_index1][q][1];
+                                  data->shape_values[dof_index2_2][q][2] = -1.0*data->shape_values[dof_index1][q][2];
+                                }
+                              // Type 3: (note we re-use k and j for convenience):
+                              const unsigned int shift_ij(j+shift_j); // here we've subbed j for i, k for j.
+                              const unsigned int dof_index3_1(cell_type3_offset1 + shift_ij);
+                              const unsigned int dof_index3_2(cell_type3_offset2 + shift_ij);
+                              const unsigned int dof_index3_3(cell_type3_offset3 + shift_ij);
+
+                              data->shape_values[dof_index3_1][q][0] = polyy[j][0]*polyz[k][0];
+                              data->shape_values[dof_index3_1][q][1] = 0.0;
+                              data->shape_values[dof_index3_1][q][2] = 0.0;
+
+                              data->shape_values[dof_index3_2][q][0] = 0.0;
+                              data->shape_values[dof_index3_2][q][1] = polyx[j][0]*polyz[k][0];
+                              data->shape_values[dof_index3_2][q][2] = 0.0;
+
+                              data->shape_values[dof_index3_3][q][0] = 0.0;
+                              data->shape_values[dof_index3_3][q][1] = 0.0;
+                              data->shape_values[dof_index3_3][q][2] = polyx[j][0]*polyy[k][0];
+
+                            }
+                        }
+                    }
+                  if (flags & update_gradients)
+                    {
+                      for (unsigned int k=0; k<degree; ++k)
+                        {
+                          const unsigned int shift_k(k*degree*degree);
+                          const unsigned int shift_j(k*degree); // Used below when subbing k for j (type 3)
+                          for (unsigned int j=0; j<degree; ++j)
+                            {
+                              const unsigned int shift_jk(j*degree + shift_k);
+                              for (unsigned int i=0; i<degree; ++i)
+                                {
+                                  const unsigned int shift_ijk(shift_jk + i);
+
+                                  // Type 1:
+                                  const unsigned int dof_index1(cell_type1_offset + shift_ijk);
+
+                                  data->shape_grads[dof_index1][q][0][0] = 4.0*polyx[i][2]*polyy[j][0]*polyz[k][0];
+                                  data->shape_grads[dof_index1][q][0][1] = 4.0*polyx[i][1]*polyy[j][1]*polyz[k][0];
+                                  data->shape_grads[dof_index1][q][0][2] = 4.0*polyx[i][1]*polyy[j][0]*polyz[k][1];
+
+                                  data->shape_grads[dof_index1][q][1][0] = data->shape_grads[dof_index1][q][0][1];
+                                  data->shape_grads[dof_index1][q][1][1] = 4.0*polyx[i][0]*polyy[j][2]*polyz[k][0];
+                                  data->shape_grads[dof_index1][q][1][2] = 4.0*polyx[i][0]*polyy[j][1]*polyz[k][1];
+
+                                  data->shape_grads[dof_index1][q][2][0] = data->shape_grads[dof_index1][q][0][2];
+                                  data->shape_grads[dof_index1][q][2][1] = data->shape_grads[dof_index1][q][1][2];
+                                  data->shape_grads[dof_index1][q][2][2] = 4.0*polyx[i][0]*polyy[j][0]*polyz[k][2];
+
+                                  // Type 2:
+                                  const unsigned int dof_index2_1(cell_type2_offset1 + shift_ijk);
+                                  const unsigned int dof_index2_2(cell_type2_offset2 + shift_ijk);
+
+                                  for (unsigned int d=0; d<dim; ++d)
+                                    {
+                                      data->shape_grads[dof_index2_1][q][0][d] = data->shape_grads[dof_index1][q][0][d];
+                                      data->shape_grads[dof_index2_1][q][1][d] = -1.0*data->shape_grads[dof_index1][q][1][d];
+                                      data->shape_grads[dof_index2_1][q][2][d] = data->shape_grads[dof_index1][q][2][d];
+
+                                      data->shape_grads[dof_index2_2][q][0][d] = data->shape_grads[dof_index1][q][0][d];
+                                      data->shape_grads[dof_index2_2][q][1][d] = -1.0*data->shape_grads[dof_index1][q][1][d];
+                                      data->shape_grads[dof_index2_2][q][2][d] = -1.0*data->shape_grads[dof_index1][q][2][d];
+                                    }
+                                }
+                              // Type 3: (note we re-use k and j for convenience):
+                              const unsigned int shift_ij(j+shift_j); // here we've subbed j for i, k for j.
+                              const unsigned int dof_index3_1(cell_type3_offset1 + shift_ij);
+                              const unsigned int dof_index3_2(cell_type3_offset2 + shift_ij);
+                              const unsigned int dof_index3_3(cell_type3_offset3 + shift_ij);
+                              for (unsigned int d1=0; d1<dim; ++d1)
+                                {
+                                  for (unsigned int d2=0; d2<dim; ++d2)
+                                    {
+                                      data->shape_grads[dof_index3_1][q][d1][d2] = 0.0;
+                                      data->shape_grads[dof_index3_2][q][d1][d2] = 0.0;
+                                      data->shape_grads[dof_index3_3][q][d1][d2] = 0.0;
+                                    }
+                                }
+                              data->shape_grads[dof_index3_1][q][0][1] = 2.0*polyy[j][1]*polyz[k][0];
+                              data->shape_grads[dof_index3_1][q][0][2] = 2.0*polyy[j][0]*polyz[k][1];
+
+                              data->shape_grads[dof_index3_2][q][1][0] = 2.0*polyx[j][1]*polyz[k][0];
+                              data->shape_grads[dof_index3_2][q][1][2] = 2.0*polyx[j][0]*polyz[k][1];
+
+                              data->shape_grads[dof_index3_3][q][2][0] = 2.0*polyx[j][1]*polyy[k][0];
+                              data->shape_grads[dof_index3_3][q][2][1] = 2.0*polyx[j][0]*polyy[k][1];
+                            }
+                        }
+                    }
+                }
+            }
+        }
+      break;
+    }
+    default:
+    {
+      Assert (false, ExcNotImplemented ());
+    }
+    }
+  return data;
+}
+
+template <int dim>
+void FE_NedelecSZ<dim>::fill_edge_values(const typename Triangulation<dim,dim>::cell_iterator &cell,
+                                         const Quadrature<dim>                                &quadrature,
+                                         const InternalData                                   &fe_data) const
+{
+  // This function handles the cell-dependent construction of the EDGE-based shape functions.
+  //
+  // Note it will handle both 2D and 3D, in 2D, the edges are faces, but we handle them here.
+  //
+  // It will fill in the missing parts of fe_data which were not possible to fill
+  // in the get_data routine, with respect to the edge-based shape functions.
+  //
+  // It should be called by the fill_fe_*_values routines in order to complete the
+  // basis set at quadrature points on the current cell for each edge.
+
+  const UpdateFlags flags(fe_data.update_each);
+  const unsigned int n_q_points = quadrature.size();
+
+  Assert(!(flags & update_values) || fe_data.shape_values.size() == this->dofs_per_cell,
+         ExcDimensionMismatch(fe_data.shape_values.size(), this->dofs_per_cell));
+  Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_q_points,
+         ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
+
+  // Useful constants:
+  const unsigned int degree (this->degree-1); // Note: constructor takes input degree + 1, so need to knock 1 off.
+
+  // Useful geometry info:
+  const unsigned int vertices_per_line (2);
+  const unsigned int lines_per_cell (GeometryInfo<dim>::lines_per_cell);
+
+  // Calculate edge orderings:
+  std::vector<std::vector<unsigned int>> edge_order(lines_per_cell,
+                                                    std::vector<unsigned int> (vertices_per_line));
+
+
+  switch (dim)
+    {
+    case 2:
+    {
+      if (flags & (update_values | update_gradients))
+        {
+
+          // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1}, e^{m}_{2}]
+          // e1 = higher global numbering of the two local vertices
+          // e2 = lower global numbering of the two local vertices
+          std::vector<int> edge_sign(lines_per_cell);
+          for (unsigned int m=0; m<lines_per_cell; ++m)
+            {
+              unsigned int v0_loc = GeometryInfo<dim>::line_to_cell_vertices(m,0);
+              unsigned int v1_loc = GeometryInfo<dim>::line_to_cell_vertices(m,1);
+              unsigned int v0_glob = cell->vertex_index(v0_loc);
+              unsigned int v1_glob = cell->vertex_index(v1_loc);
+
+              if (v0_glob > v1_glob)
+                {
+                  // Opposite to global numbering on our reference element
+                  edge_sign[m] = -1.0;
+                }
+              else
+                {
+                  // Aligns with global numbering on our reference element.
+                  edge_sign[m] = 1.0;
+                }
+            }
+
+          // Define \sigma_{m} = sigma_{e^{m}_{2}} - sigma_{e^{m}_{2}}
+          //        \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
+          //
+          // To help things, in fe_data, we have precomputed (sigma_{i} - sigma_{j})
+          // and (lambda_{i} + lambda_{j}) for 0<= i,j < lines_per_cell.
+          //
+          // There are two types:
+          // - lower order (1 per edge, m):
+          //   \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
+          //
+          // - higher order (degree per edge, m):
+          //   \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
+          //
+          //   NOTE: sigma_{m} and lambda_{m} are either a function of x OR y
+          //         and if sigma is of x, then lambda is of y, and vice versa.
+          //         This means that grad(\sigma) requires multiplication by d(sigma)/dx_{i}
+          //         for the i^th comp of grad(sigma) and similarly when taking derivatives of lambda.
+          //
+          // First handle the lowest order edges (dofs 0 to 3)
+          // 0 and 1 are the edges in the y dir. (sigma is function of y, lambda is function of x).
+          // 2 and 3 are the edges in the x dir. (sigma is function of x, lambda is function of y).
+          //
+          // More more info: see GeometryInfo for picture of the standard element.
+          //
+          // Fill edge-based points:
+          //      std::vector<std::vector< Point<dim> > > edge_points(lines_per_cell, std::vector<Point<dim>> (n_q_points));
+
+          std::vector<std::vector< double > > edge_sigma_values (fe_data.edge_sigma_values);
+          std::vector<std::vector< double > > edge_sigma_grads (fe_data.edge_sigma_grads);
+
+          std::vector<std::vector< double > > edge_lambda_values (fe_data.edge_lambda_values);
+          std::vector<std::vector< double > > edge_lambda_grads (fe_data.edge_lambda_grads_2d);
+
+          // Adjust the edge_sigma_* for the current cell:
+          for (unsigned int m=0; m<lines_per_cell; ++m)
+            {
+              std::transform (edge_sigma_values[m].begin(), edge_sigma_values[m].end(),
+                              edge_sigma_values[m].begin(),
+                              std::bind1st (std::multiplies<double> (), static_cast<double>(edge_sign[m])));
+
+              std::transform (edge_sigma_grads[m].begin(), edge_sigma_grads[m].end(),
+                              edge_sigma_grads[m].begin(),
+                              std::bind1st (std::multiplies<double> (), static_cast<double>(edge_sign[m])));
+            }
+
+          // If we want to generate shape gradients then we need second derivatives of the 1d polynomials,
+          // but only first derivatives for the shape values.
+          const unsigned int poly_length( (flags & update_gradients) ? 3 : 2);
+
+          for (unsigned int m=0; m<lines_per_cell; ++m)
+            {
+              const unsigned int shift_m (m*this->dofs_per_line);
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  // Only compute 1d polynomials if degree>0.
+                  std::vector<std::vector<double> > poly(degree, std::vector<double> (poly_length));
+                  for (unsigned int i=1; i<degree+1; ++i)
+                    {
+                      // Compute all required 1d polynomials and their derivatives,
+                      // starting at degree 2. e.g. to access L'_{3}(2x-1) use polyx[1][1].
+                      IntegratedLegendrePolynomials[i+1].value(edge_sigma_values[m][q], poly[i-1]);
+                    }
+                  if (flags & update_values)
+                    {
+                      // Lowest order edge shape functions:
+                      for (unsigned int d=0; d<dim; ++d)
+                        {
+                          fe_data.shape_values[shift_m][q][d] = 0.5*edge_sigma_grads[m][d]*edge_lambda_values[m][q];
+                        }
+                      // Higher order edge shape functions:
+                      for (unsigned int i=1; i<degree+1; ++i)
+                        {
+                          const unsigned int poly_index = i-1;
+                          const unsigned int dof_index(i + shift_m);
+                          for (unsigned int d=0; d<dim; ++d)
+                            {
+                              fe_data.shape_values[dof_index][q][d]
+                                = edge_sigma_grads[m][d]*poly[poly_index][1]*edge_lambda_values[m][q]
+                                  + poly[poly_index][0]*edge_lambda_grads[m][d];
+                            }
+                        }
+                    }
+                  if (flags & update_gradients)
+                    {
+                      // Lowest order edge shape functions:
+                      for (unsigned int d1=0; d1<dim; ++d1)
+                        {
+                          for (unsigned int d2=0; d2<dim; ++d2)
+                            {
+                              // Note: gradient is constant for a given edge.
+                              fe_data.shape_grads[shift_m][q][d1][d2]
+                                = 0.5*edge_sigma_grads[m][d1]*edge_lambda_grads[m][d2];
+                            }
+                        }
+                      // Higher order edge shape functions:
+                      for (unsigned int i=1; i<degree+1; ++i)
+                        {
+                          const unsigned int poly_index = i-1;
+                          const unsigned int dof_index(i + shift_m);
+
+                          fe_data.shape_grads[dof_index][q][0][0]
+                            = edge_sigma_grads[m][0]*edge_sigma_grads[m][0]*edge_lambda_values[m][q]*poly[poly_index][2];
+
+                          fe_data.shape_grads[dof_index][q][0][1]
+                            = (edge_sigma_grads[m][0]*edge_lambda_grads[m][1] + edge_sigma_grads[m][1]*edge_lambda_grads[m][0])*poly[poly_index][1];
+
+                          fe_data.shape_grads[dof_index][q][1][0] = fe_data.shape_grads[dof_index][q][0][1];
+
+                          fe_data.shape_grads[dof_index][q][1][1]
+                            = edge_sigma_grads[m][1]*edge_sigma_grads[m][1]*edge_lambda_values[m][q]*poly[poly_index][2];
+                        }
+                    }
+                }
+            }
+        }
+      break;
+    }
+    case 3:
+    {
+      if (flags & (update_values | update_gradients))
+        {
+          // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1}, e^{m}_{2}]
+          // e1 = higher global numbering of the two local vertices
+          // e2 = lower global numbering of the two local vertices
+          std::vector<int> edge_sign(lines_per_cell);
+          for (unsigned int m=0; m<lines_per_cell; ++m)
+            {
+              unsigned int v0_loc = GeometryInfo<dim>::line_to_cell_vertices(m,0);
+              unsigned int v1_loc = GeometryInfo<dim>::line_to_cell_vertices(m,1);
+              unsigned int v0_glob = cell->vertex_index(v0_loc);
+              unsigned int v1_glob = cell->vertex_index(v1_loc);
+
+              if (v0_glob > v1_glob)
+                {
+                  // Opposite to global numbering on our reference element
+                  edge_sign[m] = -1.0;
+                }
+              else
+                {
+                  // Aligns with global numbering on our reference element.
+                  edge_sign[m] = 1.0;
+                }
+            }
+
+          // Define \sigma_{m} = sigma_{e^{m}_{2}} - sigma_{e^{m}_{2}}
+          //        \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
+          //
+          // To help things, in fe_data, we have precomputed (sigma_{i} - sigma_{j})
+          // and (lambda_{i} + lambda_{j}) for 0<= i,j < lines_per_cell.
+          //
+          // There are two types:
+          // - lower order (1 per edge, m):
+          //   \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
+          //
+          // - higher order (degree per edge, m):
+          //   \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
+          //
+          //   NOTE: In the ref cell, sigma_{m} is a function of x OR y OR Z and lambda_{m} a function of the remaining co-ords.
+          //         for example, if sigma is of x, then lambda is of y AND z, and so on.
+          //         This means that grad(\sigma) requires multiplication by d(sigma)/dx_{i}
+          //         for the i^th comp of grad(sigma) and similarly when taking derivatives of lambda.
+          //
+          // First handle the lowest order edges (dofs 0 to 11)
+          // 0 and 1 are the edges in the y dir at z=0. (sigma is a fn of y, lambda is a fn of x & z).
+          // 2 and 3 are the edges in the x dir at z=0. (sigma is a fn of x, lambda is a fn of y & z).
+          // 4 and 5 are the edges in the y dir at z=1. (sigma is a fn of y, lambda is a fn of x & z).
+          // 6 and 7 are the edges in the x dir at z=1. (sigma is a fn of x, lambda is a fn of y & z).
+          // 8 and 9 are the edges in the z dir at y=0. (sigma is a fn of z, lambda is a fn of x & y).
+          // 10 and 11 are the edges in the z dir at y=1. (sigma is a fn of z, lambda is a fn of x & y).
+          //
+          // For more info: see GeometryInfo for picture of the standard element.
+
+          // Copy over required edge-based data:
+          std::vector<std::vector< double > > edge_sigma_values(fe_data.edge_sigma_values);
+          std::vector<std::vector< double > > edge_lambda_values(fe_data.edge_lambda_values);
+          std::vector<std::vector< double > > edge_sigma_grads(fe_data.edge_sigma_grads);
+          std::vector<std::vector<std::vector< double > > > edge_lambda_grads(fe_data.edge_lambda_grads_3d);
+          std::vector<std::vector<std::vector< double > > > edge_lambda_gradgrads_3d(fe_data.edge_lambda_gradgrads_3d);
+
+          // Adjust the edge_sigma_* for the current cell:
+          for (unsigned int m=0; m<lines_per_cell; ++m)
+            {
+              std::transform (edge_sigma_values[m].begin(), edge_sigma_values[m].end(),
+                              edge_sigma_values[m].begin(),
+                              std::bind1st (std::multiplies<double> (), static_cast<double>(edge_sign[m])));
+              std::transform (edge_sigma_grads[m].begin(), edge_sigma_grads[m].end(),
+                              edge_sigma_grads[m].begin(),
+                              std::bind1st (std::multiplies<double> (), static_cast<double>(edge_sign[m])));
+            }
+
+          // Now calculate the edge-based shape functions:
+          // If we want to generate shape gradients then we need second derivatives of the 1d polynomials,
+          // but only first derivatives for the shape values.
+          const unsigned int poly_length( (flags & update_gradients) ? 3 : 2);
+          std::vector<std::vector<double> > poly(degree, std::vector<double> (poly_length));
+          for (unsigned int m=0; m<lines_per_cell; ++m)
+            {
+              const unsigned int shift_m(m*this->dofs_per_line);
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  // precompute values of all 1d polynomials required:
+                  if (degree>0)
+                    {
+                      for (unsigned int i=0; i<degree; ++i)
+                        {
+                          IntegratedLegendrePolynomials[i+2].value(edge_sigma_values[m][q], poly[i]);
+                        }
+                    }
+                  if (flags & update_values)
+                    {
+                      // Lowest order shape functions:
+                      for (unsigned int d=0; d<dim; ++d)
+                        {
+                          fe_data.shape_values[shift_m][q][d] = 0.5*edge_sigma_grads[m][d]*edge_lambda_values[m][q];
+                        }
+                      if (degree>0)
+                        {
+                          for (unsigned int i=0; i<degree; ++i)
+                            {
+                              const unsigned int dof_index(i+1 + shift_m);
+                              for (unsigned int d=0; d<dim; ++d)
+                                {
+                                  fe_data.shape_values[dof_index][q][d] = edge_sigma_grads[m][d]*poly[i][1]*edge_lambda_values[m][q]
+                                                                          + poly[i][0]*edge_lambda_grads[m][q][d];
+                                }
+                            }
+                        }
+                    }
+                  if (flags & update_gradients)
+                    {
+                      // Lowest order shape functions:
+                      for (unsigned int d1=0; d1<dim; ++d1)
+                        {
+                          for (unsigned int d2=0; d2<dim; ++d2)
+                            {
+                              fe_data.shape_grads[shift_m][q][d1][d2]
+                                = 0.5*edge_sigma_grads[m][d1]*edge_lambda_grads[m][q][d2];
+                            }
+                        }
+                      if (degree>0)
+                        {
+                          for (unsigned int i=0; i<degree; ++i)
+                            {
+                              const unsigned int dof_index(i+1 + shift_m);
+
+                              for (unsigned int d1=0; d1<dim; ++d1)
+                                {
+                                  for (unsigned int d2=0; d2<dim; ++d2)
+                                    {
+                                      fe_data.shape_grads[dof_index][q][d1][d2] = edge_sigma_grads[m][d1]*edge_sigma_grads[m][d2]*poly[i][2]*edge_lambda_values[m][q]
+                                                                                  + ( edge_sigma_grads[m][d1]*edge_lambda_grads[m][q][d2]
+                                                                                      + edge_sigma_grads[m][d2]*edge_lambda_grads[m][q][d1] )*poly[i][1]
+                                                                                  + edge_lambda_gradgrads_3d[m][d1][d2]*poly[i][0];
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+      break;
+    }
+    default:
+    {
+      Assert (false, ExcNotImplemented ());
+    }
+    }
+}
+
+template <int dim>
+void FE_NedelecSZ<dim>::fill_face_values(const typename Triangulation<dim,dim>::cell_iterator &cell,
+                                         const Quadrature<dim>                                &quadrature,
+                                         const InternalData                                   &fe_data) const
+{
+  // This function handles the cell-dependent construction of the FACE-based shape functions.
+  //
+  // Note that it should only be called in 3D.
+  Assert(dim == 3, ExcDimensionMismatch(dim, 3));
+  //
+  // It will fill in the missing parts of fe_data which were not possible to fill
+  // in the get_data routine, with respect to face-based shape functions.
+  //
+  // It should be called by the fill_fe_*_values routines in order to complete the
+  // basis set at quadrature points on the current cell for each face.
+
+  // Useful constants:
+  const unsigned int degree (this->degree-1); // Note: constructor takes input degree + 1, so need to knock 1 off.
+
+  // Do nothing if FE degree is 0.
+  if (degree>0)
+    {
+      const UpdateFlags flags(fe_data.update_each);
+
+      if (flags & (update_values | update_gradients))
+        {
+          const unsigned int n_q_points = quadrature.size();
+
+          Assert(!(flags & update_values) || fe_data.shape_values.size() == this->dofs_per_cell,
+                 ExcDimensionMismatch(fe_data.shape_values.size(), this->dofs_per_cell));
+          Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_q_points,
+                 ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
+
+          // Useful geometry info:
+          const unsigned int vertices_per_face (GeometryInfo<dim>::vertices_per_face);
+          const unsigned int faces_per_cell (GeometryInfo<dim>::faces_per_cell);
+
+          // DoF info:
+          const unsigned int n_line_dofs = this->dofs_per_line*GeometryInfo<dim>::lines_per_cell;
+
+          // First we find the global face orientations on the current cell.
+          std::vector<std::vector<unsigned int> > face_orientation(faces_per_cell,
+                                                                   std::vector<unsigned int> (vertices_per_face));
+
+          const unsigned int vertex_opposite_on_face[GeometryInfo<3>::vertices_per_face]
+            = {3, 2, 1, 0};
+
+          const unsigned int vertices_adjacent_on_face[GeometryInfo<3>::vertices_per_face][2]
+          = {{1,2}, {0,3}, {0,3}, {1,2}};
+
+          for (unsigned int m=0; m<faces_per_cell; ++m)
+            {
+              // Find the local vertex on this face with the highest global numbering. This is f^m_0.
+              unsigned int current_max = 0;
+              unsigned int current_glob = cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(m,0));
+              for (unsigned int v=1; v<vertices_per_face; ++v)
+                {
+                  if (current_glob < cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(m,v)))
+                    {
+                      current_max = v;
+                      current_glob = cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(m,v));
+                    }
+                }
+              face_orientation[m][0] = GeometryInfo<dim>::face_to_cell_vertices(m,current_max);
+
+              // f^m_2 is the vertex opposite f^m_0.
+              face_orientation[m][2] = GeometryInfo<dim>::face_to_cell_vertices(m,vertex_opposite_on_face[current_max]);
+
+              // Finally, f^m_1 is the vertex with the greater global numbering of the remaining
+              // two local vertices. Then, f^m_3 is the other.
+              if ( cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][0]))
+                   > cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][1])) )
+                {
+                  face_orientation[m][1] = GeometryInfo<dim>::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][0]);
+                  face_orientation[m][3] = GeometryInfo<dim>::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][1]);
+                }
+              else
+                {
+                  face_orientation[m][1] = GeometryInfo<dim>::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][1]);
+                  face_orientation[m][3] = GeometryInfo<dim>::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][0]);
+                }
+            }
+
+          // Now we know the face orientation on the current cell, we can generate the parameterisation:
+          std::vector<std::vector<double> > face_xi_values(faces_per_cell,
+                                                           std::vector<double> (n_q_points));
+          std::vector<std::vector<double> > face_xi_grads(faces_per_cell,
+                                                          std::vector<double> (dim));
+          std::vector<std::vector<double> > face_eta_values(faces_per_cell,
+                                                            std::vector<double> (n_q_points));
+          std::vector<std::vector<double> > face_eta_grads(faces_per_cell,
+                                                           std::vector<double> (dim));
+
+          std::vector<std::vector<double> > face_lambda_values(faces_per_cell,
+                                                               std::vector<double> (n_q_points));
+          std::vector<std::vector<double> > face_lambda_grads(faces_per_cell,
+                                                              std::vector<double> (dim));
+          for (unsigned int m=0; m<faces_per_cell; ++m)
+            {
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  face_xi_values[m][q] = fe_data.sigma_imj_values[q][face_orientation[m][0]][face_orientation[m][1]];
+                  face_eta_values[m][q] = fe_data.sigma_imj_values[q][face_orientation[m][0]][face_orientation[m][3]];
+                  face_lambda_values[m][q] = fe_data.face_lambda_values[m][q];
+                }
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  face_xi_grads[m][d] = fe_data.sigma_imj_grads[face_orientation[m][0]][face_orientation[m][1]][d];
+                  face_eta_grads[m][d] = fe_data.sigma_imj_grads[face_orientation[m][0]][face_orientation[m][3]][d];
+
+                  face_lambda_grads[m][d] = fe_data.face_lambda_grads[m][d];
+                }
+            }
+          // Now can generate the basis
+          const unsigned int poly_length( (flags & update_gradients) ? 3 : 2);
+          std::vector<std::vector<double> > polyxi(degree, std::vector<double> (poly_length));
+          std::vector<std::vector<double> > polyeta(degree, std::vector<double> (poly_length));
+
+          // Loop through quad points:
+          for (unsigned int m=0; m<faces_per_cell; ++m)
+            {
+              const unsigned int shift_m(m*this->dofs_per_quad);
+              // Calculate the offsets for each face-based shape function:
+              //
+              // Type-1 (gradients)
+              // \phi^{F_m,1}_{ij} = \nabla( L_{i+2}(\xi_{F_{m}}) L_{j+2}(\eta_{F_{m}}) \lambda_{F_{m}} )
+              //
+              // 0 <= i,j < degree (in a group of degree*degree)
+              const unsigned int face_type1_offset(n_line_dofs + shift_m);
+              // Type-2:
+              //
+              // \phi^{F_m,2}_{ij} = ( L'_{i+2}(\xi_{F_{m}}) L_{j+2}(\eta_{F_{m}}) \nabla\xi_{F_{m}}
+              //                       - L_{i+2}(\xi_{F_{m}}) L'_{j+2}(\eta_{F_{m}}) \nabla\eta_{F_{m}} ) \lambda_{F_{m}}
+              //
+              // 0 <= i,j < degree (in a group of degree*degree)
+              const unsigned int face_type2_offset(face_type1_offset + degree*degree);
+              // Type-3:
+              //
+              // \phi^{F_m,3}_{i} = L_{i+2}(\eta_{F_{m}}) \lambda_{F_{m}} \nabla\xi_{F_{m}}
+              // \phi^{F_m,3}_{i+p} = L_{i+2}(\xi_{F_{m}}) \lambda_{F_{m}} \nabla\eta_{F_{m}}
+              //
+              // 0 <= i < degree.
+              //
+              // here we order so that all of subtype 1 comes first, then subtype 2.
+              const unsigned int face_type3_offset1(face_type2_offset + degree*degree);
+              const unsigned int face_type3_offset2(face_type3_offset1 + degree);
+
+              // Loop over all faces:
+              for (unsigned int q=0; q<n_q_points; ++q)
+                {
+                  // pre-compute values & required derivatives at this quad point:
+                  // polyxi = L_{i+2}(\xi_{F_{m}}),
+                  // polyeta = L_{j+2}(\eta_{F_{m}}),
+                  //
+                  // each polypoint[k][d], contains the dth derivative of L_{k+2} at the point \xi or \eta.
+                  // Note that this doesn't include the derivative of xi/eta via the chain rule.
+                  for (unsigned int i=0; i<degree; ++i)
+                    {
+                      // compute all required 1d polynomials:
+                      IntegratedLegendrePolynomials[i+2].value(face_xi_values[m][q], polyxi[i]);
+                      IntegratedLegendrePolynomials[i+2].value(face_eta_values[m][q], polyeta[i]);
+                    }
+                  // Now use these to compute the shape functions:
+                  if (flags & update_values)
+                    {
+                      for (unsigned int j=0; j<degree; ++j)
+                        {
+                          const unsigned int shift_j(j*degree);
+                          for (unsigned int i=0; i<degree; ++i)
+                            {
+                              const unsigned int shift_ij(shift_j + i);
+                              // Type 1:
+                              const unsigned int dof_index1(face_type1_offset + shift_ij);
+                              for (unsigned int d=0; d<dim; ++d)
+                                {
+                                  fe_data.shape_values[dof_index1][q][d] = ( face_xi_grads[m][d]*polyxi[i][1]*polyeta[j][0]
+                                                                             + face_eta_grads[m][d]*polyxi[i][0]*polyeta[j][1] )*face_lambda_values[m][q]
+                                                                           + face_lambda_grads[m][d]*polyxi[i][0]*polyeta[j][0];
+                                }
+                              // Type 2:
+                              const unsigned int dof_index2(face_type2_offset + shift_ij);
+                              for (unsigned int d=0; d<dim; ++d)
+                                {
+                                  fe_data.shape_values[dof_index2][q][d] = ( face_xi_grads[m][d]*polyxi[i][1]*polyeta[j][0]
+                                                                             - face_eta_grads[m][d]*polyxi[i][0]*polyeta[j][1] )*face_lambda_values[m][q];
+                                }
+                            }
+                          // Type 3:
+                          const unsigned int dof_index3_1(face_type3_offset1 + j);
+                          const unsigned int dof_index3_2(face_type3_offset2 + j);
+                          for (unsigned int d=0; d<dim; ++d)
+                            {
+                              fe_data.shape_values[dof_index3_1][q][d] = face_xi_grads[m][d]*polyeta[j][0]*face_lambda_values[m][q];
+
+                              fe_data.shape_values[dof_index3_2][q][d] = face_eta_grads[m][d]*polyxi[j][0]*face_lambda_values[m][q];
+                            }
+                        }
+                    }
+                  if (flags & update_gradients)
+                    {
+                      for (unsigned int j=0; j<degree; ++j)
+                        {
+                          const unsigned int shift_j(j*degree);
+                          for (unsigned int i=0; i<degree; ++i)
+                            {
+                              const unsigned int shift_ij(shift_j + i);
+                              // Type 1:
+                              const unsigned int dof_index1(face_type1_offset + shift_ij);
+                              for (unsigned int d1=0; d1<dim; ++d1)
+                                {
+                                  for (unsigned int d2=0; d2<dim; ++d2)
+                                    {
+                                      fe_data.shape_grads[dof_index1][q][d1][d2] = ( face_xi_grads[m][d1]*face_xi_grads[m][d2]*polyxi[i][2]*polyeta[j][0]
+                                                                                     + ( face_xi_grads[m][d1]*face_eta_grads[m][d2]
+                                                                                         + face_xi_grads[m][d2]*face_eta_grads[m][d1]
+                                                                                       )*polyxi[i][1]*polyeta[j][1]
+                                                                                     + face_eta_grads[m][d1]*face_eta_grads[m][d2]*polyxi[i][0]*polyeta[j][2]
+                                                                                   )*face_lambda_values[m][q]
+                                                                                   + ( face_xi_grads[m][d2]*polyxi[i][1]*polyeta[j][0]
+                                                                                       + face_eta_grads[m][d2]*polyxi[i][0]*polyeta[j][1]
+                                                                                     )*face_lambda_grads[m][d1]
+                                                                                   + ( face_xi_grads[m][d1]*polyxi[i][1]*polyeta[j][0]
+                                                                                       + face_eta_grads[m][d1]*polyxi[i][0]*polyeta[j][1]
+                                                                                     )*face_lambda_grads[m][d2];
+                                    }
+                                }
+                              // Type 2:
+                              const unsigned int dof_index2(face_type2_offset + shift_ij);
+                              for (unsigned int d1=0; d1<dim; ++d1)
+                                {
+                                  for (unsigned int d2=0; d2<dim; ++d2)
+                                    {
+                                      fe_data.shape_grads[dof_index2][q][d1][d2] = ( face_xi_grads[m][d1]*face_xi_grads[m][d2]*polyxi[i][2]*polyeta[j][0]
+                                                                                     + ( face_xi_grads[m][d1]*face_eta_grads[m][d2]
+                                                                                         - face_xi_grads[m][d2]*face_eta_grads[m][d1]
+                                                                                       )*polyxi[i][1]*polyeta[j][1]
+                                                                                     - face_eta_grads[m][d1]*face_eta_grads[m][d2]*polyxi[i][0]*polyeta[j][2]
+                                                                                   )*face_lambda_values[m][q]
+                                                                                   + ( face_xi_grads[m][d1]*polyxi[i][1]*polyeta[j][0]
+                                                                                       - face_eta_grads[m][d1]*polyxi[i][0]*polyeta[j][1]
+                                                                                     )*face_lambda_grads[m][d2];
+                                    }
+                                }
+                            }
+                          // Type 3:
+                          const unsigned int dof_index3_1(face_type3_offset1 + j);
+                          const unsigned int dof_index3_2(face_type3_offset2 + j);
+                          for (unsigned int d1=0; d1<dim; ++d1)
+                            {
+                              for (unsigned int d2=0; d2<dim; ++d2)
+                                {
+                                  fe_data.shape_grads[dof_index3_1][q][d1][d2] = face_xi_grads[m][d1]*( face_eta_grads[m][d2]*polyeta[j][1]*face_lambda_values[m][q]
+                                                                                 + face_lambda_grads[m][d2]*polyeta[j][0] );
+
+                                  fe_data.shape_grads[dof_index3_2][q][d1][d2] = face_eta_grads[m][d1]*( face_xi_grads[m][d2]*polyxi[j][1]*face_lambda_values[m][q]
+                                                                                 + face_lambda_grads[m][d2]*polyxi[j][0] );
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+      if (flags & update_hessians)
+        {
+          Assert(false, ExcNotImplemented());
+        }
+    }
+}
+
+template <int dim>
+void FE_NedelecSZ<dim>::fill_fe_values(
+  const typename Triangulation<dim,dim>::cell_iterator                         &cell,
+  const CellSimilarity::Similarity                                              /*cell_similarity*/,
+  const Quadrature<dim>                                                        &quadrature,
+  const Mapping<dim,dim>                                                       &mapping,
+  const typename Mapping<dim,dim>::InternalDataBase                            &mapping_internal,
+  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
+  const typename FiniteElement<dim,dim>::InternalDataBase                      &fe_internal,
+  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &data) const
+{
+
+  // Convert to the correct internal data class for this FE class.
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != 0,
+          ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+  // Now update the edge-based DoFs, which depend on the cell.
+  // This will fill in the missing items in the InternalData (fe_internal/fe_data)
+  // which was not filled in by get_data.
+  fill_edge_values(cell, quadrature, fe_data);
+  if (dim == 3 && this->degree>1)
+    {
+      fill_face_values(cell, quadrature, fe_data);
+    }
+
+  const UpdateFlags flags(fe_data.update_each);
+  const unsigned int n_q_points = quadrature.size();
+
+  Assert(!(flags & update_values) || fe_data.shape_values.size() == this->dofs_per_cell,
+         ExcDimensionMismatch(fe_data.shape_values.size(), this->dofs_per_cell));
+  Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_q_points,
+         ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
+
+  if (flags & update_values)
+    {
+      // Now have all shape_values stored on the reference cell.
+      // Must now transform to the physical cell.
+      std::vector<Tensor<1,dim>> transformed_shape_values(n_q_points);
+      for (unsigned int dof=0; dof<this->dofs_per_cell; ++dof)
+        {
+          const unsigned int first
+            = data.shape_function_to_row_table[dof * this->n_components() +
+                                               this->get_nonzero_components(dof).first_selected_component()];
+
+          mapping.transform(make_array_view(fe_data.shape_values[dof]),
+                            mapping_covariant,
+                            mapping_internal,
+                            make_array_view(transformed_shape_values));
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  data.shape_values(first+d,q) = transformed_shape_values[q][d];
+                }
+            }
+        }
+    }
+  if (flags & update_gradients)
+    {
+      // Now have all shape_grads stored on the reference cell.
+      // Must now transform to the physical cell.
+      std::vector<Tensor<2,dim>> input(n_q_points);
+      std::vector<Tensor<2,dim>> transformed_shape_grads(n_q_points);
+      for (unsigned int dof=0; dof<this->dofs_per_cell; ++dof)
+        {
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              input[q] = fe_data.shape_grads[dof][q];
+            }
+          mapping.transform(make_array_view(input),
+                            mapping_covariant_gradient,
+                            mapping_internal,
+                            make_array_view(transformed_shape_grads));
+
+          const unsigned int first
+            = data.shape_function_to_row_table[dof * this->n_components() +
+                                               this->get_nonzero_components(dof).first_selected_component()];
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d1=0; d1<dim; ++d1)
+                {
+                  for (unsigned int d2=0; d2<dim; ++d2)
+                    {
+                      transformed_shape_grads[q][d1] -= data.shape_values(first+d2, q)
+                                                        * mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
+                    }
+                }
+            }
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  data.shape_gradients[first+d][q] = transformed_shape_grads[q][d];
+                }
+            }
+        }
+    }
+}
+
+template <int dim>
+void FE_NedelecSZ<dim>::fill_fe_face_values (
+  const typename Triangulation<dim,dim>::cell_iterator                         &cell,
+  const unsigned int                                                            face_no,
+  const Quadrature<dim-1>                                                      &quadrature,
+  const Mapping<dim,dim>                                                       &mapping,
+  const typename Mapping<dim,dim>::InternalDataBase                            &mapping_internal,
+  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
+  const typename FiniteElement<dim,dim>::InternalDataBase                      &fe_internal,
+  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &data) const
+{
+  // Note for future improvement:
+  // We don't have the full quadrature - should use QProjector to create the 2D quadrature.
+  //
+  // For now I am effectively generating all of the shape function vals/grads, etc.
+  // On all quad points on all faces and then only using them for one face.
+  // This is obviously inefficient. I should cache the cell number and cache
+  // all of the shape_values/gradients etc and then reuse them for each face.
+
+  // convert data object to internal
+  // data for this class. fails with
+  // an exception if that is not
+  // possible
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != 0,
+          ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+  // Now update the edge-based DoFs, which depend on the cell.
+  // This will fill in the missing items in the InternalData (fe_internal/fe_data)
+  // which was not filled in by get_data.
+  fill_edge_values(cell,
+                   QProjector<dim>::project_to_all_faces(quadrature),
+                   fe_data);
+  if (dim == 3 && this->degree>1)
+    {
+      fill_face_values(cell,
+                       QProjector<dim>::project_to_all_faces(quadrature),
+                       fe_data);
+    }
+
+  const UpdateFlags flags(fe_data.update_each);
+  const unsigned int n_q_points = quadrature.size();
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = QProjector<dim>::DataSetDescriptor::face (face_no,
+                                                cell->face_orientation(face_no),
+                                                cell->face_flip(face_no),
+                                                cell->face_rotation(face_no),
+                                                n_q_points);
+
+  if (flags & update_values)
+    {
+      // Now have all shape_values stored on the reference cell.
+      // Must now transform to the physical cell.
+      std::vector<Tensor<1,dim>> transformed_shape_values(n_q_points);
+      for (unsigned int dof=0; dof<this->dofs_per_cell; ++dof)
+        {
+          mapping.transform (make_array_view(fe_data.shape_values[dof], offset, n_q_points),
+                             mapping_covariant,
+                             mapping_internal,
+                             make_array_view(transformed_shape_values));
+
+          const unsigned int first
+            = data.shape_function_to_row_table[dof * this->n_components() +
+                                               this->get_nonzero_components(dof).first_selected_component()];
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  data.shape_values(first+d,q) = transformed_shape_values[q][d];
+                }
+            }
+        }
+    }
+  if (flags & update_gradients)
+    {
+      // Now have all shape_grads stored on the reference cell.
+      // Must now transform to the physical cell.
+      std::vector<Tensor<2,dim>> transformed_shape_grads(n_q_points);
+      for (unsigned int dof=0; dof<this->dofs_per_cell; ++dof)
+        {
+          mapping.transform(make_array_view(fe_data.shape_grads[dof], offset, n_q_points),
+                            mapping_covariant_gradient,
+                            mapping_internal,
+                            make_array_view(transformed_shape_grads));
+
+          const unsigned int first
+            = data.shape_function_to_row_table[dof * this->n_components() +
+                                               this->get_nonzero_components(dof).first_selected_component()];
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d1=0; d1<dim; ++d1)
+                {
+                  for (unsigned int d2=0; d2<dim; ++d2)
+                    {
+                      transformed_shape_grads[q][d1] -= data.shape_values(first+d2, q)
+                                                        * mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
+                    }
+                }
+            }
+
+          for (unsigned int q=0; q<n_q_points; ++q)
+            {
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  data.shape_gradients[first+d][q] = transformed_shape_grads[q][d];
+                }
+            }
+        }
+    }
+}
+
+template <int dim>
+void FE_NedelecSZ<dim>::fill_fe_subface_values(
+  const typename Triangulation<dim,dim>::cell_iterator                         &/*cell*/,
+  const unsigned int                                                            /*face_no*/,
+  const unsigned int                                                            /*sub_no*/,
+  const Quadrature<dim-1>                                                      &/*quadrature*/,
+  const Mapping<dim,dim>                                                       &/*mapping*/,
+  const typename Mapping<dim,dim>::InternalDataBase                            &/*mapping_internal*/,
+  const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> &/*mapping_data*/,
+  const typename FiniteElement<dim,dim>::InternalDataBase                      &/*fe_internal*/,
+  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &/*data*/) const
+{
+  Assert (false, ExcNotImplemented ());
+}
+
+template <int dim>
+UpdateFlags
+FE_NedelecSZ<dim>::requires_update_flags(const UpdateFlags flags) const
+{
+  return update_once(flags) | update_each(flags);
+}
+
+template <int dim>
+UpdateFlags
+FE_NedelecSZ<dim>::update_once (const UpdateFlags flags) const
+{
+  const bool values_once = (mapping_type == mapping_none);
+
+  UpdateFlags out = update_default;
+  if (values_once && (flags & update_values))
+    out |= update_values;
+
+  return out;
+}
+
+template <int dim>
+UpdateFlags
+FE_NedelecSZ<dim>::update_each (const UpdateFlags flags) const
+{
+  UpdateFlags out = update_default;
+
+  if (flags & update_values)
+    out |= update_values | update_covariant_transformation;
+
+  if (flags & update_gradients)
+    out |= update_gradients | update_values |
+           update_jacobian_pushed_forward_grads |
+           update_covariant_transformation;
+
+  if (flags & update_hessians)
+    //     Assert (false, ExcNotImplemented());
+    out |= update_hessians | update_values | update_gradients |
+           update_jacobian_pushed_forward_grads |
+           update_jacobian_pushed_forward_2nd_derivatives |
+           update_covariant_transformation;
+
+  return out;
+}
+
+template <int dim>
+std::string
+FE_NedelecSZ<dim>::get_name () const
+{
+  // note that the
+  // FETools::get_fe_from_name
+  // function depends on the
+  // particular format of the string
+  // this function returns, so they
+  // have to be kept in synch
+
+  std::ostringstream namebuf;
+  namebuf << "FE_NedelecSZ<" << dim << ">(" << this->degree-1 << ")";
+
+  return namebuf.str();
+}
+
+template <int dim>
+std::unique_ptr<FiniteElement<dim,dim> >
+FE_NedelecSZ<dim>::clone () const
+{
+  return std_cxx14::make_unique<FE_NedelecSZ<dim> >(*this);
+}
+
+template <int dim>
+std::vector<unsigned int> FE_NedelecSZ<dim>::get_dpo_vector (const unsigned int degree)
+{
+  // internal function to return a vector of "dofs per object"
+  // where the objects inside the vector refer to:
+  // 0 = vertex
+  // 1 = edge
+  // 2 = face (which is a cell in 2D)
+  // 3 = cell
+  std::vector<unsigned int> dpo (dim+1);
+  dpo[0] = 0;
+  dpo[1] = degree + 1;
+  dpo[2] = 2 * degree * (degree + 1);
+  if (dim == 3)
+    {
+      dpo[3] = 3 * degree * degree * (degree + 1);
+    }
+  return dpo;
+}
+
+template <int dim>
+unsigned int FE_NedelecSZ<dim>::compute_num_dofs (const unsigned int degree) const
+{
+  // Internal function to compute the number of DoFs
+  // for a given dimension & polynomial order.
+  switch (dim)
+    {
+    case 2:
+      return 2 * (degree + 1) * (degree + 2);
+
+    case 3:
+      return 3 * (degree + 1) * (degree + 2) * (degree + 2);
+
+    default:
+    {
+      Assert (false, ExcNotImplemented ());
+      return 0;
+    }
+    }
+}
+
+template <int dim>
+void FE_NedelecSZ<dim>::create_polynomials (const unsigned int degree)
+{
+  // fill the 1d polynomials vector:
+  IntegratedLegendrePolynomials = IntegratedLegendreSZ::generate_complete_basis (degree + 1);
+}
+
+// explicit instantiations
+#include "fe_nedelec_sz.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/fe_nedelec_sz.inst.in b/source/fe/fe_nedelec_sz.inst.in
new file mode 100644 (file)
index 0000000..dd9bef6
--- /dev/null
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS)
+{
+#if deal_II_dimension != 1
+    template class FE_NedelecSZ<deal_II_dimension>;
+#endif
+}
+
diff --git a/tests/fe/fe_nedelec_sz_face_orientation.cc b/tests/fe/fe_nedelec_sz_face_orientation.cc
new file mode 100644 (file)
index 0000000..1577502
--- /dev/null
@@ -0,0 +1,223 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test FE_NedelecSZ<3> for meshes with faces with non-standard orientation.
+
+#include "../tests.h"
+#include<deal.II/base/quadrature_lib.h>
+#include<deal.II/dofs/dof_handler.h>
+#include<deal.II/fe/fe_nedelec_sz.h>
+#include<deal.II/fe/fe_nedelec.h>
+#include<deal.II/fe/fe_values.h>
+#include<deal.II/grid/grid_generator.h>
+#include<deal.II/grid/tria.h>
+#include<deal.II/lac/constraint_matrix.h>
+#include<deal.II/lac/vector.h>
+#include<deal.II/numerics/fe_field_function.h>
+#include<deal.II/numerics/vector_tools.h>
+
+void create_reference_triangulation (Triangulation<3> &tria)
+{
+  std::vector<unsigned int> repetitions (3, 1);
+
+  repetitions[0] = 2;
+  GridGenerator::subdivided_hyper_rectangle (tria, repetitions, Point<3> (-1.0, 0.0, 0.0), Point<3> (1.0, 1.0, 1.0));
+}
+
+void create_triangulation (Triangulation<3> &tria, const bool face_orientation, const bool face_flip, const bool face_rotation)
+{
+  std::vector<CellData<3> > cells (2);
+
+  cells[0].vertices[0] = 0;
+  cells[0].vertices[1] = 1;
+  cells[0].vertices[2] = 6;
+  cells[0].vertices[3] = 7;
+  cells[0].vertices[4] = 3;
+  cells[0].vertices[5] = 4;
+  cells[0].vertices[6] = 9;
+  cells[0].vertices[7] = 10;
+  cells[0].material_id = 0;
+
+  if (face_orientation)
+    {
+      if (face_flip)
+        {
+          if (face_rotation)
+            {
+              cells[1].vertices[0] = 7;
+              cells[1].vertices[1] = 8;
+              cells[1].vertices[2] = 10;
+              cells[1].vertices[3] = 11;
+              cells[1].vertices[4] = 1;
+              cells[1].vertices[5] = 2;
+              cells[1].vertices[6] = 4;
+              cells[1].vertices[7] = 5;
+            }
+
+          else
+            {
+              cells[1].vertices[0] = 10;
+              cells[1].vertices[1] = 11;
+              cells[1].vertices[2] = 4;
+              cells[1].vertices[3] = 5;
+              cells[1].vertices[4] = 7;
+              cells[1].vertices[5] = 8;
+              cells[1].vertices[6] = 1;
+              cells[1].vertices[7] = 2;
+            }
+        }
+
+      else if (face_rotation)
+        {
+          cells[1].vertices[0] = 4;
+          cells[1].vertices[1] = 5;
+          cells[1].vertices[2] = 1;
+          cells[1].vertices[3] = 2;
+          cells[1].vertices[4] = 10;
+          cells[1].vertices[5] = 11;
+          cells[1].vertices[6] = 7;
+          cells[1].vertices[7] = 8;
+        }
+
+      else
+        {
+          cells[1].vertices[0] = 2;
+          cells[1].vertices[1] = 1;
+          cells[1].vertices[2] = 5;
+          cells[1].vertices[3] = 4;
+          cells[1].vertices[4] = 8;
+          cells[1].vertices[5] = 7;
+          cells[1].vertices[6] = 11;
+          cells[1].vertices[7] = 10;
+        }
+    }
+
+  else if (face_flip)
+    {
+      if (face_rotation)
+        {
+          cells[1].vertices[0] = 8;
+          cells[1].vertices[1] = 7;
+          cells[1].vertices[2] = 2;
+          cells[1].vertices[3] = 1;
+          cells[1].vertices[4] = 11;
+          cells[1].vertices[5] = 10;
+          cells[1].vertices[6] = 5;
+          cells[1].vertices[7] = 4;
+        }
+
+      else
+        {
+          cells[1].vertices[0] = 11;
+          cells[1].vertices[1] = 10;
+          cells[1].vertices[2] = 8;
+          cells[1].vertices[3] = 7;
+          cells[1].vertices[4] = 5;
+          cells[1].vertices[5] = 4;
+          cells[1].vertices[6] = 2;
+          cells[1].vertices[7] = 1;
+        }
+    }
+
+  else if (face_rotation)
+    {
+      cells[1].vertices[0] = 5;
+      cells[1].vertices[1] = 4;
+      cells[1].vertices[2] = 11;
+      cells[1].vertices[3] = 10;
+      cells[1].vertices[4] = 2;
+      cells[1].vertices[5] = 1;
+      cells[1].vertices[6] = 8;
+      cells[1].vertices[7] = 7;
+    }
+
+  else
+    {
+      cells[1].vertices[0] = 1;
+      cells[1].vertices[1] = 2;
+      cells[1].vertices[2] = 7;
+      cells[1].vertices[3] = 8;
+      cells[1].vertices[4] = 4;
+      cells[1].vertices[5] = 5;
+      cells[1].vertices[6] = 10;
+      cells[1].vertices[7] = 11;
+    }
+
+  cells[1].material_id = 0;
+
+  std::vector<Point<3> > vertices (12);
+
+  vertices[0] = Point<3> (-1.0, 0.0, 0.0);
+  vertices[1] = Point<3> ();
+  vertices[2] = Point<3> (1.0, 0.0, 0.0);
+  vertices[3] = Point<3> (-1.0, 0.0, 1.0);
+  vertices[4] = Point<3> (0.0, 0.0, 1.0);
+  vertices[5] = Point<3> (1.0, 0.0, 1.0);
+  vertices[6] = Point<3> (-1.0, 1.0, 0.0);
+  vertices[7] = Point<3> (0.0, 1.0, 0.0);
+  vertices[8] = Point<3> (1.0, 1.0, 0.0);
+  vertices[9] = Point<3> (-1.0, 1.0, 1.0);
+  vertices[10] = Point<3> (0.0, 1.0, 1.0);
+  vertices[11] = Point<3> (1.0, 1.0, 1.0);
+  tria.create_triangulation (vertices, cells, SubCellData ());
+}
+
+void evaluate (const FiniteElement<3> &fe, const DoFHandler<3> &dof_handler)
+{
+  const FEValuesExtractors::Vector component (0);
+  const Quadrature<3> quadrature(Point<3>(0.5, 0.5, 0.5));
+  FEValues<3> fe_values (fe, quadrature, update_quadrature_points | update_values | update_gradients);
+
+  for (DoFHandler<3>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell)
+    {
+      fe_values.reinit (cell);
+
+      for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+        {
+          deallog << "DoF#" << i << ", value=[" << fe_values[component].value(i, 0) << "], curl=[" << fe_values[component].curl(i, 0) << "]" << std::endl;
+        }
+    }
+}
+
+void run (const bool face_orientation, const bool face_flip, const bool face_rotation)
+{
+//  Triangulation<3> tria_ref;
+//  create_reference_triangulation (tria_ref);
+
+  FE_NedelecSZ<3> fe (1);
+//  DoFHandler<3> dof_handler_ref (tria_ref);
+//  dof_handler_ref.distribute_dofs (fe);
+
+  Triangulation<3> tria;
+  create_triangulation (tria, face_orientation, face_flip, face_rotation);
+
+  DoFHandler<3> dof_handler (tria);
+  dof_handler.distribute_dofs (fe);
+
+  evaluate (fe, dof_handler);
+}
+
+int main()
+{
+  initlog();
+
+  for (int face_orientation = 0; face_orientation <= 1; ++face_orientation)
+    for (int face_flip = 0; face_flip <= 1; ++face_flip)
+      for (int face_rotation = 0; face_rotation <= 1; ++face_rotation)
+        {
+          deallog << face_orientation << face_flip << face_rotation << std::endl;
+          run (face_orientation, face_flip, face_rotation);
+        }
+}
diff --git a/tests/fe/fe_nedelec_sz_face_orientation.output b/tests/fe/fe_nedelec_sz_face_orientation.output
new file mode 100644 (file)
index 0000000..774a045
--- /dev/null
@@ -0,0 +1,873 @@
+
+DEAL::000
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::001
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#1, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#3, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#9, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#11, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#17, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#19, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#21, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#23, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::010
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#1, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#3, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#5, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#9, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#11, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#15, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#17, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#19, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#21, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#23, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#32, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::011
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#1, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#3, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#5, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#7, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#9, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#11, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#13, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#15, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#17, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#19, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#21, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#23, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#32, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::100
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#1, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#3, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#7, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#9, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#11, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#13, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#17, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#19, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#21, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#23, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#32, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::101
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#1, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#3, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#7, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#9, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#11, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#13, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#17, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#19, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#21, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#23, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::110
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#5, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#7, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#13, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#15, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::111
+DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000]
+DEAL::DoF#1, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000]
+DEAL::DoF#3, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000]
+DEAL::DoF#5, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000]
+DEAL::DoF#7, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000]
+DEAL::DoF#9, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000]
+DEAL::DoF#11, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000]
+DEAL::DoF#13, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000]
+DEAL::DoF#15, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000]
+DEAL::DoF#17, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000]
+DEAL::DoF#19, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000]
+DEAL::DoF#21, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000]
+DEAL::DoF#23, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#32, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000]
+DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#36, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000]
+DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#40, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000]
+DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000]
+DEAL::DoF#44, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000]
+DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000]
+DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#52, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000]
+DEAL::DoF#53, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000]
diff --git a/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc
new file mode 100644 (file)
index 0000000..84ce85c
--- /dev/null
@@ -0,0 +1,165 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// check the correctness of fe_values.shape_gradient for FE_NedelecSZ by comparing
+// the integral of all shape gradients with the flux over the boundary by the
+// divergence theorem
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_nedelec_sz.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <sstream>
+
+template <int dim>
+Tensor<1,dim> ones ()
+{
+  Tensor<1,dim> result;
+  for (unsigned int i=0; i<dim; ++i)
+    result[i] = 1.0;
+  return result;
+}
+
+template <int dim>
+void test (const Triangulation<dim> &tr,
+           const FiniteElement<dim> &fe,
+           const double tolerance)
+{
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  std::stringstream ss;
+
+  deallog << "FE=" << fe.get_name() << std::endl;
+
+  const QGauss<dim> quadrature(4);
+  FEValues<dim> fe_values (fe, quadrature, update_gradients
+                           | update_quadrature_points
+                           | update_JxW_values);
+
+  const QGauss<dim-1> face_quadrature(4);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature,
+                                    update_values
+                                    | update_quadrature_points
+                                    | update_normal_vectors
+                                    | update_JxW_values);
+
+  for (typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+       cell != dof.end();
+       ++cell)
+    {
+      fe_values.reinit (cell);
+
+      deallog << "Cell nodes:" << std::endl;
+      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+        {
+          deallog << i << ": ( ";
+          for (unsigned int d=0; d<dim; ++d)
+            deallog << cell->vertex(i)[d] << " ";
+          deallog << ")" << std::endl;
+        }
+
+      bool cell_ok = true;
+
+      for (unsigned int c=0; c<fe.n_components(); ++c)
+        {
+          FEValuesExtractors::Scalar single_component (c);
+
+          for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+            {
+              ss << "component=" << c
+                 << ", dof=" << i
+                 << std::endl;
+
+              Tensor<1,dim> bulk_integral;
+              for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+                {
+                  bulk_integral += fe_values[single_component].gradient(i,q) * fe_values.JxW(q);
+                }
+
+              Tensor<1,dim> boundary_integral;
+              for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+                {
+                  fe_face_values.reinit(cell,face);
+                  for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
+                    {
+                      boundary_integral += fe_face_values[single_component].value (i,q)
+                                           * fe_face_values.normal_vector(q)
+                                           * fe_face_values.JxW(q);
+                    }
+                }
+
+              if ((bulk_integral-boundary_integral).norm_square() > tolerance * (bulk_integral.norm() + boundary_integral.norm()))
+                {
+                  deallog << "Failed:" << std::endl;
+                  deallog << ss.str() << std::endl;
+                  deallog << "    bulk integral=" << bulk_integral << std::endl;
+                  deallog << "boundary integral=" << boundary_integral << std::endl;
+                  deallog << "Error! difference between bulk and surface integrals is greater than " << tolerance << "!\n\n" << std::endl;
+                  cell_ok = false;
+                }
+
+              ss.str("");
+            }
+        }
+
+      deallog << (cell_ok? "OK: cell bulk and boundary integrals match...\n" : "Failed divergence test...\n") << std::endl;
+    }
+}
+
+
+
+template <int dim>
+void test_hyper_ball(const double tolerance)
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_ball(tr);
+
+  static const SphericalManifold<dim> boundary;
+  tr.set_manifold (0, boundary);
+
+//  tr.refine_global(1); // taking too long,
+
+  FE_NedelecSZ<dim> fe(1);
+  test(tr, fe, tolerance);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (8);
+
+  deallog.attach(logfile);
+
+  test_hyper_ball<2>(1e-6);
+  test_hyper_ball<3>(1e-6);
+
+  deallog << "done..." << std::endl;
+}
+
diff --git a/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.output b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.output
new file mode 100644 (file)
index 0000000..6c5fcb6
--- /dev/null
@@ -0,0 +1,116 @@
+
+DEAL::FE=FE_NedelecSZ<2>(1)
+DEAL::Cell nodes:
+DEAL::0: ( -0.70710678 -0.70710678 )
+DEAL::1: ( 0.70710678 -0.70710678 )
+DEAL::2: ( -0.29289322 -0.29289322 )
+DEAL::3: ( 0.29289322 -0.29289322 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.70710678 -0.70710678 )
+DEAL::1: ( -0.29289322 -0.29289322 )
+DEAL::2: ( -0.70710678 0.70710678 )
+DEAL::3: ( -0.29289322 0.29289322 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.29289322 -0.29289322 )
+DEAL::1: ( 0.29289322 -0.29289322 )
+DEAL::2: ( -0.29289322 0.29289322 )
+DEAL::3: ( 0.29289322 0.29289322 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.70710678 -0.70710678 )
+DEAL::1: ( 0.70710678 0.70710678 )
+DEAL::2: ( 0.29289322 -0.29289322 )
+DEAL::3: ( 0.29289322 0.29289322 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.70710678 0.70710678 )
+DEAL::1: ( -0.29289322 0.29289322 )
+DEAL::2: ( 0.70710678 0.70710678 )
+DEAL::3: ( 0.29289322 0.29289322 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::FE=FE_NedelecSZ<3>(1)
+DEAL::Cell nodes:
+DEAL::0: ( -0.21132487 -0.21132487 -0.21132487 )
+DEAL::1: ( 0.21132487 -0.21132487 -0.21132487 )
+DEAL::2: ( -0.21132487 0.21132487 -0.21132487 )
+DEAL::3: ( 0.21132487 0.21132487 -0.21132487 )
+DEAL::4: ( -0.21132487 -0.21132487 0.21132487 )
+DEAL::5: ( 0.21132487 -0.21132487 0.21132487 )
+DEAL::6: ( -0.21132487 0.21132487 0.21132487 )
+DEAL::7: ( 0.21132487 0.21132487 0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.57735027 -0.57735027 -0.57735027 )
+DEAL::1: ( 0.57735027 -0.57735027 -0.57735027 )
+DEAL::2: ( -0.57735027 0.57735027 -0.57735027 )
+DEAL::3: ( 0.57735027 0.57735027 -0.57735027 )
+DEAL::4: ( -0.21132487 -0.21132487 -0.21132487 )
+DEAL::5: ( 0.21132487 -0.21132487 -0.21132487 )
+DEAL::6: ( -0.21132487 0.21132487 -0.21132487 )
+DEAL::7: ( 0.21132487 0.21132487 -0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( 0.57735027 -0.57735027 -0.57735027 )
+DEAL::1: ( 0.57735027 0.57735027 -0.57735027 )
+DEAL::2: ( 0.21132487 -0.21132487 -0.21132487 )
+DEAL::3: ( 0.21132487 0.21132487 -0.21132487 )
+DEAL::4: ( 0.57735027 -0.57735027 0.57735027 )
+DEAL::5: ( 0.57735027 0.57735027 0.57735027 )
+DEAL::6: ( 0.21132487 -0.21132487 0.21132487 )
+DEAL::7: ( 0.21132487 0.21132487 0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.57735027 -0.57735027 0.57735027 )
+DEAL::1: ( 0.57735027 -0.57735027 0.57735027 )
+DEAL::2: ( -0.21132487 -0.21132487 0.21132487 )
+DEAL::3: ( 0.21132487 -0.21132487 0.21132487 )
+DEAL::4: ( -0.57735027 0.57735027 0.57735027 )
+DEAL::5: ( 0.57735027 0.57735027 0.57735027 )
+DEAL::6: ( -0.21132487 0.21132487 0.21132487 )
+DEAL::7: ( 0.21132487 0.21132487 0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.57735027 -0.57735027 -0.57735027 )
+DEAL::1: ( -0.21132487 -0.21132487 -0.21132487 )
+DEAL::2: ( -0.57735027 0.57735027 -0.57735027 )
+DEAL::3: ( -0.21132487 0.21132487 -0.21132487 )
+DEAL::4: ( -0.57735027 -0.57735027 0.57735027 )
+DEAL::5: ( -0.21132487 -0.21132487 0.21132487 )
+DEAL::6: ( -0.57735027 0.57735027 0.57735027 )
+DEAL::7: ( -0.21132487 0.21132487 0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.57735027 -0.57735027 -0.57735027 )
+DEAL::1: ( 0.57735027 -0.57735027 -0.57735027 )
+DEAL::2: ( -0.21132487 -0.21132487 -0.21132487 )
+DEAL::3: ( 0.21132487 -0.21132487 -0.21132487 )
+DEAL::4: ( -0.57735027 -0.57735027 0.57735027 )
+DEAL::5: ( 0.57735027 -0.57735027 0.57735027 )
+DEAL::6: ( -0.21132487 -0.21132487 0.21132487 )
+DEAL::7: ( 0.21132487 -0.21132487 0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::Cell nodes:
+DEAL::0: ( -0.57735027 0.57735027 -0.57735027 )
+DEAL::1: ( -0.21132487 0.21132487 -0.21132487 )
+DEAL::2: ( 0.57735027 0.57735027 -0.57735027 )
+DEAL::3: ( 0.21132487 0.21132487 -0.21132487 )
+DEAL::4: ( -0.57735027 0.57735027 0.57735027 )
+DEAL::5: ( -0.21132487 0.21132487 0.21132487 )
+DEAL::6: ( 0.57735027 0.57735027 0.57735027 )
+DEAL::7: ( 0.21132487 0.21132487 0.21132487 )
+DEAL::OK: cell bulk and boundary integrals match...
+
+DEAL::done...
diff --git a/tests/fe/fe_nedelec_sz_non_rect_2d.cc b/tests/fe/fe_nedelec_sz_non_rect_2d.cc
new file mode 100644 (file)
index 0000000..028b4c9
--- /dev/null
@@ -0,0 +1,322 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+//
+// By Ross Kynch, Alexander Grayver
+//
+// Test to confirm that FE_NedelecSZ works on deformed elements in 2D when using
+// Dirichlet boundary conditions (n x E = f, type). This is handled by the
+// function project_boundary_values_curl_conforming_l2().
+//
+// This test solves the real valued curl-curl equation in 2D:
+//
+// curl(curl(E)) + E = Js
+//
+// where the solution is:
+// E = (0, x^2)
+//
+// so, Js = (0,-2) + E.
+//
+// i.e the solution should be "exact" at p=2.
+//
+// The domain is a distorted quad after 1 refinement. The distortion is
+// performed by GridTools::distort_random.
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/fe/fe_nedelec_sz.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <iostream>
+#include <sstream>
+
+
+using namespace dealii;
+
+namespace polytest
+{
+  template <int dim>
+  class SimplePolynomial : public Function<dim>
+  {
+  public:
+    SimplePolynomial() : Function<dim>() {}
+    void vector_value_list (const std::vector<Point<dim> > &points,
+                            std::vector<Vector<double> >   &values) const;
+
+    void rhs_value_list (const std::vector<Point<dim> > &points,
+                         std::vector<Vector<double> >   &values) const;
+
+  };
+  template <int dim>
+  void SimplePolynomial<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                                 std::vector<Vector<double> >   &values) const
+  {
+    Assert (dim == 2, ExcNotImplemented());
+    Assert (values.size() == points.size(),
+            ExcDimensionMismatch(values.size(), points.size()));
+    for (unsigned int i=0; i<points.size(); ++i)
+      {
+        const Point<dim> &p = points[i];
+        // non-zero curl-curl:
+        values[i][0] = 0.0;
+        values[i][1] = p[0]*p[0];
+
+      }
+
+  }
+
+  template <int dim>
+  void SimplePolynomial<dim>::rhs_value_list (const std::vector<Point<dim> > &points,
+                                              std::vector<Vector<double> >   &values) const
+  {
+    Assert (dim == 2, ExcNotImplemented());
+    Assert (values.size() == points.size(),
+            ExcDimensionMismatch(values.size(), points.size()));
+
+    for (unsigned int i=0; i<points.size(); ++i)
+      {
+        const Point<dim> &p = points[i];
+        // non-zero curl-curl:
+        values[i][0] = 0.0;
+        values[i][1] = -2.0 + p[0]*p[0];
+      }
+
+  }
+
+  template <int dim>
+  class polytest
+  {
+  public:
+    polytest(unsigned int degree);
+    ~polytest ();
+
+    void run();
+
+  private:
+    void setup_system();
+    void assemble_system();
+    void solve();
+
+    void output_error();
+
+    unsigned int p_order;
+    unsigned int quad_order;
+
+    Triangulation<dim> tria;
+    DoFHandler<dim> dof_handler;
+    FE_NedelecSZ<dim> fe;
+    ConstraintMatrix constraints;
+    SparsityPattern sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double> solution;
+    Vector<double> system_rhs;
+
+  };
+
+  template <int dim>
+  polytest<dim>::polytest (unsigned int degree)
+    :
+    p_order(degree),
+    quad_order(2*degree + 3),
+    dof_handler(tria),
+    fe (degree)
+  {
+  }
+  template <int dim>
+  polytest<dim>::~polytest ()
+  {
+    dof_handler.clear ();
+  }
+
+  template <int dim>
+  void polytest<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs (fe);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+
+
+    constraints.clear ();
+    SimplePolynomial<dim> boundary_function;
+
+    VectorTools::project_boundary_values_curl_conforming_l2 (dof_handler,
+                                                             0,
+                                                             boundary_function,
+                                                             0,
+                                                             constraints);
+    constraints.close ();
+    DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    c_sparsity,
+                                    constraints,false);
+
+    sparsity_pattern.copy_from(c_sparsity);
+    system_matrix.reinit (sparsity_pattern);
+  }
+
+  template <int dim>
+  void polytest<dim>::assemble_system()
+  {
+    const QGauss<dim> test_quad(quad_order);
+    FEValues<dim> fe_values_test (fe, test_quad,
+                                  update_values    |  update_gradients |
+                                  update_quadrature_points  |  update_JxW_values);
+
+    const QGauss<dim> quadrature_formula(quad_order);
+    const unsigned int n_q_points = quadrature_formula.size();
+
+    const QGauss<dim-1> face_quadrature_formula(quad_order);
+    const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+
+    FEFaceValues<dim> fe_face_values(fe, face_quadrature_formula,
+                                     update_values | update_quadrature_points |
+                                     update_normal_vectors | update_JxW_values);
+
+    const FEValuesExtractors::Vector vec(0);
+
+    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+
+
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    SimplePolynomial<dim>      right_hand_side;
+    std::vector<Vector<double> > rhs_value_list(n_q_points,
+                                                Vector<double>(fe.n_components()));
+
+    typename DoFHandler<dim>::active_cell_iterator cell, endc;
+    endc = dof_handler.end();
+    cell = dof_handler.begin_active();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values_test.reinit(cell);
+        fe_values.reinit (cell);
+        cell_matrix=0;
+        cell_rhs=0;
+
+        right_hand_side.rhs_value_list(fe_values.get_quadrature_points(), rhs_value_list);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            Tensor<1,dim> rhs_value;
+            for (unsigned int d=0; d<dim; ++d)
+              {
+                rhs_value[d] = rhs_value_list[q](d);
+              }
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              {
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  {
+                    cell_matrix(i,j)+= (
+                                         fe_values[vec].curl(i,q)*fe_values[vec].curl(j,q)
+                                         + fe_values[vec].value(i,q)*fe_values[vec].value(j,q)
+                                       )*fe_values.JxW(q);
+                  }
+                cell_rhs(j) += rhs_value*fe_values[vec].value(j,q)*fe_values.JxW(q);
+              }
+          }
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global(cell_matrix,
+                                               cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix, system_rhs);
+      }
+  }
+
+  template <int dim>
+  void polytest<dim>::solve()
+  {
+    SolverControl      solver_control (1000, 1e-12);
+    SolverCG<>         solver (solver_control);
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+    solver.solve (system_matrix, solution, system_rhs,
+                  preconditioner);
+    constraints.distribute (solution);
+  }
+
+  template <int dim>
+  void polytest<dim>::output_error()
+  {
+    SimplePolynomial<dim> exact_solution;
+    Vector<double> diff_per_cell(tria.n_active_cells());
+    VectorTools::integrate_difference(dof_handler, solution, exact_solution,
+                                      diff_per_cell, QGauss<dim>(quad_order), VectorTools::L2_norm);
+    const double L2_error = diff_per_cell.l2_norm();
+
+    deallog << "p=" << p_order << " L2_error: " << L2_error << std::endl;
+  }
+
+  template <int dim>
+  void polytest<dim>::run()
+  {
+    GridGenerator::hyper_cube (tria, -1.2, 1.3);
+    // REFINE
+    tria.refine_global(1);
+
+    // DISTORT ALL:
+    GridTools::distort_random (0.2, tria, false);
+
+    setup_system();
+    assemble_system();
+    solve();
+    output_error();
+  }
+}
+int main ()
+{
+  const unsigned int dim(2);
+
+  initlog();
+
+  for (unsigned int p=0; p<3; ++p)
+    {
+      polytest::polytest<dim> poly(p);
+      poly.run();
+    }
+}
+
diff --git a/tests/fe/fe_nedelec_sz_non_rect_2d.output b/tests/fe/fe_nedelec_sz_non_rect_2d.output
new file mode 100644 (file)
index 0000000..80eb253
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL:cg::Starting value 0.747607
+DEAL:cg::Convergence step 4 value 6.71226e-19
+DEAL::p=0 L2_error: 0.864275
+DEAL:cg::Starting value 2.34680
+DEAL:cg::Convergence step 13 value 3.46263e-13
+DEAL::p=1 L2_error: 0.0650158
+DEAL:cg::Starting value 2.34814
+DEAL:cg::Convergence step 18 value 3.36532e-13
+DEAL::p=2 L2_error: 1.40992e-13
diff --git a/tests/fe/fe_nedelec_sz_non_rect_face.cc b/tests/fe/fe_nedelec_sz_non_rect_face.cc
new file mode 100644 (file)
index 0000000..d7317b0
--- /dev/null
@@ -0,0 +1,468 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+//
+// By Ross Kynch
+//
+// Test to confirm that FE_NedelecSZ works on non-rectangular faces in 3D.
+// The previous project_boundary_values_curl_conforming function only
+// produced the correct constraints on rectangular faces.
+//
+// The new function project_boundary_values_curl_conforming_l2 fixes this
+// and should work on any quad face, assuming the underlying mesh conforms to the
+// standard orientation (i.e. spheres may still cause problems).
+//
+// This test solves the real valued curl-curl equation in 3D:
+//
+// curl(curl(E)) + E = Js
+//
+// where the solution is:
+// E = (x^2 , y^2 + , z^2).
+//
+// so:
+// Js = E.
+//
+// The domain is a distorted cube which will contain non-rectangular faces.
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/convergence_table.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/fe/fe_nedelec_sz.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <iostream>
+#include <sstream>
+
+namespace Maxwell
+{
+  using namespace dealii;
+
+  // Dirichlet BCs / exact solution:.
+  template <int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    ExactSolution();
+
+    virtual void vector_value_list (const std::vector<Point<dim> > &points,
+                                    std::vector<Vector<double> >   &values) const;
+
+
+    void curl_value_list(const std::vector<Point<dim> > &points,
+                         std::vector<Vector<double> > &value_list);
+  };
+
+  template <int dim>
+  ExactSolution<dim>::ExactSolution()
+    :
+    Function<dim> (dim)
+  {}
+  template <int dim>
+  void ExactSolution<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                              std::vector<Vector<double> > &value_list) const
+  {
+    Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size()));
+    const unsigned int n_points = points.size();
+
+    for (unsigned int i=0; i<n_points; ++i)
+      {
+        const Point<dim> &p = points[i];
+
+        /* quadratic: */
+        value_list[i](0) = p(0)*p(0);
+        value_list[i](1) = p(1)*p(1);
+        value_list[i](2) = p(2)*p(2);
+
+      }
+  }
+  // Additional functions to create Neumann conditions, zero in this case.
+  template <int dim>
+  void ExactSolution<dim>::curl_value_list(const std::vector<Point<dim> > &points,
+                                           std::vector<Vector<double> > &value_list)
+  {
+    Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size()));
+    const unsigned int n_points = points.size();
+
+
+    double exponent;
+    for (unsigned int i=0; i<n_points; ++i)
+      {
+        const Point<dim> &p = points[i];
+        // Real:
+        value_list[i](0) = 0.0;
+        value_list[i](1) = 0.0;
+        value_list[i](2) = 0.0;
+      }
+  }
+  // END Dirichlet BC
+
+  // MAIN MAXWELL CLASS
+  template <int dim>
+  class MaxwellProblem
+  {
+  public:
+    MaxwellProblem (const unsigned int order);
+    ~MaxwellProblem ();
+    void run ();
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void process_solution(const unsigned int cycle);
+
+    double calcErrorHcurlNorm();
+
+    Triangulation<dim>   triangulation;
+    MappingQ<dim>        mapping;
+    DoFHandler<dim>      dof_handler;
+    FE_NedelecSZ<dim>            fe;
+    ConstraintMatrix     constraints;
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    ConvergenceTable       convergence_table;
+
+    ExactSolution<dim> exact_solution;
+
+    unsigned int p_order;
+    unsigned int quad_order;
+  };
+
+  template <int dim>
+  MaxwellProblem<dim>::MaxwellProblem (const unsigned int order)
+    :
+    mapping (1),
+    dof_handler (triangulation),
+    fe(order),
+    exact_solution()
+  {
+    p_order = order;
+    quad_order = 2*(p_order+1);
+  }
+
+  template <int dim>
+  MaxwellProblem<dim>::~MaxwellProblem ()
+  {
+    dof_handler.clear ();
+  }
+  template <int dim>
+  double MaxwellProblem<dim>::calcErrorHcurlNorm()
+  {
+    QGauss<dim>  quadrature_formula(quad_order);
+    const unsigned int n_q_points = quadrature_formula.size();
+
+    FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+
+    // Extractor
+    const FEValuesExtractors::Vector E_re(0);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    // storage for exact sol:
+    std::vector<Vector<double> > exactsol_list(n_q_points, Vector<double>(fe.n_components()));
+    std::vector<Vector<double> > exactcurlsol_list(n_q_points, Vector<double>(fe.n_components()));
+    Tensor<1,dim>  exactsol;
+    Tensor<1,dim>  exactcurlsol;
+
+    // storage for computed sol:
+    std::vector<Tensor<1,dim> > sol(n_q_points);
+    Tensor<1,dim> curlsol;
+
+    double h_curl_norm=0.0;
+
+    unsigned int block_index_i;
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+
+        // Store exact values of E and curlE:
+        exact_solution.vector_value_list(fe_values.get_quadrature_points(), exactsol_list);
+        exact_solution.curl_value_list(fe_values.get_quadrature_points(), exactcurlsol_list);
+
+
+        // Store computed values at quad points:
+        fe_values[E_re].get_function_values(solution, sol);
+
+        // Calc values of curlE from fe solution:
+        cell->get_dof_indices (local_dof_indices);
+        // Loop over quad points to calculate solution:
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            // Split exact solution into real/imaginary parts:
+            for (unsigned int component=0; component<dim; component++)
+              {
+                exactsol[component] = exactsol_list[q_point][component];
+                exactcurlsol[component] = exactcurlsol_list[q_point][component];
+              }
+            // Loop over DoFs to calculate curl of solution @ quad point
+            curlsol=0.0;
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                // Construct local curl value @ quad point
+                curlsol += solution(local_dof_indices[i])*fe_values[E_re].curl(i,q_point);
+              }
+            // Integrate difference at each point:
+            h_curl_norm += ( (exactsol-sol[q_point])*(exactsol-sol[q_point])
+                             + (exactcurlsol-curlsol)*(exactcurlsol-curlsol)
+                           )*fe_values.JxW(q_point);
+          }
+      }
+    return sqrt(h_curl_norm);
+  }
+
+  template <int dim>
+  void MaxwellProblem<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+    DoFRenumbering::block_wise (dof_handler);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+    constraints.clear ();
+
+
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+    // FE_Nedelec boundary condition.
+    VectorTools::project_boundary_values_curl_conforming_l2 (dof_handler, 0, exact_solution, 0, constraints);
+
+
+    constraints.close ();
+
+    DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    c_sparsity,
+                                    constraints,false);
+
+    sparsity_pattern.copy_from(c_sparsity);
+    system_matrix.reinit (sparsity_pattern);
+
+  }
+  template <int dim>
+  void MaxwellProblem<dim>::assemble_system ()
+  {
+    QGauss<dim>  quadrature_formula(quad_order);
+    QGauss<dim-1> face_quadrature_formula(quad_order);
+
+    const unsigned int n_q_points    = quadrature_formula.size();
+    const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+    FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+
+    FEFaceValues<dim> fe_face_values(mapping, fe, face_quadrature_formula,
+                                     update_values | update_quadrature_points |
+                                     update_normal_vectors | update_JxW_values);
+
+    // Extractor
+    const FEValuesExtractors::Vector E_re(0);
+
+    // Local cell storage:
+    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double> cell_rhs (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    //RHS storage:
+    std::vector<Vector<double> > rhs_value_list(n_q_points, Vector<double>(fe.n_components()));
+    Tensor<1,dim> rhs_value_vector;
+
+    // Neumann storage
+    std::vector<Vector<double> > neumann_value_list(n_face_q_points, Vector<double>(fe.n_components()));
+    std::vector<Tensor<1,dim> > normal_vector_list(fe_face_values.get_all_normal_vectors());
+    Tensor<1,dim> neumann_value_vector;
+    Tensor<1,dim> neumann_value;
+    Tensor<1,dim> normal_vector;
+
+    // loop over all cells:
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit(cell);
+        cell_matrix = 0;
+        cell_rhs = 0;
+
+        // Calc RHS values:
+        exact_solution.vector_value_list(fe_values.get_quadrature_points(), rhs_value_list);
+
+        // Loop over all element quad points:
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            // store rhs value at this q point & turn into tensor
+            for (unsigned int component=0; component<dim; component++)
+              {
+                rhs_value_vector[component] = rhs_value_list[q_point](component);
+              }
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                // Construct local matrix:
+                for (unsigned int j=0; j<dofs_per_cell; ++j)
+                  {
+                    cell_matrix(i,j) += ( (fe_values[E_re].curl(i,q_point)*fe_values[E_re].curl(j,q_point))
+                                          + fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)
+                                        )*fe_values.JxW(q_point);
+
+                  }
+                // construct local RHS:
+                cell_rhs(i) += rhs_value_vector*fe_values[E_re].value(i,q_point)*fe_values.JxW(q_point);
+
+              }
+          }
+
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global(cell_matrix,
+                                               cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix, system_rhs);
+      }
+  }
+  template <int dim>
+  void MaxwellProblem<dim>::solve ()
+  {
+    /* Direct */
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+
+    A_direct.vmult (solution, system_rhs);
+    constraints.distribute (solution);
+
+  }
+  template <int dim>
+  void MaxwellProblem<dim>::process_solution(const unsigned int cycle)
+  {
+
+    Vector<double> diff_per_cell(triangulation.n_active_cells());
+
+
+    VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
+                                      diff_per_cell, QGauss<dim>(quad_order+2),
+                                      VectorTools::L2_norm);
+
+
+    const double L2_error = diff_per_cell.l2_norm();
+
+    double Hcurl_error = calcErrorHcurlNorm();
+
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", triangulation.n_active_cells());
+    convergence_table.add_value("dofs", dof_handler.n_dofs());
+    convergence_table.add_value("L2 Error", L2_error);
+    convergence_table.add_value("H(curl) Error", Hcurl_error);
+
+  }
+
+  template <int dim>
+  void MaxwellProblem<dim>::run ()
+  {
+    unsigned int numcycles=(p_order>1)?2:3;
+
+    for (unsigned int cycle=0; cycle<numcycles; ++cycle)
+      {
+        if (cycle == 0)
+          {
+            /* Cube mesh */
+            GridGenerator::hyper_cube (triangulation, 0, 1);
+            GridTools::distort_random (0.2, triangulation, false);
+
+
+            /* Testing hanging nodes
+             *        triangulation.refine_global(1);
+             *        cell = triangulation.begin_active();
+             *        cell->set_refine_flag(); // create single refined cell
+             *        triangulation.execute_coarsening_and_refinement();
+             */
+
+          }
+        else
+          {
+            triangulation.refine_global (1);
+          }
+
+        setup_system ();
+        assemble_system ();
+        solve ();
+        process_solution (cycle);
+      }
+
+    deallog << "p = " << p_order << std::endl;
+
+    //convergence_table.set_precision("L2 Error",4);
+    convergence_table.set_scientific("L2 Error",true);
+
+    //convergence_table.set_precision("H(curl) Error",4);
+    convergence_table.set_scientific("H(curl) Error",true);
+
+    convergence_table.write_text(deallog.get_file_stream());
+
+  }
+}
+// END MAXWELL CLASS
+int main ()
+{
+  using namespace Maxwell;
+
+  initlog();
+
+  for (unsigned int p=0; p<3; ++p)
+    {
+      MaxwellProblem<3> (p).run();
+    }
+
+  return 0;
+}
diff --git a/tests/fe/fe_nedelec_sz_non_rect_face.with_umfpack=on.output b/tests/fe/fe_nedelec_sz_non_rect_face.with_umfpack=on.output
new file mode 100644 (file)
index 0000000..0a838e3
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL::p = 0
+cycle cells dofs  L2 Error  H(curl) Error 
+0     1     12   5.8761e-01 5.7154e-01    
+1     8     54   3.1310e-01 3.1132e-01    
+2     64    300  1.5881e-01 1.5860e-01    
+DEAL::p = 1
+cycle cells dofs  L2 Error  H(curl) Error 
+0     1     54   1.3757e-01 1.4097e-01    
+1     8     300  3.4151e-02 3.5926e-02    
+2     64    1944 8.5308e-03 9.4069e-03    
+DEAL::p = 2
+cycle cells dofs  L2 Error  H(curl) Error 
+0     1     144  5.3517e-16 7.5744e-16    
+1     8     882  6.0713e-15 6.5454e-15    

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.