From 2ef7c00be09cd9bb663ba977f3796dcb0e1b601c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 8 Apr 1998 11:27:14 +0000 Subject: [PATCH] Add new test case to demonstrate the error reduction for a simple problem. git-svn-id: https://svn.dealii.org/trunk@158 0785d39b-7218-0410-832d-ea1e28bc413d --- .../Attic/examples/convergence/convergence.cc | 306 ++++++++++++++++++ tests/big-tests/convergence/convergence.cc | 306 ++++++++++++++++++ 2 files changed, 612 insertions(+) create mode 100644 deal.II/deal.II/Attic/examples/convergence/convergence.cc create mode 100644 tests/big-tests/convergence/convergence.cc diff --git a/deal.II/deal.II/Attic/examples/convergence/convergence.cc b/deal.II/deal.II/Attic/examples/convergence/convergence.cc new file mode 100644 index 0000000000..331afc5baf --- /dev/null +++ b/deal.II/deal.II/Attic/examples/convergence/convergence.cc @@ -0,0 +1,306 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include +extern "C" { +# include +} + + + + +template +class PoissonEquation : public Equation { + public: + PoissonEquation (const Function &rhs) : + Equation(1), + right_hand_side (rhs) {}; + + virtual void assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + protected: + const Function &right_hand_side; +}; + + + + + + +template +class PoissonProblem : public ProblemBase { + public: + PoissonProblem (); + + void clear (); + virtual void create_new (); + virtual void run (unsigned int level); + + protected: + Triangulation *tria; + DoFHandler *dof; + + Function *rhs; + Function *boundary_values; +}; + + + + + +/** + Right hand side constructed such that the exact solution is + $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc. + */ +template +class RHSPoly : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + +template +class Solution : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + + +template +double RHSPoly::operator () (const Point &p) const { + double ret_val = 0; + for (unsigned int i=0; i +double Solution::operator () (const Point &p) const { + double ret_val = 1; + for (unsigned int i=0; i::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point +void PoissonEquation::assemble (dFMatrix &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + +template +void PoissonEquation::assemble (dVector &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + + + + + + + +template +PoissonProblem::PoissonProblem () : + tria(0), dof(0), rhs(0), boundary_values(0) {}; + + + + +template +void PoissonProblem::clear () { + if (tria != 0) { + delete tria; + tria = 0; + }; + + if (dof != 0) { + delete dof; + dof = 0; + }; + + if (rhs != 0) + { + delete rhs; + rhs = 0; + }; + + if (boundary_values != 0) + { + delete boundary_values; + boundary_values = 0; + }; + + ProblemBase::clear (); +}; + + + + +template +void PoissonProblem::create_new () { + clear (); + + tria = new Triangulation(); + dof = new DoFHandler (tria); + set_tria_and_dof (tria, dof); +}; + + + + + + +template +void PoissonProblem::run (const unsigned int level) { + create_new (); + + cout << "Refinement level = " << level + << endl; + + cout << " Making grid... "; + tria->create_hypercube (); + tria->refine_global (level); + cout << tria->n_active_cells() << " active cells." << endl; + + rhs = new RHSPoly(); + boundary_values = new ZeroFunction (); + + + FELinear fe; + PoissonEquation equation (*rhs); + QGauss3 quadrature; + + cout << " Distributing dofs... "; + dof->distribute_dofs (fe); + cout << dof->n_dofs() << " degrees of freedom." << endl; + + cout << " Assembling matrices..." << endl; + FEValues::UpdateStruct update_flags; + update_flags.q_points = update_flags.gradients = true; + update_flags.jacobians = update_flags.JxW_values = true; + + ProblemBase::DirichletBC dirichlet_bc; + dirichlet_bc[0] = boundary_values; + assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + + cout << " Solving..." << endl; + solve (); + + Solution sol; + dVector l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell; + QGauss4 q; + + cout << " Calculating L1 error... "; + integrate_difference (sol, l1_error_per_cell, q, fe, L1_norm); + cout << l1_error_per_cell.l1_norm() << endl; + + cout << " Calculating L2 error... "; + integrate_difference (sol, l2_error_per_cell, q, fe, L2_norm); + cout << l2_error_per_cell.l2_norm() << endl; + + cout << " Calculating L-infinity error... "; + integrate_difference (sol, linfty_error_per_cell, q, fe, Linfty_norm); + cout << linfty_error_per_cell.linfty_norm() << endl; + + dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; + dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof); + dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); + dof->distribute_cell_to_dof_vector (linfty_error_per_cell, linfty_error_per_dof); + + string filename = "gnuplot."; + filename += ('0'+level); + cout << " Writing error plots to <" << filename << ">..." << endl; + + DataOut out; + ofstream gnuplot(filename.c_str()); + fill_data (out); + out.add_data_vector (l1_error_per_dof, "L1-Error"); + out.add_data_vector (l2_error_per_dof, "L2-Error"); + out.add_data_vector (linfty_error_per_dof, "L3-Error"); + out.write_gnuplot (gnuplot); + gnuplot.close (); + + cout << endl; +}; + + + + + +int main () { + PoissonProblem<2> problem; + + for (unsigned int level=1; level<5; ++level) + problem.run (level); + + return 0; +}; diff --git a/tests/big-tests/convergence/convergence.cc b/tests/big-tests/convergence/convergence.cc new file mode 100644 index 0000000000..331afc5baf --- /dev/null +++ b/tests/big-tests/convergence/convergence.cc @@ -0,0 +1,306 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include +extern "C" { +# include +} + + + + +template +class PoissonEquation : public Equation { + public: + PoissonEquation (const Function &rhs) : + Equation(1), + right_hand_side (rhs) {}; + + virtual void assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + protected: + const Function &right_hand_side; +}; + + + + + + +template +class PoissonProblem : public ProblemBase { + public: + PoissonProblem (); + + void clear (); + virtual void create_new (); + virtual void run (unsigned int level); + + protected: + Triangulation *tria; + DoFHandler *dof; + + Function *rhs; + Function *boundary_values; +}; + + + + + +/** + Right hand side constructed such that the exact solution is + $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc. + */ +template +class RHSPoly : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + +template +class Solution : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + + +template +double RHSPoly::operator () (const Point &p) const { + double ret_val = 0; + for (unsigned int i=0; i +double Solution::operator () (const Point &p) const { + double ret_val = 1; + for (unsigned int i=0; i::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point +void PoissonEquation::assemble (dFMatrix &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + +template +void PoissonEquation::assemble (dVector &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + + + + + + + +template +PoissonProblem::PoissonProblem () : + tria(0), dof(0), rhs(0), boundary_values(0) {}; + + + + +template +void PoissonProblem::clear () { + if (tria != 0) { + delete tria; + tria = 0; + }; + + if (dof != 0) { + delete dof; + dof = 0; + }; + + if (rhs != 0) + { + delete rhs; + rhs = 0; + }; + + if (boundary_values != 0) + { + delete boundary_values; + boundary_values = 0; + }; + + ProblemBase::clear (); +}; + + + + +template +void PoissonProblem::create_new () { + clear (); + + tria = new Triangulation(); + dof = new DoFHandler (tria); + set_tria_and_dof (tria, dof); +}; + + + + + + +template +void PoissonProblem::run (const unsigned int level) { + create_new (); + + cout << "Refinement level = " << level + << endl; + + cout << " Making grid... "; + tria->create_hypercube (); + tria->refine_global (level); + cout << tria->n_active_cells() << " active cells." << endl; + + rhs = new RHSPoly(); + boundary_values = new ZeroFunction (); + + + FELinear fe; + PoissonEquation equation (*rhs); + QGauss3 quadrature; + + cout << " Distributing dofs... "; + dof->distribute_dofs (fe); + cout << dof->n_dofs() << " degrees of freedom." << endl; + + cout << " Assembling matrices..." << endl; + FEValues::UpdateStruct update_flags; + update_flags.q_points = update_flags.gradients = true; + update_flags.jacobians = update_flags.JxW_values = true; + + ProblemBase::DirichletBC dirichlet_bc; + dirichlet_bc[0] = boundary_values; + assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + + cout << " Solving..." << endl; + solve (); + + Solution sol; + dVector l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell; + QGauss4 q; + + cout << " Calculating L1 error... "; + integrate_difference (sol, l1_error_per_cell, q, fe, L1_norm); + cout << l1_error_per_cell.l1_norm() << endl; + + cout << " Calculating L2 error... "; + integrate_difference (sol, l2_error_per_cell, q, fe, L2_norm); + cout << l2_error_per_cell.l2_norm() << endl; + + cout << " Calculating L-infinity error... "; + integrate_difference (sol, linfty_error_per_cell, q, fe, Linfty_norm); + cout << linfty_error_per_cell.linfty_norm() << endl; + + dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; + dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof); + dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); + dof->distribute_cell_to_dof_vector (linfty_error_per_cell, linfty_error_per_dof); + + string filename = "gnuplot."; + filename += ('0'+level); + cout << " Writing error plots to <" << filename << ">..." << endl; + + DataOut out; + ofstream gnuplot(filename.c_str()); + fill_data (out); + out.add_data_vector (l1_error_per_dof, "L1-Error"); + out.add_data_vector (l2_error_per_dof, "L2-Error"); + out.add_data_vector (linfty_error_per_dof, "L3-Error"); + out.write_gnuplot (gnuplot); + gnuplot.close (); + + cout << endl; +}; + + + + + +int main () { + PoissonProblem<2> problem; + + for (unsigned int level=1; level<5; ++level) + problem.run (level); + + return 0; +}; -- 2.39.5