Implement W1infty norm and seminorm in VectorTools::integrate_difference.

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

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 8a4f303a513d7edce221df3e0fbd55a77bdcc2e2..f52970946e14e7f3da8d8bd958a553f81ae7e0f5 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -47,6 +47,13 @@ used for boundary indicators.
 
 
 <ol>
+<li>
+New: step-15 has been replaced by a program that demonstrates the
+solution of nonlinear problem (the minimal surface equation) using
+Newton's method.
+<br>
+(Sven Wetterauer, 2012/06/03)
+
 
 <li>
 New: step-48 demonstrates the solution of a nonlinear wave equation
@@ -201,6 +208,12 @@ enabled due to a missing include file in file
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: Computing the $W^{1,\infty}$ norm and seminorm in
+VectorTools::integrate_difference was not implemented. This is now
+fixed.
+<br>
+(Wolfgang Bangerth 2012/06/02)
+
 <li> Fixed: A problem in MappingQ::transform_real_to_unit_cell
 that sometimes led the algorithm in this function to abort.
 <br>

diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h
index 72b02f2b7bcc7465144683a9e42e5820130bcb5a..941af268f5646e60d577d56ae182dbfac213ccf3 100644 (file)
--- a/deal.II/include/deal.II/numerics/vectors.templates.h
+++ b/deal.II/include/deal.II/numerics/vectors.templates.h
@@ -5192,9 +5192,9 @@ namespace VectorTools
                                                   0.0);
                       diff = std::sqrt(diff);
                       break;
+
                 case W1infty_seminorm:
                 case W1infty_norm:
-                      Assert(false, ExcNotImplemented());
                       std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
                       for (unsigned int k=0; k<n_components; ++k)
                         for (unsigned int q=0; q<n_q_points; ++q)
@@ -5207,8 +5207,16 @@ namespace VectorTools
                             psi_scalar[q] = std::max(psi_scalar[q],t);
                           }
 
-                      for (unsigned int i=0;i<psi_scalar.size();++i)
-                        diff = std::max (diff, psi_scalar[i]);
+                      // compute seminorm
+                      {
+                        double t = 0;
+                        for (unsigned int i=0;i<psi_scalar.size();++i)
+                          t = std::max (t, psi_scalar[i]);
+
+                        // then add seminorm to norm if that had previously been
+                        // computed
+                        diff += t;
+                      }
                       break;
                 default:
                       break;

diff --git a/tests/hp/integrate_difference_02.cc b/tests/hp/integrate_difference_02.cc
new file mode 100644 (file)
index 0000000..f8df00e
--- /dev/null
+++ b/tests/hp/integrate_difference_02.cc
@@ -0,0 +1,156 @@
+//----------------------------  integrate_difference.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005, 2006, 2008, 2012 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  integrate_difference.cc  ---------------------------
+
+
+// call VectorTools::integrate_difference with fe's distributed in the
+// same random way as in hp/random and on a function that is d-linear
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vectors.h>
+
+#include <fstream>
+
+
+
+template <int dim>
+void test ()
+{
+  deallog << "dim=" << dim << std::endl;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global (2);
+  tria.begin_active()->set_refine_flag ();
+  tria.execute_coarsening_and_refinement ();
+  tria.refine_global (4-dim);
+
+  hp::FECollection<dim> fe_collection;
+  hp::QCollection<dim> q_collection;
+  for (unsigned int i=1; i<=4; ++i)
+    {
+      fe_collection.push_back(FE_Q<dim> (i));
+      q_collection.push_back (QGauss<dim> (i+2));
+    }
+
+
+  hp::DoFHandler<dim> dof_handler(tria);
+
+  for (typename hp::DoFHandler<dim>::active_cell_iterator
+        cell = dof_handler.begin_active();
+       cell != dof_handler.end(); ++cell)
+    cell->set_active_fe_index (rand() % fe_collection.size());
+
+  dof_handler.distribute_dofs(fe_collection);
+
+                                  // interpolate a linear function
+  Vector<double> vec (dof_handler.n_dofs());
+  VectorTools::interpolate (dof_handler,
+                           Functions::Monomial<dim>
+                           (dim == 1 ?
+                            Point<dim>(1.) :
+                            (dim == 2 ?
+                             Point<dim>(1.,0.) :
+                             Point<dim>(1.,0.,0.))),
+                           vec);
+
+  Vector<float> diff (tria.n_active_cells());
+
+                                  // L1 norm. the function is u(x)=x, so its
+                                  // L1 norm should be equal to 1/2
+  {
+    VectorTools::integrate_difference (dof_handler,
+                                      vec,
+                                      ZeroFunction<dim>(),
+                                      diff,
+                                      q_collection,
+                                      VectorTools::L1_norm);
+    deallog << "L1, diff=" << diff.l1_norm() << std::endl;
+  }
+
+                                  // H1 seminorm. the function is u(x)=x, so
+                                  // its H1 seminorm should be equal to 1/2
+  {
+    VectorTools::integrate_difference (dof_handler,
+                                      vec,
+                                      ZeroFunction<dim>(),
+                                      diff,
+                                      q_collection,
+                                      VectorTools::H1_seminorm);
+    deallog << "H1 seminorm, diff=" << diff.l2_norm() << std::endl;
+  }
+
+                                  // W1infty seminorm. the function is
+                                  // u(x)=x, so the norm must be equal to 1
+                                  // on every cell
+  {
+    VectorTools::integrate_difference (dof_handler,
+                                      vec,
+                                      ZeroFunction<dim>(),
+                                      diff,
+                                      q_collection,
+                                      VectorTools::W1infty_seminorm);
+    deallog << "W1infty semi, diff=" << diff.linfty_norm() << std::endl;
+                                    // also ensure that we indeed get the
+                                    // same value on every cell
+    diff. add(-1);
+    Assert (diff.l2_norm() == 0, ExcInternalError());
+  }
+
+                                  // W1infty norm. the Linfty norm is one, so
+                                  // the W1infty norm must be two. but not on
+                                  // every cell
+  {
+    VectorTools::integrate_difference (dof_handler,
+                                      vec,
+                                      ZeroFunction<dim>(),
+                                      diff,
+                                      q_collection,
+                                      VectorTools::W1infty_norm);
+    deallog << "W1infty, diff=" << diff.linfty_norm() << std::endl;
+    diff.add(-2);
+    Assert (diff.l1_norm() > 0.5, ExcInternalError());
+  }
+}
+
+
+int main ()
+{
+  std::ofstream logfile("integrate_difference_02/output");
+  logfile.precision(2);
+  deallog << std::setprecision(2);
+
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+
+  deallog << "OK" << std::endl;
+}

diff --git a/tests/hp/integrate_difference_02/cmp/generic b/tests/hp/integrate_difference_02/cmp/generic
new file mode 100644 (file)
index 0000000..b5969a8
--- /dev/null
+++ b/tests/hp/integrate_difference_02/cmp/generic
@@ -0,0 +1,17 @@
+
+DEAL::dim=1
+DEAL::L1, diff=0.50
+DEAL::H1 seminorm, diff=1.0
+DEAL::W1infty semi, diff=1.0
+DEAL::W1infty, diff=2.0
+DEAL::dim=2
+DEAL::L1, diff=0.50
+DEAL::H1 seminorm, diff=1.0
+DEAL::W1infty semi, diff=1.0
+DEAL::W1infty, diff=2.0
+DEAL::dim=3
+DEAL::L1, diff=0.50
+DEAL::H1 seminorm, diff=1.0
+DEAL::W1infty semi, diff=1.0
+DEAL::W1infty, diff=2.0
+DEAL::OK