From 181fdf206e9d2969654308ad49af086163e4155c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 11 May 1998 15:51:55 +0000 Subject: [PATCH] Example update: use ParameterHandler class and other changes. git-svn-id: https://svn.dealii.org/trunk@275 0785d39b-7218-0410-832d-ea1e28bc413d --- .../Attic/examples/error-estimation/Makefile | 6 +- .../examples/error-estimation/ee.gauss.prm | 6 + .../examples/error-estimation/ee.singular.prm | 6 + .../error-estimation/error-estimation.cc | 249 ++++++++++++------ .../Attic/examples/error-estimation/make_ps | 39 ++- .../examples/error-estimation/strip_comments | 2 +- tests/big-tests/error-estimation/Makefile | 6 +- tests/big-tests/error-estimation/ee.gauss.prm | 6 + .../error-estimation/ee.singular.prm | 6 + .../error-estimation/error-estimation.cc | 249 ++++++++++++------ tests/big-tests/error-estimation/make_ps | 39 ++- .../big-tests/error-estimation/strip_comments | 2 +- 12 files changed, 438 insertions(+), 178 deletions(-) create mode 100644 deal.II/deal.II/Attic/examples/error-estimation/ee.gauss.prm create mode 100644 deal.II/deal.II/Attic/examples/error-estimation/ee.singular.prm create mode 100644 tests/big-tests/error-estimation/ee.gauss.prm create mode 100644 tests/big-tests/error-estimation/ee.singular.prm diff --git a/deal.II/deal.II/Attic/examples/error-estimation/Makefile b/deal.II/deal.II/Attic/examples/error-estimation/Makefile index 873179dc03..3f035617cb 100644 --- a/deal.II/deal.II/Attic/examples/error-estimation/Makefile +++ b/deal.II/deal.II/Attic/examples/error-estimation/Makefile @@ -19,16 +19,16 @@ user-libs = ../../../mia/control.o # To run the program, use "make run"; to give parameters to the program, # give the parameters to the following variable: -run-parameters = +run-parameters = ee.gauss.prm # To execute additional action apart from running the program, fill # in this list: -additional-run-action = gnuplot make_ps ; ./strip_comments +additional-run-action = $(target) ee.singular.prm ; gnuplot make_ps ; ./strip_comments # To specify which files are to be deleted by "make clean" (apart from # the usual ones: object files, executables, backups, etc), fill in the # following list -delete-files = gnuplot* *inp *.ps *.eps history.* +delete-files = data-*/* diff --git a/deal.II/deal.II/Attic/examples/error-estimation/ee.gauss.prm b/deal.II/deal.II/Attic/examples/error-estimation/ee.gauss.prm new file mode 100644 index 0000000000..588bb9a327 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/error-estimation/ee.gauss.prm @@ -0,0 +1,6 @@ +set Test case = Gauss shape +set Initial refinement = 2 +set Refinement criterion = { global | true error | estimated error } +set Refinement fraction = 0.3 +set Output base filename = data-gauss/ +set Output format = ucd diff --git a/deal.II/deal.II/Attic/examples/error-estimation/ee.singular.prm b/deal.II/deal.II/Attic/examples/error-estimation/ee.singular.prm new file mode 100644 index 0000000000..d906f3f768 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/error-estimation/ee.singular.prm @@ -0,0 +1,6 @@ +set Test case = Singular +set Initial refinement = 1 +set Refinement criterion = { global | true error | estimated error } +set Refinement fraction = 0.2 +set Output base filename = data-singular/ +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 4566b60ab5..2b42b1b95a 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 @@ -9,20 +9,18 @@ #include #include #include +#include #include #include #include #include #include - #include #include #include #include -extern "C" { -# include -} +#include @@ -54,7 +52,8 @@ class PoissonEquation : public Equation { template -class PoissonProblem : public ProblemBase { +class PoissonProblem : public ProblemBase, + public MultipleParameterLoop::UserClass { public: enum RefineMode { global, true_error, error_estimator @@ -63,19 +62,20 @@ class PoissonProblem : public ProblemBase { PoissonProblem (); void clear (); - void create_new (); - void run (const unsigned int start_level, - const RefineMode refine_mode); - void print_history (const RefineMode refine_mode) const; + void create_new (const unsigned int); + void declare_parameters (ParameterHandler &prm); + void run (ParameterHandler &prm); + void print_history (const ParameterHandler &prm, + const RefineMode refine_mode) const; protected: Triangulation *tria; DoFHandler *dof; Function *rhs; - Function *boundary_values; + Function *solution_function; - HyperBallBoundary boundary; + Boundary *boundary; vector l2_error, linfty_error; vector h1_error, estimated_error; @@ -88,56 +88,91 @@ class PoissonProblem : public ProblemBase { -/** - Right hand side constructed such that the exact solution is - $x*y*exp(-(x**2+y**2)*40)$. - */ template -class RHSPoly : public Function { +class Solution { public: - /** - * Return the value of the function - * at the given point. - */ - virtual double operator () (const Point &p) const; + + class GaussShape : public Function { + public: + virtual double operator () (const Point &p) const; + virtual Point gradient (const Point &p) const; + }; + + class Singular : public Function { + public: + virtual double operator () (const Point &p) const; + virtual Point gradient (const Point &p) const; + }; }; + template -class Solution : public Function { +class RHS { 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. + * Right hand side constructed such that + * the exact solution is + * $x*y*exp(-(x**2+y**2)*40)$. */ - virtual Point gradient (const Point &p) const; -}; - + class GaussShape : public Function { + public: + virtual double operator () (const Point &p) const; + }; + /** + * Right hand side constructed such that + * the exact solution is + * $r^{2/3} sin(2\phi)$. + */ + class Singular : public Function { + public: + virtual double operator () (const Point &p) const; + }; +}; -double RHSPoly<2>::operator () (const Point<2> &p) const { - return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); -}; -double Solution<2>::operator () (const Point<2> &p) const { +double Solution<2>::GaussShape::operator () (const Point<2> &p) const { return p(0)*p(1)*exp(-40*p.square()); }; -Point<2> Solution<2>::gradient (const Point<2> &p) const { +Point<2> Solution<2>::GaussShape::gradient (const Point<2> &p) const { return Point<2> ((1-80.*p(0)*p(0))*p(1)*exp(-40*p.square()), (1-80.*p(1)*p(1))*p(0)*exp(-40*p.square())); }; + + +double Solution<2>::Singular::operator () (const Point<2> &p) const { + return pow(p.square(), 1./3.); +}; + + +Point<2> Solution<2>::Singular::gradient (const Point<2> &p) const { + return 2./3.*pow(p.square(), -2./3.) * p; +}; + + + +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.); +}; + + + + @@ -189,7 +224,7 @@ void PoissonEquation::assemble (dVector &, template PoissonProblem::PoissonProblem () : tria(0), dof(0), rhs(0), - boundary_values(0), boundary() {}; + solution_function(0), boundary(0) {}; @@ -212,12 +247,18 @@ void PoissonProblem::clear () { rhs = 0; }; - if (boundary_values != 0) + if (solution_function != 0) { - delete boundary_values; - boundary_values = 0; + delete solution_function; + solution_function = 0; }; + if (boundary != 0) + { + delete boundary; + boundary = 0; + }; + l2_error.clear (); linfty_error.clear (); h1_error.clear (); @@ -231,27 +272,52 @@ void PoissonProblem::clear () { template -void PoissonProblem::create_new () { +void PoissonProblem::create_new (const unsigned int) { clear (); - + tria = new Triangulation(); dof = new DoFHandler (tria); set_tria_and_dof (tria, dof); + boundary = new HyperBallBoundary (); }; +template +void PoissonProblem::declare_parameters (ParameterHandler &prm) { + prm.declare_entry ("Test case", "Gauss shape", "Gauss shape\\|Singular"); + prm.declare_entry ("Initial refinement", "2", + ParameterHandler::RegularExpressions::Integer); + prm.declare_entry ("Refinement criterion", "estimated error", + "global\\|true error\\|estimated error"); + prm.declare_entry ("Refinement fraction", "0.3", + ParameterHandler::RegularExpressions::Double); + prm.declare_entry ("Output base filename", ""); + prm.declare_entry ("Output format", "ucd" + "ucd\\|gnuplot"); +}; -template -void PoissonProblem::run (const unsigned int start_level, - const RefineMode refine_mode) { - create_new (); +template +void PoissonProblem::run (ParameterHandler &prm) { cout << "=======================================" << "=======================================" << endl + << "===== Test case: " << prm.get ("Test case") << endl << "===== Doing computation with refinement criterion: "; + RefineMode refine_mode; + if (prm.get("Refinement criterion")=="global") + refine_mode = global; + else + if (prm.get("Refinement criterion")=="true error") + refine_mode = true_error; + else + if (prm.get("Refinement criterion")=="estimated error") + refine_mode = error_estimator; + else + return; + switch (refine_mode) { case global: @@ -264,16 +330,27 @@ void PoissonProblem::run (const unsigned int start_level, cout << "error estimator"; break; }; + + const unsigned int start_level = prm.get_integer("Initial refinement"); cout << endl << "=======================================" << "=======================================" << endl; cout << "Making initial grid... " << endl; - tria->set_boundary (&boundary); + tria->set_boundary (boundary); tria->create_hyper_ball (); tria->refine_global (start_level); - rhs = new RHSPoly(); - boundary_values = new Solution (); + if (prm.get("Test case")=="Gauss shape") + rhs = new RHS::GaussShape(); + else + if (prm.get("Test case")=="Singular") + rhs = new RHS::Singular(); + + if (prm.get("Test case")=="Gauss shape") + solution_function = new Solution::GaussShape (); + else + if (prm.get("Test case")=="Singular") + solution_function = new Solution::Singular (); FELinear fe; @@ -297,37 +374,38 @@ void PoissonProblem::run (const unsigned int start_level, update_jacobians | update_JxW_values); ProblemBase::FunctionMap dirichlet_bc; - dirichlet_bc[0] = boundary_values; + dirichlet_bc[0] = solution_function; assemble (equation, quadrature, fe, update_flags, dirichlet_bc); cout << " Solving..." << endl; solve (); - Solution sol; dVector l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell; dVector estimated_error_per_cell; dVector laplacian_per_cell; QGauss3 q; cout << " Calculating L2 error... "; - integrate_difference (sol, l2_error_per_cell, q, fe, L2_norm); + integrate_difference (*solution_function, 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); + integrate_difference (*solution_function, 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 error... "; - integrate_difference (sol, h1_error_per_cell, q, fe, H1_norm); + integrate_difference (*solution_function, 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()); cout << " Estimating H1 error... "; KellyErrorEstimator ee; QSimpson eq; - ee.estimate_error (*dof, eq, fe, boundary, + ee.estimate_error (*dof, eq, fe, *boundary, KellyErrorEstimator::FunctionMap(), solution, estimated_error_per_cell); @@ -364,8 +442,7 @@ void PoissonProblem::run (const unsigned int start_level, out.add_data_vector (h1_error_per_dof, "H1-Error"); out.add_data_vector (estimated_error_per_dof, "Estimated Error"); out.add_data_vector (laplacian_per_dof, "Second derivative pointwise"); -// string filename = "gnuplot."; - string filename = "ee."; + String filename = prm.get ("Output base filename"); switch (refine_mode) { case global: @@ -380,15 +457,22 @@ void PoissonProblem::run (const unsigned int start_level, }; filename += ('0'+(start_level+refine_step)/10); filename += ('0'+(start_level+refine_step)%10); - filename += ".inp"; //* + + if (prm.get("Output format")=="ucd") + filename += ".inp"; + else + if (prm.get("Output format")=="gnuplot") + filename += ".gnuplot"; + cout << " Writing error plots to <" << filename << ">..." << endl; - -// ofstream gnuplot(filename.c_str()); -// out.write_gnuplot (gnuplot); -// gnuplot.close (); - ofstream ucd(filename.c_str()); - out.write_ucd (ucd); - ucd.close(); + ofstream outfile(filename); + if (prm.get("Output format")=="ucd") + out.write_ucd (outfile); + else + if (prm.get("Output format")=="gnuplot") + out.write_gnuplot (outfile); + + outfile.close(); cout << " Refining triangulation..."; switch (refine_mode) @@ -397,11 +481,13 @@ void PoissonProblem::run (const unsigned int start_level, tria->refine_global (1); break; case true_error: - tria->refine_fixed_fraction (h1_error_per_cell, 0.5); + tria->refine_fixed_number (h1_error_per_cell, + prm.get_double("Refinement fraction")); tria->execute_refinement (); break; case error_estimator: - tria->refine_fixed_number (estimated_error_per_cell, 0.2); + tria->refine_fixed_number (estimated_error_per_cell, + prm.get_double("Refinement fraction")); tria->execute_refinement (); break; }; @@ -409,14 +495,16 @@ void PoissonProblem::run (const unsigned int start_level, ++refine_step; }; - print_history (refine_mode); + print_history (prm, refine_mode); cout << endl << endl << endl; }; template -void PoissonProblem::print_history (const RefineMode refine_mode) const { - string filename("history."); +void PoissonProblem::print_history (const ParameterHandler &prm, + const RefineMode refine_mode) const { + String filename(prm.get("Output base filename")); + filename += "history."; switch (refine_mode) { case global: @@ -433,7 +521,7 @@ void PoissonProblem::print_history (const RefineMode refine_mode) const { cout << endl << "Printing convergence history to <" << filename << ">..." << endl; - ofstream out(filename.c_str()); + ofstream out(filename); out << "# n_dofs l2_error linfty_error " << "h1_error estimated_error " << "laplacian" @@ -476,13 +564,20 @@ void PoissonProblem::print_history (const RefineMode refine_mode) const { -int main () { - PoissonProblem<2> problem; +int main (int argc, char **argv) { + if (argc!=2) + { + cout << "Usage: error-estimation parameterfile" << endl << endl; + return 1; + }; - problem.run (2, PoissonProblem<2>::global); - problem.run (2, PoissonProblem<2>::true_error); - problem.run (2, PoissonProblem<2>::error_estimator); + PoissonProblem<2> poisson; + MultipleParameterLoop input_data; + poisson.declare_parameters(input_data); + input_data.read_input (argv[1]); + input_data.loop (poisson); + return 0; }; 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 df9f68edaf..ef9c54cf76 100644 --- a/deal.II/deal.II/Attic/examples/error-estimation/make_ps +++ b/deal.II/deal.II/Attic/examples/error-estimation/make_ps @@ -6,21 +6,44 @@ set logscale xy set term postscript eps -set output "history.global.eps" +set output "data-gauss/history.global.eps" -plot "history.global.gnuplot" using 1:2 title "L2 error","history.global.gnuplot" using 1:3 title "Linfty error","history.global.gnuplot" using 1:4 title "H1 error","history.global.gnuplot" using 1:5 title "Estimated H1 error" +plot "data-gauss/history.global.gnuplot" using 1:2 title "L2 error","data-gauss/history.global.gnuplot" using 1:3 title "Linfty error","data-gauss/history.global.gnuplot" using 1:4 title "H1 error","data-gauss/history.global.gnuplot" using 1:5 title "Estimated H1 error" -set output "history.true_error.eps" +set output "data-gauss/history.true_error.eps" -plot "history.true_error.gnuplot" using 1:2 title "L2 error","history.true_error.gnuplot" using 1:3 title "Linfty error","history.true_error.gnuplot" using 1:4 title "H1 error","history.true_error.gnuplot" using 1:5 title "Estimated H1 error" +plot "data-gauss/history.true_error.gnuplot" using 1:2 title "L2 error","data-gauss/history.true_error.gnuplot" using 1:3 title "Linfty error","data-gauss/history.true_error.gnuplot" using 1:4 title "H1 error","data-gauss/history.true_error.gnuplot" using 1:5 title "Estimated H1 error" -set output "history.estimated_error.eps" +set output "data-gauss/history.estimated_error.eps" -plot "history.estimated_error.gnuplot" using 1:2 title "L2 error","history.estimated_error.gnuplot" using 1:3 title "Linfty error","history.estimated_error.gnuplot" using 1:4 title "H1 error","history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" +plot "data-gauss/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-gauss/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-gauss/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-gauss/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" -set output "history.compare.eps" -plot "history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "history.true_error.gnuplot" using 1:2 title "ref. by true error -- L2 error", "history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "history.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" +set output "data-gauss/history.compare.eps" +plot "data-gauss/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-gauss/history.true_error.gnuplot" using 1:2 title "ref. by true error -- L2 error", "data-gauss/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "data-gauss/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-gauss/history.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "data-gauss/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" + + + + + +set output "data-singular/history.global.eps" + +plot "data-singular/history.global.gnuplot" using 1:2 title "L2 error","data-singular/history.global.gnuplot" using 1:3 title "Linfty error","data-singular/history.global.gnuplot" using 1:4 title "H1 error","data-singular/history.global.gnuplot" using 1:5 title "Estimated H1 error" + + +set output "data-singular/history.true_error.eps" + +plot "data-singular/history.true_error.gnuplot" using 1:2 title "L2 error","data-singular/history.true_error.gnuplot" using 1:3 title "Linfty error","data-singular/history.true_error.gnuplot" using 1:4 title "H1 error","data-singular/history.true_error.gnuplot" using 1:5 title "Estimated H1 error" + + + +set output "data-singular/history.estimated_error.eps" + +plot "data-singular/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-singular/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-singular/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-singular/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 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.true_error.gnuplot" using 1:2 title "ref. by true error -- 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.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "data-singular/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" diff --git a/deal.II/deal.II/Attic/examples/error-estimation/strip_comments b/deal.II/deal.II/Attic/examples/error-estimation/strip_comments index 5d2acbe004..779b6b16c7 100755 --- a/deal.II/deal.II/Attic/examples/error-estimation/strip_comments +++ b/deal.II/deal.II/Attic/examples/error-estimation/strip_comments @@ -1 +1 @@ -perl -pi -e 's/^#.*$\\n//g' *.inp +perl -pi -e 's/^#.*$\\n//g' data-*/*.inp diff --git a/tests/big-tests/error-estimation/Makefile b/tests/big-tests/error-estimation/Makefile index 873179dc03..3f035617cb 100644 --- a/tests/big-tests/error-estimation/Makefile +++ b/tests/big-tests/error-estimation/Makefile @@ -19,16 +19,16 @@ user-libs = ../../../mia/control.o # To run the program, use "make run"; to give parameters to the program, # give the parameters to the following variable: -run-parameters = +run-parameters = ee.gauss.prm # To execute additional action apart from running the program, fill # in this list: -additional-run-action = gnuplot make_ps ; ./strip_comments +additional-run-action = $(target) ee.singular.prm ; gnuplot make_ps ; ./strip_comments # To specify which files are to be deleted by "make clean" (apart from # the usual ones: object files, executables, backups, etc), fill in the # following list -delete-files = gnuplot* *inp *.ps *.eps history.* +delete-files = data-*/* diff --git a/tests/big-tests/error-estimation/ee.gauss.prm b/tests/big-tests/error-estimation/ee.gauss.prm new file mode 100644 index 0000000000..588bb9a327 --- /dev/null +++ b/tests/big-tests/error-estimation/ee.gauss.prm @@ -0,0 +1,6 @@ +set Test case = Gauss shape +set Initial refinement = 2 +set Refinement criterion = { global | true error | estimated error } +set Refinement fraction = 0.3 +set Output base filename = data-gauss/ +set Output format = ucd diff --git a/tests/big-tests/error-estimation/ee.singular.prm b/tests/big-tests/error-estimation/ee.singular.prm new file mode 100644 index 0000000000..d906f3f768 --- /dev/null +++ b/tests/big-tests/error-estimation/ee.singular.prm @@ -0,0 +1,6 @@ +set Test case = Singular +set Initial refinement = 1 +set Refinement criterion = { global | true error | estimated error } +set Refinement fraction = 0.2 +set Output base filename = data-singular/ +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 4566b60ab5..2b42b1b95a 100644 --- a/tests/big-tests/error-estimation/error-estimation.cc +++ b/tests/big-tests/error-estimation/error-estimation.cc @@ -9,20 +9,18 @@ #include #include #include +#include #include #include #include #include #include - #include #include #include #include -extern "C" { -# include -} +#include @@ -54,7 +52,8 @@ class PoissonEquation : public Equation { template -class PoissonProblem : public ProblemBase { +class PoissonProblem : public ProblemBase, + public MultipleParameterLoop::UserClass { public: enum RefineMode { global, true_error, error_estimator @@ -63,19 +62,20 @@ class PoissonProblem : public ProblemBase { PoissonProblem (); void clear (); - void create_new (); - void run (const unsigned int start_level, - const RefineMode refine_mode); - void print_history (const RefineMode refine_mode) const; + void create_new (const unsigned int); + void declare_parameters (ParameterHandler &prm); + void run (ParameterHandler &prm); + void print_history (const ParameterHandler &prm, + const RefineMode refine_mode) const; protected: Triangulation *tria; DoFHandler *dof; Function *rhs; - Function *boundary_values; + Function *solution_function; - HyperBallBoundary boundary; + Boundary *boundary; vector l2_error, linfty_error; vector h1_error, estimated_error; @@ -88,56 +88,91 @@ class PoissonProblem : public ProblemBase { -/** - Right hand side constructed such that the exact solution is - $x*y*exp(-(x**2+y**2)*40)$. - */ template -class RHSPoly : public Function { +class Solution { public: - /** - * Return the value of the function - * at the given point. - */ - virtual double operator () (const Point &p) const; + + class GaussShape : public Function { + public: + virtual double operator () (const Point &p) const; + virtual Point gradient (const Point &p) const; + }; + + class Singular : public Function { + public: + virtual double operator () (const Point &p) const; + virtual Point gradient (const Point &p) const; + }; }; + template -class Solution : public Function { +class RHS { 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. + * Right hand side constructed such that + * the exact solution is + * $x*y*exp(-(x**2+y**2)*40)$. */ - virtual Point gradient (const Point &p) const; -}; - + class GaussShape : public Function { + public: + virtual double operator () (const Point &p) const; + }; + /** + * Right hand side constructed such that + * the exact solution is + * $r^{2/3} sin(2\phi)$. + */ + class Singular : public Function { + public: + virtual double operator () (const Point &p) const; + }; +}; -double RHSPoly<2>::operator () (const Point<2> &p) const { - return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); -}; -double Solution<2>::operator () (const Point<2> &p) const { +double Solution<2>::GaussShape::operator () (const Point<2> &p) const { return p(0)*p(1)*exp(-40*p.square()); }; -Point<2> Solution<2>::gradient (const Point<2> &p) const { +Point<2> Solution<2>::GaussShape::gradient (const Point<2> &p) const { return Point<2> ((1-80.*p(0)*p(0))*p(1)*exp(-40*p.square()), (1-80.*p(1)*p(1))*p(0)*exp(-40*p.square())); }; + + +double Solution<2>::Singular::operator () (const Point<2> &p) const { + return pow(p.square(), 1./3.); +}; + + +Point<2> Solution<2>::Singular::gradient (const Point<2> &p) const { + return 2./3.*pow(p.square(), -2./3.) * p; +}; + + + +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.); +}; + + + + @@ -189,7 +224,7 @@ void PoissonEquation::assemble (dVector &, template PoissonProblem::PoissonProblem () : tria(0), dof(0), rhs(0), - boundary_values(0), boundary() {}; + solution_function(0), boundary(0) {}; @@ -212,12 +247,18 @@ void PoissonProblem::clear () { rhs = 0; }; - if (boundary_values != 0) + if (solution_function != 0) { - delete boundary_values; - boundary_values = 0; + delete solution_function; + solution_function = 0; }; + if (boundary != 0) + { + delete boundary; + boundary = 0; + }; + l2_error.clear (); linfty_error.clear (); h1_error.clear (); @@ -231,27 +272,52 @@ void PoissonProblem::clear () { template -void PoissonProblem::create_new () { +void PoissonProblem::create_new (const unsigned int) { clear (); - + tria = new Triangulation(); dof = new DoFHandler (tria); set_tria_and_dof (tria, dof); + boundary = new HyperBallBoundary (); }; +template +void PoissonProblem::declare_parameters (ParameterHandler &prm) { + prm.declare_entry ("Test case", "Gauss shape", "Gauss shape\\|Singular"); + prm.declare_entry ("Initial refinement", "2", + ParameterHandler::RegularExpressions::Integer); + prm.declare_entry ("Refinement criterion", "estimated error", + "global\\|true error\\|estimated error"); + prm.declare_entry ("Refinement fraction", "0.3", + ParameterHandler::RegularExpressions::Double); + prm.declare_entry ("Output base filename", ""); + prm.declare_entry ("Output format", "ucd" + "ucd\\|gnuplot"); +}; -template -void PoissonProblem::run (const unsigned int start_level, - const RefineMode refine_mode) { - create_new (); +template +void PoissonProblem::run (ParameterHandler &prm) { cout << "=======================================" << "=======================================" << endl + << "===== Test case: " << prm.get ("Test case") << endl << "===== Doing computation with refinement criterion: "; + RefineMode refine_mode; + if (prm.get("Refinement criterion")=="global") + refine_mode = global; + else + if (prm.get("Refinement criterion")=="true error") + refine_mode = true_error; + else + if (prm.get("Refinement criterion")=="estimated error") + refine_mode = error_estimator; + else + return; + switch (refine_mode) { case global: @@ -264,16 +330,27 @@ void PoissonProblem::run (const unsigned int start_level, cout << "error estimator"; break; }; + + const unsigned int start_level = prm.get_integer("Initial refinement"); cout << endl << "=======================================" << "=======================================" << endl; cout << "Making initial grid... " << endl; - tria->set_boundary (&boundary); + tria->set_boundary (boundary); tria->create_hyper_ball (); tria->refine_global (start_level); - rhs = new RHSPoly(); - boundary_values = new Solution (); + if (prm.get("Test case")=="Gauss shape") + rhs = new RHS::GaussShape(); + else + if (prm.get("Test case")=="Singular") + rhs = new RHS::Singular(); + + if (prm.get("Test case")=="Gauss shape") + solution_function = new Solution::GaussShape (); + else + if (prm.get("Test case")=="Singular") + solution_function = new Solution::Singular (); FELinear fe; @@ -297,37 +374,38 @@ void PoissonProblem::run (const unsigned int start_level, update_jacobians | update_JxW_values); ProblemBase::FunctionMap dirichlet_bc; - dirichlet_bc[0] = boundary_values; + dirichlet_bc[0] = solution_function; assemble (equation, quadrature, fe, update_flags, dirichlet_bc); cout << " Solving..." << endl; solve (); - Solution sol; dVector l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell; dVector estimated_error_per_cell; dVector laplacian_per_cell; QGauss3 q; cout << " Calculating L2 error... "; - integrate_difference (sol, l2_error_per_cell, q, fe, L2_norm); + integrate_difference (*solution_function, 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); + integrate_difference (*solution_function, 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 error... "; - integrate_difference (sol, h1_error_per_cell, q, fe, H1_norm); + integrate_difference (*solution_function, 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()); cout << " Estimating H1 error... "; KellyErrorEstimator ee; QSimpson eq; - ee.estimate_error (*dof, eq, fe, boundary, + ee.estimate_error (*dof, eq, fe, *boundary, KellyErrorEstimator::FunctionMap(), solution, estimated_error_per_cell); @@ -364,8 +442,7 @@ void PoissonProblem::run (const unsigned int start_level, out.add_data_vector (h1_error_per_dof, "H1-Error"); out.add_data_vector (estimated_error_per_dof, "Estimated Error"); out.add_data_vector (laplacian_per_dof, "Second derivative pointwise"); -// string filename = "gnuplot."; - string filename = "ee."; + String filename = prm.get ("Output base filename"); switch (refine_mode) { case global: @@ -380,15 +457,22 @@ void PoissonProblem::run (const unsigned int start_level, }; filename += ('0'+(start_level+refine_step)/10); filename += ('0'+(start_level+refine_step)%10); - filename += ".inp"; //* + + if (prm.get("Output format")=="ucd") + filename += ".inp"; + else + if (prm.get("Output format")=="gnuplot") + filename += ".gnuplot"; + cout << " Writing error plots to <" << filename << ">..." << endl; - -// ofstream gnuplot(filename.c_str()); -// out.write_gnuplot (gnuplot); -// gnuplot.close (); - ofstream ucd(filename.c_str()); - out.write_ucd (ucd); - ucd.close(); + ofstream outfile(filename); + if (prm.get("Output format")=="ucd") + out.write_ucd (outfile); + else + if (prm.get("Output format")=="gnuplot") + out.write_gnuplot (outfile); + + outfile.close(); cout << " Refining triangulation..."; switch (refine_mode) @@ -397,11 +481,13 @@ void PoissonProblem::run (const unsigned int start_level, tria->refine_global (1); break; case true_error: - tria->refine_fixed_fraction (h1_error_per_cell, 0.5); + tria->refine_fixed_number (h1_error_per_cell, + prm.get_double("Refinement fraction")); tria->execute_refinement (); break; case error_estimator: - tria->refine_fixed_number (estimated_error_per_cell, 0.2); + tria->refine_fixed_number (estimated_error_per_cell, + prm.get_double("Refinement fraction")); tria->execute_refinement (); break; }; @@ -409,14 +495,16 @@ void PoissonProblem::run (const unsigned int start_level, ++refine_step; }; - print_history (refine_mode); + print_history (prm, refine_mode); cout << endl << endl << endl; }; template -void PoissonProblem::print_history (const RefineMode refine_mode) const { - string filename("history."); +void PoissonProblem::print_history (const ParameterHandler &prm, + const RefineMode refine_mode) const { + String filename(prm.get("Output base filename")); + filename += "history."; switch (refine_mode) { case global: @@ -433,7 +521,7 @@ void PoissonProblem::print_history (const RefineMode refine_mode) const { cout << endl << "Printing convergence history to <" << filename << ">..." << endl; - ofstream out(filename.c_str()); + ofstream out(filename); out << "# n_dofs l2_error linfty_error " << "h1_error estimated_error " << "laplacian" @@ -476,13 +564,20 @@ void PoissonProblem::print_history (const RefineMode refine_mode) const { -int main () { - PoissonProblem<2> problem; +int main (int argc, char **argv) { + if (argc!=2) + { + cout << "Usage: error-estimation parameterfile" << endl << endl; + return 1; + }; - problem.run (2, PoissonProblem<2>::global); - problem.run (2, PoissonProblem<2>::true_error); - problem.run (2, PoissonProblem<2>::error_estimator); + PoissonProblem<2> poisson; + MultipleParameterLoop input_data; + poisson.declare_parameters(input_data); + input_data.read_input (argv[1]); + input_data.loop (poisson); + return 0; }; diff --git a/tests/big-tests/error-estimation/make_ps b/tests/big-tests/error-estimation/make_ps index df9f68edaf..ef9c54cf76 100644 --- a/tests/big-tests/error-estimation/make_ps +++ b/tests/big-tests/error-estimation/make_ps @@ -6,21 +6,44 @@ set logscale xy set term postscript eps -set output "history.global.eps" +set output "data-gauss/history.global.eps" -plot "history.global.gnuplot" using 1:2 title "L2 error","history.global.gnuplot" using 1:3 title "Linfty error","history.global.gnuplot" using 1:4 title "H1 error","history.global.gnuplot" using 1:5 title "Estimated H1 error" +plot "data-gauss/history.global.gnuplot" using 1:2 title "L2 error","data-gauss/history.global.gnuplot" using 1:3 title "Linfty error","data-gauss/history.global.gnuplot" using 1:4 title "H1 error","data-gauss/history.global.gnuplot" using 1:5 title "Estimated H1 error" -set output "history.true_error.eps" +set output "data-gauss/history.true_error.eps" -plot "history.true_error.gnuplot" using 1:2 title "L2 error","history.true_error.gnuplot" using 1:3 title "Linfty error","history.true_error.gnuplot" using 1:4 title "H1 error","history.true_error.gnuplot" using 1:5 title "Estimated H1 error" +plot "data-gauss/history.true_error.gnuplot" using 1:2 title "L2 error","data-gauss/history.true_error.gnuplot" using 1:3 title "Linfty error","data-gauss/history.true_error.gnuplot" using 1:4 title "H1 error","data-gauss/history.true_error.gnuplot" using 1:5 title "Estimated H1 error" -set output "history.estimated_error.eps" +set output "data-gauss/history.estimated_error.eps" -plot "history.estimated_error.gnuplot" using 1:2 title "L2 error","history.estimated_error.gnuplot" using 1:3 title "Linfty error","history.estimated_error.gnuplot" using 1:4 title "H1 error","history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" +plot "data-gauss/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-gauss/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-gauss/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-gauss/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 error" -set output "history.compare.eps" -plot "history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "history.true_error.gnuplot" using 1:2 title "ref. by true error -- L2 error", "history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "history.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" +set output "data-gauss/history.compare.eps" +plot "data-gauss/history.global.gnuplot" using 1:2 title "global refinement -- L2 error", "data-gauss/history.true_error.gnuplot" using 1:2 title "ref. by true error -- L2 error", "data-gauss/history.estimated_error.gnuplot" using 1:2 title "ref. by estimated error -- L2 error", "data-gauss/history.global.gnuplot" using 1:4 title "global refinement -- H1 error", "data-gauss/history.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "data-gauss/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" + + + + + +set output "data-singular/history.global.eps" + +plot "data-singular/history.global.gnuplot" using 1:2 title "L2 error","data-singular/history.global.gnuplot" using 1:3 title "Linfty error","data-singular/history.global.gnuplot" using 1:4 title "H1 error","data-singular/history.global.gnuplot" using 1:5 title "Estimated H1 error" + + +set output "data-singular/history.true_error.eps" + +plot "data-singular/history.true_error.gnuplot" using 1:2 title "L2 error","data-singular/history.true_error.gnuplot" using 1:3 title "Linfty error","data-singular/history.true_error.gnuplot" using 1:4 title "H1 error","data-singular/history.true_error.gnuplot" using 1:5 title "Estimated H1 error" + + + +set output "data-singular/history.estimated_error.eps" + +plot "data-singular/history.estimated_error.gnuplot" using 1:2 title "L2 error","data-singular/history.estimated_error.gnuplot" using 1:3 title "Linfty error","data-singular/history.estimated_error.gnuplot" using 1:4 title "H1 error","data-singular/history.estimated_error.gnuplot" using 1:5 title "Estimated H1 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.true_error.gnuplot" using 1:2 title "ref. by true error -- 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.true_error.gnuplot" using 1:4 title "ref. by true error -- H1 error", "data-singular/history.estimated_error.gnuplot" using 1:4 title "ref. by estimated error -- H1 error" diff --git a/tests/big-tests/error-estimation/strip_comments b/tests/big-tests/error-estimation/strip_comments index 5d2acbe004..779b6b16c7 100755 --- a/tests/big-tests/error-estimation/strip_comments +++ b/tests/big-tests/error-estimation/strip_comments @@ -1 +1 @@ -perl -pi -e 's/^#.*$\\n//g' *.inp +perl -pi -e 's/^#.*$\\n//g' data-*/*.inp -- 2.39.5