]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added project_boundary_values_curl_conforming_l2() to VectorTools for 197/head
authorRoss Kynch <w0ss4g3@gmail.com>
Tue, 14 Oct 2014 13:54:57 +0000 (14:54 +0100)
committerRoss Kynch <w0ss4g3@gmail.com>
Tue, 24 Feb 2015 13:54:54 +0000 (13:54 +0000)
use with FE_Nedelec elements in 2D & 3D.

This includes a test code included in tests/deal.II to check it works correctly.

This new function uses an L2-projection on edges and
faces (only edges in 2D, both for 3D) to compute
the constaints and factors in the residual from the edge
projection when computing the face constraints. The original
function project_boundary_values_curl_conforming() has
problems with non-rectangular faces, this one does not. The
function assumes that the face orientations are standard (e.g. a
sphere may cause problems due to mismatched face orientations).

The test shows that project_boundary_curl_conforming_l2() works on a
distorted cube (non-rectangular faces) for a simple test case
(polynomial function).

include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/vector_tools_boundary.inst.in
tests/deal.II/nedelec_non_rect_face.cc [new file with mode: 0644]
tests/deal.II/nedelec_non_rect_face.output [new file with mode: 0644]

index 1c86ec796e3cce7fc6b653ff633889ff933360ec..0274d29689bfee32d2cb528c0c7936712bdcf904 100644 (file)
@@ -1110,6 +1110,110 @@ namespace VectorTools
                                                 ConstraintMatrix &constraints,
                                                 const hp::MappingCollection<dim, dim> &mapping_collection = hp::StaticMappingQ1<dim>::mapping_collection);
 
+  /**
+   * This function is an updated version of the project_boundary_values_curl_conforming
+   * function. The intention is to fix a problem when using the previous function in
+   * conjunction with non-rectangular geometries (i.e. elements with non-rectangular faces).
+   * The L2-projection method used has been taken from the paper "Electromagnetic scattering
+   * simulation using an H (curl) conforming hp finite element method in three dimensions"
+   * by PD Ledger, K Morgan and O Hassan ( Int. J.  Num. Meth. Fluids, Volume 53, Issue 8, pages 1267–1296).
+   *
+   * This function will compute constraints that correspond to Dirichlet boundary conditions of the form
+   * $\vec{n}\times\vec{E}=\vec{n}\times\vec{F}$
+   * i.e. the tangential components of $\vec{E}$ and $f$ shall coincide.
+   *
+   * <h4>Computing constraints</h4>
+   *
+   * To compute the constraints we use a projection method based upon the paper mentioned
+   * above. In 2D this is done in a single stage for the edge-based shape functions, regardless
+   * of the order of the finite element. In 3D this is done in two stages, edges first and then
+   * faces.
+   *
+   * For each cell, each edge, $e$, is projected by solving the linear system $Ax=b$ where $x$
+   * is the vector of contraints on degrees of freedom on the edge and
+   *
+   * $A_{ij} = \int_{e} (\vec{s}_{i}\cdot\vec{t})(\vec{s}_{j}\cdot\vec{t}) dS$
+   *
+   * $b_{i} = \int_{e} (\vec{s}_{i}\cdot\vec{t})(\vec{F}\cdot\vec{t}) dS$
+   *
+   * with $\vec{s}_{i}$ the $i^{th}$ shape function and $\vec{t}$ the tangent vector.
+   *
+   * Once all edge constraints, $x$, have been computed, we may compute the face constraints
+   * in a similar fashion, taking into account the residuals from the edges.
+   *
+   * For each face on the cell, $f$, we solve the linear system By=c where $y$ is the vector of
+   * constraints on degrees of freedom on the face and
+   *
+   * B_{ij} = \int_{f} (\vec{n} \cross \vec{s}_{i}) \cdot (\vec{n} \cross \vec{s}_{j}) dS
+   *
+   * $c_{i} = \int_{f} (\vec{n} \cross \vec{r}) \cdot (\vec{n} \cross \vec{s}_i} dS
+   *
+   * and $\vec{r} = \vec{F} - \sum_{e \in f} \sum{i \in e} \x_{i}\vec{s}_i}$, the edge residual.
+   *
+   * The resulting constraints are then given in the solutions $x$ and $y$.
+   *
+   * If the ConstraintMatrix @p constraints contained values or other
+   * constraints before, the new ones are added or the old ones overwritten,
+   * if a node of the boundary part to be used was already in the list of
+   * constraints. This is handled by using inhomogeneous constraints. Please
+   * note that when combining adaptive meshes and this kind of constraints, the
+   * Dirichlet conditions should be set first, and then completed by hanging
+   * node constraints, in order to make sure that the discretization remains
+   * consistent. See the discussion on conflicting constraints in the
+   * module on @ref constraints
+   *
+   * <h4>Arguments to this function></h4>
+   *
+   * This function is explicitly for use with FE_Nedelec elements, or with FESystem
+   * elements which contain FE_Nedelec elements. It will throw an exception if called
+   * with any other finite element. The user must ensure that FESystem elements are
+   * correctly setup when using this function as this check not possible in this case.
+   *
+   * The second argument of this function denotes the first vector component of the
+   * finite element which corresponds to the vector function that you wish to constrain.
+   * For example, if we are solving Maxwell's equations in 3D and have components
+   * $(E_x,E_y,E_z,B_x,B_y,B_z)$ and we want the boundary conditions
+   * $\vec{n}\times\vec{B}=\vec{n}\times\vec{f}$, then @p first_vector_component would
+   * be 3. The @p boundary_function must return 6 components in this example, with the first 3
+   * corresponding to $\vec{E}$ and the second 3 corresponding to $\vec{B}$.
+   * Vectors are implicitly assumed to have exactly <code>dim</code> components
+   * that are ordered in the same way as we usually order the coordinate directions,
+   * i.e. $x$-, $y$-, and finally $z$-component.
+   *
+   * The parameter @p boundary_component corresponds to the number @p boundary_indicator
+   * of the face. numbers::internal_face_boundary_id is an illegal value, since it is
+   * reserved for interior faces.
+   *
+   * The last argument is denoted to compute the normal vector $\vec{n}$ at the
+   * boundary points.
+   *
+   *
+   * @ingroup constraints
+   *
+   * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   */
+  template <int dim>
+  void project_boundary_values_curl_conforming_l2 (const DoFHandler<dim> &dof_handler,
+                                                   const unsigned int first_vector_component,
+                                                   const Function<dim,double> &boundary_function,
+                                                   const types::boundary_id boundary_component,
+                                                   ConstraintMatrix &constraints,
+                                                   const Mapping<dim> &mapping = StaticMappingQ1<dim>::mapping);
+
+
+  /**
+   * hp-namespace version of project_boundary_values_curl_conforming_l2 (above).
+   *
+   * @ingroup constraints
+   */
+  template <int dim>
+  void project_boundary_values_curl_conforming_l2 (const hp::DoFHandler<dim> &dof_handler,
+                                                   const unsigned int first_vector_component,
+                                                   const Function<dim,double> &boundary_function,
+                                                   const types::boundary_id boundary_component,
+                                                   ConstraintMatrix &constraints,
+                                                   const hp::MappingCollection<dim, dim> &mapping_collection = hp::StaticMappingQ1<dim>::mapping_collection);
+
 
   /**
    * Compute constraints that correspond to boundary conditions of the form
index d5db44f78f745ca5cfc1c20e8521788e6f178dcf..7c5f9ee516832c4b2362fd314558c2679f2c79d8 100644 (file)
@@ -3801,6 +3801,926 @@ namespace VectorTools
   }
 
 
+  namespace internals
+  {
+    template<typename cell_iterator>
+    void
+    compute_edge_projection_l2 (const cell_iterator &cell,
+                                const unsigned int face,
+                                const unsigned int line,
+                                hp::FEValues<3> &hp_fe_values,
+                                const Function<3> &boundary_function,
+                                const unsigned int first_vector_component,
+                                std::vector<double> &dof_values,
+                                std::vector<bool> &dofs_processed)
+    {
+      // This function computes the L2-projection of the given
+      // boundary function on 3D edges and returns the constraints
+      // associated with the edge functions for the given cell.
+      //
+      // In the context of this function, by associated DoFs we mean:
+      // the DoFs corresponding to the group of components making up the vector
+      // with first component first_vector_component (length dim).
+      const unsigned int dim = 3;
+      const unsigned int spacedim = 3;
+      const FiniteElement<dim> &fe = cell->get_fe ();
+
+      // reinit for this cell, face and line.
+      hp_fe_values.reinit
+      (cell,
+       (cell->active_fe_index () * GeometryInfo<dim>::faces_per_cell + face)
+       * GeometryInfo<dim>::lines_per_face + line);
+
+      // Initialize the required objects.
+      const FEValues<dim> &
+      fe_values = hp_fe_values.get_present_fe_values ();
+      // Store degree as fe.degree-1
+      // For nedelec elements FE_Nedelec<dim> (0) returns fe.degree = 1.
+      const unsigned int degree = fe.degree - 1;
+
+      const std::vector<Point<dim> > &
+      quadrature_points = fe_values.get_quadrature_points ();
+      std::vector<Vector<double> > values (fe_values.n_quadrature_points,
+                                           Vector<double> (fe.n_components ()));
+
+      // Get boundary function values
+      // at quadrature points.
+      boundary_function.vector_value_list (quadrature_points, values);
+
+      // Find the group of vector components (dim of them,
+      // starting at first_vector_component) are within an FESystem.
+      //
+      // If not using FESystem then must be using FE_Nedelec,
+      // which has one base element and one copy of it (with 3 components).
+      std::pair<unsigned int, unsigned int> base_indices (0, 0);
+      if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) != 0)
+        {
+          unsigned int fe_index = 0;
+          unsigned int fe_index_old = 0;
+          unsigned int i = 0;
+
+          // Find base element:
+          // base_indices.first
+          //
+          // Then select which copy of that base element
+          // [ each copy is of length fe.base_element(base_indices.first).n_components() ]
+          // corresponds to first_vector_component:
+          // base_index.second
+          for (; i < fe.n_base_elements (); ++i)
+            {
+              fe_index_old = fe_index;
+              fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components ();
+
+              if (fe_index > first_vector_component)
+                break;
+            }
+
+          base_indices.first = i;
+          base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components ();
+        }
+
+      // Find DoFs we want to constrain:
+      // There are fe.dofs_per_line DoFs associated with the
+      // given line on the given face on the given cell.
+      //
+      // Want to know which of these DoFs (there are degree+1 of interest)
+      // are associated with the components given by first_vector_component.
+      // Then we can make a map from the associated line DoFs to the face DoFs.
+      //
+      // For a single FE_Nedelec<3> element this is simple:
+      //    We know the ordering of local DoFs goes
+      //    lines -> faces -> cells
+      //
+      // For a set of FESystem<3> elements we need to pick out the matching base element and
+      // the index within this ordering.
+      //
+      // We call the map associated_edge_dof_to_face_dof
+      std::vector<unsigned int> associated_edge_dof_to_face_dof (degree + 1);
+
+      // Lowest DoF in the base element allowed for this edge:
+      const unsigned int lower_bound =
+        fe.base_element(base_indices.first).face_to_cell_index(line * (degree + 1), face);
+      // Highest DoF in the base element allowed for this edge:
+      const unsigned int upper_bound =
+        fe.base_element(base_indices.first).face_to_cell_index((line + 1) * (degree + 1) - 1, face);
+
+      unsigned int associated_edge_dof_index = 0;
+      //       for (unsigned int face_idx = 0; face_idx < fe.dofs_per_face; ++face_idx)
+      for (unsigned int line_idx = 0; line_idx < fe.dofs_per_line; ++line_idx)
+        {
+          // Assuming DoFs on a face are numbered in order by lines then faces.
+          // i.e. line 0 has degree+1 dofs numbered 0,..,degree
+          //      line 1 has degree+1 dofs numbered (degree+1),..,2*(degree+1) and so on.
+          const unsigned int face_idx = line*fe.dofs_per_line + line_idx;
+          // Note, assuming that the edge orientations are "standard"
+          //       i.e. cell->line_orientation(line) = true.
+          const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face);
+
+          // Check this cell_idx belongs to the correct base_element, component and line:
+          if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+               && (fe.system_to_base_index (cell_idx).first == base_indices)
+               && (lower_bound <= fe.system_to_base_index (cell_idx).second)
+               && (fe.system_to_base_index (cell_idx).second <= upper_bound))
+              || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
+                  && (line * (degree + 1) <= face_idx)
+                  && (face_idx <= (line + 1) * (degree + 1) - 1)))
+            {
+              associated_edge_dof_to_face_dof[associated_edge_dof_index] = face_idx;
+              ++associated_edge_dof_index;
+            }
+        }
+      // Sanity check:
+      const unsigned int associated_edge_dofs = associated_edge_dof_index;
+      Assert (associated_edge_dofs == degree + 1,
+              ExcMessage ("Error: Unexpected number of 3D edge DoFs"));
+
+      // Matrix and RHS vectors to store linear system:
+      // We have (degree+1) basis functions for an edge
+      FullMatrix<double> edge_matrix (degree + 1,degree + 1);
+      FullMatrix<double> edge_matrix_inv (degree + 1,degree + 1);
+      Vector<double> edge_rhs (degree + 1);
+      Vector<double> edge_solution (degree + 1);
+
+      const FEValuesExtractors::Vector vec (first_vector_component);
+
+      // coordinate directions of
+      // the edges of the face.
+      const unsigned int
+      edge_coordinate_direction
+      [GeometryInfo<dim>::faces_per_cell]
+      [GeometryInfo<dim>::lines_per_face]
+      = { { 2, 2, 1, 1 },
+        { 2, 2, 1, 1 },
+        { 0, 0, 2, 2 },
+        { 0, 0, 2, 2 },
+        { 1, 1, 0, 0 },
+        { 1, 1, 0, 0 }
+      };
+
+      const double tol = 0.5 * cell->face (face)->line (line)->diameter () / fe.degree;
+      const std::vector<Point<dim> > &
+      reference_quadrature_points = fe_values.get_quadrature ().get_points ();
+
+      // Project the boundary function onto the shape functions for this edge
+      // and set up a linear system of equations to get the values for the DoFs
+      // associated with this edge.
+      for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
+        {
+          // Compute the tangential
+          // of the edge at
+          // the quadrature point.
+          Point<dim> shifted_reference_point_1 = reference_quadrature_points[q_point];
+          Point<dim> shifted_reference_point_2 = reference_quadrature_points[q_point];
+
+          shifted_reference_point_1 (edge_coordinate_direction[face][line]) += tol;
+          shifted_reference_point_2 (edge_coordinate_direction[face][line]) -= tol;
+          Tensor<1,dim> tangential
+            = (0.5 *
+               (fe_values.get_mapping ()
+                .transform_unit_to_real_cell (cell,
+                                              shifted_reference_point_1)
+                -
+                fe_values.get_mapping ()
+                .transform_unit_to_real_cell (cell,
+                                              shifted_reference_point_2))
+               / tol);
+          tangential
+          /= tangential.norm ();
+
+          // Compute the entires of the linear system
+          // Note the system is symmetric so we could only compute the lower/upper triangle.
+          //
+          // The matrix entries are
+          // \int_{edge} (tangential*edge_shape_function_i)*(tangential*edge_shape_function_j) dS
+          //
+          // The RHS entries are:
+          // \int_{edge} (tangential*boundary_value)*(tangential*edge_shape_function_i) dS.
+          for (unsigned int j = 0; j < associated_edge_dofs; ++j)
+            {
+              const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j];
+              const unsigned int j_cell_idx = fe.face_to_cell_index (j_face_idx, face);
+              for (unsigned int i = 0; i < associated_edge_dofs; ++i)
+                {
+                  const unsigned int i_face_idx = associated_edge_dof_to_face_dof[i];
+                  const unsigned int i_cell_idx = fe.face_to_cell_index (i_face_idx, face);
+
+                  edge_matrix(i,j)
+                  += fe_values.JxW (q_point)
+                     * (fe_values[vec].value (i_cell_idx, q_point) * tangential)
+                     * (fe_values[vec].value (j_cell_idx, q_point) * tangential);
+                }
+              // Compute the RHS entries:
+              edge_rhs(j) += fe_values.JxW (q_point)
+                             * (values[q_point] (first_vector_component) * tangential [0]
+                                + values[q_point] (first_vector_component + 1) * tangential [1]
+                                + values[q_point] (first_vector_component + 2) * tangential [2])
+                             * (fe_values[vec].value (j_cell_idx, q_point) * tangential);
+            }
+        }
+
+      // Invert linear system
+      edge_matrix_inv.invert(edge_matrix);
+      edge_matrix_inv.vmult(edge_solution,edge_rhs);
+
+      // Store computed DoFs
+      for (unsigned int i = 0; i < associated_edge_dofs; ++i)
+        {
+          dof_values[associated_edge_dof_to_face_dof[i]] = edge_solution(i);
+          dofs_processed[associated_edge_dof_to_face_dof[i]] = true;
+        }
+    }
+
+
+    template<int dim, typename cell_iterator>
+    void
+    compute_edge_projection_l2 (const cell_iterator &,
+                                const unsigned int,
+                                const unsigned int,
+                                hp::FEValues<dim> &,
+                                const Function<dim> &,
+                                const unsigned int,
+                                std::vector<double> &,
+                                std::vector<bool> &)
+    {
+      // dummy implementation of above function
+      // for all other dimensions
+      Assert (false, ExcInternalError ());
+    }
+
+    template<int dim, typename cell_iterator>
+    void
+    compute_face_projection_curl_conforming_l2 (const cell_iterator &cell,
+                                                const unsigned int face,
+                                                hp::FEFaceValues<dim> &hp_fe_face_values,
+                                                const Function<dim> &boundary_function,
+                                                const unsigned int first_vector_component,
+                                                std::vector<double> &dof_values,
+                                                std::vector<bool> &dofs_processed)
+    {
+      // This function computes the L2-projection of the boundary
+      // function on the interior of faces only. In 3D, this should only be
+      // called after first calling compute_edge_projection_l2, as it relies on
+      // edge constraints which are found.
+
+      // In the context of this function, by associated DoFs we mean:
+      // the DoFs corresponding to the group of components making up the vector
+      // with first component first_vector_component (with total components dim).
+
+      // Copy to the standard FEFaceValues object:
+      hp_fe_face_values.reinit (cell, face);
+      const FEFaceValues<dim> &fe_face_values = hp_fe_face_values.get_present_fe_values();
+
+      // Initialize the required objects.
+      const FiniteElement<dim> &fe = cell->get_fe ();
+      const std::vector<Point<dim> > &
+      quadrature_points = fe_face_values.get_quadrature_points ();
+      const unsigned int degree = fe.degree - 1;
+
+      std::vector<Vector<double> >
+      values (fe_face_values.n_quadrature_points, Vector<double> (fe.n_components ()));
+
+      // Get boundary function values at quadrature points.
+      boundary_function.vector_value_list (quadrature_points, values);
+
+      // Find where the group of vector components (dim of them,
+      // starting at first_vector_component) are within an FESystem.
+      //
+      // If not using FESystem then must be using FE_Nedelec,
+      // which has one base element and one copy of it (with 3 components).
+      std::pair<unsigned int, unsigned int> base_indices (0, 0);
+      if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) != 0)
+        {
+          unsigned int fe_index = 0;
+          unsigned int fe_index_old = 0;
+          unsigned int i = 0;
+
+          // Find base element:
+          // base_indices.first
+          //
+          // Then select which copy of that base element
+          // [ each copy is of length fe.base_element(base_indices.first).n_components() ]
+          // corresponds to first_vector_component:
+          // base_index.second
+          for (; i < fe.n_base_elements (); ++i)
+            {
+              fe_index_old = fe_index;
+              fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components ();
+
+              if (fe_index > first_vector_component)
+                break;
+            }
+          base_indices.first = i;
+          base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components ();
+        }
+
+      switch (dim)
+        {
+        case 2:
+          // NOTE: This is very similar to compute_edge_projection as used in 3D,
+          //       and contains a lot of overlap with that function.
+        {
+
+          // Find the DoFs we want to constrain. There are degree+1 in total.
+          // Create a map from these to the face index
+          // Note:
+          //    - for a single FE_Nedelec<2> element this is
+          //      simply 0 to fe.dofs_per_face
+          //    - for FESystem<2> this just requires matching the
+          //      base element, fe.system_to_base_index.first.first
+          //      and the copy of the base element we're interested
+          //      in, fe.system_to_base_index.first.second
+          std::vector<unsigned int> associated_edge_dof_to_face_dof (degree + 1);
+
+          unsigned int associated_edge_dof_index = 0;
+          for (unsigned int face_idx = 0; face_idx < fe.dofs_per_face; ++face_idx)
+            {
+              const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face);
+              if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+                   && (fe.system_to_base_index (cell_idx).first == base_indices))
+                  || (dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0))
+                {
+                  associated_edge_dof_to_face_dof[associated_edge_dof_index] = face_idx;
+                  ++associated_edge_dof_index;
+                }
+            }
+          // Sanity check:
+          const unsigned int associated_edge_dofs = associated_edge_dof_index;
+          Assert (associated_edge_dofs == degree + 1,
+                  ExcMessage ("Error: Unexpected number of 2D edge DoFs"));
+
+          // Matrix and RHS vectors to store:
+          // We have (degree+1) edge basis functions
+          FullMatrix<double> edge_matrix (degree + 1,degree + 1);
+          FullMatrix<double> edge_matrix_inv (degree + 1,degree + 1);
+          Vector<double> edge_rhs (degree + 1);
+          Vector<double> edge_solution (degree + 1);
+
+          const FEValuesExtractors::Vector vec (first_vector_component);
+
+          // coordinate directions of the face.
+          const unsigned int
+          face_coordinate_direction[GeometryInfo<dim>::faces_per_cell]
+            = { 1, 1, 0, 0 };
+
+          const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
+          const std::vector<Point<dim> > &
+          reference_quadrature_points = fe_face_values.get_quadrature_points ();
+
+          // Project the boundary function onto the shape functions for this edge
+          // and set up a linear system of equations to get the values for the DoFs
+          // associated with this edge.
+          for (unsigned int q_point = 0;
+               q_point < fe_face_values.n_quadrature_points; ++q_point)
+            {
+              // Compute the tangent of the face
+              // at the quadrature point.
+              Point<dim> shifted_reference_point_1
+                = reference_quadrature_points[q_point];
+              Point<dim> shifted_reference_point_2
+                = reference_quadrature_points[q_point];
+
+              shifted_reference_point_1 (face_coordinate_direction[face])
+              += tol;
+              shifted_reference_point_2 (face_coordinate_direction[face])
+              -= tol;
+              Tensor<1,dim> tangential
+                = (fe_face_values.get_mapping ()
+                   .transform_unit_to_real_cell (cell,
+                                                 shifted_reference_point_1)
+                   -
+                   fe_face_values.get_mapping ()
+                   .transform_unit_to_real_cell (cell,
+                                                 shifted_reference_point_2))
+                  / tol;
+              tangential
+              /= tangential.norm ();
+
+              // Compute the entires of the linear system
+              // Note the system is symmetric so we could only compute the lower/upper triangle.
+              //
+              // The matrix entries are
+              // \int_{edge} (tangential * edge_shape_function_i) * (tangential * edge_shape_function_j) dS
+              //
+              // The RHS entries are:
+              // \int_{edge} (tangential* boundary_value) * (tangential * edge_shape_function_i) dS.
+              for (unsigned int j = 0; j < associated_edge_dofs; ++j)
+                {
+                  const unsigned int j_face_idx = associated_edge_dof_to_face_dof[j];
+                  for (unsigned int i = 0; i < associated_edge_dofs; ++i)
+                    {
+                      const unsigned int i_face_idx = associated_edge_dof_to_face_dof[i];
+                      edge_matrix(i,j)
+                      += fe_face_values.JxW (q_point)
+                         * (fe_face_values[vec].value (fe.face_to_cell_index (i_face_idx, face), q_point) * tangential)
+                         * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential);
+                    }
+                  edge_rhs(j)
+                  += fe_face_values.JxW (q_point)
+                     * (values[q_point] (first_vector_component) * tangential [0]
+                        + values[q_point] (first_vector_component + 1) * tangential [1])
+                     * (fe_face_values[vec].value (fe.face_to_cell_index (j_face_idx, face), q_point) * tangential);
+                }
+            }
+
+          // Invert linear system
+          edge_matrix_inv.invert(edge_matrix);
+          edge_matrix_inv.vmult(edge_solution,edge_rhs);
+
+          // Store computed DoFs
+          for (unsigned int associated_edge_dof_index = 0; associated_edge_dof_index < associated_edge_dofs; ++associated_edge_dof_index)
+            {
+              dof_values[associated_edge_dof_to_face_dof[associated_edge_dof_index]] = edge_solution (associated_edge_dof_index);
+              dofs_processed[associated_edge_dof_to_face_dof[associated_edge_dof_index]] = true;
+            }
+          break;
+        }
+
+        case 3:
+        {
+          const FEValuesExtractors::Vector vec (first_vector_component);
+
+          // First group DoFs associated with edges which we already know.
+          // Sort these into groups of dofs (0 -> degree+1 of them) by each edge.
+          // This will help when computing the residual for the face projections.
+          //
+          // This matches with the search done in compute_edge_projection.
+          const unsigned int lines_per_face = GeometryInfo<dim>::lines_per_face;
+          std::vector<std::vector<unsigned int>>
+                                              associated_edge_dof_to_face_dof(lines_per_face, std::vector<unsigned int> (degree + 1));
+          std::vector<unsigned int> associated_edge_dofs (lines_per_face);
+
+          for (unsigned int line = 0; line < lines_per_face; ++line)
+            {
+              // Lowest DoF in the base element allowed for this edge:
+              const unsigned int lower_bound =
+                fe.base_element(base_indices.first).face_to_cell_index(line * (degree + 1), face);
+              // Highest DoF in the base element allowed for this edge:
+              const unsigned int upper_bound =
+                fe.base_element(base_indices.first).face_to_cell_index((line + 1) * (degree + 1) - 1, face);
+              unsigned int associated_edge_dof_index = 0;
+              for (unsigned int line_idx = 0; line_idx < fe.dofs_per_line; ++line_idx)
+                {
+                  const unsigned int face_idx = line*fe.dofs_per_line + line_idx;
+                  const unsigned int cell_idx = fe.face_to_cell_index(face_idx, face);
+                  // Check this cell_idx belongs to the correct base_element, component and line:
+                  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+                       && (fe.system_to_base_index (cell_idx).first == base_indices)
+                       && (lower_bound <= fe.system_to_base_index (cell_idx).second)
+                       && (fe.system_to_base_index (cell_idx).second <= upper_bound))
+                      || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
+                          && (line * (degree + 1) <= face_idx)
+                          && (face_idx <= (line + 1) * (degree + 1) - 1)))
+                    {
+                      associated_edge_dof_to_face_dof[line][associated_edge_dof_index] = face_idx;
+                      ++associated_edge_dof_index;
+                    }
+                }
+              // Sanity check:
+              associated_edge_dofs[line] = associated_edge_dof_index;
+              Assert (associated_edge_dofs[line] == degree + 1,
+                      ExcMessage ("Error: Unexpected number of 3D edge DoFs"));
+            }
+
+          // Next find the face DoFs associated with the vector components
+          // we're interested in. There are 2*degree*(degree+1) DoFs associated
+          // with each face (not including edges!).
+          //
+          // Create a map mapping from the consecutively numbered associated_dofs
+          // to the face DoF (which can be transferred to a local cell index).
+          //
+          // For FE_Nedelec<3> we just need to have a face numbering greater than
+          // the number of edge DoFs (=lines_per_face*(degree+1).
+          //
+          // For FE_System<3> we need to base the base_indices (base element and
+          // copy within that base element) and ensure we're above the number of
+          // edge DoFs within that base element.
+          std::vector<unsigned int> associated_face_dof_to_face_dof (2*degree*(degree + 1));
+
+          // Skip the edge DoFs, so we start at lines_per_face*(fe.dofs_per_line).
+          unsigned int associated_face_dof_index = 0;
+          for (unsigned int face_idx = lines_per_face*(fe.dofs_per_line); face_idx < fe.dofs_per_face; ++face_idx)
+            {
+              const unsigned int cell_idx = fe.face_to_cell_index (face_idx, face);
+              if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+                   && (fe.system_to_base_index (cell_idx).first == base_indices))
+                  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)))
+                {
+                  associated_face_dof_to_face_dof[associated_face_dof_index] = face_idx;
+                  ++associated_face_dof_index;
+                }
+            }
+          // Sanity check:
+          const unsigned int associated_face_dofs = associated_face_dof_index;
+          Assert (associated_face_dofs == 2*degree*(degree + 1),
+                  ExcMessage ("Error: Unexpected number of 3D face DoFs"));
+
+          // Storage for the linear system.
+          // There are 2*degree*(degree+1) DoFs associated with a face in 3D.
+          // Note this doesn't include the DoFs associated with edges on that face.
+          FullMatrix<double> face_matrix (2*degree*(degree + 1));
+          FullMatrix<double> face_matrix_inv (2*degree*(degree + 1));
+          Vector<double> face_rhs (2*degree*(degree + 1));
+          Vector<double> face_solution (2*degree*(degree + 1));
+
+          // Project the boundary function onto the shape functions for this face
+          // and set up a linear system of equations to get the values for the DoFs
+          // associated with this face. We also must include the residuals from the
+          // shape funcations associated with edges.
+          Tensor<1, dim> tmp;
+          Tensor<1, dim> normal_vector,
+                 cross_product_i,
+                 cross_product_j,
+                 cross_product_rhs;
+
+          // Store all normal vectors at quad points:
+          std::vector<Point<dim>> normal_vector_list(fe_face_values.get_normal_vectors());
+
+          // Loop to construct face linear system.
+          for (unsigned int q_point = 0;
+               q_point < fe_face_values.n_quadrature_points; ++q_point)
+            {
+              // First calculate the residual from the edge functions
+              // store the result in tmp.
+              //
+              // Edge_residual =
+              //        boundary_value - (
+              //            \sum_(edges on face)
+              //                 \sum_(DoFs on edge) edge_dof_value*edge_shape_function
+              //                   )
+              for (unsigned int d = 0; d < dim; ++d)
+                {
+                  tmp[d]=0.0;
+                }
+              for (unsigned int line = 0; line < lines_per_face; ++line)
+                {
+                  for (unsigned int associated_edge_dof = 0; associated_edge_dof < associated_edge_dofs[line]; ++associated_edge_dof)
+                    {
+                      const unsigned int face_idx = associated_edge_dof_to_face_dof[line][associated_edge_dof];
+                      const unsigned int cell_idx = fe.face_to_cell_index(face_idx, face);
+                      tmp -= dof_values[face_idx]*fe_face_values[vec].value(cell_idx, q_point);
+                    }
+                }
+
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  tmp[d] += values[q_point] (first_vector_component + d);
+                }
+
+              // Tensor of normal vector on the face at q_point;
+              for (unsigned int d = 0; d < dim; ++d)
+                {
+                  normal_vector[d] = normal_vector_list[q_point](d);
+                }
+              // Now compute the linear system:
+              // On a face:
+              // The matrix entries are:
+              // \int_{face} (n x face_shape_function_i) \cdot ( n x face_shape_function_j) dS
+              //
+              // The RHS entries are:
+              // \int_{face} (n x (Edge_residual) \cdot (n x face_shape_function_i) dS
+
+              for (unsigned int j = 0; j < associated_face_dofs; ++j)
+                {
+                  const unsigned int j_face_idx = associated_face_dof_to_face_dof[j];
+                  const unsigned int cell_j = fe.face_to_cell_index (j_face_idx, face);
+
+                  cross_product(cross_product_j,
+                                normal_vector,
+                                fe_face_values[vec].value(cell_j, q_point));
+
+                  for (unsigned int i = 0; i < associated_face_dofs; ++i)
+                    {
+                      const unsigned int i_face_idx = associated_face_dof_to_face_dof[i];
+                      const unsigned int cell_i = fe.face_to_cell_index (i_face_idx, face);
+                      cross_product(cross_product_i,
+                                    normal_vector,
+                                    fe_face_values[vec].value(cell_i, q_point));
+
+                      face_matrix(i,j)
+                      += fe_face_values.JxW (q_point)
+                         *cross_product_i
+                         *cross_product_j;
+
+                    }
+                  // compute rhs
+                  cross_product(cross_product_rhs,
+                                normal_vector,
+                                tmp);
+                  face_rhs(j)
+                  += fe_face_values.JxW (q_point)
+                     *cross_product_rhs
+                     *cross_product_j;
+
+                }
+            }
+
+          // Solve lienar system:
+          face_matrix_inv.invert(face_matrix);
+          face_matrix_inv.vmult(face_solution, face_rhs);
+
+
+          // Store computed DoFs:
+          for (unsigned int associated_face_dof = 0; associated_face_dof < associated_face_dofs; ++associated_face_dof)
+            {
+              dof_values[associated_face_dof_to_face_dof[associated_face_dof]] = face_solution(associated_face_dof);
+              dofs_processed[associated_face_dof_to_face_dof[associated_face_dof]] = true;
+            }
+          break;
+        }
+        default:
+          Assert (false, ExcNotImplemented ());
+        }
+    }
+
+    template <int dim, class DH>
+    void
+    compute_project_boundary_values_curl_conforming_l2 (const DH &dof_handler,
+                                                        const unsigned int first_vector_component,
+                                                        const Function<dim> &boundary_function,
+                                                        const types::boundary_id boundary_component,
+                                                        ConstraintMatrix &constraints,
+                                                        const hp::MappingCollection<dim, dim> &mapping_collection)
+    {
+      // L2-projection based interpolation formed in one (in 2D) or two (in 3D) steps.
+      //
+      // In 2D we only need to constrain edge DoFs.
+      //
+      // In 3D we need to constrain both edge and face DoFs. This is done in two parts.
+      //
+      // For edges, since the face shape functions are zero here ("bubble functions"),
+      // we project the tangential component of the boundary function and compute
+      // the L2-projection. This returns the values for the DoFs associated with each edge shape
+      // function. In 3D, this is computed by internals::compute_edge_projection_l2, in 2D,
+      // it is handled by compute_face_projection_curl_conforming_l2.
+      //
+      // For faces we compute the residual of the boundary function which is satisfied
+      // by the edge shape functions alone. Which can then be used to calculate the
+      // remaining face DoF values via a projection which leads to a linear system to
+      // solve. This is handled by compute_face_projection_curl_conforming_l2
+      //
+      // For details see (for example) section 4.2:
+      // Electromagnetic scattering simulation using an H (curl) conforming hp finite element
+      // method in three dimensions, PD Ledger, K Morgan, O Hassan,
+      // Int. J.  Num. Meth. Fluids, Volume 53, Issue 8, pages 1267–1296, 20 March 2007:
+      // http://onlinelibrary.wiley.com/doi/10.1002/fld.1223/abstract
+
+      // Create hp FEcollection, dof_handler can be either hp or standard type.
+      // From here on we can treat it like a hp-namespace object.
+      const hp::FECollection<dim> fe_collection (dof_handler.get_fe ());
+
+      // Create face quadrature collection
+      hp::QCollection<dim - 1> face_quadrature_collection;
+      for (unsigned int i = 0; i < fe_collection.size (); ++i)
+        {
+          const QGauss<dim - 1>  reference_face_quadrature (2 * fe_collection[i].degree + 1);
+          face_quadrature_collection.push_back(reference_face_quadrature);
+        }
+
+      hp::FEFaceValues<dim> fe_face_values (mapping_collection, fe_collection,
+                                            face_quadrature_collection,
+                                            update_values |
+                                            update_quadrature_points |
+                                            update_normal_vectors |
+                                            update_JxW_values);
+
+      // Storage for dof values found and whether they have been processed:
+      std::vector<bool> dofs_processed;
+      std::vector<double> dof_values;
+      std::vector<types::global_dof_index> face_dof_indices;
+      typename DH::active_cell_iterator cell = dof_handler.begin_active ();
+
+      switch (dim)
+        {
+        case 2:
+        {
+          for (; cell != dof_handler.end (); ++cell)
+            {
+              if (cell->at_boundary () && cell->is_locally_owned ())
+                {
+                  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+                    {
+                      if (cell->face (face)->boundary_indicator () == boundary_component)
+                        {
+                          // If the FE is an FE_Nothing object there is no work to do
+                          if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
+                            {
+                              return;
+                            }
+
+                          // This is only implemented for FE_Nedelec elements.
+                          // If the FE is a FESystem we cannot check this.
+                          if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
+                            {
+                              typedef FiniteElement<dim> FEL;
+                              AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
+                                           typename FEL::ExcInterpolationNotImplemented ());
+
+                            }
+
+                          const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face;
+
+                          dofs_processed.resize (dofs_per_face);
+                          dof_values.resize (dofs_per_face);
+
+                          for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                            {
+                              dof_values[dof] = 0.0;
+                              dofs_processed[dof] = false;
+                            }
+
+                          // Compute the projection of the boundary function on the edge.
+                          // In 2D this is all that's required.
+                          compute_face_projection_curl_conforming_l2 (cell, face, fe_face_values,
+                                                                      boundary_function,
+                                                                      first_vector_component,
+                                                                      dof_values, dofs_processed);
+
+                          // store the local->global map:
+                          face_dof_indices.resize(dofs_per_face);
+                          cell->face (face)->get_dof_indices (face_dof_indices,
+                                                              cell->active_fe_index ());
+
+                          // Add the computed constraints to the constraint matrix,
+                          // assuming the degree of freedom is not already constrained.
+                          for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                            {
+                              if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
+                                  && !(constraints.is_constrained (face_dof_indices[dof])))
+                                {
+                                  constraints.add_line (face_dof_indices[dof]);
+                                  if (std::abs (dof_values[dof]) > 1e-13)
+                                    {
+                                      constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+          break;
+        }
+
+        case 3:
+        {
+          hp::QCollection<dim> edge_quadrature_collection;
+
+          // Create equivalent of FEEdgeValues:
+          for (unsigned int i = 0; i < fe_collection.size (); ++i)
+            {
+              const QGauss<dim-2> reference_edge_quadrature (2 * fe_collection[i].degree + 1);
+              for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+                {
+                  for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_face; ++line)
+                    {
+                      edge_quadrature_collection.push_back
+                      (QProjector<dim>::project_to_face
+                       (QProjector<dim - 1>::project_to_face
+                        (reference_edge_quadrature, line), face));
+                    }
+                }
+            }
+
+          hp::FEValues<dim> fe_edge_values (mapping_collection, fe_collection,
+                                            edge_quadrature_collection,
+                                            update_jacobians |
+                                            update_JxW_values |
+                                            update_quadrature_points |
+                                            update_values);
+
+          for (; cell != dof_handler.end (); ++cell)
+            {
+              if (cell->at_boundary () && cell->is_locally_owned ())
+                {
+                  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
+                    {
+                      if (cell->face (face)->boundary_indicator () == boundary_component)
+                        {
+                          // If the FE is an FE_Nothing object there is no work to do
+                          if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
+                            {
+                              return;
+                            }
+
+                          // This is only implemented for FE_Nedelec elements.
+                          // If the FE is a FESystem we cannot check this.
+                          if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
+                            {
+                              typedef FiniteElement<dim> FEL;
+
+                              AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
+                                           typename FEL::ExcInterpolationNotImplemented ());
+                            }
+
+                          const unsigned int superdegree = cell->get_fe ().degree;
+                          const unsigned int degree = superdegree - 1;
+                          const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face;
+
+                          dofs_processed.resize (dofs_per_face);
+                          dof_values.resize (dofs_per_face);
+                          for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                            {
+                              dof_values[dof] = 0.0;
+                              dofs_processed[dof] = false;
+                            }
+
+                          // First compute the projection on the edges.
+                          for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line)
+                            {
+                              compute_edge_projection_l2 (cell, face, line, fe_edge_values,
+                                                          boundary_function,
+                                                          first_vector_component,
+                                                          dof_values, dofs_processed);
+                            }
+
+                          // If there are higher order shape functions, then we
+                          // still need to compute the face projection
+                          if (degree > 0)
+                            {
+                              compute_face_projection_curl_conforming_l2 (cell, face, fe_face_values,
+                                                                          boundary_function,
+                                                                          first_vector_component,
+                                                                          dof_values,
+                                                                          dofs_processed);
+                            }
+
+                          // Store the computed values in the global vector.
+                          face_dof_indices.resize(dofs_per_face);
+                          cell->face (face)->get_dof_indices (face_dof_indices,
+                                                              cell->active_fe_index ());
+
+                          for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+                            {
+                              if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
+                                  && !(constraints.is_constrained (face_dof_indices[dof])))
+                                {
+                                  constraints.add_line (face_dof_indices[dof]);
+
+                                  if (std::abs (dof_values[dof]) > 1e-13)
+                                    {
+                                      constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+          break;
+        }
+        default:
+          Assert (false, ExcNotImplemented ());
+        }
+    }
+
+  }
+
+
+  template <int dim>
+  void
+  project_boundary_values_curl_conforming_l2 (const DoFHandler<dim> &dof_handler,
+                                              const unsigned int first_vector_component,
+                                              const Function<dim> &boundary_function,
+                                              const types::boundary_id boundary_component,
+                                              ConstraintMatrix &constraints,
+                                              const Mapping<dim> &mapping)
+  {
+    // non-hp version - calls the internal
+    // compute_project_boundary_values_curl_conforming_l2() function
+    // above after recasting the mapping.
+
+    hp::MappingCollection<dim> mapping_collection (mapping);
+    internals::
+    compute_project_boundary_values_curl_conforming_l2(dof_handler,
+                                                       first_vector_component,
+                                                       boundary_function,
+                                                       boundary_component,
+                                                       constraints,
+                                                       mapping_collection);
+  }
+
+  template <int dim>
+  void
+  project_boundary_values_curl_conforming_l2 (const hp::DoFHandler<dim> &dof_handler,
+                                              const unsigned int first_vector_component,
+                                              const Function<dim> &boundary_function,
+                                              const types::boundary_id boundary_component,
+                                              ConstraintMatrix &constraints,
+                                              const hp::MappingCollection<dim, dim> &mapping_collection)
+  {
+    // hp version - calls the internal
+    // compute_project_boundary_values_curl_conforming_l2() function above.
+    internals::
+    compute_project_boundary_values_curl_conforming_l2(dof_handler,
+                                                       first_vector_component,
+                                                       boundary_function,
+                                                       boundary_component,
+                                                       constraints,
+                                                       mapping_collection);
+  }
+
+
+
   namespace internals
   {
     // This function computes the projection of the boundary function on the
index 2b006e4727550fd4fe53f10bc30f132099044633..dec32b7e2a1b4693f4b0b15153c6e33d60708fb4 100644 (file)
@@ -141,6 +141,22 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
          const types::boundary_id,
          ConstraintMatrix&,
          const hp::MappingCollection<deal_II_dimension>&);
+      template
+        void project_boundary_values_curl_conforming_l2<deal_II_dimension>
+        (const DoFHandler<deal_II_dimension>&,
+         const unsigned int,
+         const Function<deal_II_dimension>&,
+         const types::boundary_id,
+         ConstraintMatrix&,
+         const Mapping<deal_II_dimension>&);
+      template
+        void project_boundary_values_curl_conforming_l2<deal_II_dimension>
+        (const hp::DoFHandler<deal_II_dimension>&,
+         const unsigned int,
+         const Function<deal_II_dimension>&,
+         const types::boundary_id,
+         ConstraintMatrix&,
+         const hp::MappingCollection<deal_II_dimension>&);
       template
         void project_boundary_values_div_conforming<deal_II_dimension>
         (const DoFHandler<deal_II_dimension>&,
diff --git a/tests/deal.II/nedelec_non_rect_face.cc b/tests/deal.II/nedelec_non_rect_face.cc
new file mode 100644 (file)
index 0000000..283ad77
--- /dev/null
@@ -0,0 +1,472 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2014 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_Nedelec 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/convergence_table.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+namespace Maxwell
+{
+  using namespace dealii;
+
+  // Dirichlet BCs / exact solution:.
+  template<int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    ExactSolution();
+
+    virtual void vector_value_list (const std::vector<Point<dim> > &points,
+                                    std::vector<Vector<double> >   &values) const;
+
+
+    void curl_value_list(const std::vector<Point<dim> > &points,
+                         std::vector<Vector<double> > &value_list);
+  };
+
+  template<int dim>
+  ExactSolution<dim>::ExactSolution()
+    :
+    Function<dim> (dim)
+  {}
+  template <int dim>
+  void ExactSolution<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                              std::vector<Vector<double> > &value_list) const
+  {
+    Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size()));
+    const unsigned int n_points = points.size();
+
+    for (unsigned int i=0; i<n_points; ++i)
+      {
+        const Point<dim> &p = points[i];
+
+        /* quadratic: */
+        value_list[i](0) = p(0)*p(0);
+        value_list[i](1) = p(1)*p(1);
+        value_list[i](2) = p(2)*p(2);
+
+      }
+  }
+  // Additional functions to create Neumann conditions, zero in this case.
+  template <int dim>
+  void ExactSolution<dim>::curl_value_list(const std::vector<Point<dim> > &points,
+                                           std::vector<Vector<double> > &value_list)
+  {
+    Assert(value_list.size() == points.size(), ExcDimensionMismatch(value_list.size(), points.size()));
+    const unsigned int n_points = points.size();
+
+
+    double exponent;
+    for (unsigned int i=0; i<n_points; ++i)
+      {
+        const Point<dim> &p = points[i];
+        // Real:
+        value_list[i](0) = 0.0;
+        value_list[i](1) = 0.0;
+        value_list[i](2) = 0.0;
+      }
+  }
+  // END Dirichlet BC
+
+  // MAIN MAXWELL CLASS
+  template <int dim>
+  class MaxwellProblem
+  {
+  public:
+    MaxwellProblem (const unsigned int order);
+    ~MaxwellProblem ();
+    void run ();
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void process_solution(const unsigned int cycle);
+
+    double calcErrorHcurlNorm();
+
+    Triangulation<dim>   triangulation;
+    MappingQ<dim>        mapping;
+    DoFHandler<dim>      dof_handler;
+    FE_Nedelec<dim>            fe;
+    ConstraintMatrix     constraints;
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    ConvergenceTable       convergence_table;
+
+    ExactSolution<dim> exact_solution;
+
+    unsigned int p_order;
+    unsigned int quad_order;
+  };
+
+  template <int dim>
+  MaxwellProblem<dim>::MaxwellProblem (const unsigned int order)
+    :
+    mapping (1),
+    dof_handler (triangulation),
+    fe(order),
+    exact_solution()
+  {
+    p_order = order;
+    quad_order = 2*(p_order+1);
+  }
+
+  template <int dim>
+  MaxwellProblem<dim>::~MaxwellProblem ()
+  {
+    dof_handler.clear ();
+  }
+  template<int dim>
+  double MaxwellProblem<dim>::calcErrorHcurlNorm()
+  {
+    QGauss<dim>  quadrature_formula(quad_order);
+    const unsigned int n_q_points = quadrature_formula.size();
+
+    FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+
+    // Extractor
+    const FEValuesExtractors::Vector E_re(0);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    // storage for exact sol:
+    std::vector<Vector<double> > exactsol_list(n_q_points, Vector<double>(fe.n_components()));
+    std::vector<Vector<double> > exactcurlsol_list(n_q_points, Vector<double>(fe.n_components()));
+    Tensor<1,dim>  exactsol;
+    Tensor<1,dim>  exactcurlsol;
+
+    // storage for computed sol:
+    std::vector<Tensor<1,dim> > sol(n_q_points);
+    Tensor<1,dim> curlsol(n_q_points);
+
+    double h_curl_norm=0.0;
+
+    unsigned int block_index_i;
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+
+        // Store exact values of E and curlE:
+        exact_solution.vector_value_list(fe_values.get_quadrature_points(), exactsol_list);
+        exact_solution.curl_value_list(fe_values.get_quadrature_points(), exactcurlsol_list);
+
+
+        // Store computed values at quad points:
+        fe_values[E_re].get_function_values(solution, sol);
+
+        // Calc values of curlE from fe solution:
+        cell->get_dof_indices (local_dof_indices);
+        // Loop over quad points to calculate solution:
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            // Split exact solution into real/imaginary parts:
+            for (unsigned int component=0; component<dim; component++)
+              {
+                exactsol[component] = exactsol_list[q_point][component];
+                exactcurlsol[component] = exactcurlsol_list[q_point][component];
+              }
+            // Loop over DoFs to calculate curl of solution @ quad point
+            curlsol=0.0;
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                // Construct local curl value @ quad point
+                curlsol += solution(local_dof_indices[i])*fe_values[E_re].curl(i,q_point);
+              }
+            // Integrate difference at each point:
+            h_curl_norm += ( (exactsol-sol[q_point])*(exactsol-sol[q_point])
+                             + (exactcurlsol-curlsol)*(exactcurlsol-curlsol)
+                           )*fe_values.JxW(q_point);
+          }
+      }
+    return sqrt(h_curl_norm);
+  }
+
+  template <int dim>
+  void MaxwellProblem<dim>::setup_system ()
+  {
+    dof_handler.distribute_dofs (fe);
+    DoFRenumbering::block_wise (dof_handler);
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+    constraints.clear ();
+
+
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+    // FE_Nedelec boundary condition.
+    VectorTools::project_boundary_values_curl_conforming_l2 (dof_handler, 0, exact_solution, 0, constraints);
+
+
+    constraints.close ();
+
+    CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    c_sparsity,
+                                    constraints,false);
+
+    sparsity_pattern.copy_from(c_sparsity);
+    system_matrix.reinit (sparsity_pattern);
+
+  }
+  template <int dim>
+  void MaxwellProblem<dim>::assemble_system ()
+  {
+    QGauss<dim>  quadrature_formula(quad_order);
+    QGauss<dim-1> face_quadrature_formula(quad_order);
+
+    const unsigned int n_q_points    = quadrature_formula.size();
+    const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+    FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_quadrature_points  |  update_JxW_values);
+
+    FEFaceValues<dim> fe_face_values(mapping, fe, face_quadrature_formula,
+                                     update_values | update_quadrature_points |
+                                     update_normal_vectors | update_JxW_values);
+
+    // Extractor
+    const FEValuesExtractors::Vector E_re(0);
+
+    // Local cell storage:
+    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double> cell_rhs (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    //RHS storage:
+    std::vector<Vector<double> > rhs_value_list(n_q_points, Vector<double>(fe.n_components()));
+    Tensor<1,dim> rhs_value_vector(dim);
+
+    // Neumann storage
+    std::vector<Vector<double> > neumann_value_list(n_face_q_points, Vector<double>(fe.n_components()));
+    std::vector<Point<dim>> normal_vector_list(fe_face_values.get_normal_vectors());
+    Tensor<1,dim> neumann_value_vector(dim);
+    Tensor<1,dim> neumann_value(dim);
+    Tensor<1,dim> normal_vector;
+
+    // loop over all cells:
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit(cell);
+        cell_matrix = 0;
+        cell_rhs = 0;
+
+        // Calc RHS values:
+        exact_solution.vector_value_list(fe_values.get_quadrature_points(), rhs_value_list);
+
+        // Loop over all element quad points:
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          {
+            // store rhs value at this q point & turn into tensor
+            for (unsigned int component=0; component<dim; component++)
+              {
+                rhs_value_vector[component] = rhs_value_list[q_point](component);
+              }
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                // Construct local matrix:
+                for (unsigned int j=0; j<dofs_per_cell; ++j)
+                  {
+                    cell_matrix(i,j) += ( (fe_values[E_re].curl(i,q_point)*fe_values[E_re].curl(j,q_point))
+                                          + fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)
+                                        )*fe_values.JxW(q_point);
+
+                  }
+                // construct local RHS:
+                cell_rhs(i) += rhs_value_vector*fe_values[E_re].value(i,q_point)*fe_values.JxW(q_point);
+
+              }
+          }
+
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global(cell_matrix,
+                                               cell_rhs,
+                                               local_dof_indices,
+                                               system_matrix, system_rhs);
+      }
+  }
+  template <int dim>
+  void MaxwellProblem<dim>::solve ()
+  {
+    /* Direct */
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+
+    A_direct.vmult (solution, system_rhs);
+    constraints.distribute (solution);
+
+  }
+  template<int dim>
+  void MaxwellProblem<dim>::process_solution(const unsigned int cycle)
+  {
+
+    Vector<double> diff_per_cell(triangulation.n_active_cells());
+
+
+    VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
+                                      diff_per_cell, QGauss<dim>(quad_order+2),
+                                      VectorTools::L2_norm);
+
+
+    const double L2_error = diff_per_cell.l2_norm();
+
+    double Hcurl_error = calcErrorHcurlNorm();
+
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", triangulation.n_active_cells());
+    convergence_table.add_value("dofs", dof_handler.n_dofs());
+    convergence_table.add_value("L2 Error", L2_error);
+    convergence_table.add_value("H(curl) Error", Hcurl_error);
+
+  }
+
+  template <int dim>
+  void MaxwellProblem<dim>::run ()
+  {
+
+    unsigned int numcycles=3;
+
+    for (unsigned int cycle=0; cycle<numcycles; ++cycle)
+      {
+        if (cycle == 0)
+          {
+            /* Cube mesh */
+            GridGenerator::hyper_cube (triangulation, 0, 1);
+            GridTools::distort_random (0.2, triangulation, false);
+
+
+            /* Testing hanging nodes
+             *        triangulation.refine_global(1);
+             *        cell = triangulation.begin_active();
+             *        cell->set_refine_flag(); // create single refined cell
+             *        triangulation.execute_coarsening_and_refinement();
+             */
+
+          }
+        else
+          {
+            triangulation.refine_global (1);
+          }
+
+        setup_system ();
+        assemble_system ();
+        solve ();
+        process_solution (cycle);
+      }
+
+    deallog << "p = " << p_order << std::endl;
+    
+    //convergence_table.set_precision("L2 Error",4);
+    convergence_table.set_scientific("L2 Error",true);
+      
+    //convergence_table.set_precision("H(curl) Error",4);
+    convergence_table.set_scientific("H(curl) Error",true);
+
+    convergence_table.write_text(deallog.get_file_stream());
+
+  }
+}
+// END MAXWELL CLASS
+int main ()
+{
+  using namespace Maxwell;
+  
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+  for (unsigned int p=0; p<4; ++p)
+    {
+      MaxwellProblem<3> (p).run();
+    }
+
+  return 0;
+}
diff --git a/tests/deal.II/nedelec_non_rect_face.output b/tests/deal.II/nedelec_non_rect_face.output
new file mode 100644 (file)
index 0000000..aee54e3
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::p = 0
+cycle cells dofs  L2 Error  H(curl) Error 
+0     1     12   6.2187e-01 6.0931e-01    
+1     8     54   3.3123e-01 3.2983e-01    
+2     64    300  1.6777e-01 1.6760e-01    
+DEAL::p = 1
+cycle cells dofs  L2 Error  H(curl) Error 
+0     1     54   1.5549e-01 1.6230e-01    
+1     8     300  3.8832e-02 4.1960e-02    
+2     64    1944 9.7185e-03 1.1143e-02    
+DEAL::p = 2
+cycle cells dofs  L2 Error  H(curl) Error 
+0     1     144  6.0548e-16 1.2581e-15    
+1     8     882  4.7208e-15 5.6992e-15    
+2     64    6084 3.7516e-14 4.1547e-14    
+DEAL::p = 3
+cycle cells dofs   L2 Error  H(curl) Error 
+0     1     300   9.4037e-16 1.8758e-15    
+1     8     1944  1.5731e-14 1.6809e-14    
+2     64    13872 8.7253e-14 9.1275e-14    

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.