From bae481ec94acdad663f00bef46df5bd9987c7ed6 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 16 Dec 2010 20:45:03 +0000
Subject: [PATCH] Make this whole thing more similar to step-4.

git-svn-id: https://svn.dealii.org/trunk@22978 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-38/step-38.cc | 379 ++++++++++++++--------------
 1 file changed, 184 insertions(+), 195 deletions(-)

diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc
index f514c2bd0b..005b4fa875 100755
--- a/deal.II/examples/step-38/step-38.cc
+++ b/deal.II/examples/step-38/step-38.cc
@@ -1,25 +1,19 @@
-#include <base/polynomial_space.h>
-#include <base/parsed_function.h>
-#include <base/smartpointer.h>
-#include <base/convergence_table.h>
+
+
 #include <base/quadrature_lib.h>
-#include <base/quadrature_selector.h>
-#include <base/utilities.h>
+#include <base/function.h>
 
 #include <lac/full_matrix.h>
 #include <lac/vector.h>
 #include <lac/solver_control.h>
-#include <lac/solver_gmres.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
-#include <lac/constraint_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <lac/compressed_sparsity_pattern.h>
 
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_accessor.h>
-#include <grid/grid_out.h>
 #include <grid/tria_boundary_lib.h>
 #include <grid/grid_generator.h>
 #include <grid/grid_tools.h>
@@ -37,33 +31,26 @@
 #include <fstream>
 #include <iostream>
 
-#include <sys/types.h>
-#include <sys/stat.h>
-
-#define deal_II_dimension 3
-
-using std::cout;
-using std::endl;
 using namespace dealii;
 
 
 template <int dim>
 class Solution  : public Function<dim>
 {
-public:
-  Solution () : Function<dim>() {}
+  public:
+    Solution () : Function<dim>() {}
   
-  virtual double value (const Point<dim>   &p,
-			const unsigned int  component = 0) const;
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
   
-  virtual Tensor<1,dim> gradient (const Point<dim>   &p,
-				  const unsigned int  component = 0) const;
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+				    const unsigned int  component = 0) const;
 
 };
  
 template <int dim>
 double Solution<dim>::value (const Point<dim> &p,
-			  const unsigned int) const 
+			     const unsigned int) const 
 {
   return sin(numbers::PI * p(0))*cos(numbers::PI * p(1))*exp(p(2)); 
 }
@@ -86,11 +73,11 @@ Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
 template <int dim>
 class RightHandSide : public Function<dim>
 {
-public:
-  RightHandSide () : Function<dim>() {}
+  public:
+    RightHandSide () : Function<dim>() {}
   
-  virtual double value (const Point<dim>   &p,
-			const unsigned int  component = 0) const;
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component = 0) const;
 };
 
 template <int dim>
@@ -101,7 +88,7 @@ double RightHandSide<dim>::value (const Point<dim> &p,
   
   double dPi = numbers::PI;
 
-  // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
+				   // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
 
   Tensor<2,dim> hessian;
 
@@ -146,92 +133,115 @@ template <int dim>
 class LaplaceBeltrami 
 {
   public:
-  LaplaceBeltrami (Triangulation<dim-1,dim> *tria, Function<dim> &func_data, 
-		   unsigned int fe_degree = 1, unsigned int mapping_degree = 1,
-		   Function<dim> *pExact = 0);
-  // arguments are:
-  // triangulation
-  // right-hand side
-  // fe_degree for solution
-  // fe_degree for mapping
-  // exact solution is known
+    LaplaceBeltrami ();
+				     // arguments are:
+				     // triangulation
+				     // right-hand side
+				     // fe_degree for solution
+				     // fe_degree for mapping
+				     // exact solution is known
   
-  ~LaplaceBeltrami ();
-  void run ();
-  double compute_error (VectorTools::NormType norm_type) const;
+    ~LaplaceBeltrami ();
+    void run ();
   
   private:
-  
-  void setup_system ();
-  void assemble_system ();
-  void solve ();
-  void output_results () const;
+
+    void make_mesh ();
+    void setup_system ();
+    void assemble_system ();
+    void solve ();
+    void output_results () const;
+    void compute_error () const;
   
   
   
-  Triangulation<dim-1,dim>   *pTria;
-  FE_Q<dim-1,dim>            fe;
-  DoFHandler<dim-1,dim>      dh;
-  MappingQ<dim-1, dim>       mapping;
+    Triangulation<dim-1,dim>   triangulation;
+    FE_Q<dim-1,dim>            fe;
+    DoFHandler<dim-1,dim>      dof_handler;
+    MappingQ<dim-1, dim>       mapping;
 
-  ConstraintMatrix     matrix_constraints;
+    ConstraintMatrix     matrix_constraints;
 
-  SparsityPattern      sparsity_pattern;
-  SparseMatrix<double> system_matrix;
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
   
-  Vector<double>       solution;
-  Vector<double>       system_rhs;
-
-  Function<dim> &rhs_func; // function data
-
-  Function<dim> *pExact; // exact solution if provided
-
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
 };
 
 
 template <int dim>
-LaplaceBeltrami<dim>::LaplaceBeltrami (Triangulation<dim-1,dim> *tria,
-				       Function<dim> &func_data,
-				       unsigned int fe_degree,
-				       unsigned int mapping_degree,
-				       Function<dim> *pExact
-				       )
+LaplaceBeltrami<dim>::LaplaceBeltrami ()
 		:
-  fe (fe_degree),
-  dh(*tria),
-  mapping(mapping_degree),
-  rhs_func(func_data)
-{
-  pTria = tria;
-  this->pExact = pExact;
+		fe (2),
+		dof_handler(triangulation),
+		mapping(2)
+{}
 
-}
 
 template <int dim>
 LaplaceBeltrami<dim>::~LaplaceBeltrami () 
 {
-  dh.clear ();
+  dof_handler.clear ();
 }
 
+
+template <int dim>
+void LaplaceBeltrami<dim>::make_mesh ()
+{
+  std::map<typename Triangulation<dim-1,dim>::cell_iterator,
+           typename Triangulation<dim,dim>::face_iterator>
+  surface_to_volume_mapping;
+
+  HyperBallBoundary<dim> boundary_description;
+  Triangulation<dim> volume_mesh;
+  GridGenerator::half_hyper_ball(volume_mesh);
+  
+  volume_mesh.set_boundary (1, boundary_description);
+  volume_mesh.set_boundary (0, boundary_description);
+  volume_mesh.refine_global (1);
+  
+  static HyperBallBoundary<dim-1,dim> surface_description;
+  triangulation.set_boundary (1, surface_description);
+  triangulation.set_boundary (0, surface_description);
+  
+  std::set<unsigned char> boundary_ids;
+  boundary_ids.insert(0);
+  
+  GridTools::extract_boundary_mesh (volume_mesh, triangulation,
+				    surface_to_volume_mapping,
+				    boundary_ids);
+  triangulation.refine_global(3);
+
+  std::cout << "Surface mesh has " << triangulation.n_active_cells()
+	    << " cells."
+	    << std::endl;
+}
+
+
 template <int dim>
 void LaplaceBeltrami<dim>::setup_system ()
 {
-  dh.distribute_dofs (fe);
+  dof_handler.distribute_dofs (fe);
 
+  std::cout << "Surface mesh has " << dof_handler.n_dofs()
+	    << " degrees of freedom."
+	    << std::endl;
+  
   matrix_constraints.clear ();
   matrix_constraints.close ();
 
-  CompressedSparsityPattern csp (dh.n_dofs(),
-                                 dh.n_dofs());
+  CompressedSparsityPattern csp (dof_handler.n_dofs(),
+                                 dof_handler.n_dofs());
 
-  DoFTools::make_sparsity_pattern (dh, csp);
+  DoFTools::make_sparsity_pattern (dof_handler, csp);
   matrix_constraints.condense (csp);
 
   sparsity_pattern.copy_from (csp);
 
   system_matrix.reinit (sparsity_pattern);
-  solution.reinit (dh.n_dofs());
-  system_rhs.reinit (dh.n_dofs());
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
 }
 
 template <int dim>
@@ -257,9 +267,11 @@ void LaplaceBeltrami<dim>::assemble_system ()
   std::vector< double > rhs_values(n_q_points);
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
+  const RightHandSide<dim> rhs;
+  
   typename DoFHandler<dim-1,dim>::active_cell_iterator
-    cell = dh.begin_active(),
-    endc = dh.end();
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
 
   for (; cell!=endc; ++cell)
     {
@@ -268,59 +280,59 @@ void LaplaceBeltrami<dim>::assemble_system ()
 
       fe_values.reinit (cell);
 
-      rhs_func.value_list (fe_values.get_quadrature_points(), rhs_values); 
+      rhs.value_list (fe_values.get_quadrature_points(), rhs_values); 
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-      {	  
-	for (unsigned int j=0; j<dofs_per_cell; ++j) 
-	{	      
-	  for (unsigned int q_point=0; q_point<n_q_points;
-	       ++q_point)
-	  {
-	    cell_matrix(i,j) 
-	      += fe_values.shape_grad(i,q_point) *
-	      fe_values.shape_grad(j,q_point) *
-	      fe_values.JxW(q_point);
-	  }
+	{	  
+	  for (unsigned int j=0; j<dofs_per_cell; ++j) 
+	    {	      
+	      for (unsigned int q_point=0; q_point<n_q_points;
+		   ++q_point)
+		{
+		  cell_matrix(i,j) 
+		    += fe_values.shape_grad(i,q_point) *
+		    fe_values.shape_grad(j,q_point) *
+		    fe_values.JxW(q_point);
+		}
+	    }
 	}
-      }
 
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-      {
+	{
 	  
-	for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-	  cell_rhs(i) += fe_values.shape_value(i,q_point) *
-	    rhs_values[q_point]*
-	    fe_values.JxW(q_point);
-      }
+	  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	    cell_rhs(i) += fe_values.shape_value(i,q_point) *
+			   rhs_values[q_point]*
+			   fe_values.JxW(q_point);
+	}
 
       cell->get_dof_indices (local_dof_indices);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-      {
-	for (unsigned int j=0; j<dofs_per_cell; ++j)
-	  system_matrix.add (local_dof_indices[i],
-			     local_dof_indices[j],
-			     cell_matrix(i,j));
+	{
+	  for (unsigned int j=0; j<dofs_per_cell; ++j)
+	    system_matrix.add (local_dof_indices[i],
+			       local_dof_indices[j],
+			       cell_matrix(i,j));
 	
-	system_rhs(local_dof_indices[i]) += cell_rhs(i);
-      }
+	  system_rhs(local_dof_indices[i]) += cell_rhs(i);
+	}
     }
 
 
-   std::map<unsigned int,double> bdy_values; 
-   VectorTools::interpolate_boundary_values (mapping,dh,
-					     0,
-					     *pExact,
-					     bdy_values
-					     );
+  std::map<unsigned int,double> bdy_values; 
+  VectorTools::interpolate_boundary_values (mapping,dof_handler,
+					    0,
+					    Solution<dim>(),
+					    bdy_values
+  );
 
-   MatrixTools::apply_boundary_values (bdy_values,
-				       system_matrix,
-				       solution,
-				       system_rhs,false);
+  MatrixTools::apply_boundary_values (bdy_values,
+				      system_matrix,
+				      solution,
+				      system_rhs,false);
 
-  // condense matrices
+				   // condense matrices
   matrix_constraints.condense (system_matrix);
   matrix_constraints.condense (system_rhs);
 }
@@ -351,7 +363,7 @@ void LaplaceBeltrami<dim>::output_results () const
   std::ofstream output (filename.c_str());
 
   DataOut<dim-1,DoFHandler<dim-1,dim> > data_out;
-  data_out.attach_dof_handler (dh);
+  data_out.attach_dof_handler (dof_handler);
 
  
 
@@ -364,93 +376,70 @@ void LaplaceBeltrami<dim>::output_results () const
 
 
 template <int dim>
-void LaplaceBeltrami<dim>::run () 
-{
-
-  setup_system();
-  assemble_system();
-  solve();
-  output_results();
-
+void LaplaceBeltrami<dim>::compute_error () const
+{  
+  Vector<float> difference_per_cell (triangulation.n_active_cells());
+  VectorTools::integrate_difference (mapping, dof_handler, solution,
+				     Solution<dim>(),
+				     difference_per_cell,
+				     QGauss<dim-1>(2*fe.degree+1),
+				     VectorTools::H1_norm);
+  
+  
+  std::cout << "H1 error = "
+	    << difference_per_cell.l2_norm()
+	    << std::endl;
 }
 
-//################################################################################//
 
-template <int dim>
-double LaplaceBeltrami<dim>::compute_error (VectorTools::NormType norm_type) const
-{  
-  Assert(pExact != 0, ExcInternalError());
-  
-  Vector<float> difference_per_cell (pTria->n_active_cells());
-  VectorTools::integrate_difference (mapping, dh, solution,
-				     *pExact, difference_per_cell,
-				     QGauss<(dim-1)>(2*fe.degree+1),
-				     norm_type);
 
-  return difference_per_cell.l2_norm();
+template <int dim>
+void LaplaceBeltrami<dim>::run () 
+{
+  make_mesh ();
+  setup_system ();
+  assemble_system ();
+  solve ();
+  output_results ();
+  compute_error ();
 }
 
 
 
-
 int main ( int argc, char **argv )
 {
-  std::cout<<std::endl;
-  std::cout<<"================================"<<std::endl;
-  std::cout<<"          LB PROBLEM "<<std::endl;
-  std::cout<<"================================"<<std::endl;
-  std::cout<<std::endl;
-  
-  
-  Triangulation<deal_II_dimension-1,deal_II_dimension> tria;
-  
-  // create a mesh consisting on the boundary of the half sphere... thx Seba
-  std::map< Triangulation<deal_II_dimension-1,deal_II_dimension>::cell_iterator,
-    Triangulation<deal_II_dimension,deal_II_dimension>::face_iterator> surface_to_volume_mapping;
-
-  HyperBallBoundary<deal_II_dimension> boundary_description;
-  Triangulation<deal_II_dimension> volume_mesh;
-  GridGenerator::half_hyper_ball(volume_mesh);
-  
-  volume_mesh.set_boundary (1, boundary_description);
-  volume_mesh.set_boundary (0, boundary_description);
-  volume_mesh.refine_global (1);
-  
-  static HyperBallBoundary<deal_II_dimension-1,deal_II_dimension> surface_description;
-  tria.set_boundary (1, surface_description);
-  tria.set_boundary (0, surface_description);
-  
-  std::set<unsigned char> boundary_ids;
-  boundary_ids.insert(0);
-  
-  GridTools::extract_boundary_mesh (volume_mesh, tria,
-				    surface_to_volume_mapping,
-				    boundary_ids);
- 
-
-  tria.refine_global(4);
-      
+  try
+    {
+      deallog.depth_console (0);
 
-  RightHandSide<deal_II_dimension> rhs;
-  Solution<deal_II_dimension> exact;
-  
-  LaplaceBeltrami<deal_II_dimension> laplace_beltrami_2d(&tria,rhs,2,2,&exact);
+      LaplaceBeltrami<3> laplace_beltrami_2d;
+      laplace_beltrami_2d.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+		<< exc.what() << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+		<< "Aborting!" << std::endl
+		<< "----------------------------------------------------"
+		<< std::endl;
+      return 1;
+    }
   
-  laplace_beltrami_2d.run();
   
-  std::cout<<laplace_beltrami_2d.compute_error(VectorTools::H1_norm)<<std::endl;
-
-  tria.set_boundary(0);
   return 0;
 }
-
-
-
-
-
-
-
-
-
-
-
-- 
2.39.5