From: buerg <buerg@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 10 Jan 2014 22:35:44 +0000 (+0000)
Subject: Compute hessians for FE_Nedelec correctly.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=104ecad246dbb84f87b2db616493d84660ea49b6;p=dealii-svn.git

Compute hessians for FE_Nedelec correctly.

git-svn-id: https://svn.dealii.org/trunk@32191 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h
index 40b4c742c5..55898960f7 100644
--- a/deal.II/include/deal.II/fe/fe.h
+++ b/deal.II/include/deal.II/fe/fe.h
@@ -931,65 +931,6 @@ public:
   std::pair<unsigned int, unsigned int>
   face_system_to_component_index (const unsigned int index) const;
 
-  /**
-   * Given an index in the natural ordering of indices on a face, return the
-   * index of the same degree of freedom on the cell.
-   *
-   * To explain the concept, consider the case where we would like to know
-   * whether a degree of freedom on a face, for example as part of an FESystem
-   * element, is primitive. Unfortunately, the
-   * is_primitive() function in the FiniteElement class takes a cell index, so
-   * we would need to find the cell index of the shape function that
-   * corresponds to the present face index. This function does that.
-   *
-   * Code implementing this would then look like this:
-   * @code
-   * for (i=0; i<dofs_per_face; ++i)
-   *  if (fe.is_primitive(fe.face_to_equivalent_cell_index(i, some_face_no)))
-   *   ... do whatever
-   * @endcode
-   * The function takes additional arguments that account for the fact that
-   * actual faces can be in their standard ordering with respect to the cell
-   * under consideration, or can be flipped, oriented, etc.
-   *
-   * @param face_dof_index The index of the degree of freedom on a face.
-   *   This index must be between zero and dofs_per_face.
-   * @param face The number of the face this degree of freedom lives on.
-   *   This number must be between zero and GeometryInfo::faces_per_cell.
-   * @param face_orientation One part of the description of the orientation
-   *   of the face. See @ref GlossFaceOrientation .
-   * @param face_flip One part of the description of the orientation
-   *   of the face. See @ref GlossFaceOrientation .
-   * @param face_rotation One part of the description of the orientation
-   *   of the face. See @ref GlossFaceOrientation .
-   * @return The index of this degree of freedom within the set
-   *   of degrees of freedom on the entire cell. The returned value
-   *   will be between zero and dofs_per_cell.
-   *
-   * @note This function exists in this class because that is where it
-   * was first implemented. However, it can't really work in the most
-   * general case without knowing what element we have. The reason is that
-   * when a face is flipped or rotated, we also need to know whether we
-   * need to swap the degrees of freedom on this face, or whether they
-   * are immune from this. For this, consider the situation of a $Q_3$
-   * element in 2d. If face_flip is true, then we need to consider
-   * the two degrees of freedom on the edge in reverse order. On the other
-   * hand, if the element were a $Q_1^2$, then because the two degrees of
-   * freedom on this edge belong to different vector components, they
-   * should not be considered in reverse order. What all of this shows is
-   * that the function can't work if there are more than one degree of
-   * freedom per line or quad, and that in these cases the function will
-   * throw an exception pointing out that this functionality will need
-   * to be provided by a derived class that knows what degrees of freedom
-   * actually represent.
-   */
-  virtual
-  unsigned int face_to_cell_index (const unsigned int face_dof_index,
-                                   const unsigned int face,
-                                   const bool face_orientation = true,
-                                   const bool face_flip        = false,
-                                   const bool face_rotation    = false) const;
-
   /**
    * For faces with non-standard face_orientation in 3D, the dofs on faces
    * (quads) have to be permuted in order to be combined with the correct
@@ -1749,7 +1690,8 @@ protected:
                     const unsigned int                       offset,
                     typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
                     InternalDataBase                        &fe_internal,
-                    FEValuesData<dim,spacedim>                       &data) const;
+                    FEValuesData<dim,spacedim>                       &data,
+                    const MappingType                        mapping_type = mapping_covariant) const;
 
   /**
    * Given the pattern of nonzero components for each shape function, compute
diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc
index df61dbc73a..a9f9e0aeca 100644
--- a/deal.II/source/fe/fe.cc
+++ b/deal.II/source/fe/fe.cc
@@ -544,96 +544,6 @@ block_mask (const ComponentMask &component_mask) const
 
 
 
-template <int dim, int spacedim>
-unsigned int
-FiniteElement<dim,spacedim>::
-face_to_cell_index (const unsigned int face_index,
-                    const unsigned int face,
-                    const bool face_orientation,
-                    const bool face_flip,
-                    const bool face_rotation) const
-{
-  Assert (face_index < this->dofs_per_face,
-          ExcIndexRange(face_index, 0, this->dofs_per_face));
-  Assert (face < GeometryInfo<dim>::faces_per_cell,
-          ExcIndexRange(face, 0, GeometryInfo<dim>::faces_per_cell));
-
-//TODO: we could presumably solve the 3d case below using the
-// adjust_quad_dof_index_for_face_orientation_table field. for the
-// 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
-// since that array is empty (presumably because we thought that
-// there are no flipped edges in 2d, but these can happen in
-// DoFTools::make_periodicity_constraints, for example). so we
-// would need to either fill this field, or rely on derived classes
-// implementing this function, as we currently do
-
-  // see the function's documentation for an explanation of this
-  // assertion -- in essence, derived classes have to implement
-  // an overloaded version of this function if we are to use any
-  // other than standard orientation
-  if ((face_orientation != true) || (face_flip != false) || (face_rotation != false))
-    Assert ((this->dofs_per_line <= 1) && (this->dofs_per_quad <= 1),
-            ExcMessage ("The function in this base class can not handle this case. "
-                        "Rather, the derived class you are using must provide "
-                        "an overloaded version but apparently hasn't done so. See "
-                        "the documentation of this function for more information."));
-
-  // we need to distinguish between DoFs on vertices, lines and in 3d quads.
-  // do so in a sequence of if-else statements
-  if (face_index < this->first_face_line_index)
-    // DoF is on a vertex
-    {
-      // get the number of the vertex on the face that corresponds to this DoF,
-      // along with the number of the DoF on this vertex
-      const unsigned int face_vertex         = face_index / this->dofs_per_vertex;
-      const unsigned int dof_index_on_vertex = face_index % this->dofs_per_vertex;
-
-      // then get the number of this vertex on the cell and translate
-      // this to a DoF number on the cell
-      return (GeometryInfo<dim>::face_to_cell_vertices(face, face_vertex,
-                                                       face_orientation,
-                                                       face_flip,
-                                                       face_rotation)
-              * this->dofs_per_vertex
-              +
-              dof_index_on_vertex);
-    }
-  else if (face_index < this->first_face_quad_index)
-    // DoF is on a face
-    {
-      // do the same kind of translation as before. we need to only consider
-      // DoFs on the lines, i.e., ignoring those on the vertices
-      const unsigned int index = face_index - this->first_face_line_index;
-
-      const unsigned int face_line         = index / this->dofs_per_line;
-      const unsigned int dof_index_on_line = index % this->dofs_per_line;
-
-      return (this->first_line_index
-              + GeometryInfo<dim>::face_to_cell_lines(face, face_line,
-                                                      face_orientation,
-                                                      face_flip,
-                                                      face_rotation)
-              * this->dofs_per_line
-              +
-              dof_index_on_line);
-    }
-  else
-    // DoF is on a quad
-    {
-      Assert (dim >= 3, ExcInternalError());
-
-      // ignore vertex and line dofs
-      const unsigned int index = face_index - this->first_face_quad_index;
-
-      return (this->first_quad_index
-              + face * this->dofs_per_quad
-              + index);
-    }
-}
-
-
-
-
 template <int dim, int spacedim>
 unsigned int
 FiniteElement<dim,spacedim>::adjust_quad_dof_index_for_face_orientation (const unsigned int index,
@@ -1195,7 +1105,8 @@ FiniteElement<1,2>::compute_2nd (
   const unsigned int,
   Mapping<1,2>::InternalDataBase &,
   InternalDataBase &,
-  FEValuesData<1,2> &) const
+  FEValuesData<1,2> &,
+  MappingType) const
 {
 
   Assert(false, ExcNotImplemented());
@@ -1210,7 +1121,8 @@ FiniteElement<1,3>::compute_2nd (
   const unsigned int,
   Mapping<1,3>::InternalDataBase &,
   InternalDataBase &,
-  FEValuesData<1,3> &) const
+  FEValuesData<1,3> &,
+  MappingType) const
 {
 
   Assert(false, ExcNotImplemented());
@@ -1226,7 +1138,8 @@ FiniteElement<2,3>::compute_2nd (
   const unsigned int,
   Mapping<2,3>::InternalDataBase &,
   InternalDataBase &,
-  FEValuesData<2,3> &) const
+  FEValuesData<2,3> &,
+  MappingType) const
 {
 
   Assert(false, ExcNotImplemented());
@@ -1242,7 +1155,8 @@ FiniteElement<dim,spacedim>::compute_2nd (
   const unsigned int offset,
   typename Mapping<dim,spacedim>::InternalDataBase &mapping_internal,
   InternalDataBase                     &fe_internal,
-  FEValuesData<dim,spacedim>                    &data) const
+  FEValuesData<dim,spacedim>                    &data,
+  MappingType mapping_type) const
 {
   Assert ((fe_internal.update_each | fe_internal.update_once)
           & update_hessians,
@@ -1368,7 +1282,7 @@ FiniteElement<dim,spacedim>::compute_2nd (
         // cell
         for (unsigned int d=0; d<spacedim; ++d)
           {
-            mapping.transform (diff_quot[d], diff_quot2, mapping_internal, mapping_covariant);
+            mapping.transform (diff_quot[d], diff_quot2, mapping_internal, mapping_type);
 
             for (unsigned int q=0; q<n_q_points; ++q)
               for (unsigned int d1=0; d1<spacedim; ++d1)
diff --git a/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/source/fe/fe_poly_tensor.cc
index 9f6460c363..84a20f2810 100644
--- a/deal.II/source/fe/fe_poly_tensor.cc
+++ b/deal.II/source/fe/fe_poly_tensor.cc
@@ -585,9 +585,27 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_values (
     }
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
-    this->compute_2nd (mapping, cell,
-                       typename QProjector<dim>::DataSetDescriptor().cell(),
-                       mapping_data, fe_data, data);
+    if (mapping_type == mapping_nedelec)
+    {
+      this->compute_2nd (mapping, cell,
+                         typename QProjector<dim>::DataSetDescriptor().cell(),
+                         mapping_data, fe_data, data, mapping_covariant_gradient);
+      
+      for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+      {
+        const unsigned int first = data.shape_function_to_row_table[i * this->n_components() +
+                                                                    this->get_nonzero_components(i).first_selected_component()];
+        
+        for (unsigned int k = 0; k < n_q_points; ++k)
+          for (unsigned int d = 0; d < dim; ++d)
+            data.shape_hessians[first + d][k] *= sign_change[i];
+      }
+    }
+    
+    else
+      this->compute_2nd (mapping, cell,
+                         typename QProjector<dim>::DataSetDescriptor().cell(),
+                         mapping_data, fe_data, data);
 }
 
 
@@ -791,7 +809,23 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_face_values (
     }
 
   if (flags & update_hessians)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+    if (mapping_type == mapping_nedelec)
+    {
+      this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data, mapping_covariant_gradient);
+      
+      for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+      {
+        const unsigned int first = data.shape_function_to_row_table[i * this->n_components() +
+                                                                    this->get_nonzero_components(i).first_selected_component()];
+        
+        for (unsigned int k = 0; k < n_q_points; ++k)
+          for (unsigned int d = 0; d < dim; ++d)
+            data.shape_hessians[first + d][k] *= sign_change[i];
+      }
+    }
+    
+    else
+      this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
 }
 
 
@@ -989,7 +1023,7 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
 
               for (unsigned int k = 0; k < n_q_points; ++k)
                 for (unsigned int d = 0; d < dim; ++d)
-                  data.shape_gradients[first + d][k] = shape_grads1[k][d];
+                  data.shape_gradients[first + d][k] = sign_change[i] * shape_grads1[k][d];
 
               break;
             }
@@ -1001,7 +1035,23 @@ FE_PolyTensor<POLY,dim,spacedim>::fill_fe_subface_values (
     }
 
   if (flags & update_hessians)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+    if (mapping_type == mapping_nedelec)
+    {
+      this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data, mapping_covariant_gradient);
+      
+      for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+      {
+        const unsigned int first = data.shape_function_to_row_table[i * this->n_components() +
+                                                                    this->get_nonzero_components(i).first_selected_component()];
+        
+        for (unsigned int k = 0; k < n_q_points; ++k)
+          for (unsigned int d = 0; d < dim; ++d)
+            data.shape_hessians[first + d][k] *= sign_change[i];
+      }
+    }
+    
+    else
+      this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
 }