--- /dev/null
+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)
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
#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>
&& (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)))
{
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;
&& (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)))
{
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;
// 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 ());
}
// 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 ());
}
fe_dg_vector.cc
fe_face.cc
fe_nedelec.cc
+ fe_nedelec_sz.cc
fe_nothing.cc
fe_poly.cc
fe_poly_tensor.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+ }
+}
--- /dev/null
+
+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]
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
+
--- /dev/null
+
+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...
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+ }
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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