From 5483b9468a79a0db9164fba59231f8d9ba2f5095 Mon Sep 17 00:00:00 2001 From: pauletti Date: Sat, 18 Dec 2010 22:55:01 +0000 Subject: [PATCH] Modified the code a little bit to make it similar to step-4 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 | 7 + deal.II/examples/step-38/step-38.cc | 247 +++++++++++-------------- 2 files changed, 112 insertions(+), 142 deletions(-) diff --git a/deal.II/examples/step-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox index 623bf96a9b..ea30923933 100644 --- a/deal.II/examples/step-38/doc/intro.dox +++ b/deal.II/examples/step-38/doc/intro.dox @@ -120,5 +120,12 @@ and

Implementation

Testcase

+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 diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index e2985c233e..efa7d178ed 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -10,10 +10,21 @@ /* 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 #include - +#include +#include +#include +#include +#include +#include #include #include #include @@ -21,20 +32,10 @@ #include #include #include - -#include -#include -#include -#include -#include -#include - #include #include #include - #include - #include #include #include @@ -44,6 +45,49 @@ using namespace dealii; + /** + @sect3{The LaplaceBeltramiProblem class template} + + This class is extremely similar to the + LaplaceProblem 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 MappingQ appears. + */ +template +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 triangulation; + FE_Q fe; + DoFHandler dof_handler; + MappingQ mapping; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; template class Solution : public Function @@ -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 class RightHandSide : public Function { @@ -96,7 +140,6 @@ class RightHandSide : public Function 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 -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 triangulation; - FE_Q fe; - DoFHandler dof_handler; - MappingQ mapping; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; -}; - - -template -LaplaceBeltrami::LaplaceBeltrami () +template +LaplaceBeltramiProblem::LaplaceBeltramiProblem (const unsigned degree) : - fe (2), + fe (degree), dof_handler(triangulation), - mapping(2) + mapping(degree) {} -template -void LaplaceBeltrami::make_mesh () +template +void LaplaceBeltramiProblem::make_grid_and_dofs () { - HyperBallBoundary boundary_description; - Triangulation 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 volume_mesh; + GridGenerator::half_hyper_ball(volume_mesh); - static HyperBallBoundary surface_description; - triangulation.set_boundary (1, surface_description); + static HyperBallBoundary surface_description; triangulation.set_boundary (0, surface_description); std::set boundary_ids; @@ -207,49 +207,42 @@ void LaplaceBeltrami::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 -void LaplaceBeltrami::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 -void LaplaceBeltrami::assemble_system () +template +void LaplaceBeltramiProblem::assemble_system () { system_matrix = 0; system_rhs = 0; - QGauss quadrature_formula(2); + QGauss quadrature_formula(2); - FEValues fe_values (mapping, fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); + FEValues 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::assemble_system () std::vector< double > rhs_values(n_q_points); std::vector local_dof_indices (dofs_per_cell); - const RightHandSide rhs; + const RightHandSide rhs; - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - for (; cell!=endc; ++cell) + for (typename DoFHandler::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::assemble_system () VectorTools::interpolate_boundary_values (mapping, dof_handler, 0, - Solution(), + Solution(), boundary_values); MatrixTools::apply_boundary_values (boundary_values, system_matrix, @@ -315,8 +306,8 @@ void LaplaceBeltrami::assemble_system () } -template -void LaplaceBeltrami::solve () +template +void LaplaceBeltramiProblem::solve () { SolverControl solver_control (solution.size(), 1e-7); SolverCG<> cg (solver_control); @@ -330,13 +321,13 @@ void LaplaceBeltrami::solve () -template -void LaplaceBeltrami::output_results () const +template +void LaplaceBeltramiProblem::output_results () const { - DataOut > data_out; + DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution", - DataOut >::type_dof_data); + DataOut >::type_dof_data); data_out.build_patches (mapping, mapping.get_degree()); @@ -346,14 +337,14 @@ void LaplaceBeltrami::output_results () const -template -void LaplaceBeltrami::compute_error () const +template +void LaplaceBeltramiProblem::compute_error () const { Vector difference_per_cell (triangulation.n_active_cells()); VectorTools::integrate_difference (mapping, dof_handler, solution, - Solution(), + Solution(), difference_per_cell, - QGauss(2*fe.degree+1), + QGauss(2*fe.degree+1), VectorTools::H1_norm); std::cout << "H1 error = " @@ -363,11 +354,10 @@ void LaplaceBeltrami::compute_error () const -template -void LaplaceBeltrami::run () +template +void LaplaceBeltramiProblem::run () { - make_mesh (); - setup_system (); + make_grid_and_dofs(); assemble_system (); solve (); output_results (); @@ -378,38 +368,11 @@ void LaplaceBeltrami::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; } -- 2.39.5