rewrite normal vector computation

diff --git a/include/deal.II/numerics/error_estimator.templates.h b/include/deal.II/numerics/error_estimator.templates.h
index 483cd2c5f0ab98d68256dc05e0a0ceb4c4803588..82771dec1ffbf00ff7b9cb05139e24a9c11ad1e0 100644 (file)
--- a/include/deal.II/numerics/error_estimator.templates.h
+++ b/include/deal.II/numerics/error_estimator.templates.h
@@ -151,6 +151,11 @@ namespace internal
        */
       std::vector<Tensor<1,spacedim> > normal_vectors;
 
+      /**
+      * Normal vectors of the opposing face.
+      */
+      std::vector<Tensor<1,spacedim> > neighbor_normal_vectors;
+
       /**
        * Two arrays needed for the values of coefficients in the jumps, if
        * they are given.
@@ -254,6 +259,7 @@ namespace internal
                     (face_quadratures.max_n_quadrature_points(),
                      std::vector<Tensor<1,spacedim,number> > (fe.n_components()))),
       normal_vectors (face_quadratures.max_n_quadrature_points()),
+      neighbor_normal_vectors (face_quadratures.max_n_quadrature_points()),
       coefficient_values1 (face_quadratures.max_n_quadrature_points()),
       coefficient_values (face_quadratures.max_n_quadrature_points(),
                           dealii::Vector<double> (fe.n_components())),
@@ -275,6 +281,7 @@ namespace internal
       const unsigned int n_components = finite_element.n_components();
 
       normal_vectors.resize(n_q_points);
+      neighbor_normal_vectors.resize(n_q_points);
       coefficient_values1.resize(n_q_points);
       coefficient_values.resize(n_q_points);
       JxW_values.resize(n_q_points);
@@ -342,8 +349,7 @@ 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,
-     const std::vector<Tensor<1,DoFHandlerType::space_dimension> > & other_normals)
+     dealii::hp::FEFaceValues<DoFHandlerType::dimension, DoFHandlerType::space_dimension> &fe_face_values_cell)
     {
       const unsigned int n_q_points         = parallel_data.psi[0].size(),
                          n_components       = parallel_data.finite_element.n_components(),
@@ -378,16 +384,19 @@ 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]);
+                += (parallel_data.neighbor_psi[n][p][component] *
+                    parallel_data.neighbor_normal_vectors[p]);
         }
+
       // if a coefficient was given: use that to scale the jump in the
       // gradient
       if (parallel_data.coefficients != 0)
@@ -459,6 +468,8 @@ namespace internal
         }
 
 
+
+
       // now phi contains the following:
       // - for an internal face, phi=[a du/dn]
       // - for a neumann boundary face, phi=a du/dn-g
@@ -671,6 +682,10 @@ namespace internal
               .get_function_gradients (*solutions[n],
                                        parallel_data.neighbor_psi[n]);
             }
+
+          parallel_data.neighbor_normal_vectors =
+            fe_face_values_neighbor.get_present_fe_values().get_all_normal_vectors();
+
         }
       else
         {
@@ -682,8 +697,8 @@ 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,
-                             fe_face_values_neighbor.get_present_fe_values().get_all_normal_vectors());
+                             fe_face_values_cell);
+
       for (unsigned int i = 0; i < local_face_integrals[face].size(); i++)
         local_face_integrals[face][i] *= factor;
     }
@@ -775,9 +790,11 @@ namespace internal
             .get_function_gradients (*solutions[n], parallel_data.neighbor_psi[n]);
 
           // call generic evaluate function
+          parallel_data.neighbor_normal_vectors =
+            fe_subface_values.get_present_fe_values().get_all_normal_vectors();
+
           local_face_integrals[neighbor_child->face(neighbor_neighbor)] =
-            integrate_over_face (parallel_data, face, fe_face_values,
-                                 fe_subface_values.get_present_fe_values().get_all_normal_vectors());
+            integrate_over_face (parallel_data, face, fe_face_values);
           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;
         }