Fixed integrate difference, and updated step-38 to reflect the new functionality.

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

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 0ea8e105ad8a191177546960fcf67f8030ebfb07..c08db3a93b7648e6205d632d00a354a605523482 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -94,6 +94,23 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
 <h3>General</h3>
 
 <ol>
+  <li><p>New: When computing errors using
+  VectorTools::integrate_difference in codimension one problems, if
+  you specified a norm that requires the computation of the gradients,
+  you'd get in trouble, because on codimension one manifolds we only
+  have information on the <em>tangential</em> gradient. This is now
+  fixed, by discarding the normal component of the provided function,
+  before computing the difference with the finite element function.
+  <br>
+  (Luca Heltai, 2010/12/05)
+  </p></li>
+  
+  <li><p>New: We now have a new tutorial program (step-38), written by
+  Andrea Bonito, to compute the Laplace Beltrami operator on
+  a half sphere.
+  <br> (Luca Heltai, 2010/12/03)
+  </p></li>
+
   <li><p>New: Added a new constructor for PETScWrappers::Vector that takes
   an existing PETSc Vec that was created by the user and is only wrapped for
   usage.

diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc
index 33516965c61c155ccb36e30390d425b3b5f1c583..f514c2bd0ba6a661c9f82d254c277b9dc0316c98 100755 (executable)
--- a/deal.II/examples/step-38/step-38.cc
+++ b/deal.II/examples/step-38/step-38.cc
@@ -80,17 +80,7 @@ Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
   return_value[1] = -dPi *sin(dPi * p(0))*sin(dPi * p(1))*exp(p(2));
   return_value[2] = sin(dPi * p(0))*cos(dPi * p(1))*exp(p(2));
   
-  // tangential gradient: nabla u - (nabla u nu)nu
-  Point<dim> normal;
-  double dLength;  
-
-  dLength = sqrt(p(0)*p(0)+p(1)*p(1)+p(2)*p(2));
-  
-  normal[0] = p(0)/dLength;
-  normal[1] = p(1)/dLength;
-  normal[2] = p(2)/dLength;
-  
-  return return_value - (return_value*normal)*normal;
+  return return_value;
 }
 
 template <int dim>

diff --git a/deal.II/include/deal.II/numerics/vectors.h b/deal.II/include/deal.II/numerics/vectors.h
index aa5103da74227137cb7064e568edb7808ea1e39b..fb3e9e4b8448d14b18afa4dedf1c3bbc4b9b5e6e 100644 (file)
--- a/deal.II/include/deal.II/numerics/vectors.h
+++ b/deal.II/include/deal.II/numerics/vectors.h
@@ -304,7 +304,15 @@ class ConstraintMatrix;
  *   using the Vector::linfty_norm() function.
  *
  *   For the global <i>H</i><sup>1</sup> norm and seminorm, the same rule applies as for the
- *   <i>L</i><sup>2</sup> norm: compute the <i>l</i><sub>2</sub> norm of the cell error vector.
+ *   <i>L</i><sup>2</sup> norm: compute the <i>l</i><sub>2</sub> norm
+ *   of the cell error vector.
+ *
+ *   Note that, in the codimension one case, if you ask for a norm
+ *   that requires the computation of a gradient, then the provided
+ *   function is automatically projected along the curve, and the
+ *   difference is only computed on the tangential part of the
+ *   gradient, since no information is available, on the finite
+ *   dimensional one, on the normal component of the gradient. 
  * </ul>
  *
  * All functions use the finite element given to the DoFHandler object the last

diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h
index ec4ac460eb73fdd61ed5ba7e3061ba9454f9d476..f79ead7a06e288afd5bad213f112141fd8ea5a8d 100644 (file)
--- a/deal.II/include/deal.II/numerics/vectors.templates.h
+++ b/deal.II/include/deal.II/numerics/vectors.templates.h
@@ -4191,11 +4191,14 @@ namespace internal
           case dealii::VectorTools::W1p_seminorm:
           case dealii::VectorTools::W1infty_seminorm:
                 update_flags |= UpdateFlags (update_gradients);
+               if(spacedim == dim+1) update_flags |= UpdateFlags (update_normal_vectors);
+               
                 break;
           case dealii::VectorTools::H1_norm:
           case dealii::VectorTools::W1p_norm:
           case dealii::VectorTools::W1infty_norm:
                 update_flags |= UpdateFlags (update_gradients);
+               if(spacedim == dim+1) update_flags |= UpdateFlags (update_normal_vectors);
                                                  // no break!
           default:
                 update_flags |= UpdateFlags (update_values);
@@ -4332,13 +4335,35 @@ namespace internal
                 }
 
                                                // then subtract finite element
-                                               // function_grads
+                                               // function_grads. We
+                                               // need to be careful
+                                               // in the codimension
+                                               // one case, since
+                                               // there we only have
+                                               // tangential gradients
+                                               // in the finite
+                                               // element function,
+                                               // not the full
+                                               // gradient. This is
+                                               // taken care of, by
+                                               // subtracting the
+                                               // normal component of
+                                               // the gradient from
+                                               // the exact function.  
               fe_values.get_function_grads (fe_function, function_grads);
-              for (unsigned int k=0; k<n_components; ++k)
-                for (unsigned int q=0; q<n_q_points; ++q)
-                  psi_grads[q][k] -= function_grads[q][k];
+             if(update_flags & update_normal_vectors)
+               for (unsigned int k=0; k<n_components; ++k)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   psi_grads[q][k] -= (function_grads[q][k] +
+                                       (psi_grads[q][k]* // (f.n) n
+                                        fe_values.normal_vector(q))*
+                                       fe_values.normal_vector(q));
+             else
+               for (unsigned int k=0; k<n_components; ++k)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   psi_grads[q][k] -= function_grads[q][k];
             }
-
+         
           switch (norm)
             {
               case dealii::VectorTools::mean: