From 4af251ebd802148f25231de93e9ec063d7195b1f Mon Sep 17 00:00:00 2001 From: heltai Date: Sun, 5 Dec 2010 15:12:56 +0000 Subject: [PATCH] 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 --- deal.II/doc/news/changes.h | 17 +++++++++ deal.II/examples/step-38/step-38.cc | 12 +------ deal.II/include/deal.II/numerics/vectors.h | 10 +++++- .../deal.II/numerics/vectors.templates.h | 35 ++++++++++++++++--- 4 files changed, 57 insertions(+), 17 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0ea8e105ad..c08db3a93b 100644 --- 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.

General

    +

  1. 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 tangential gradient. This is now + fixed, by discarding the normal component of the provided function, + before computing the difference with the finite element function. +
    + (Luca Heltai, 2010/12/05) +

    2. + +

  2. New: We now have a new tutorial program (step-38), written by + Andrea Bonito, to compute the Laplace Beltrami operator on + a half sphere. +
    (Luca Heltai, 2010/12/03) +

    3. +

  3. 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 33516965c6..f514c2bd0b 100755 --- 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::gradient (const Point &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 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 diff --git a/deal.II/include/deal.II/numerics/vectors.h b/deal.II/include/deal.II/numerics/vectors.h index aa5103da74..fb3e9e4b84 100644 --- 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 H1 norm and seminorm, the same rule applies as for the - * L2 norm: compute the l2 norm of the cell error vector. + * L2 norm: compute the l2 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. * * * 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 ec4ac460eb..f79ead7a06 100644 --- 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