From 5775dc74d7771613614ce03e648fccac439f7cdb Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 16 Jul 2013 19:45:42 +0000
Subject: [PATCH] Fix a bug in VectorTools::project.

git-svn-id: https://svn.dealii.org/trunk@30020 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/news/changes.h                    |   7 +
 .../deal.II/numerics/vector_tools.templates.h |   2 +-
 tests/deal.II/project_01_curved_boundary.cc   | 171 ++++++++++++++++++
 .../project_01_curved_boundary/cmp/generic    | 165 +++++++++++++++++
 4 files changed, 344 insertions(+), 1 deletion(-)
 create mode 100644 tests/deal.II/project_01_curved_boundary.cc
 create mode 100644 tests/deal.II/project_01_curved_boundary/cmp/generic

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index d8442a4bc5..fae3f4bc26 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -164,6 +164,13 @@ this function.
 
 <ol>
 
+<li>Fixed: VectorTools::project has an option to first project onto the
+boundary. However, the implementation of this option ignored the mapping
+that is provided to the function. This is now fixed.
+<br>
+(Wolfgang Bangerth, 2013/07/16)
+</li>
+
 <li>Improved: The WorkStream class used throughout deal.II is now using
 thread local variables and initializes temporary variables on the thread
 that uses them, leading to better cache locality.
diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h
index 6e90a698fa..89452ceffc 100644
--- a/deal.II/include/deal.II/numerics/vector_tools.templates.h
+++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h
@@ -575,7 +575,7 @@ namespace VectorTools
           typename FunctionMap<spacedim>::type boundary_functions;
           for (types::boundary_id c=0; c<numbers::internal_face_boundary_id; ++c)
             boundary_functions[c] = &function;
-          project_boundary_values (dof, boundary_functions, q_boundary,
+          project_boundary_values (mapping, dof, boundary_functions, q_boundary,
                                    boundary_values);
         }
 
diff --git a/tests/deal.II/project_01_curved_boundary.cc b/tests/deal.II/project_01_curved_boundary.cc
new file mode 100644
index 0000000000..9df1f03723
--- /dev/null
+++ b/tests/deal.II/project_01_curved_boundary.cc
@@ -0,0 +1,171 @@
+//----------------------------  project_01_curved_boundary.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2006, 2013 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.
+//
+//----------------------------  project_01_curved_boundary.cc  ---------------------------
+
+// check VectorTools::project for Vector<double> arguments
+
+
+#include "../tests.h"
+#include <deal.II/base/function.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/grid_tools.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <fstream>
+
+
+// define a function that is zero within the square
+//   [-sqrt(2)/2,sqrt(2)/2]^d
+// which is the domain produced by GridGenerator::hyper_sphere
+// when using a linear mapping. we will want to see a solution
+// that is not zero when using a higher order mapping and a Q2
+// element
+template <int dim>
+class F : public Function<dim>
+{
+  public:
+    virtual double value (const Point<dim> &p,
+                          const unsigned int = 0) const
+      {
+	// compute the linfty norm of p
+	double m = 0;
+	for (unsigned int d=0; d<dim; ++d)
+	  m = std::max(m, std::fabs(p[d]));
+
+	// let the value increase linearly away from the square
+	return std::max (0., m-std::sqrt(2.)/2);
+      }
+};
+
+
+template<int dim>
+void test()
+{
+  HyperBallBoundary<dim> boundary;
+
+  // create 2 triangulations with the
+  // same coarse grid, and refine
+  // them differently
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball (tria);
+  tria.set_boundary (0, boundary);
+  
+
+  FE_Q<dim> fe(1);
+  DoFHandler<dim> dh (tria);
+  dh.distribute_dofs (fe);
+
+  Vector<double> v(dh.n_dofs());
+
+  ConstraintMatrix cm;
+  cm.close ();
+
+  // use the implicit Q1 mapping. this will yield a zero solution
+  {
+    VectorTools::project (dh, cm, QGauss<dim>(3), F<dim>(),
+			  v);
+    deallog << v.l2_norm() << std::endl;
+    Assert (v.l2_norm() == 0, ExcInternalError());
+  }
+
+  // use an explicit Q1 mapping. this will yield a zero solution
+  {
+    VectorTools::project (MappingQ1<dim>(),
+			  dh, cm, QGauss<dim>(3), F<dim>(),
+			  v);
+    deallog << v.l2_norm() << std::endl;
+    Assert (v.l2_norm() == 0, ExcInternalError());
+  }
+
+  // use an explicit Q2 mapping. this will yield a nonzero solution if it's a
+  // straight projection since some of the quadrature points are lying outside
+  // the area where the function is zero
+  {
+    VectorTools::project (MappingQ<dim>(2),
+			  dh, cm, QGauss<dim>(3), F<dim>(),
+			  v);
+    deallog << v.l2_norm() << std::endl;
+    Assert (v.l2_norm() != 0, ExcInternalError());
+  }
+
+  // use an explicit Q2 mapping but enforce zero boundary values. this will
+  // yield a nonzero solution since some of the quadrature points are lying
+  // outside the area where the function is zero, even though the values of
+  // the DoFs at the boundary are in fact zero (they are interpolated only at
+  // points where the function is zero)
+  {
+    VectorTools::project (MappingQ<dim>(2),
+			  dh, cm, QGauss<dim>(3), F<dim>(),
+			  v, true);
+    deallog << v.l2_norm() << std::endl;
+    Assert (v.l2_norm() != 0, ExcInternalError());
+  }
+
+  // use an explicit Q2 mapping and project to the boundary first. this will
+  // yield a nonzero solution since some of the quadrature points are lying
+  // outside the area where the function is zero. furthermore, the values on
+  // the boundary will be nonzero since we project, even though the
+  // *interpolation* of boundary values onto the trace of the Q1 space is zero
+  {
+    VectorTools::project (MappingQ<dim>(2),
+			  dh, cm, QGauss<dim>(3), F<dim>(),
+			  v, false,
+			  QGauss<dim-1>(2), true);
+    deallog << v.l2_norm() << std::endl;
+    Assert (v.l2_norm() != 0, ExcInternalError());
+    for (typename DoFHandler<dim>::active_cell_iterator cell=dh.begin_active();
+	 cell != dh.end(); ++cell)
+      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+	deallog << cell->vertex(i) << ' ' << v(cell->vertex_dof_index(i,0))
+		<< std::endl;
+  }
+
+  
+  // same as above, but use a projection with a QTrapez formula. this happens
+  // to evaluate the function only at points where it is zero, and
+  // consequently the values at the boundary should be zero
+  {
+    VectorTools::project (MappingQ<dim>(2),
+			  dh, cm, QGauss<dim>(3), F<dim>(),
+			  v, false,
+			  QTrapez<dim-1>(), true);
+    deallog << v.l2_norm() << std::endl;
+    Assert (v.l2_norm() != 0, ExcInternalError());
+    for (typename DoFHandler<dim>::active_cell_iterator cell=dh.begin_active();
+	 cell != dh.end(); ++cell)
+      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+	deallog << cell->vertex(i) << ' ' << v(cell->vertex_dof_index(i,0))
+		<< std::endl;
+  }
+}
+
+
+int main()
+{
+  std::ofstream logfile ("project_01_curved_boundary/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
+
diff --git a/tests/deal.II/project_01_curved_boundary/cmp/generic b/tests/deal.II/project_01_curved_boundary/cmp/generic
new file mode 100644
index 0000000000..218888d980
--- /dev/null
+++ b/tests/deal.II/project_01_curved_boundary/cmp/generic
@@ -0,0 +1,165 @@
+
+DEAL::0
+DEAL::0
+DEAL::0.227756
+DEAL::0.0359680
+DEAL::0.421875
+DEAL::-0.707107 -0.707107 0.195262
+DEAL::0.707107 -0.707107 0.195262
+DEAL::-0.292893 -0.292893 -0.0797954
+DEAL::0.292893 -0.292893 -0.0797954
+DEAL::-0.707107 -0.707107 0.195262
+DEAL::-0.292893 -0.292893 -0.0797954
+DEAL::-0.707107 0.707107 0.195262
+DEAL::-0.292893 0.292893 -0.0797954
+DEAL::-0.292893 -0.292893 -0.0797954
+DEAL::0.292893 -0.292893 -0.0797954
+DEAL::-0.292893 0.292893 -0.0797954
+DEAL::0.292893 0.292893 -0.0797954
+DEAL::0.707107 -0.707107 0.195262
+DEAL::0.707107 0.707107 0.195262
+DEAL::0.292893 -0.292893 -0.0797954
+DEAL::0.292893 0.292893 -0.0797954
+DEAL::-0.707107 0.707107 0.195262
+DEAL::-0.292893 0.292893 -0.0797954
+DEAL::0.707107 0.707107 0.195262
+DEAL::0.292893 0.292893 -0.0797954
+DEAL::0.0359680
+DEAL::-0.707107 -0.707107 0
+DEAL::0.707107 -0.707107 0
+DEAL::-0.292893 -0.292893 0.0179840
+DEAL::0.292893 -0.292893 0.0179840
+DEAL::-0.707107 -0.707107 0
+DEAL::-0.292893 -0.292893 0.0179840
+DEAL::-0.707107 0.707107 0
+DEAL::-0.292893 0.292893 0.0179840
+DEAL::-0.292893 -0.292893 0.0179840
+DEAL::0.292893 -0.292893 0.0179840
+DEAL::-0.292893 0.292893 0.0179840
+DEAL::0.292893 0.292893 0.0179840
+DEAL::0.707107 -0.707107 0
+DEAL::0.707107 0.707107 0
+DEAL::0.292893 -0.292893 0.0179840
+DEAL::0.292893 0.292893 0.0179840
+DEAL::-0.707107 0.707107 0
+DEAL::-0.292893 0.292893 0.0179840
+DEAL::0.707107 0.707107 0
+DEAL::0.292893 0.292893 0.0179840
+DEAL::0
+DEAL::0
+DEAL::0.166690
+DEAL::0.0709393
+DEAL::0.440214
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 -0.577350 0.115757
+DEAL::0.577350 -0.577350 -0.577350 0.115757
+DEAL::-0.577350 0.577350 -0.577350 0.115757
+DEAL::0.577350 0.577350 -0.577350 0.115757
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.577350 -0.577350 -0.577350 0.115757
+DEAL::0.577350 0.577350 -0.577350 0.115757
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.577350 -0.577350 0.577350 0.115757
+DEAL::0.577350 0.577350 0.577350 0.115757
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 0.577350 0.115757
+DEAL::0.577350 -0.577350 0.577350 0.115757
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.577350 0.577350 0.577350 0.115757
+DEAL::0.577350 0.577350 0.577350 0.115757
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 -0.577350 0.115757
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.577350 0.577350 -0.577350 0.115757
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::-0.577350 -0.577350 0.577350 0.115757
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.577350 0.577350 0.577350 0.115757
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::-0.577350 -0.577350 -0.577350 0.115757
+DEAL::0.577350 -0.577350 -0.577350 0.115757
+DEAL::-0.211325 -0.211325 -0.211325 -0.104038
+DEAL::0.211325 -0.211325 -0.211325 -0.104038
+DEAL::-0.577350 -0.577350 0.577350 0.115757
+DEAL::0.577350 -0.577350 0.577350 0.115757
+DEAL::-0.211325 -0.211325 0.211325 -0.104038
+DEAL::0.211325 -0.211325 0.211325 -0.104038
+DEAL::-0.577350 0.577350 -0.577350 0.115757
+DEAL::-0.211325 0.211325 -0.211325 -0.104038
+DEAL::0.577350 0.577350 -0.577350 0.115757
+DEAL::0.211325 0.211325 -0.211325 -0.104038
+DEAL::-0.577350 0.577350 0.577350 0.115757
+DEAL::-0.211325 0.211325 0.211325 -0.104038
+DEAL::0.577350 0.577350 0.577350 0.115757
+DEAL::0.211325 0.211325 0.211325 -0.104038
+DEAL::0.0709393
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 -0.577350 0
+DEAL::0.577350 -0.577350 -0.577350 0
+DEAL::-0.577350 0.577350 -0.577350 0
+DEAL::0.577350 0.577350 -0.577350 0
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.577350 -0.577350 -0.577350 0
+DEAL::0.577350 0.577350 -0.577350 0
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.577350 -0.577350 0.577350 0
+DEAL::0.577350 0.577350 0.577350 0
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 0.577350 0
+DEAL::0.577350 -0.577350 0.577350 0
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.577350 0.577350 0.577350 0
+DEAL::0.577350 0.577350 0.577350 0
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 -0.577350 0
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.577350 0.577350 -0.577350 0
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::-0.577350 -0.577350 0.577350 0
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.577350 0.577350 0.577350 0
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::-0.577350 -0.577350 -0.577350 0
+DEAL::0.577350 -0.577350 -0.577350 0
+DEAL::-0.211325 -0.211325 -0.211325 0.0250808
+DEAL::0.211325 -0.211325 -0.211325 0.0250808
+DEAL::-0.577350 -0.577350 0.577350 0
+DEAL::0.577350 -0.577350 0.577350 0
+DEAL::-0.211325 -0.211325 0.211325 0.0250808
+DEAL::0.211325 -0.211325 0.211325 0.0250808
+DEAL::-0.577350 0.577350 -0.577350 0
+DEAL::-0.211325 0.211325 -0.211325 0.0250808
+DEAL::0.577350 0.577350 -0.577350 0
+DEAL::0.211325 0.211325 -0.211325 0.0250808
+DEAL::-0.577350 0.577350 0.577350 0
+DEAL::-0.211325 0.211325 0.211325 0.0250808
+DEAL::0.577350 0.577350 0.577350 0
+DEAL::0.211325 0.211325 0.211325 0.0250808
-- 
2.39.5