From 39b01ca86810d6925d5841f728c38a1537a07669 Mon Sep 17 00:00:00 2001
From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 29 Apr 2006 13:24:03 +0000
Subject: [PATCH] bug fix in 2nd derivatives on faces

git-svn-id: https://svn.dealii.org/branches/Branch-5-2@12939 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/deal.II/source/fe/fe.cc | 28 +++++++++++++++-------------
 1 file changed, 15 insertions(+), 13 deletions(-)

diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc
index 81eaf2ff58..7bb316829b 100644
--- a/deal.II/deal.II/source/fe/fe.cc
+++ b/deal.II/deal.II/source/fe/fe.cc
@@ -660,15 +660,14 @@ FiniteElement<dim>::memory_consumption () const
 
 template <int dim>
 void
-FiniteElement<dim>::
-compute_2nd (const Mapping<dim>                   &mapping,
-	     const typename Triangulation<dim>::cell_iterator &cell,
-	     const unsigned int,
-	     typename Mapping<dim>::InternalDataBase &mapping_internal,
-	     InternalDataBase                     &fe_internal,
-	     FEValuesData<dim>                    &data) const
-{
-//TODO[GK]: This function presently has a flaw: it ignores the offset parameter. This means that if it is called for faces or subfaces, it computes the data for _all_ faces or subfaces, even though we only need to have this for one of them. given that computing second derivatives is expensive, this is significant. What should be done is to honor the offset parameter and of course only compute what we actually need
+FiniteElement<dim>::compute_2nd (
+  const Mapping<dim>                   &mapping,
+  const typename Triangulation<dim>::cell_iterator &cell,
+  const unsigned int offset,
+  typename Mapping<dim>::InternalDataBase &mapping_internal,
+  InternalDataBase                     &fe_internal,
+  FEValuesData<dim>                    &data) const
+{
   Assert ((fe_internal.update_each | fe_internal.update_once)
 	  & update_second_derivatives,
 	  ExcInternalError());
@@ -689,6 +688,9 @@ compute_2nd (const Mapping<dim>                   &mapping,
     {
       fe_internal.differences[d]->reinit(cell);
       fe_internal.differences[d+dim]->reinit(cell);
+      Assert(offset <= fe_internal.differences[d]->n_quadrature_points - n_q_points,
+	     ExcIndexRange(offset, 0, fe_internal.differences[d]->n_quadrature_points
+			   - n_q_points));
     }
 
 				   // collection of difference
@@ -728,8 +730,8 @@ compute_2nd (const Mapping<dim>                   &mapping,
               Tensor<1,dim> right, left;
               if (is_primitive(shape_index))
                 {
-                  right = fe_internal.differences[d1]->shape_grad(shape_index, q);
-                  left  = fe_internal.differences[d1+dim]->shape_grad(shape_index, q);
+                  right = fe_internal.differences[d1]->shape_grad(shape_index, q+offset);
+                  left  = fe_internal.differences[d1+dim]->shape_grad(shape_index, q+offset);
                 }
               else
                 {
@@ -765,9 +767,9 @@ compute_2nd (const Mapping<dim>                   &mapping,
                           ExcInternalError());
 
                   right = fe_internal.differences[d1]
-                          ->shape_grad_component(shape_index, q, component);
+                          ->shape_grad_component(shape_index, q+offset, component);
                   left  = fe_internal.differences[d1+dim]
-                          ->shape_grad_component(shape_index, q, component);
+                          ->shape_grad_component(shape_index, q+offset, component);
                 };
               
                                                // compute the second
-- 
2.39.5