From 3eec369159c80b0a81f379dd75fc83ff0d2f526e Mon Sep 17 00:00:00 2001
From: Timo Heister <timo.heister@gmail.com>
Date: Sat, 16 Jan 2016 17:46:28 -0600
Subject: [PATCH] WIP kelly fix

---
 .../numerics/error_estimator.templates.h      | 37 +++++++++++++------
 1 file changed, 26 insertions(+), 11 deletions(-)

diff --git a/include/deal.II/numerics/error_estimator.templates.h b/include/deal.II/numerics/error_estimator.templates.h
index facca6b8fe..483cd2c5f0 100644
--- a/include/deal.II/numerics/error_estimator.templates.h
+++ b/include/deal.II/numerics/error_estimator.templates.h
@@ -234,11 +234,13 @@ namespace internal
       fe_face_values_neighbor (mapping,
                                finite_element,
                                face_quadratures,
-                               update_gradients),
+                               update_gradients|
+                               update_normal_vectors),
       fe_subface_values (mapping,
                          finite_element,
                          face_quadratures,
-                         update_gradients),
+                         update_gradients|
+                         update_normal_vectors),
       phi (n_solution_vectors,
            std::vector<std::vector<number> >
            (face_quadratures.max_n_quadrature_points(),
@@ -340,7 +342,8 @@ namespace internal
     integrate_over_face
     (ParallelData<DoFHandlerType,number>                 &parallel_data,
      const typename DoFHandlerType::face_iterator        &face,
-     dealii::hp::FEFaceValues<DoFHandlerType::dimension, DoFHandlerType::space_dimension> &fe_face_values_cell)
+     dealii::hp::FEFaceValues<DoFHandlerType::dimension, DoFHandlerType::space_dimension> &fe_face_values_cell,
+     const std::vector<Tensor<1,DoFHandlerType::space_dimension> > & other_normals)
     {
       const unsigned int n_q_points         = parallel_data.psi[0].size(),
                          n_components       = parallel_data.finite_element.n_components(),
@@ -348,11 +351,11 @@ namespace internal
 
       if (face->at_boundary() == false)
         {
-          // compute the jump in the gradients
-          for (unsigned int n=0; n<n_solution_vectors; ++n)
-            for (unsigned int component=0; component<n_components; ++component)
-              for (unsigned int p=0; p<n_q_points; ++p)
-                parallel_data.psi[n][p][component] -= parallel_data.neighbor_psi[n][p][component];
+//          // compute the jump in the gradients
+//          for (unsigned int n=0; n<n_solution_vectors; ++n)
+//            for (unsigned int component=0; component<n_components; ++component)
+//              for (unsigned int p=0; p<n_q_points; ++p)
+//                parallel_data.psi[n][p][component] -= parallel_data.neighbor_psi[n][p][component];
         }
 
       // now psi contains the following:
@@ -375,7 +378,16 @@ namespace internal
             parallel_data.phi[n][point][component]
               = (parallel_data.psi[n][point][component] *
                  parallel_data.normal_vectors[point]);
-
+      if (face->at_boundary() == false)
+        {
+          // compute the jump in the gradients
+          for (unsigned int n=0; n<n_solution_vectors; ++n)
+            for (unsigned int component=0; component<n_components; ++component)
+              for (unsigned int p=0; p<n_q_points; ++p)
+                parallel_data.phi[n][p][component]
+                  -= (parallel_data.neighbor_psi[n][p][component] *
+                      other_normals[p]);
+        }
       // if a coefficient was given: use that to scale the jump in the
       // gradient
       if (parallel_data.coefficients != 0)
@@ -669,7 +681,9 @@ namespace internal
 
       // now go to the generic function that does all the other things
       local_face_integrals[face] =
-        integrate_over_face (parallel_data, face, fe_face_values_cell);
+        integrate_over_face (parallel_data, face,
+                             fe_face_values_cell,
+                             fe_face_values_neighbor.get_present_fe_values().get_all_normal_vectors());
       for (unsigned int i = 0; i < local_face_integrals[face].size(); i++)
         local_face_integrals[face][i] *= factor;
     }
@@ -762,7 +776,8 @@ namespace internal
 
           // call generic evaluate function
           local_face_integrals[neighbor_child->face(neighbor_neighbor)] =
-            integrate_over_face (parallel_data, face, fe_face_values);
+            integrate_over_face (parallel_data, face, fe_face_values,
+                                 fe_subface_values.get_present_fe_values().get_all_normal_vectors());
           for (unsigned int i = 0; i < local_face_integrals[neighbor_child->face(neighbor_neighbor)].size(); i++)
             local_face_integrals[neighbor_child->face(neighbor_neighbor)][i] *= factor;
         }
-- 
2.39.5