From 7526f84a7396032f69afb073881f84c9d3a8d1a3 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 1 Sep 2009 13:49:44 +0000
Subject: [PATCH] The Laplacian is a symmetric operation, so only store the
 symmetric part of it. Produce even 3D output in results.

git-svn-id: https://svn.dealii.org/trunk@19361 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-37/Makefile        |  1 +
 deal.II/examples/step-37/doc/results.dox | 41 +++++++++++++++++++-
 deal.II/examples/step-37/step-37.cc      | 49 +++++++++++++++++++-----
 3 files changed, 79 insertions(+), 12 deletions(-)

diff --git a/deal.II/examples/step-37/Makefile b/deal.II/examples/step-37/Makefile
index feef6229a0..ebd829a479 100644
--- a/deal.II/examples/step-37/Makefile
+++ b/deal.II/examples/step-37/Makefile
@@ -59,6 +59,7 @@ include $D/common/Make.global_options
 # You may need to augment the lists of libraries when compiling your
 # program for other dimensions, or when using third party libraries
 libs.g   = $(lib-deal2-2d.g) \
+	   $(lib-deal2-3d.g) \
 	   $(lib-lac.g)      \
            $(lib-base.g)
 libs.o   = $(lib-deal2-2d.o) \
diff --git a/deal.II/examples/step-37/doc/results.dox b/deal.II/examples/step-37/doc/results.dox
index dfc70b7f4a..01670dfeb4 100644
--- a/deal.II/examples/step-37/doc/results.dox
+++ b/deal.II/examples/step-37/doc/results.dox
@@ -4,7 +4,8 @@ Since this example solves the same problem as @ref step_5 "step-5", we
 refer to the graphical output there. The only difference between the two is 
 the solver and the implementation of the matrix-vector products. 
 
-The output produced by this program is the following:
+When we run this program in 2D for quadratic elements, we get the following 
+output:
 @code
 Cycle 0
 Number of degrees of freedom: 337
@@ -44,4 +45,40 @@ Convergence in 9 CG iterations.
 @endcode
 
 As in step-16, we see that the number of CG iterations stays constant with 
-increasing number of degrees of freedom.
\ No newline at end of file
+increasing number of degrees of freedom. Not much changes if we run the 
+program in three dimensions:
+@code
+Cycle 0
+Number of degrees of freedom: 517
+System matrix memory consumption: 0.09649 MBytes.
+Multigrid objects memory consumption: 0.1389 MBytes.
+Convergence in 7 CG iterations.
+
+Cycle 1
+Number of degrees of freedom: 3817
+System matrix memory consumption: 0.6333 MBytes.
+Multigrid objects memory consumption: 0.8396 MBytes.
+Convergence in 8 CG iterations.
+
+Cycle 2
+Number of degrees of freedom: 29521
+System matrix memory consumption: 4.882 MBytes.
+Multigrid objects memory consumption: 6.273 MBytes.
+Convergence in 9 CG iterations.
+
+Cycle 3
+Number of degrees of freedom: 232609
+System matrix memory consumption: 38.68 MBytes.
+Multigrid objects memory consumption: 49.26 MBytes.
+Convergence in 10 CG iterations.
+
+Cycle 4
+Number of degrees of freedom: 1847617
+System matrix memory consumption: 308.4 MBytes.
+Multigrid objects memory consumption: 391.5 MBytes.
+Convergence in 11 CG iterations.
+@endcode
+
+There is a slight increase in the number of CG iterations, but also for the
+largest size with almost two million unknowns, the iteration count is still
+very modest.
\ No newline at end of file
diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc
index 83d9fb172e..8a1e398dd6 100644
--- a/deal.II/examples/step-37/step-37.cc
+++ b/deal.II/examples/step-37/step-37.cc
@@ -469,7 +469,7 @@ public:
   LaplaceOperator<dim,number>&
   operator = (const Tensor<2,dim> &tensor);
 
-  number transformation[dim][dim];
+  number transformation[dim*(dim+1)/2];
 };
 
 template<int dim,typename number>
@@ -490,20 +490,49 @@ void LaplaceOperator<dim,number>::transform (number* result) const
     temp_result[d] = result[d];
   for (unsigned int d=0; d<dim; ++d)
     {
-      number output = 0;
-      for (unsigned int e=0; e<dim; ++e)
-	output += transformation[d][e] * temp_result[e];
+      number output = transformation[d]*temp_result[d];
+      if (dim == 2)
+	output += transformation[2]*temp_result[1-d];
+      else if (dim == 3)
+	{
+	  if (d==0)
+	    output += transformation[3]*temp_result[1] + transformation[4]*temp_result[2];
+	  else if (d==1)
+	    output += transformation[3]*temp_result[0] + transformation[5]*temp_result[2];
+	  else
+	    output += transformation[4]*temp_result[0] + transformation[5]*temp_result[1];
+	}
       result[d] = output;
     }
 }
 
 template <int dim, typename number>
 LaplaceOperator<dim,number>&
-LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor)
+LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor) 
 {
-  for (unsigned int d=0;d<dim;++d)
-    for (unsigned int e=0;e<dim;++e)
-      transformation[d][e] = tensor[d][e];
+  if (dim == 2)
+    {
+      transformation[0] = tensor[0][0];
+      transformation[1] = tensor[1][1];
+      transformation[2] = tensor[0][1];
+      Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15,
+	      ExcInternalError());
+    }
+  else if (dim == 3)
+    {
+      transformation[0] = tensor[0][0];
+      transformation[1] = tensor[1][1];
+      transformation[2] = tensor[2][2];
+      transformation[3] = tensor[0][1];
+      transformation[4] = tensor[0][2];
+      transformation[5] = tensor[1][2];
+      Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15,
+	      ExcInternalError());
+      Assert (std::fabs(tensor[2][0]-tensor[0][2])<1e-15,
+	      ExcInternalError());
+      Assert (std::fabs(tensor[2][1]-tensor[1][2])<1e-15,
+	      ExcInternalError());
+    }
   return *this;
 }
 
@@ -851,7 +880,7 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 template <int dim>
 void LaplaceProblem<dim>::run () 
 {
-  for (unsigned int cycle=0; cycle<6; ++cycle)
+  for (unsigned int cycle=0; cycle<8-dim; ++cycle)
     {
       std::cout << "Cycle " << cycle << std::endl;
 
@@ -861,7 +890,7 @@ void LaplaceProblem<dim>::run ()
 	  GridGenerator::hyper_ball(triangulation);
 	  static const HyperBallBoundary<dim> boundary;
 	  triangulation.set_boundary (0, boundary);
-	  triangulation.refine_global (0);
+	  triangulation.refine_global (3-dim);
 	}
       triangulation.refine_global (1);
       setup_system ();
-- 
2.39.5