]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Modified the code a little bit to make it similar to step-4
authorpauletti <pauletti@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 18 Dec 2010 22:55:01 +0000 (22:55 +0000)
committerpauletti <pauletti@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 18 Dec 2010 22:55:01 +0000 (22:55 +0000)
Remove a few redundant lines of code

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

deal.II/examples/step-38/doc/intro.dox
deal.II/examples/step-38/step-38.cc

index 623bf96a9bb265fc544c1f0f8edc95d6f427288d..ea3092393310e3e69e71dd9cabe6786fbe4ab5f2 100644 (file)
@@ -120,5 +120,12 @@ and
 <h3>Implementation</h3>
 
 <h3>Testcase</h3>
+In general when you want to test numerically the accuracy and/or 
+order of convergence of an algorithm you need to provide an exact
+solution. 
+The usual trick is to pick a function that we want to be the solution applied 
+the differential operator to define a forcing term for the right hand side.
+This is what we do in this example.
+
 
 Mapping objects
index e2985c233e00e74963bc2ffac0a921b6f73c3eff..efa7d178ed74ea3833a0b10e9b9f92dc891306d1 100644 (file)
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
+                                // @sect3{Include files}
 
+                                // The first few (many?) include
+                                // files have already been used in
+                                // example 4, so we will
+                                // not explain their meaning here
+                                // again.
 #include <base/quadrature_lib.h>
 #include <base/function.h>
-
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
 #include <lac/full_matrix.h>
 #include <lac/vector.h>
 #include <lac/solver_control.h>
 #include <lac/precondition.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/tria_boundary_lib.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_tools.h>
-
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
-
 #include <fe/fe_values.h>
-
 #include <numerics/data_out.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 
 using namespace dealii;
 
+                                /**
+                                   @sect3{The <code>LaplaceBeltramiProblem</code> class template}
+
+                                   This class is extremely similar to the
+                                   <code>LaplaceProblem</code> class as in 
+                                   example 4.
+                                   One difference is that now some members
+                                   will be defined with two template parameters
+                                   one for the dimension of the mesh @p dim,
+                                   and the other for the dimension of
+                                   the embedding space @p spacedim.
+                                   Now <code>MappingQ</code> appears.
+                                */
+template <int spacedim>
+class LaplaceBeltramiProblem 
+{
+  private:
+    static const unsigned int dim = spacedim-1;
+        
+  public:
+    LaplaceBeltramiProblem (const unsigned degree = 2);
+    void run ();
+  
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results () const;
+    void compute_error () const;
+  
+  
+  
+    Triangulation<dim,spacedim>   triangulation;
+    FE_Q<dim,spacedim>            fe;
+    DoFHandler<dim,spacedim>      dof_handler;
+    MappingQ<dim, spacedim>       mapping;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+  
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
 
 template <int dim>
 class Solution  : public Function<dim>
@@ -85,7 +129,7 @@ Solution<3>::gradient (const Point<3>   &p,
   return return_value;
 }
 
-
+// LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
 template <int dim>
 class RightHandSide : public Function<dim>
 {
@@ -96,7 +140,6 @@ class RightHandSide : public Function<dim>
                          const unsigned int  component = 0) const;
 };
 
-
 template <>
 double
 RightHandSide<3>::value (const Point<3> &p,
@@ -104,7 +147,7 @@ RightHandSide<3>::value (const Point<3> &p,
 {
   using numbers::PI;
 
-                                  // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
+                                  
 
   Tensor<2,3> hessian;
 
@@ -126,80 +169,37 @@ RightHandSide<3>::value (const Point<3> &p,
   gradient[1] = - PI * sin(PI*p(0))*sin(PI*p(1))*exp(p(2));
   gradient[2] = sin(PI*p(0))*cos(PI*p(1))*exp(p(2));
 
-  double curvature;
-  Point<3> normal;
-  double dLength;  
-
-  curvature = 3-1; // dim-1
-  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;
+  Point<3> normal = p;
+  normal /= p.norm();
   
+   
   return (-trace(hessian) +
          (hessian * normal) * normal +
-         (gradient * normal) * curvature);
+         (gradient * normal) * 2.);
 }
 
 
 
 
 
-
-template <int dim>
-class LaplaceBeltrami 
-{
-  public:
-    LaplaceBeltrami ();
-    void run ();
-  
-  private:
-    void make_mesh ();
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void output_results () const;
-    void compute_error () const;
-  
-  
-  
-    Triangulation<dim-1,dim>   triangulation;
-    FE_Q<dim-1,dim>            fe;
-    DoFHandler<dim-1,dim>      dof_handler;
-    MappingQ<dim-1, dim>       mapping;
-
-    SparsityPattern      sparsity_pattern;
-    SparseMatrix<double> system_matrix;
-  
-    Vector<double>       solution;
-    Vector<double>       system_rhs;
-};
-
-
-template <int dim>
-LaplaceBeltrami<dim>::LaplaceBeltrami ()
+template <int spacedim>
+LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem (const unsigned degree)
                :
-               fe (2),
+               fe (degree),
                dof_handler(triangulation),
-               mapping(2)
+               mapping(degree)
 {}
 
 
 
-template <int dim>
-void LaplaceBeltrami<dim>::make_mesh ()
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs ()
 {
-  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);
+  Triangulation<spacedim> volume_mesh;
+  GridGenerator::half_hyper_ball(volume_mesh);
   
-  static HyperBallBoundary<dim-1,dim> surface_description;
-  triangulation.set_boundary (1, surface_description);
+  static HyperBallBoundary<dim,spacedim> surface_description;
   triangulation.set_boundary (0, surface_description);
   
   std::set<unsigned char> boundary_ids;
@@ -207,49 +207,42 @@ void LaplaceBeltrami<dim>::make_mesh ()
   
   GridTools::extract_boundary_mesh (volume_mesh, triangulation,
                                    boundary_ids);
-  triangulation.refine_global(3);
+  triangulation.refine_global(4);
 
   std::cout << "Surface mesh has " << triangulation.n_active_cells()
            << " cells."
            << std::endl;
-}
-
 
-template <int dim>
-void LaplaceBeltrami<dim>::setup_system ()
-{
   dof_handler.distribute_dofs (fe);
 
   std::cout << "Surface mesh has " << dof_handler.n_dofs()
            << " degrees of freedom."
            << std::endl;
   
-  CompressedSparsityPattern csp (dof_handler.n_dofs(),
-                                 dof_handler.n_dofs());
-
+  CompressedSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs());
   DoFTools::make_sparsity_pattern (dof_handler, csp);
-
   sparsity_pattern.copy_from (csp);
 
   system_matrix.reinit (sparsity_pattern);
+  
   solution.reinit (dof_handler.n_dofs());
   system_rhs.reinit (dof_handler.n_dofs());
 }
 
 
-template <int dim>
-void LaplaceBeltrami<dim>::assemble_system () 
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::assemble_system () 
 {
   system_matrix = 0;
   system_rhs = 0;
   
-  QGauss<dim-1>  quadrature_formula(2);
+  QGauss<dim>  quadrature_formula(2);
 
-  FEValues<dim-1,dim> fe_values (mapping, fe, quadrature_formula, 
-                                update_values              |
-                                update_gradients           |
-                                update_quadrature_points   |
-                                update_JxW_values);
+  FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula, 
+                                   update_values              |
+                                   update_gradients           |
+                                   update_quadrature_points   |
+                                   update_JxW_values);
 
   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
   const unsigned int   n_q_points    = quadrature_formula.size();
@@ -260,13 +253,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;
+  const RightHandSide<spacedim> rhs;
   
-  typename DoFHandler<dim-1,dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-
-  for (; cell!=endc; ++cell)
+  for (typename DoFHandler<dim,spacedim>::active_cell_iterator
+        cell = dof_handler.begin_active(),
+        endc = dof_handler.end(); cell!=endc; ++cell)
     {
       cell_matrix = 0;
       cell_rhs = 0;
@@ -306,7 +297,7 @@ void LaplaceBeltrami<dim>::assemble_system ()
   VectorTools::interpolate_boundary_values (mapping,
                                            dof_handler,
                                            0,
-                                           Solution<dim>(),
+                                           Solution<spacedim>(),
                                            boundary_values);
   MatrixTools::apply_boundary_values (boundary_values,
                                      system_matrix,
@@ -315,8 +306,8 @@ void LaplaceBeltrami<dim>::assemble_system ()
 }
 
 
-template <int dim>
-void LaplaceBeltrami<dim>::solve () 
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::solve () 
 {
   SolverControl solver_control (solution.size(), 1e-7);
   SolverCG<>    cg (solver_control);
@@ -330,13 +321,13 @@ void LaplaceBeltrami<dim>::solve ()
 
 
 
-template <int dim>
-void LaplaceBeltrami<dim>::output_results () const
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::output_results () const
 {
-  DataOut<dim-1,DoFHandler<dim-1,dim> > data_out;
+  DataOut<dim,DoFHandler<dim,spacedim> > data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution, "solution",
-                           DataOut<dim-1,DoFHandler<dim-1,dim> >::type_dof_data);
+                           DataOut<dim,DoFHandler<dim,spacedim> >::type_dof_data);
   data_out.build_patches (mapping,
                          mapping.get_degree());
 
@@ -346,14 +337,14 @@ void LaplaceBeltrami<dim>::output_results () const
 
 
 
-template <int dim>
-void LaplaceBeltrami<dim>::compute_error () const
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::compute_error () const
 {  
   Vector<float> difference_per_cell (triangulation.n_active_cells());
   VectorTools::integrate_difference (mapping, dof_handler, solution,
-                                    Solution<dim>(),
+                                    Solution<spacedim>(),
                                     difference_per_cell,
-                                    QGauss<dim-1>(2*fe.degree+1),
+                                    QGauss<dim>(2*fe.degree+1),
                                     VectorTools::H1_norm);  
   
   std::cout << "H1 error = "
@@ -363,11 +354,10 @@ void LaplaceBeltrami<dim>::compute_error () const
 
 
 
-template <int dim>
-void LaplaceBeltrami<dim>::run () 
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::run () 
 {
-  make_mesh ();
-  setup_system ();
+  make_grid_and_dofs();
   assemble_system ();
   solve ();
   output_results ();
@@ -378,38 +368,11 @@ void LaplaceBeltrami<dim>::run ()
 
 int main ( int argc, char **argv )
 {
-  try
-    {
-      deallog.depth_console (0);
-
-      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;
-    }
-  
-  
+  deallog.depth_console (0);
+  {
+    LaplaceBeltramiProblem<3> laplace_beltrami_2d;
+    laplace_beltrami_2d.run();
+  }
+    
   return 0;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.