From: Wolfgang Bangerth Date: Tue, 12 May 1998 13:34:11 +0000 (+0000) Subject: Extend example by a function with a kink (discontinuity of the gradient). X-Git-Tag: v8.0.0~22973 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=99f49d6a247f6774e75deaf275f45ca1c979d8bd;p=dealii.git Extend example by a function with a kink (discontinuity of the gradient). git-svn-id: https://svn.dealii.org/trunk@281 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/error-estimation/ee.kink.prm b/deal.II/deal.II/Attic/examples/error-estimation/ee.kink.prm new file mode 100644 index 0000000000..5a3c5447d8 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/error-estimation/ee.kink.prm @@ -0,0 +1,7 @@ +set Test case = Kink +set Initial refinement = 1 +set Refinement criterion = { global | estimated error } +set Refinement fraction = 0.1 +set Maximum cells = 30000 +set Output base filename = data-kink/ +set Output format = ucd 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 index 18c19ef5d5..6f055f91e2 100644 --- a/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc +++ b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc @@ -30,7 +30,16 @@ class PoissonEquation : public Equation { public: PoissonEquation (const Function &rhs) : Equation(1), - right_hand_side (rhs) {}; + right_hand_side (rhs), + coefficient (default_coefficient), + use_coefficient(false) {}; + + PoissonEquation (const Function &rhs, + const Function &coefficient ) : + Equation(1), + right_hand_side (rhs), + coefficient (coefficient), + use_coefficient(true) {}; virtual void assemble (dFMatrix &cell_matrix, dVector &rhs, @@ -43,24 +52,29 @@ class PoissonEquation : public Equation { const FEValues &fe_values, const Triangulation::cell_iterator &cell) const; protected: + const bool use_coefficient; const Function &right_hand_side; + const Function &coefficient; + + static const ConstantFunction default_coefficient; }; +const ConstantFunction<2> PoissonEquation<2>::default_coefficient(1); + template -class PoissonProblem : public ProblemBase, - public MultipleParameterLoop::UserClass { +class PoissonProblem : public ProblemBase, public MultipleParameterLoop::UserClass { public: enum RefineMode { global, true_error, error_estimator }; PoissonProblem (); - + void clear (); void create_new (const unsigned int); void declare_parameters (ParameterHandler &prm); @@ -74,7 +88,8 @@ class PoissonProblem : public ProblemBase, Function *rhs; Function *solution_function; - + Function *coefficient; + Boundary *boundary; vector l2_error, linfty_error; @@ -101,6 +116,17 @@ class Solution { virtual double operator () (const Point &p) const; virtual Point gradient (const Point &p) const; }; + + class Kink : public Function { + public: + class Coefficient : public Function { + public: + virtual double operator () (const Point &p) const; + }; + + virtual double operator () (const Point &p) const; + virtual Point gradient (const Point &p) const; + }; }; @@ -123,11 +149,22 @@ class RHS { /** * Right hand side constructed such that * the exact solution is - * $r^{2/3} sin(2\phi)$. + * $r^{2/3}$. */ class Singular : public Function { public: virtual double operator () (const Point &p) const; + }; + + /** + * Right hand side constructed such that + * the exact solution is + * $(1+4\theta(f))*f$ with + * $f=y-x**2$. + */ + class Kink : public Function { + public: + virtual double operator () (const Point &p) const; }; }; @@ -158,17 +195,46 @@ Point<2> Solution<2>::Singular::gradient (const Point<2> &p) const { + +inline double theta(const double x) { + return (x>0 ? 1 : 0); +}; + + +double Solution<2>::Kink::operator () (const Point<2> &p) const { + const double s = p(1)-p(0)*p(0); + return (1+4*theta(s))*s; +}; + + +Point<2> Solution<2>::Kink::gradient (const Point<2> &p) const { + const double s = p(1)-p(0)*p(0); + return (1+4*theta(s))*Point<2>(-2*p(0),1); +}; + + +double Solution<2>::Kink::Coefficient::operator () (const Point<2> &p) const { + const double s = p(1)-p(0)*p(0); + return 1./(1.+4.*theta(s)); +}; + + + double RHS<2>::GaussShape::operator () (const Point<2> &p) const { return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); }; - double RHS<2>::Singular::operator () (const Point<2> &p) const { return -4./9. * pow(p.square(), -2./3.); }; +double RHS<2>::Kink::operator () (const Point<2> &) const { + return 2; +}; + + @@ -180,17 +246,23 @@ 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::assemble (dVector &, template PoissonProblem::PoissonProblem () : tria(0), dof(0), rhs(0), - solution_function(0), boundary(0) {}; + solution_function(0), coefficient(0), + boundary(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 (solution_function != 0) - { - delete solution_function; - solution_function = 0; - }; - - if (boundary != 0) - { - delete boundary; - boundary = 0; - }; + if (tria != 0) { delete tria; tria = 0; }; + if (dof != 0) { delete dof; dof = 0; }; + if (rhs != 0) { delete rhs; rhs = 0; }; + if (solution_function != 0) { delete solution_function; solution_function = 0; }; + if (coefficient != 0) { delete coefficient; coefficient = 0; }; + if (boundary != 0) { delete boundary; boundary = 0; }; l2_error.clear (); linfty_error.clear (); @@ -283,7 +335,8 @@ void PoissonProblem::create_new (const unsigned int) { template void PoissonProblem::declare_parameters (ParameterHandler &prm) { - prm.declare_entry ("Test case", "Gauss shape", "Gauss shape\\|Singular"); + prm.declare_entry ("Test case", "Gauss shape", + "Gauss shape\\|Singular\\|Kink"); prm.declare_entry ("Initial refinement", "2", ParameterHandler::RegularExpressions::Integer); prm.declare_entry ("Refinement criterion", "estimated error", @@ -345,17 +398,29 @@ void PoissonProblem::run (ParameterHandler &prm) { else if (prm.get("Test case")=="Singular") rhs = new RHS::Singular(); + else + if (prm.get("Test case")=="Kink") + rhs = new RHS::Kink(); if (prm.get("Test case")=="Gauss shape") solution_function = new Solution::GaussShape (); else if (prm.get("Test case")=="Singular") solution_function = new Solution::Singular (); + else + if (prm.get("Test case")=="Kink") + solution_function = new Solution::Kink (); - FELinear fe; - PoissonEquation equation (*rhs); - QGauss3 quadrature; + FELinear fe; + QGauss3 quadrature; + PoissonEquation *equation; + + static Solution::Kink::Coefficient kink_coefficient; + if (prm.get("Test case")=="Kink") + equation = new PoissonEquation(*rhs, kink_coefficient); + else + equation = new PoissonEquation(*rhs); unsigned int refine_step = 0; const unsigned int max_cells = prm.get_integer("Maximum cells"); @@ -376,7 +441,7 @@ void PoissonProblem::run (ParameterHandler &prm) { ProblemBase::FunctionMap dirichlet_bc; dirichlet_bc[0] = solution_function; - assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + assemble (*equation, quadrature, fe, update_flags, dirichlet_bc); cout << " Solving..." << endl; solve (); @@ -485,6 +550,8 @@ void PoissonProblem::run (ParameterHandler &prm) { print_history (prm, refine_mode); cout << endl << endl << endl; + + delete equation; }; diff --git a/deal.II/deal.II/Attic/examples/error-estimation/make_ps b/deal.II/deal.II/Attic/examples/error-estimation/make_ps index 38d05f569a..06340e21a4 100644 --- a/deal.II/deal.II/Attic/examples/error-estimation/make_ps +++ b/deal.II/deal.II/Attic/examples/error-estimation/make_ps @@ -41,3 +41,22 @@ plot "data-singular/history.estimated_error.gnuplot" using 1:2 title "L2 error", set output "data-singular/history.compare.eps" plot "data-singular/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-singular/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "data-singular/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-singular/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" + + + + + + + +set output "data-kink/history.global.eps" + +plot "data-kink/history.global.gnuplot" using 1:2 title "L2 error","data-kink/history.global.gnuplot" using 1:3 title "Linfty error","data-kink/history.global.gnuplot" using 1:4 title "H1 error","data-kink/history.global.gnuplot" using 1:5 title "Estimated H1 error" + + +set output "data-kink/history.estimated_error.eps" + +plot "data-kink/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-kink/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-kink/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-kink/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" + + +set output "data-kink/history.compare.eps" +plot "data-kink/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-kink/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "data-kink/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-kink/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" diff --git a/tests/big-tests/error-estimation/ee.kink.prm b/tests/big-tests/error-estimation/ee.kink.prm new file mode 100644 index 0000000000..5a3c5447d8 --- /dev/null +++ b/tests/big-tests/error-estimation/ee.kink.prm @@ -0,0 +1,7 @@ +set Test case = Kink +set Initial refinement = 1 +set Refinement criterion = { global | estimated error } +set Refinement fraction = 0.1 +set Maximum cells = 30000 +set Output base filename = data-kink/ +set Output format = ucd diff --git a/tests/big-tests/error-estimation/error-estimation.cc b/tests/big-tests/error-estimation/error-estimation.cc index 18c19ef5d5..6f055f91e2 100644 --- a/tests/big-tests/error-estimation/error-estimation.cc +++ b/tests/big-tests/error-estimation/error-estimation.cc @@ -30,7 +30,16 @@ class PoissonEquation : public Equation { public: PoissonEquation (const Function &rhs) : Equation(1), - right_hand_side (rhs) {}; + right_hand_side (rhs), + coefficient (default_coefficient), + use_coefficient(false) {}; + + PoissonEquation (const Function &rhs, + const Function &coefficient ) : + Equation(1), + right_hand_side (rhs), + coefficient (coefficient), + use_coefficient(true) {}; virtual void assemble (dFMatrix &cell_matrix, dVector &rhs, @@ -43,24 +52,29 @@ class PoissonEquation : public Equation { const FEValues &fe_values, const Triangulation::cell_iterator &cell) const; protected: + const bool use_coefficient; const Function &right_hand_side; + const Function &coefficient; + + static const ConstantFunction default_coefficient; }; +const ConstantFunction<2> PoissonEquation<2>::default_coefficient(1); + template -class PoissonProblem : public ProblemBase, - public MultipleParameterLoop::UserClass { +class PoissonProblem : public ProblemBase, public MultipleParameterLoop::UserClass { public: enum RefineMode { global, true_error, error_estimator }; PoissonProblem (); - + void clear (); void create_new (const unsigned int); void declare_parameters (ParameterHandler &prm); @@ -74,7 +88,8 @@ class PoissonProblem : public ProblemBase, Function *rhs; Function *solution_function; - + Function *coefficient; + Boundary *boundary; vector l2_error, linfty_error; @@ -101,6 +116,17 @@ class Solution { virtual double operator () (const Point &p) const; virtual Point gradient (const Point &p) const; }; + + class Kink : public Function { + public: + class Coefficient : public Function { + public: + virtual double operator () (const Point &p) const; + }; + + virtual double operator () (const Point &p) const; + virtual Point gradient (const Point &p) const; + }; }; @@ -123,11 +149,22 @@ class RHS { /** * Right hand side constructed such that * the exact solution is - * $r^{2/3} sin(2\phi)$. + * $r^{2/3}$. */ class Singular : public Function { public: virtual double operator () (const Point &p) const; + }; + + /** + * Right hand side constructed such that + * the exact solution is + * $(1+4\theta(f))*f$ with + * $f=y-x**2$. + */ + class Kink : public Function { + public: + virtual double operator () (const Point &p) const; }; }; @@ -158,17 +195,46 @@ Point<2> Solution<2>::Singular::gradient (const Point<2> &p) const { + +inline double theta(const double x) { + return (x>0 ? 1 : 0); +}; + + +double Solution<2>::Kink::operator () (const Point<2> &p) const { + const double s = p(1)-p(0)*p(0); + return (1+4*theta(s))*s; +}; + + +Point<2> Solution<2>::Kink::gradient (const Point<2> &p) const { + const double s = p(1)-p(0)*p(0); + return (1+4*theta(s))*Point<2>(-2*p(0),1); +}; + + +double Solution<2>::Kink::Coefficient::operator () (const Point<2> &p) const { + const double s = p(1)-p(0)*p(0); + return 1./(1.+4.*theta(s)); +}; + + + double RHS<2>::GaussShape::operator () (const Point<2> &p) const { return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); }; - double RHS<2>::Singular::operator () (const Point<2> &p) const { return -4./9. * pow(p.square(), -2./3.); }; +double RHS<2>::Kink::operator () (const Point<2> &) const { + return 2; +}; + + @@ -180,17 +246,23 @@ 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::assemble (dVector &, template PoissonProblem::PoissonProblem () : tria(0), dof(0), rhs(0), - solution_function(0), boundary(0) {}; + solution_function(0), coefficient(0), + boundary(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 (solution_function != 0) - { - delete solution_function; - solution_function = 0; - }; - - if (boundary != 0) - { - delete boundary; - boundary = 0; - }; + if (tria != 0) { delete tria; tria = 0; }; + if (dof != 0) { delete dof; dof = 0; }; + if (rhs != 0) { delete rhs; rhs = 0; }; + if (solution_function != 0) { delete solution_function; solution_function = 0; }; + if (coefficient != 0) { delete coefficient; coefficient = 0; }; + if (boundary != 0) { delete boundary; boundary = 0; }; l2_error.clear (); linfty_error.clear (); @@ -283,7 +335,8 @@ void PoissonProblem::create_new (const unsigned int) { template void PoissonProblem::declare_parameters (ParameterHandler &prm) { - prm.declare_entry ("Test case", "Gauss shape", "Gauss shape\\|Singular"); + prm.declare_entry ("Test case", "Gauss shape", + "Gauss shape\\|Singular\\|Kink"); prm.declare_entry ("Initial refinement", "2", ParameterHandler::RegularExpressions::Integer); prm.declare_entry ("Refinement criterion", "estimated error", @@ -345,17 +398,29 @@ void PoissonProblem::run (ParameterHandler &prm) { else if (prm.get("Test case")=="Singular") rhs = new RHS::Singular(); + else + if (prm.get("Test case")=="Kink") + rhs = new RHS::Kink(); if (prm.get("Test case")=="Gauss shape") solution_function = new Solution::GaussShape (); else if (prm.get("Test case")=="Singular") solution_function = new Solution::Singular (); + else + if (prm.get("Test case")=="Kink") + solution_function = new Solution::Kink (); - FELinear fe; - PoissonEquation equation (*rhs); - QGauss3 quadrature; + FELinear fe; + QGauss3 quadrature; + PoissonEquation *equation; + + static Solution::Kink::Coefficient kink_coefficient; + if (prm.get("Test case")=="Kink") + equation = new PoissonEquation(*rhs, kink_coefficient); + else + equation = new PoissonEquation(*rhs); unsigned int refine_step = 0; const unsigned int max_cells = prm.get_integer("Maximum cells"); @@ -376,7 +441,7 @@ void PoissonProblem::run (ParameterHandler &prm) { ProblemBase::FunctionMap dirichlet_bc; dirichlet_bc[0] = solution_function; - assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + assemble (*equation, quadrature, fe, update_flags, dirichlet_bc); cout << " Solving..." << endl; solve (); @@ -485,6 +550,8 @@ void PoissonProblem::run (ParameterHandler &prm) { print_history (prm, refine_mode); cout << endl << endl << endl; + + delete equation; }; diff --git a/tests/big-tests/error-estimation/make_ps b/tests/big-tests/error-estimation/make_ps index 38d05f569a..06340e21a4 100644 --- a/tests/big-tests/error-estimation/make_ps +++ b/tests/big-tests/error-estimation/make_ps @@ -41,3 +41,22 @@ plot "data-singular/history.estimated_error.gnuplot" using 1:2 title "L2 error", set output "data-singular/history.compare.eps" plot "data-singular/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-singular/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "data-singular/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-singular/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" + + + + + + + +set output "data-kink/history.global.eps" + +plot "data-kink/history.global.gnuplot" using 1:2 title "L2 error","data-kink/history.global.gnuplot" using 1:3 title "Linfty error","data-kink/history.global.gnuplot" using 1:4 title "H1 error","data-kink/history.global.gnuplot" using 1:5 title "Estimated H1 error" + + +set output "data-kink/history.estimated_error.eps" + +plot "data-kink/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-kink/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-kink/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-kink/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" + + +set output "data-kink/history.compare.eps" +plot "data-kink/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-kink/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "data-kink/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-kink/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error"