From 62084c5d9c468c81d3324fbdb7465ee64de99763 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 5 Aug 2011 19:44:57 +0000 Subject: [PATCH] Produce output that shows the solution as a whole and inside the cells as well. Also compute the L2 error. git-svn-id: https://svn.dealii.org/trunk@24023 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-47/step-47.cc | 163 +++++++++++++++++++++++++--- 1 file changed, 148 insertions(+), 15 deletions(-) diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index 420b63728a..981905e091 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -74,9 +74,6 @@ class LaplaceProblem void run (); private: - double level_set (const Point &p) const; - Tensor<1,dim> grad_level_set (const Point &p) const; - bool interface_intersects_cell (const typename Triangulation::cell_iterator &cell) const; std::pair > compute_quadrature(const Quadrature &plain_quadrature, const typename hp::DoFHandler::active_cell_iterator &cell, const std::vector &level_set_values); void append_quadrature(const Quadrature &plain_quadrature, @@ -89,6 +86,7 @@ class LaplaceProblem void solve (); void refine_grid (); void output_results (const unsigned int cycle) const; + void compute_error () const; Triangulation triangulation; @@ -183,8 +181,7 @@ LaplaceProblem::~LaplaceProblem () template double -LaplaceProblem:: -level_set (const Point &p) const +level_set (const Point &p) { return p.norm() - 0.5; } @@ -193,8 +190,7 @@ level_set (const Point &p) const template Tensor<1,dim> -LaplaceProblem:: -grad_level_set (const Point &p) const +grad_level_set (const Point &p) { return p / p.norm(); } @@ -336,11 +332,6 @@ void LaplaceProblem::assemble_system () cell = dof_handler.begin_active(), endc = dof_handler.end(); - std::vector level_set_values; - level_set_values.push_back(1.); - level_set_values.push_back(1.); - level_set_values.push_back(1.); - level_set_values.push_back(-3.); for (; cell!=endc; ++cell) { const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; @@ -691,7 +682,6 @@ LaplaceProblem::compute_quadrature (const Quadrature &plain_quadrature {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}} }; - std::cout << "Pos : " << Pos << std::endl; for (unsigned int subcell = 0; subcell<5; subcell++) { //std::cout << "subcell : " << subcell << std::endl; @@ -939,6 +929,92 @@ void LaplaceProblem::refine_grid () +template +class Postprocessor : public DataPostprocessor +{ + public: + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; + + virtual std::vector get_names () const; + + virtual unsigned int n_output_variables() const; + + virtual + std::vector + get_data_component_interpretation () const; + + virtual UpdateFlags get_needed_update_flags () const; +}; + + +template +std::vector +Postprocessor::get_names() const +{ + std::vector solution_names (1, "total_solution"); + return solution_names; +} + + +template +unsigned int +Postprocessor::n_output_variables() const +{ + return get_names().size(); +} + + +template +std::vector +Postprocessor:: +get_data_component_interpretation () const +{ + std::vector + interpretation (1, + DataComponentInterpretation::component_is_scalar); + return interpretation; +} + + +template +UpdateFlags +Postprocessor::get_needed_update_flags() const +{ + return update_values | update_q_points; +} + + +template +void +Postprocessor:: +compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const +{ + const unsigned int n_quadrature_points = uh.size(); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + Assert (uh[0].size() == 2, ExcInternalError()); + Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError()); + + for (unsigned int q=0; q void LaplaceProblem::output_results (const unsigned int cycle) const { @@ -950,17 +1026,72 @@ void LaplaceProblem::output_results (const unsigned int cycle) const std::ofstream output (filename.c_str()); + Postprocessor postprocessor; DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); - data_out.build_patches (); + data_out.add_data_vector (solution, postprocessor); + data_out.build_patches (5); data_out.write_vtk (output); } +template +void LaplaceProblem::compute_error () const +{ + hp::QCollection q_collection; + q_collection.push_back (QGauss(2)); + q_collection.push_back (QIterated(QGauss<1>(2), 4)); + + hp::FEValues hp_fe_values (fe_collection, q_collection, + update_values | update_q_points | update_JxW_values); + + double l2_error_square = 0; + + std::vector > solution_values; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + hp_fe_values.reinit (cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values (); + + solution_values.resize (fe_values.n_quadrature_points, + Vector(2)); + fe_values.get_function_values (solution, + solution_values); + + for (unsigned int q=0; q::run () assemble_system (); solve (); + compute_error (); output_results (cycle); } } -- 2.39.5