From 995b7124820bb87595633a4e80b574874aba3b6c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 24 Apr 1998 08:12:43 +0000 Subject: [PATCH] Add new test case (an interesting one). git-svn-id: https://svn.dealii.org/trunk@189 0785d39b-7218-0410-832d-ea1e28bc413d --- .../error-estimation/error-estimation.cc | 382 ++++++++++++++++++ .../error-estimation/error-estimation.cc | 382 ++++++++++++++++++ 2 files changed, 764 insertions(+) create mode 100644 deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc create mode 100644 tests/big-tests/error-estimation/error-estimation.cc diff --git a/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc new file mode 100644 index 0000000000..5ea362273f --- /dev/null +++ b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc @@ -0,0 +1,382 @@ +/* $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 (); + void create_new (); + void run (unsigned int level); + void print_history () const; + + protected: + Triangulation *tria; + DoFHandler *dof; + + Function *rhs; + Function *boundary_values; + + vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; + vector n_dofs; +}; + + + + + +/** + Right hand side constructed such that the exact solution is + $x*y*exp(-(x**2+y**2)*10)$. + */ +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; + /** + * Return the gradient of the function + * at the given point. + */ + virtual Point gradient (const Point &p) const; +}; + + + + +double RHSPoly<2>::operator () (const Point<2> &p) const { + return (120.-400.*p.square())*p(0)*p(1)*exp(-10.*p.square()); +}; + + + +double Solution<2>::operator () (const Point<2> &p) const { + return p(0)*p(1)*exp(-10*p.square()); +}; + + +Point<2> Solution<2>::gradient (const Point<2> &p) const { + return Point<2> ((1-20.*p(0)*p(0))*p(1)*exp(-10*p.square()), + (1-20.*p(1)*p(1))*p(0)*exp(-10*p.square())); +}; + + + + + + +void PoissonEquation<2>::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 (-1,+1); + tria->refine_global (level); + cout << tria->n_active_cells() << " active cells." << endl; + + rhs = new RHSPoly(); + boundary_values = new Solution (); + + + FELinear fe; + PoissonEquation equation (*rhs); + QGauss3 quadrature; + + cout << " Distributing dofs... "; + dof->distribute_dofs (fe); + cout << dof->n_dofs() << " degrees of freedom." << endl; + n_dofs.push_back (dof->n_dofs()); + + cout << " Assembling matrices..." << endl; + UpdateFields update_flags = UpdateFields(update_q_points | update_gradients | + update_jacobians | update_JxW_values); + + 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; + dVector h1_seminorm_error_per_cell, h1_error_per_cell; + QGauss3 q; + + cout << " Calculating L1 error... "; + integrate_difference (sol, l1_error_per_cell, q, fe, L1_norm); + cout << l1_error_per_cell.l1_norm() << endl; + l1_error.push_back (l1_error_per_cell.l1_norm()); + + cout << " Calculating L2 error... "; + integrate_difference (sol, l2_error_per_cell, q, fe, L2_norm); + cout << l2_error_per_cell.l2_norm() << endl; + l2_error.push_back (l2_error_per_cell.l2_norm()); + + cout << " Calculating L-infinity error... "; + integrate_difference (sol, linfty_error_per_cell, q, fe, Linfty_norm); + cout << linfty_error_per_cell.linfty_norm() << endl; + linfty_error.push_back (linfty_error_per_cell.linfty_norm()); + + cout << " Calculating H1-seminorm error... "; + integrate_difference (sol, h1_seminorm_error_per_cell, q, fe, H1_seminorm); + cout << h1_seminorm_error_per_cell.l2_norm() << endl; + h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); + + cout << " Calculating H1 error... "; + integrate_difference (sol, h1_error_per_cell, q, fe, H1_norm); + cout << h1_error_per_cell.l2_norm() << endl; + h1_error.push_back (h1_error_per_cell.l2_norm()); + + if (level<=5) + { + dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; + dVector h1_seminorm_error_per_dof, h1_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); + dof->distribute_cell_to_dof_vector (h1_seminorm_error_per_cell, h1_seminorm_error_per_dof); + dof->distribute_cell_to_dof_vector (h1_error_per_cell, h1_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, "Linfty-Error"); + out.add_data_vector (h1_seminorm_error_per_dof, "H1-seminorm-Error"); + out.add_data_vector (h1_error_per_dof, "H1-Error"); + out.write_gnuplot (gnuplot); + gnuplot.close (); + } + else + cout << " Not writing error as grid." << endl; + + cout << endl; +}; + + +template +void PoissonProblem::print_history () const { + ofstream out("gnuplot.history"); + out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" + << endl; + for (unsigned int i=0; ih/2:" << endl; + cout << " L1 error : " << 1./average_l1 << endl + << " L2 error : " << 1./average_l2 << endl + << " Linfty error : " << 1./average_linfty << endl + << " H1 seminorm error: " << 1./average_h1_semi << endl + << " H1 error : " << 1./average_h1 << endl; +}; + + + + +int main () { + PoissonProblem<2> problem; + + for (unsigned int level=3; level<7; ++level) + problem.run (level); + + cout << endl << "Printing convergence history to ..." << endl; + problem.print_history (); + + return 0; +}; diff --git a/tests/big-tests/error-estimation/error-estimation.cc b/tests/big-tests/error-estimation/error-estimation.cc new file mode 100644 index 0000000000..5ea362273f --- /dev/null +++ b/tests/big-tests/error-estimation/error-estimation.cc @@ -0,0 +1,382 @@ +/* $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 (); + void create_new (); + void run (unsigned int level); + void print_history () const; + + protected: + Triangulation *tria; + DoFHandler *dof; + + Function *rhs; + Function *boundary_values; + + vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; + vector n_dofs; +}; + + + + + +/** + Right hand side constructed such that the exact solution is + $x*y*exp(-(x**2+y**2)*10)$. + */ +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; + /** + * Return the gradient of the function + * at the given point. + */ + virtual Point gradient (const Point &p) const; +}; + + + + +double RHSPoly<2>::operator () (const Point<2> &p) const { + return (120.-400.*p.square())*p(0)*p(1)*exp(-10.*p.square()); +}; + + + +double Solution<2>::operator () (const Point<2> &p) const { + return p(0)*p(1)*exp(-10*p.square()); +}; + + +Point<2> Solution<2>::gradient (const Point<2> &p) const { + return Point<2> ((1-20.*p(0)*p(0))*p(1)*exp(-10*p.square()), + (1-20.*p(1)*p(1))*p(0)*exp(-10*p.square())); +}; + + + + + + +void PoissonEquation<2>::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 (-1,+1); + tria->refine_global (level); + cout << tria->n_active_cells() << " active cells." << endl; + + rhs = new RHSPoly(); + boundary_values = new Solution (); + + + FELinear fe; + PoissonEquation equation (*rhs); + QGauss3 quadrature; + + cout << " Distributing dofs... "; + dof->distribute_dofs (fe); + cout << dof->n_dofs() << " degrees of freedom." << endl; + n_dofs.push_back (dof->n_dofs()); + + cout << " Assembling matrices..." << endl; + UpdateFields update_flags = UpdateFields(update_q_points | update_gradients | + update_jacobians | update_JxW_values); + + 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; + dVector h1_seminorm_error_per_cell, h1_error_per_cell; + QGauss3 q; + + cout << " Calculating L1 error... "; + integrate_difference (sol, l1_error_per_cell, q, fe, L1_norm); + cout << l1_error_per_cell.l1_norm() << endl; + l1_error.push_back (l1_error_per_cell.l1_norm()); + + cout << " Calculating L2 error... "; + integrate_difference (sol, l2_error_per_cell, q, fe, L2_norm); + cout << l2_error_per_cell.l2_norm() << endl; + l2_error.push_back (l2_error_per_cell.l2_norm()); + + cout << " Calculating L-infinity error... "; + integrate_difference (sol, linfty_error_per_cell, q, fe, Linfty_norm); + cout << linfty_error_per_cell.linfty_norm() << endl; + linfty_error.push_back (linfty_error_per_cell.linfty_norm()); + + cout << " Calculating H1-seminorm error... "; + integrate_difference (sol, h1_seminorm_error_per_cell, q, fe, H1_seminorm); + cout << h1_seminorm_error_per_cell.l2_norm() << endl; + h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); + + cout << " Calculating H1 error... "; + integrate_difference (sol, h1_error_per_cell, q, fe, H1_norm); + cout << h1_error_per_cell.l2_norm() << endl; + h1_error.push_back (h1_error_per_cell.l2_norm()); + + if (level<=5) + { + dVector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; + dVector h1_seminorm_error_per_dof, h1_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); + dof->distribute_cell_to_dof_vector (h1_seminorm_error_per_cell, h1_seminorm_error_per_dof); + dof->distribute_cell_to_dof_vector (h1_error_per_cell, h1_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, "Linfty-Error"); + out.add_data_vector (h1_seminorm_error_per_dof, "H1-seminorm-Error"); + out.add_data_vector (h1_error_per_dof, "H1-Error"); + out.write_gnuplot (gnuplot); + gnuplot.close (); + } + else + cout << " Not writing error as grid." << endl; + + cout << endl; +}; + + +template +void PoissonProblem::print_history () const { + ofstream out("gnuplot.history"); + out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" + << endl; + for (unsigned int i=0; ih/2:" << endl; + cout << " L1 error : " << 1./average_l1 << endl + << " L2 error : " << 1./average_l2 << endl + << " Linfty error : " << 1./average_linfty << endl + << " H1 seminorm error: " << 1./average_h1_semi << endl + << " H1 error : " << 1./average_h1 << endl; +}; + + + + +int main () { + PoissonProblem<2> problem; + + for (unsigned int level=3; level<7; ++level) + problem.run (level); + + cout << endl << "Printing convergence history to ..." << endl; + problem.print_history (); + + return 0; +}; -- 2.39.5