From: bangerth Date: Thu, 16 Dec 2010 20:45:03 +0000 (+0000) Subject: Make this whole thing more similar to step-4. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bae481ec94acdad663f00bef46df5bd9987c7ed6;p=dealii-svn.git Make this whole thing more similar to step-4. git-svn-id: https://svn.dealii.org/trunk@22978 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 -#include -#include -#include + + #include -#include -#include +#include #include #include #include -#include #include #include -#include #include #include #include #include #include -#include #include #include #include @@ -37,33 +31,26 @@ #include #include -#include -#include - -#define deal_II_dimension 3 - -using std::cout; -using std::endl; using namespace dealii; template class Solution : public Function { -public: - Solution () : Function() {} + public: + Solution () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; }; template double Solution::value (const Point &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::gradient (const Point &p, template class RightHandSide : public Function { -public: - RightHandSide () : Function() {} + public: + RightHandSide () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; }; template @@ -101,7 +88,7 @@ double RightHandSide::value (const Point &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 class LaplaceBeltrami { public: - LaplaceBeltrami (Triangulation *tria, Function &func_data, - unsigned int fe_degree = 1, unsigned int mapping_degree = 1, - Function *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 *pTria; - FE_Q fe; - DoFHandler dh; - MappingQ mapping; + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + MappingQ mapping; - ConstraintMatrix matrix_constraints; + ConstraintMatrix matrix_constraints; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; - Vector solution; - Vector system_rhs; - - Function &rhs_func; // function data - - Function *pExact; // exact solution if provided - + Vector solution; + Vector system_rhs; }; template -LaplaceBeltrami::LaplaceBeltrami (Triangulation *tria, - Function &func_data, - unsigned int fe_degree, - unsigned int mapping_degree, - Function *pExact - ) +LaplaceBeltrami::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 LaplaceBeltrami::~LaplaceBeltrami () { - dh.clear (); + dof_handler.clear (); } + +template +void LaplaceBeltrami::make_mesh () +{ + std::map::cell_iterator, + typename Triangulation::face_iterator> + surface_to_volume_mapping; + + 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); + + static HyperBallBoundary surface_description; + triangulation.set_boundary (1, surface_description); + triangulation.set_boundary (0, surface_description); + + std::set 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 void LaplaceBeltrami::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 @@ -257,9 +267,11 @@ void LaplaceBeltrami::assemble_system () std::vector< double > rhs_values(n_q_points); std::vector local_dof_indices (dofs_per_cell); + const RightHandSide rhs; + typename DoFHandler::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::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; iget_dof_indices (local_dof_indices); for (unsigned int i=0; i bdy_values; - VectorTools::interpolate_boundary_values (mapping,dh, - 0, - *pExact, - bdy_values - ); + std::map bdy_values; + VectorTools::interpolate_boundary_values (mapping,dof_handler, + 0, + Solution(), + 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::output_results () const std::ofstream output (filename.c_str()); DataOut > data_out; - data_out.attach_dof_handler (dh); + data_out.attach_dof_handler (dof_handler); @@ -364,93 +376,70 @@ void LaplaceBeltrami::output_results () const template -void LaplaceBeltrami::run () -{ - - setup_system(); - assemble_system(); - solve(); - output_results(); - +void LaplaceBeltrami::compute_error () const +{ + Vector difference_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (mapping, dof_handler, solution, + Solution(), + difference_per_cell, + QGauss(2*fe.degree+1), + VectorTools::H1_norm); + + + std::cout << "H1 error = " + << difference_per_cell.l2_norm() + << std::endl; } -//################################################################################// -template -double LaplaceBeltrami::compute_error (VectorTools::NormType norm_type) const -{ - Assert(pExact != 0, ExcInternalError()); - - Vector 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 +void LaplaceBeltrami::run () +{ + make_mesh (); + setup_system (); + assemble_system (); + solve (); + output_results (); + compute_error (); } - int main ( int argc, char **argv ) { - std::cout< tria; - - // create a mesh consisting on the boundary of the half sphere... thx Seba - std::map< Triangulation::cell_iterator, - Triangulation::face_iterator> surface_to_volume_mapping; - - 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); - - static HyperBallBoundary surface_description; - tria.set_boundary (1, surface_description); - tria.set_boundary (0, surface_description); - - std::set 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 rhs; - Solution exact; - - LaplaceBeltrami 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<