From: Ross Kynch Date: Sat, 28 Apr 2018 20:48:37 +0000 (+0100) Subject: Added FE_NedelecSZ, a H(curl)-conforming element which overcomes the sign-conflict... X-Git-Tag: v9.1.0-rc1~1168^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=52200077938a86328be3abe82bdd198f66a46112;p=dealii.git Added FE_NedelecSZ, a H(curl)-conforming element which overcomes the sign-conflict problem for non-standard edge/face orientations. This work includes contributions (tests and bug fixes) from Alexander Gravyer. --- diff --git a/doc/news/changes/major/20180428RossKynch b/doc/news/changes/major/20180428RossKynch new file mode 100644 index 0000000000..0d86c7d762 --- /dev/null +++ b/doc/news/changes/major/20180428RossKynch @@ -0,0 +1,8 @@ +New: FE_NedelecSZ is a new H(curl)-conforming element which overcomes the sign conflict problem +which can be encountered for Nedelec elements on hexahedral meshes. This is an implementation of +the method described in the PhD thesis of Sabine Zaglmayr and should be considered over the +FE_Nedelec element where meshes with non-standard orientation are used. Note that not all +functionality (e.g. Hessians), has been implemented for this element. +
+(Ross Kynch, 2018/04/28) + diff --git a/include/deal.II/fe/fe_nedelec_sz.h b/include/deal.II/fe/fe_nedelec_sz.h new file mode 100644 index 0000000000..5575de3395 --- /dev/null +++ b/include/deal.II/fe/fe_nedelec_sz.h @@ -0,0 +1,427 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__fe_nedelec_sz_h +#define dealii__fe_nedelec_sz_h + +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/*!@addtogroup fe */ +/*@{*/ + +/** +* This class represents an implementation of the +* Hcurl-conforming Nédélec element described in the PhD +* thesis of S. Zaglmayr, High Order Finite Element Methods +* for Electromagnetic Field Computation, 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: Resolving the sign conflict +* problem for hp–hexahedral Nédélec elements with application to eddy +* current problems, 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 +class FE_NedelecSZ : public FiniteElement +{ +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 > clone() const; + + /** + * This element is vector-valued so this function will + * throw an exception. + */ + virtual double shape_value (const unsigned int i, + const Point &p) const; + + /** + * Not implemented. + */ + virtual double shape_value_component (const unsigned int i, + const Point &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 &p) const; + + /** + * Not implemented. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &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 &p) const; + + /** + * Not implemented. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Given flags, 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 flags, 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::InternalDataBase> + get_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &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::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fedata, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &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::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fedata, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &data) const; + + /** + * Not implemented. + */ + virtual void + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fedata, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &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: + *
    + *
  • \lambda_{i} - trilinear function, equal to one at the i-th vertex + * and zero at all other vertices.
  • + *
  • \sigma_{i} - linear functional associated with the i-th vertex.
  • + *
+ * + * 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::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 > > 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 > > 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 > > 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 > > 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 > 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 > 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 > 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 > 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 > > 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 > > 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 > 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 > 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 get_dpo_vector (const unsigned int degree); + + /** + * Internal storage for all required integrated Legendre polynomials. + */ + std::vector > 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::cell_iterator &cell, + const Quadrature &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::cell_iterator &cell, + const Quadrature &quadrature, + const InternalData &fedata) const; +}; + + + +/*@}*/ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 008e03b53e..14a1aa8eaf 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -59,6 +59,7 @@ #include #include #include +#include #include #include #include @@ -4787,7 +4788,8 @@ namespace VectorTools && (fe.system_to_base_index (cell_idx).first == base_indices) && (lower_bound <= fe.system_to_base_index (cell_idx).second) && (fe.system_to_base_index (cell_idx).second <= upper_bound)) - || ((dynamic_cast*> (&fe) != nullptr) + || (((dynamic_cast*> (&fe) != nullptr) + || (dynamic_cast*> (&fe) != nullptr)) && (line * (degree + 1) <= face_idx) && (face_idx <= (line + 1) * (degree + 1) - 1))) { @@ -5005,7 +5007,8 @@ namespace VectorTools const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face); if (((dynamic_cast*> (&fe) != nullptr) && (fe.system_to_base_index (cell_idx).first == base_indices)) - || (dynamic_cast*> (&fe) != nullptr)) + || (dynamic_cast*> (&fe) != nullptr) + || (dynamic_cast*> (&fe) != nullptr)) { associated_edge_dof_to_face_dof[associated_edge_dof_index] = face_idx; ++associated_edge_dof_index; @@ -5122,7 +5125,8 @@ namespace VectorTools && (fe.system_to_base_index (cell_idx).first == base_indices) && (lower_bound <= fe.system_to_base_index (cell_idx).second) && (fe.system_to_base_index (cell_idx).second <= upper_bound)) - || ((dynamic_cast*> (&fe) != nullptr) + || (((dynamic_cast*> (&fe) != nullptr) + || (dynamic_cast*> (&fe) != nullptr)) && (line * (degree + 1) <= face_idx) && (face_idx <= (line + 1) * (degree + 1) - 1))) { @@ -5159,7 +5163,8 @@ namespace VectorTools const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face); if (((dynamic_cast*> (&fe) != nullptr) && (fe.system_to_base_index (cell_idx).first == base_indices)) - || ((dynamic_cast*> (&fe) != nullptr))) + || (dynamic_cast*> (&fe) != nullptr) + || (dynamic_cast*> (&fe) != nullptr)) { associated_face_dof_to_face_dof[associated_face_dof_index] = face_idx; ++associated_face_dof_index; @@ -5357,7 +5362,8 @@ namespace VectorTools // If the FE is a FESystem we cannot check this. if (dynamic_cast*> (&cell->get_fe ()) == nullptr) { - AssertThrow (dynamic_cast*> (&cell->get_fe ()) != nullptr, + AssertThrow ((dynamic_cast*> (&cell->get_fe ()) != nullptr) + || (dynamic_cast*> (&cell->get_fe ()) != nullptr), typename FiniteElement::ExcInterpolationNotImplemented ()); } @@ -5453,7 +5459,8 @@ namespace VectorTools // If the FE is a FESystem we cannot check this. if (dynamic_cast*> (&cell->get_fe ()) == nullptr) { - AssertThrow (dynamic_cast*> (&cell->get_fe ()) != nullptr, + AssertThrow ((dynamic_cast*> (&cell->get_fe ()) != nullptr) + || (dynamic_cast*> (&cell->get_fe ()) != nullptr), typename FiniteElement::ExcInterpolationNotImplemented ()); } diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index f8fd611253..8ec617c984 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -30,6 +30,7 @@ SET(_unity_include_src fe_dg_vector.cc fe_face.cc fe_nedelec.cc + fe_nedelec_sz.cc fe_nothing.cc fe_poly.cc fe_poly_tensor.cc @@ -94,6 +95,7 @@ SET(_inst fe_face.inst.in fe.inst.in fe_nedelec.inst.in + fe_nedelec_sz.inst.in fe_nothing.inst.in fe_poly.inst.in fe_poly_tensor.inst.in diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc new file mode 100644 index 0000000000..ffe25cff29 --- /dev/null +++ b/source/fe/fe_nedelec_sz.cc @@ -0,0 +1,1839 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +// Constructor: +template +FE_NedelecSZ::FE_NedelecSZ (const unsigned int degree) + : + FiniteElement + (FiniteElementData (get_dpo_vector (degree), dim, degree + 1, + FiniteElementData::Hcurl), + std::vector (compute_num_dofs(degree), true), + std::vector (compute_num_dofs(degree), + std::vector (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; compn_components() ; ++comp) + { + this->component_to_base_table[comp].first.second = comp; + } + + // Generate the 1-D polynomial basis. + create_polynomials(degree); +} + +// Shape functions: +template +double +FE_NedelecSZ::shape_value (const unsigned int /*i*/, + const Point &/*p*/) const +{ + typedef FiniteElement FEE; + Assert(false, typename FEE::ExcFENotPrimitive()); + return 0.; +} + + + +template +double +FE_NedelecSZ::shape_value_component (const unsigned int /*i*/, + const Point &/*p*/, + const unsigned int /*component*/) const +{ + // Not implemented yet: + Assert (false, ExcNotImplemented ()); + return 0.; +} + + + +template +Tensor<1,dim> +FE_NedelecSZ::shape_grad (const unsigned int /*i*/, + const Point &/*p*/) const +{ + typedef FiniteElement FEE; + Assert(false, typename FEE::ExcFENotPrimitive()); + return Tensor<1,dim>(); +} + + + +template +Tensor<1,dim> +FE_NedelecSZ::shape_grad_component (const unsigned int /*i*/, + const Point &/*p*/, + const unsigned int /*component*/) const +{ + Assert (false, ExcNotImplemented ()); + return Tensor<1,dim> (); +} + + + +template +Tensor<2,dim> +FE_NedelecSZ::shape_grad_grad (const unsigned int /*i*/, + const Point &/*p*/) const +{ + typedef FiniteElement FEE; + Assert(false, typename FEE::ExcFENotPrimitive()); + return Tensor<2,dim>(); +} + + + +template +Tensor<2,dim> +FE_NedelecSZ::shape_grad_grad_component (const unsigned int /*i*/, + const Point &/*p*/, + const unsigned int /*component*/) const +{ + Assert (false, ExcNotImplemented ()); + return Tensor<2,dim>(); +} + + + +template +std::unique_ptr::InternalDataBase> +FE_NedelecSZ::get_data ( + const UpdateFlags update_flags, + const Mapping &/*mapping*/, + const Quadrature &quadrature, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &/*output_data*/) const +{ + std::unique_ptr 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::vertices_per_cell; + const unsigned int lines_per_cell = GeometryInfo::lines_per_cell; + const unsigned int faces_per_cell = GeometryInfo::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> (vertices_per_cell, + std::vector (vertices_per_cell))); + + data->sigma_imj_grads.resize( + vertices_per_cell, + std::vector> (vertices_per_cell, + std::vector (dim))); + + // Resize shape function arrays according to update flags: + if (flags & update_values) + { + data->shape_values.resize( + this->dofs_per_cell, std::vector> (n_q_points)); + } + + if (flags & update_gradients) + { + data->shape_grads.resize( + this->dofs_per_cell, std::vector > (n_q_points)); + } + // Not implementing second derivatives yet: + if (flags & update_hessians) + { + Assert(false, ExcNotImplemented()); + } + + std::vector< Point > 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> sigma(n_q_points, std::vector (lines_per_cell)); + std::vector> lambda(n_q_points, std::vector (lines_per_cell)); + + for (unsigned int q=0; qsigma_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; isigma_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; medge_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 (n_q_points)); + data->edge_lambda_grads_2d.resize(lines_per_cell, std::vector (dim)); + for (unsigned int m=0; m::line_to_cell_vertices(m,1)); + const unsigned int e2 (GeometryInfo::line_to_cell_vertices(m,0)); + for (unsigned int q=0; qedge_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 > cell_points (n_q_points); + for (unsigned int q=0; q > polyx(degree, std::vector (poly_length)); + std::vector > polyy(degree, std::vector (poly_length)); + for (unsigned int i=0; ishape_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; jshape_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> sigma(n_q_points, std::vector (lines_per_cell)); + std::vector> lambda(n_q_points, std::vector (lines_per_cell)); + for (unsigned int q=0; qsigma_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; isigma_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; medge_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; qedge_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; dedge_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::line_to_cell_vertices(m,1)); + const unsigned int e2 (GeometryInfo::line_to_cell_vertices(m,0)); + + for (unsigned int q=0; qedge_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; qedge_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; medge_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; qface_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 > cell_points (n_q_points); + for (unsigned int q=0; q > polyx(degree, std::vector (poly_length)); + std::vector > polyy(degree, std::vector (poly_length)); + std::vector > polyz(degree, std::vector (poly_length)); + for (unsigned int i=0; ishape_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; kshape_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; dshape_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; d1shape_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 +void FE_NedelecSZ::fill_edge_values(const typename Triangulation::cell_iterator &cell, + const Quadrature &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::lines_per_cell); + + // Calculate edge orderings: + std::vector> edge_order(lines_per_cell, + std::vector (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 edge_sign(lines_per_cell); + for (unsigned int m=0; m::line_to_cell_vertices(m,0); + unsigned int v1_loc = GeometryInfo::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 > > edge_points(lines_per_cell, std::vector> (n_q_points)); + + std::vector > edge_sigma_values (fe_data.edge_sigma_values); + std::vector > edge_sigma_grads (fe_data.edge_sigma_grads); + + std::vector > edge_lambda_values (fe_data.edge_lambda_values); + std::vector > edge_lambda_grads (fe_data.edge_lambda_grads_2d); + + // Adjust the edge_sigma_* for the current cell: + for (unsigned int m=0; m (), static_cast(edge_sign[m]))); + + std::transform (edge_sigma_grads[m].begin(), edge_sigma_grads[m].end(), + edge_sigma_grads[m].begin(), + std::bind1st (std::multiplies (), static_cast(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; mdofs_per_line); + for (unsigned int q=0; q0. + std::vector > poly(degree, std::vector (poly_length)); + for (unsigned int i=1; i edge_sign(lines_per_cell); + for (unsigned int m=0; m::line_to_cell_vertices(m,0); + unsigned int v1_loc = GeometryInfo::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 > edge_sigma_values(fe_data.edge_sigma_values); + std::vector > edge_lambda_values(fe_data.edge_lambda_values); + std::vector > edge_sigma_grads(fe_data.edge_sigma_grads); + std::vector > > edge_lambda_grads(fe_data.edge_lambda_grads_3d); + std::vector > > edge_lambda_gradgrads_3d(fe_data.edge_lambda_gradgrads_3d); + + // Adjust the edge_sigma_* for the current cell: + for (unsigned int m=0; m (), static_cast(edge_sign[m]))); + std::transform (edge_sigma_grads[m].begin(), edge_sigma_grads[m].end(), + edge_sigma_grads[m].begin(), + std::bind1st (std::multiplies (), static_cast(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 > poly(degree, std::vector (poly_length)); + for (unsigned int m=0; mdofs_per_line); + for (unsigned int q=0; q0) + { + for (unsigned int i=0; i0) + { + for (unsigned int i=0; i0) + { + for (unsigned int i=0; i +void FE_NedelecSZ::fill_face_values(const typename Triangulation::cell_iterator &cell, + const Quadrature &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::vertices_per_face); + const unsigned int faces_per_cell (GeometryInfo::faces_per_cell); + + // DoF info: + const unsigned int n_line_dofs = this->dofs_per_line*GeometryInfo::lines_per_cell; + + // First we find the global face orientations on the current cell. + std::vector > face_orientation(faces_per_cell, + std::vector (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; mvertex_index(GeometryInfo::face_to_cell_vertices(m,0)); + for (unsigned int v=1; vvertex_index(GeometryInfo::face_to_cell_vertices(m,v))) + { + current_max = v; + current_glob = cell->vertex_index(GeometryInfo::face_to_cell_vertices(m,v)); + } + } + face_orientation[m][0] = GeometryInfo::face_to_cell_vertices(m,current_max); + + // f^m_2 is the vertex opposite f^m_0. + face_orientation[m][2] = GeometryInfo::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::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][0])) + > cell->vertex_index(GeometryInfo::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][1])) ) + { + face_orientation[m][1] = GeometryInfo::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][0]); + face_orientation[m][3] = GeometryInfo::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][1]); + } + else + { + face_orientation[m][1] = GeometryInfo::face_to_cell_vertices(m,vertices_adjacent_on_face[current_max][1]); + face_orientation[m][3] = GeometryInfo::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 > face_xi_values(faces_per_cell, + std::vector (n_q_points)); + std::vector > face_xi_grads(faces_per_cell, + std::vector (dim)); + std::vector > face_eta_values(faces_per_cell, + std::vector (n_q_points)); + std::vector > face_eta_grads(faces_per_cell, + std::vector (dim)); + + std::vector > face_lambda_values(faces_per_cell, + std::vector (n_q_points)); + std::vector > face_lambda_grads(faces_per_cell, + std::vector (dim)); + for (unsigned int m=0; m > polyxi(degree, std::vector (poly_length)); + std::vector > polyeta(degree, std::vector (poly_length)); + + // Loop through quad points: + for (unsigned int m=0; mdofs_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 +void FE_NedelecSZ::fill_fe_values( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity /*cell_similarity*/, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &data) const +{ + + // Convert to the correct internal data class for this FE class. + Assert (dynamic_cast (&fe_internal) != 0, + ExcInternalError()); + const InternalData &fe_data = static_cast (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> transformed_shape_values(n_q_points); + for (unsigned int dof=0; dofdofs_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> input(n_q_points); + std::vector> transformed_shape_grads(n_q_points); + for (unsigned int dof=0; dofdofs_per_cell; ++dof) + { + for (unsigned int q=0; qn_components() + + this->get_nonzero_components(dof).first_selected_component()]; + + for (unsigned int q=0; q +void FE_NedelecSZ::fill_fe_face_values ( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &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 (&fe_internal) != 0, + ExcInternalError()); + const InternalData &fe_data = static_cast (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::project_to_all_faces(quadrature), + fe_data); + if (dim == 3 && this->degree>1) + { + fill_face_values(cell, + QProjector::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::DataSetDescriptor offset + = QProjector::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> transformed_shape_values(n_q_points); + for (unsigned int dof=0; dofdofs_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> transformed_shape_grads(n_q_points); + for (unsigned int dof=0; dofdofs_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 +void FE_NedelecSZ::fill_fe_subface_values( + const typename Triangulation::cell_iterator &/*cell*/, + const unsigned int /*face_no*/, + const unsigned int /*sub_no*/, + const Quadrature &/*quadrature*/, + const Mapping &/*mapping*/, + const typename Mapping::InternalDataBase &/*mapping_internal*/, + const dealii::internal::FEValuesImplementation::MappingRelatedData &/*mapping_data*/, + const typename FiniteElement::InternalDataBase &/*fe_internal*/, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData &/*data*/) const +{ + Assert (false, ExcNotImplemented ()); +} + +template +UpdateFlags +FE_NedelecSZ::requires_update_flags(const UpdateFlags flags) const +{ + return update_once(flags) | update_each(flags); +} + +template +UpdateFlags +FE_NedelecSZ::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 +UpdateFlags +FE_NedelecSZ::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 +std::string +FE_NedelecSZ::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 +std::unique_ptr > +FE_NedelecSZ::clone () const +{ + return std_cxx14::make_unique >(*this); +} + +template +std::vector FE_NedelecSZ::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 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 +unsigned int FE_NedelecSZ::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 +void FE_NedelecSZ::create_polynomials (const unsigned int degree) +{ + // fill the 1d polynomials vector: + IntegratedLegendrePolynomials = IntegratedLegendreSZ::generate_complete_basis (degree + 1); +} + +// explicit instantiations +#include "fe_nedelec_sz.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_nedelec_sz.inst.in b/source/fe/fe_nedelec_sz.inst.in new file mode 100644 index 0000000000..dd9bef66ce --- /dev/null +++ b/source/fe/fe_nedelec_sz.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS) +{ +#if deal_II_dimension != 1 + template class FE_NedelecSZ; +#endif +} + diff --git a/tests/fe/fe_nedelec_sz_face_orientation.cc b/tests/fe/fe_nedelec_sz_face_orientation.cc new file mode 100644 index 0000000000..1577502d46 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_face_orientation.cc @@ -0,0 +1,223 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test FE_NedelecSZ<3> for meshes with faces with non-standard orientation. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +void create_reference_triangulation (Triangulation<3> &tria) +{ + std::vector 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 > 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 > vertices (12); + + vertices[0] = Point<3> (-1.0, 0.0, 0.0); + vertices[1] = Point<3> (); + vertices[2] = Point<3> (1.0, 0.0, 0.0); + vertices[3] = Point<3> (-1.0, 0.0, 1.0); + vertices[4] = Point<3> (0.0, 0.0, 1.0); + vertices[5] = Point<3> (1.0, 0.0, 1.0); + vertices[6] = Point<3> (-1.0, 1.0, 0.0); + vertices[7] = Point<3> (0.0, 1.0, 0.0); + vertices[8] = Point<3> (1.0, 1.0, 0.0); + vertices[9] = Point<3> (-1.0, 1.0, 1.0); + vertices[10] = Point<3> (0.0, 1.0, 1.0); + vertices[11] = Point<3> (1.0, 1.0, 1.0); + tria.create_triangulation (vertices, cells, SubCellData ()); +} + +void evaluate (const FiniteElement<3> &fe, const DoFHandler<3> &dof_handler) +{ + const FEValuesExtractors::Vector component (0); + const Quadrature<3> quadrature(Point<3>(0.5, 0.5, 0.5)); + FEValues<3> fe_values (fe, quadrature, update_quadrature_points | update_values | update_gradients); + + for (DoFHandler<3>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) + { + fe_values.reinit (cell); + + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + { + deallog << "DoF#" << i << ", value=[" << fe_values[component].value(i, 0) << "], curl=[" << fe_values[component].curl(i, 0) << "]" << std::endl; + } + } +} + +void run (const bool face_orientation, const bool face_flip, const bool face_rotation) +{ +// Triangulation<3> tria_ref; +// create_reference_triangulation (tria_ref); + + FE_NedelecSZ<3> fe (1); +// DoFHandler<3> dof_handler_ref (tria_ref); +// dof_handler_ref.distribute_dofs (fe); + + Triangulation<3> tria; + create_triangulation (tria, face_orientation, face_flip, face_rotation); + + DoFHandler<3> dof_handler (tria); + dof_handler.distribute_dofs (fe); + + evaluate (fe, dof_handler); +} + +int main() +{ + initlog(); + + for (int face_orientation = 0; face_orientation <= 1; ++face_orientation) + for (int face_flip = 0; face_flip <= 1; ++face_flip) + for (int face_rotation = 0; face_rotation <= 1; ++face_rotation) + { + deallog << face_orientation << face_flip << face_rotation << std::endl; + run (face_orientation, face_flip, face_rotation); + } +} diff --git a/tests/fe/fe_nedelec_sz_face_orientation.output b/tests/fe/fe_nedelec_sz_face_orientation.output new file mode 100644 index 0000000000..774a045685 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_face_orientation.output @@ -0,0 +1,873 @@ + +DEAL::000 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::001 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#1, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#3, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#9, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#11, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#17, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#19, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#21, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#23, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::010 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#1, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#3, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#5, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#9, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#11, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#15, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#17, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#19, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#21, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#23, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#32, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::011 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#1, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#3, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#5, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#7, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#9, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#11, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#13, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#15, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#17, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#19, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#21, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#23, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#32, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::100 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#1, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#3, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#7, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#9, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#11, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#13, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#17, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#19, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#21, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#23, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#28, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#32, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::101 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#1, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#3, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#7, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#9, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#11, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#13, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#17, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#19, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#21, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#23, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::110 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#5, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#7, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#13, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#15, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::111 +DEAL::DoF#0, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#1, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#3, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#5, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#9, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#11, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#15, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#17, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#19, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#21, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#23, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#35, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#39, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#43, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#47, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#0, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 -0.500000 0.00000] +DEAL::DoF#1, value=[0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#2, value=[0.00000 0.00000 -0.250000], curl=[-0.500000 0.500000 0.00000] +DEAL::DoF#3, value=[-0.250000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#4, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 0.500000] +DEAL::DoF#5, value=[0.00000 -0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#6, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 0.500000] +DEAL::DoF#7, value=[0.00000 -0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#8, value=[0.00000 0.00000 -0.250000], curl=[0.500000 -0.500000 0.00000] +DEAL::DoF#9, value=[0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#10, value=[0.00000 0.00000 -0.250000], curl=[0.500000 0.500000 0.00000] +DEAL::DoF#11, value=[-0.250000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#12, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.500000 -0.500000] +DEAL::DoF#13, value=[0.00000 0.250000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#14, value=[-0.250000 0.00000 0.00000], curl=[0.00000 -0.500000 -0.500000] +DEAL::DoF#15, value=[0.00000 0.250000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#16, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 0.500000] +DEAL::DoF#17, value=[0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#18, value=[0.00000 -0.250000 0.00000], curl=[-0.500000 0.00000 -0.500000] +DEAL::DoF#19, value=[-0.250000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#20, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 0.500000] +DEAL::DoF#21, value=[0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#22, value=[0.00000 -0.250000 0.00000], curl=[0.500000 0.00000 -0.500000] +DEAL::DoF#23, value=[-0.250000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#24, value=[-0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#25, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#26, value=[0.00000 0.00000 -0.500000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#27, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#28, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#29, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#30, value=[0.00000 0.00000 -0.500000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#31, value=[0.00000 -0.500000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#32, value=[0.00000 0.00000 -0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#33, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#34, value=[-0.500000 0.00000 0.00000], curl=[0.00000 1.00000 0.00000] +DEAL::DoF#35, value=[0.00000 -0.500000 0.00000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#36, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#37, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#38, value=[-0.500000 0.00000 0.00000], curl=[0.00000 -1.00000 0.00000] +DEAL::DoF#39, value=[0.00000 -0.500000 0.00000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#40, value=[0.00000 0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#41, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#42, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 1.00000] +DEAL::DoF#43, value=[0.00000 0.00000 -0.500000], curl=[-1.00000 0.00000 0.00000] +DEAL::DoF#44, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#45, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#46, value=[-0.500000 0.00000 0.00000], curl=[0.00000 0.00000 -1.00000] +DEAL::DoF#47, value=[0.00000 0.00000 -0.500000], curl=[1.00000 0.00000 0.00000] +DEAL::DoF#48, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#49, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#50, value=[0.00000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#51, value=[0.250000 0.00000 0.00000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#52, value=[0.00000 0.00000 0.250000], curl=[0.00000 0.00000 0.00000] +DEAL::DoF#53, value=[0.00000 -0.250000 0.00000], curl=[0.00000 0.00000 0.00000] diff --git a/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc new file mode 100644 index 0000000000..84ce85cb03 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.cc @@ -0,0 +1,165 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check the correctness of fe_values.shape_gradient for FE_NedelecSZ by comparing +// the integral of all shape gradients with the flux over the boundary by the +// divergence theorem + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +template +Tensor<1,dim> ones () +{ + Tensor<1,dim> result; + for (unsigned int i=0; i +void test (const Triangulation &tr, + const FiniteElement &fe, + const double tolerance) +{ + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + std::stringstream ss; + + deallog << "FE=" << fe.get_name() << std::endl; + + const QGauss quadrature(4); + FEValues fe_values (fe, quadrature, update_gradients + | update_quadrature_points + | update_JxW_values); + + const QGauss face_quadrature(4); + FEFaceValues fe_face_values (fe, face_quadrature, + update_values + | update_quadrature_points + | update_normal_vectors + | update_JxW_values); + + for (typename DoFHandler::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::vertices_per_cell; ++i) + { + deallog << i << ": ( "; + for (unsigned int d=0; dvertex(i)[d] << " "; + deallog << ")" << std::endl; + } + + bool cell_ok = true; + + for (unsigned int c=0; c +void test_hyper_ball(const double tolerance) +{ + Triangulation tr; + GridGenerator::hyper_ball(tr); + + static const SphericalManifold boundary; + tr.set_manifold (0, boundary); + +// tr.refine_global(1); // taking too long, + + FE_NedelecSZ fe(1); + test(tr, fe, tolerance); +} + + +int main() +{ + std::ofstream logfile ("output"); + deallog << std::setprecision (8); + + deallog.attach(logfile); + + test_hyper_ball<2>(1e-6); + test_hyper_ball<3>(1e-6); + + deallog << "done..." << std::endl; +} + diff --git a/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.output b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.output new file mode 100644 index 0000000000..6c5fcb6f42 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_gradient_divergence_theorem.output @@ -0,0 +1,116 @@ + +DEAL::FE=FE_NedelecSZ<2>(1) +DEAL::Cell nodes: +DEAL::0: ( -0.70710678 -0.70710678 ) +DEAL::1: ( 0.70710678 -0.70710678 ) +DEAL::2: ( -0.29289322 -0.29289322 ) +DEAL::3: ( 0.29289322 -0.29289322 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.70710678 -0.70710678 ) +DEAL::1: ( -0.29289322 -0.29289322 ) +DEAL::2: ( -0.70710678 0.70710678 ) +DEAL::3: ( -0.29289322 0.29289322 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.29289322 -0.29289322 ) +DEAL::1: ( 0.29289322 -0.29289322 ) +DEAL::2: ( -0.29289322 0.29289322 ) +DEAL::3: ( 0.29289322 0.29289322 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( 0.70710678 -0.70710678 ) +DEAL::1: ( 0.70710678 0.70710678 ) +DEAL::2: ( 0.29289322 -0.29289322 ) +DEAL::3: ( 0.29289322 0.29289322 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.70710678 0.70710678 ) +DEAL::1: ( -0.29289322 0.29289322 ) +DEAL::2: ( 0.70710678 0.70710678 ) +DEAL::3: ( 0.29289322 0.29289322 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::FE=FE_NedelecSZ<3>(1) +DEAL::Cell nodes: +DEAL::0: ( -0.21132487 -0.21132487 -0.21132487 ) +DEAL::1: ( 0.21132487 -0.21132487 -0.21132487 ) +DEAL::2: ( -0.21132487 0.21132487 -0.21132487 ) +DEAL::3: ( 0.21132487 0.21132487 -0.21132487 ) +DEAL::4: ( -0.21132487 -0.21132487 0.21132487 ) +DEAL::5: ( 0.21132487 -0.21132487 0.21132487 ) +DEAL::6: ( -0.21132487 0.21132487 0.21132487 ) +DEAL::7: ( 0.21132487 0.21132487 0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.57735027 -0.57735027 -0.57735027 ) +DEAL::1: ( 0.57735027 -0.57735027 -0.57735027 ) +DEAL::2: ( -0.57735027 0.57735027 -0.57735027 ) +DEAL::3: ( 0.57735027 0.57735027 -0.57735027 ) +DEAL::4: ( -0.21132487 -0.21132487 -0.21132487 ) +DEAL::5: ( 0.21132487 -0.21132487 -0.21132487 ) +DEAL::6: ( -0.21132487 0.21132487 -0.21132487 ) +DEAL::7: ( 0.21132487 0.21132487 -0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( 0.57735027 -0.57735027 -0.57735027 ) +DEAL::1: ( 0.57735027 0.57735027 -0.57735027 ) +DEAL::2: ( 0.21132487 -0.21132487 -0.21132487 ) +DEAL::3: ( 0.21132487 0.21132487 -0.21132487 ) +DEAL::4: ( 0.57735027 -0.57735027 0.57735027 ) +DEAL::5: ( 0.57735027 0.57735027 0.57735027 ) +DEAL::6: ( 0.21132487 -0.21132487 0.21132487 ) +DEAL::7: ( 0.21132487 0.21132487 0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.57735027 -0.57735027 0.57735027 ) +DEAL::1: ( 0.57735027 -0.57735027 0.57735027 ) +DEAL::2: ( -0.21132487 -0.21132487 0.21132487 ) +DEAL::3: ( 0.21132487 -0.21132487 0.21132487 ) +DEAL::4: ( -0.57735027 0.57735027 0.57735027 ) +DEAL::5: ( 0.57735027 0.57735027 0.57735027 ) +DEAL::6: ( -0.21132487 0.21132487 0.21132487 ) +DEAL::7: ( 0.21132487 0.21132487 0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.57735027 -0.57735027 -0.57735027 ) +DEAL::1: ( -0.21132487 -0.21132487 -0.21132487 ) +DEAL::2: ( -0.57735027 0.57735027 -0.57735027 ) +DEAL::3: ( -0.21132487 0.21132487 -0.21132487 ) +DEAL::4: ( -0.57735027 -0.57735027 0.57735027 ) +DEAL::5: ( -0.21132487 -0.21132487 0.21132487 ) +DEAL::6: ( -0.57735027 0.57735027 0.57735027 ) +DEAL::7: ( -0.21132487 0.21132487 0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.57735027 -0.57735027 -0.57735027 ) +DEAL::1: ( 0.57735027 -0.57735027 -0.57735027 ) +DEAL::2: ( -0.21132487 -0.21132487 -0.21132487 ) +DEAL::3: ( 0.21132487 -0.21132487 -0.21132487 ) +DEAL::4: ( -0.57735027 -0.57735027 0.57735027 ) +DEAL::5: ( 0.57735027 -0.57735027 0.57735027 ) +DEAL::6: ( -0.21132487 -0.21132487 0.21132487 ) +DEAL::7: ( 0.21132487 -0.21132487 0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::Cell nodes: +DEAL::0: ( -0.57735027 0.57735027 -0.57735027 ) +DEAL::1: ( -0.21132487 0.21132487 -0.21132487 ) +DEAL::2: ( 0.57735027 0.57735027 -0.57735027 ) +DEAL::3: ( 0.21132487 0.21132487 -0.21132487 ) +DEAL::4: ( -0.57735027 0.57735027 0.57735027 ) +DEAL::5: ( -0.21132487 0.21132487 0.21132487 ) +DEAL::6: ( 0.57735027 0.57735027 0.57735027 ) +DEAL::7: ( 0.21132487 0.21132487 0.21132487 ) +DEAL::OK: cell bulk and boundary integrals match... + +DEAL::done... diff --git a/tests/fe/fe_nedelec_sz_non_rect_2d.cc b/tests/fe/fe_nedelec_sz_non_rect_2d.cc new file mode 100644 index 0000000000..028b4c949c --- /dev/null +++ b/tests/fe/fe_nedelec_sz_non_rect_2d.cc @@ -0,0 +1,322 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +// +// By Ross Kynch, Alexander Grayver +// +// Test to confirm that FE_NedelecSZ works on deformed elements in 2D when using +// Dirichlet boundary conditions (n x E = f, type). This is handled by the +// function project_boundary_values_curl_conforming_l2(). +// +// This test solves the real valued curl-curl equation in 2D: +// +// curl(curl(E)) + E = Js +// +// where the solution is: +// E = (0, x^2) +// +// so, Js = (0,-2) + E. +// +// i.e the solution should be "exact" at p=2. +// +// The domain is a distorted quad after 1 refinement. The distortion is +// performed by GridTools::distort_random. + +#include "../tests.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + + +using namespace dealii; + +namespace polytest +{ + template + class SimplePolynomial : public Function + { + public: + SimplePolynomial() : Function() {} + void vector_value_list (const std::vector > &points, + std::vector > &values) const; + + void rhs_value_list (const std::vector > &points, + std::vector > &values) const; + + }; + template + void SimplePolynomial::vector_value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (dim == 2, ExcNotImplemented()); + Assert (values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + for (unsigned int i=0; i &p = points[i]; + // non-zero curl-curl: + values[i][0] = 0.0; + values[i][1] = p[0]*p[0]; + + } + + } + + template + void SimplePolynomial::rhs_value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (dim == 2, ExcNotImplemented()); + Assert (values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i=0; i &p = points[i]; + // non-zero curl-curl: + values[i][0] = 0.0; + values[i][1] = -2.0 + p[0]*p[0]; + } + + } + + template + 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 tria; + DoFHandler dof_handler; + FE_NedelecSZ fe; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + }; + + template + polytest::polytest (unsigned int degree) + : + p_order(degree), + quad_order(2*degree + 3), + dof_handler(tria), + fe (degree) + { + } + template + polytest::~polytest () + { + dof_handler.clear (); + } + + template + void polytest::setup_system() + { + dof_handler.distribute_dofs (fe); + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + + constraints.clear (); + SimplePolynomial 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 + void polytest::assemble_system() + { + const QGauss test_quad(quad_order); + FEValues fe_values_test (fe, test_quad, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const QGauss quadrature_formula(quad_order); + const unsigned int n_q_points = quadrature_formula.size(); + + const QGauss face_quadrature_formula(quad_order); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + SimplePolynomial right_hand_side; + std::vector > rhs_value_list(n_q_points, + Vector(fe.n_components())); + + typename DoFHandler::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 rhs_value; + for (unsigned int d=0; dget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } + + template + void polytest::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 + void polytest::output_error() + { + SimplePolynomial exact_solution; + Vector diff_per_cell(tria.n_active_cells()); + VectorTools::integrate_difference(dof_handler, solution, exact_solution, + diff_per_cell, QGauss(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 + void polytest::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 poly(p); + poly.run(); + } +} + diff --git a/tests/fe/fe_nedelec_sz_non_rect_2d.output b/tests/fe/fe_nedelec_sz_non_rect_2d.output new file mode 100644 index 0000000000..80eb253a95 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_non_rect_2d.output @@ -0,0 +1,10 @@ + +DEAL:cg::Starting value 0.747607 +DEAL:cg::Convergence step 4 value 6.71226e-19 +DEAL::p=0 L2_error: 0.864275 +DEAL:cg::Starting value 2.34680 +DEAL:cg::Convergence step 13 value 3.46263e-13 +DEAL::p=1 L2_error: 0.0650158 +DEAL:cg::Starting value 2.34814 +DEAL:cg::Convergence step 18 value 3.36532e-13 +DEAL::p=2 L2_error: 1.40992e-13 diff --git a/tests/fe/fe_nedelec_sz_non_rect_face.cc b/tests/fe/fe_nedelec_sz_non_rect_face.cc new file mode 100644 index 0000000000..d7317b04e5 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_non_rect_face.cc @@ -0,0 +1,468 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// +// By Ross Kynch +// +// Test to confirm that FE_NedelecSZ works on non-rectangular faces in 3D. +// The previous project_boundary_values_curl_conforming function only +// produced the correct constraints on rectangular faces. +// +// The new function project_boundary_values_curl_conforming_l2 fixes this +// and should work on any quad face, assuming the underlying mesh conforms to the +// standard orientation (i.e. spheres may still cause problems). +// +// This test solves the real valued curl-curl equation in 3D: +// +// curl(curl(E)) + E = Js +// +// where the solution is: +// E = (x^2 , y^2 + , z^2). +// +// so: +// Js = E. +// +// The domain is a distorted cube which will contain non-rectangular faces. + +#include "../tests.h" +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace Maxwell +{ + using namespace dealii; + + // Dirichlet BCs / exact solution:. + template + class ExactSolution : public Function + { + public: + ExactSolution(); + + virtual void vector_value_list (const std::vector > &points, + std::vector > &values) const; + + + void curl_value_list(const std::vector > &points, + std::vector > &value_list); + }; + + template + ExactSolution::ExactSolution() + : + Function (dim) + {} + template + void ExactSolution::vector_value_list (const std::vector > &points, + std::vector > &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 &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 + void ExactSolution::curl_value_list(const std::vector > &points, + std::vector > &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 &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 + 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 triangulation; + MappingQ mapping; + DoFHandler dof_handler; + FE_NedelecSZ fe; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + ConvergenceTable convergence_table; + + ExactSolution exact_solution; + + unsigned int p_order; + unsigned int quad_order; + }; + + template + MaxwellProblem::MaxwellProblem (const unsigned int order) + : + mapping (1), + dof_handler (triangulation), + fe(order), + exact_solution() + { + p_order = order; + quad_order = 2*(p_order+1); + } + + template + MaxwellProblem::~MaxwellProblem () + { + dof_handler.clear (); + } + template + double MaxwellProblem::calcErrorHcurlNorm() + { + QGauss quadrature_formula(quad_order); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues 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 local_dof_indices (dofs_per_cell); + + // storage for exact sol: + std::vector > exactsol_list(n_q_points, Vector(fe.n_components())); + std::vector > exactcurlsol_list(n_q_points, Vector(fe.n_components())); + Tensor<1,dim> exactsol; + Tensor<1,dim> exactcurlsol; + + // storage for computed sol: + std::vector > sol(n_q_points); + Tensor<1,dim> curlsol; + + double h_curl_norm=0.0; + + unsigned int block_index_i; + + typename DoFHandler::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 + void MaxwellProblem::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 + void MaxwellProblem::assemble_system () + { + QGauss quadrature_formula(quad_order); + QGauss 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 fe_values (mapping, fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + //RHS storage: + std::vector > rhs_value_list(n_q_points, Vector(fe.n_components())); + Tensor<1,dim> rhs_value_vector; + + // Neumann storage + std::vector > neumann_value_list(n_face_q_points, Vector(fe.n_components())); + std::vector > 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::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_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } + template + void MaxwellProblem::solve () + { + /* Direct */ + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + + A_direct.vmult (solution, system_rhs); + constraints.distribute (solution); + + } + template + void MaxwellProblem::process_solution(const unsigned int cycle) + { + + Vector diff_per_cell(triangulation.n_active_cells()); + + + VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution, + diff_per_cell, QGauss(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 + void MaxwellProblem::run () + { + unsigned int numcycles=(p_order>1)?2:3; + + for (unsigned int cycle=0; cycleset_refine_flag(); // create single refined cell + * triangulation.execute_coarsening_and_refinement(); + */ + + } + else + { + triangulation.refine_global (1); + } + + setup_system (); + assemble_system (); + solve (); + process_solution (cycle); + } + + deallog << "p = " << p_order << std::endl; + + //convergence_table.set_precision("L2 Error",4); + convergence_table.set_scientific("L2 Error",true); + + //convergence_table.set_precision("H(curl) Error",4); + convergence_table.set_scientific("H(curl) Error",true); + + convergence_table.write_text(deallog.get_file_stream()); + + } +} +// END MAXWELL CLASS +int main () +{ + using namespace Maxwell; + + initlog(); + + for (unsigned int p=0; p<3; ++p) + { + MaxwellProblem<3> (p).run(); + } + + return 0; +} diff --git a/tests/fe/fe_nedelec_sz_non_rect_face.with_umfpack=on.output b/tests/fe/fe_nedelec_sz_non_rect_face.with_umfpack=on.output new file mode 100644 index 0000000000..0a838e30e1 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_non_rect_face.with_umfpack=on.output @@ -0,0 +1,15 @@ + +DEAL::p = 0 +cycle cells dofs L2 Error H(curl) Error +0 1 12 5.8761e-01 5.7154e-01 +1 8 54 3.1310e-01 3.1132e-01 +2 64 300 1.5881e-01 1.5860e-01 +DEAL::p = 1 +cycle cells dofs L2 Error H(curl) Error +0 1 54 1.3757e-01 1.4097e-01 +1 8 300 3.4151e-02 3.5926e-02 +2 64 1944 8.5308e-03 9.4069e-03 +DEAL::p = 2 +cycle cells dofs L2 Error H(curl) Error +0 1 144 5.3517e-16 7.5744e-16 +1 8 882 6.0713e-15 6.5454e-15