From: Wolfgang Bangerth Date: Fri, 3 Jul 1998 09:11:25 +0000 (+0000) Subject: Update example to also use higher order elements. X-Git-Tag: v8.0.0~22820 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=829816d401196dd7b18677d10dd546d590605dbc;p=dealii.git Update example to also use higher order elements. git-svn-id: https://svn.dealii.org/trunk@434 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/convergence/Makefile b/deal.II/deal.II/Attic/examples/convergence/Makefile index 5361794214..d3c0974d42 100644 --- a/deal.II/deal.II/Attic/examples/convergence/Makefile +++ b/deal.II/deal.II/Attic/examples/convergence/Makefile @@ -11,7 +11,7 @@ target = convergence # the archives ./Obj.a and ./Obj.g.a. By default, the debug version # is used to link. It you don't like that, change the following # variable to "off" -debug-mode = on +debug-mode = off # If you want your program to be linked with extra object or library # files, specify them here: diff --git a/deal.II/deal.II/Attic/examples/convergence/convergence.cc b/deal.II/deal.II/Attic/examples/convergence/convergence.cc index bad4d474f6..0934316fa7 100644 --- a/deal.II/deal.II/Attic/examples/convergence/convergence.cc +++ b/deal.II/deal.II/Attic/examples/convergence/convergence.cc @@ -55,12 +55,12 @@ class PoissonEquation : public Equation { template class PoissonProblem : public ProblemBase { public: - PoissonProblem (); + PoissonProblem (unsigned int order); void clear (); void create_new (); - void run (unsigned int level); - void print_history () const; + int run (unsigned int level); + void print_history (string filename) const; protected: Triangulation *tria; @@ -71,6 +71,8 @@ class PoissonProblem : public ProblemBase { vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; vector n_dofs; + + unsigned int order; }; @@ -79,7 +81,7 @@ class PoissonProblem : public ProblemBase { /** Right hand side constructed such that the exact solution is - $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc. + $sin(2 pi x) + sin(2 pi y)$ */ template class RHSPoly : public Function { @@ -112,18 +114,28 @@ class Solution : public Function { double RHSPoly<2>::operator () (const Point<2> &p) const { - return 2*(p(0)*(1-p(0)) + p(1)*(1-p(1))); + const double x = p(0), + y = p(1); + const double pi= 3.1415926536; + return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y)); }; double Solution<2>::operator () (const Point<2> &p) const { - return p(0)*(1-p(0))*p(1)*(1-p(1)); + const double x = p(0), + y = p(1); + const double pi= 3.1415926536; + return sin(2*pi*x)+sin(2*pi*y); }; Point<2> Solution<2>::gradient (const Point<2> &p) const { - return Point<2> ((1-2*p(0))*p(1)*(1-p(1)), (1-2*p(1))*p(0)*(1-p(0))); + const double x = p(0), + y = p(1); + const double pi= 3.1415926536; + return Point<2> (2*pi*cos(2*pi*x), + 2*pi*cos(2*pi*y)); }; @@ -175,8 +187,9 @@ void PoissonEquation::assemble (dVector &, template -PoissonProblem::PoissonProblem () : - tria(0), dof(0), rhs(0), boundary_values(0) {}; +PoissonProblem::PoissonProblem (unsigned int order) : + tria(0), dof(0), rhs(0), + boundary_values(0), order(order) {}; @@ -223,30 +236,59 @@ void PoissonProblem::create_new () { -#include - template -void PoissonProblem::run (const unsigned int level) { +int PoissonProblem::run (const unsigned int level) { create_new (); cout << "Refinement level = " << level + << ", using elements of order " << order << endl; cout << " Making grid... "; - tria->create_hypercube (); + tria->create_hyper_ball (); + HyperBallBoundary boundary_description; + tria->set_boundary (&boundary_description); + tria->begin_active()->set_refine_flag(); + (++(++(tria->begin_active())))->set_refine_flag(); + tria->execute_refinement (); tria->refine_global (level); cout << tria->n_active_cells() << " active cells." << endl; rhs = new RHSPoly(); - boundary_values = new ZeroFunction (); + boundary_values = new Solution (); - - FELinear fe; - PoissonEquation equation (*rhs); - QGauss3 quadrature; + + FiniteElement *fe; + PoissonEquation equation (*rhs); + Quadrature *quadrature; + Quadrature *boundary_quadrature; + switch (order) { + case 1: + fe = new FELinear(); + quadrature = new QGauss3(); + boundary_quadrature = new QGauss2(); + break; + case 2: + fe = new FEQuadraticSub(); + quadrature = new QGauss4(); + boundary_quadrature = new QGauss3(); + break; + case 3: + fe = new FECubicSub(); + quadrature = new QGauss5(); + boundary_quadrature = new QGauss4(); + break; + case 4: + fe = new FEQuarticSub(); + quadrature = new QGauss6(); + boundary_quadrature = new QGauss5(); + break; + default: + return 100000; + }; cout << " Distributing dofs... "; - dof->distribute_dofs (fe); + dof->distribute_dofs (*fe); cout << dof->n_dofs() << " degrees of freedom." << endl; n_dofs.push_back (dof->n_dofs()); @@ -256,7 +298,7 @@ void PoissonProblem::run (const unsigned int level) { ProblemBase::FunctionMap dirichlet_bc; dirichlet_bc[0] = boundary_values; - assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + assemble (equation, *quadrature, *fe, update_flags, dirichlet_bc); cout << " Solving..." << endl; solve (); @@ -264,13 +306,12 @@ void PoissonProblem::run (const unsigned int level) { 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... "; VectorTools::integrate_difference (*dof_handler, solution, sol, l1_error_per_cell, - q, fe, L1_norm); + *quadrature, *fe, L1_norm); cout << l1_error_per_cell.l1_norm() << endl; l1_error.push_back (l1_error_per_cell.l1_norm()); @@ -278,7 +319,7 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, l2_error_per_cell, - q, fe, L2_norm); + *quadrature, *fe, L2_norm); cout << l2_error_per_cell.l2_norm() << endl; l2_error.push_back (l2_error_per_cell.l2_norm()); @@ -286,7 +327,7 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, linfty_error_per_cell, - q, fe, Linfty_norm); + *quadrature, *fe, Linfty_norm); cout << linfty_error_per_cell.linfty_norm() << endl; linfty_error.push_back (linfty_error_per_cell.linfty_norm()); @@ -294,7 +335,7 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, h1_seminorm_error_per_cell, - q, fe, H1_seminorm); + *quadrature, *fe, H1_seminorm); cout << h1_seminorm_error_per_cell.l2_norm() << endl; h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); @@ -302,45 +343,72 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, h1_error_per_cell, - q, fe, H1_norm); + *quadrature, *fe, H1_norm); cout << h1_error_per_cell.l2_norm() << endl; h1_error.push_back (h1_error_per_cell.l2_norm()); - if (level<=5) + if (dof->n_dofs()<=5000) { 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 (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."; + dVector projected_solution; + ConstraintMatrix constraints; + constraints.close (); + VectorTools::project (*dof, constraints, *fe, + StraightBoundary(), *quadrature, + sol, projected_solution, false, + *boundary_quadrature); + cout << " Calculating L2 error of projected solution... "; + VectorTools::integrate_difference (*dof_handler, + projected_solution, sol, + l2_error_per_cell, + *quadrature, *fe, L2_norm); + cout << l2_error_per_cell.l2_norm() << endl; + + + string filename; + filename = ('0'+order); + filename += "."; filename += ('0'+level); + filename += ".ucd"; cout << " Writing error plots to <" << filename << ">..." << endl; DataOut out; - ofstream gnuplot(filename.c_str()); + ofstream o(filename.c_str()); fill_data (out); + out.add_data_vector (projected_solution, "projected u"); 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 (); + out.write_ucd (o); + o.close (); } else cout << " Not writing error as grid." << endl; cout << endl; + + delete fe; + delete quadrature; + delete boundary_quadrature; + + return dof->n_dofs(); }; template -void PoissonProblem::print_history () const { - ofstream out("gnuplot.history"); +void PoissonProblem::print_history (string filename) const { + ofstream out(filename.c_str()); out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" << endl; for (unsigned int i=0; i::print_history () const { average_h1_semi /= (l1_error.size()-1); average_h1 /= (l1_error.size()-1); + cout << "==========================================================\n"; cout << "Average error reduction rates for h->h/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; + cout << "==========================================================\n"; + cout << "==========================================================\n"; }; int main () { - PoissonProblem<2> problem; - - for (unsigned int level=1; level<8; ++level) - problem.run (level); - - cout << endl << "Printing convergence history to ..." << endl; - problem.print_history (); - + for (unsigned int order=1; order<5; ++order) + { + PoissonProblem<2> problem (order); + + unsigned int level=0; + unsigned int n_dofs; + do + n_dofs = problem.run (level++); + while (n_dofs<25000); + + string filename; + switch (order) + { + case 1: + filename = "linear"; + break; + case 2: + filename = "quadratic"; + break; + case 3: + filename = "cubic"; + break; + case 4: + filename = "quartic"; + break; + }; + filename += ".history"; + + cout << endl << "Printing convergence history to <" + << filename << ">..." << endl; + problem.print_history (filename); + cout << endl << endl << endl; + }; + return 0; }; diff --git a/deal.II/deal.II/Attic/examples/convergence/make_ps b/deal.II/deal.II/Attic/examples/convergence/make_ps index ea538afd10..502062d53d 100644 --- a/deal.II/deal.II/Attic/examples/convergence/make_ps +++ b/deal.II/deal.II/Attic/examples/convergence/make_ps @@ -4,6 +4,45 @@ set data style linespoints set logscale xy set term postscript eps -set output "convergence.eps" +set output "linear.eps" + +plot "linear.history" using 1:2 title "L1 error","linear.history" using 1:3 title "L2 error","linear.history" using 1:4 title "Linfty error","linear.history" using 1:5 title "H1 seminorm error","linear.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "quadratic.eps" + +plot "quadratic.history" using 1:2 title "L1 error","quadratic.history" using 1:3 title "L2 error","quadratic.history" using 1:4 title "Linfty error","quadratic.history" using 1:5 title "H1 seminorm error","quadratic.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "cubic.eps" + +plot "cubic.history" using 1:2 title "L1 error","cubic.history" using 1:3 title "L2 error","cubic.history" using 1:4 title "Linfty error","cubic.history" using 1:5 title "H1 seminorm error","cubic.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "quartic.eps" + +plot "quartic.history" using 1:2 title "L1 error","quartic.history" using 1:3 title "L2 error","quartic.history" using 1:4 title "Linfty error","quartic.history" using 1:5 title "H1 seminorm error","quartic.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "l2error.eps" +set ylabel "L2-error" + +plot "linear.history" using 1:3 title "Linear elements", "quadratic.history" using 1:3 title "Quadratic elements", "cubic.history" using 1:3 title "Cubic elements", "quartic.history" using 1:3 title "Quartic elements" + + + +set term postscript eps +set output "h1error.eps" +set ylabel "H1-error" + +plot "linear.history" using 1:6 title "Linear elements", "quadratic.history" using 1:6 title "Quadratic elements", "cubic.history" using 1:6 title "Cubic elements", "quartic.history" using 1:6 title "Quartic elements" + -plot "gnuplot.history" using 1:2 title "L1 error","gnuplot.history" using 1:3 title "L2 error","gnuplot.history" using 1:4 title "Linfty error","gnuplot.history" using 1:5 title "H1 seminorm error","gnuplot.history" using 1:6 title "H1 error" diff --git a/tests/big-tests/convergence/Makefile b/tests/big-tests/convergence/Makefile index 5361794214..d3c0974d42 100644 --- a/tests/big-tests/convergence/Makefile +++ b/tests/big-tests/convergence/Makefile @@ -11,7 +11,7 @@ target = convergence # the archives ./Obj.a and ./Obj.g.a. By default, the debug version # is used to link. It you don't like that, change the following # variable to "off" -debug-mode = on +debug-mode = off # If you want your program to be linked with extra object or library # files, specify them here: diff --git a/tests/big-tests/convergence/convergence.cc b/tests/big-tests/convergence/convergence.cc index bad4d474f6..0934316fa7 100644 --- a/tests/big-tests/convergence/convergence.cc +++ b/tests/big-tests/convergence/convergence.cc @@ -55,12 +55,12 @@ class PoissonEquation : public Equation { template class PoissonProblem : public ProblemBase { public: - PoissonProblem (); + PoissonProblem (unsigned int order); void clear (); void create_new (); - void run (unsigned int level); - void print_history () const; + int run (unsigned int level); + void print_history (string filename) const; protected: Triangulation *tria; @@ -71,6 +71,8 @@ class PoissonProblem : public ProblemBase { vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; vector n_dofs; + + unsigned int order; }; @@ -79,7 +81,7 @@ class PoissonProblem : public ProblemBase { /** Right hand side constructed such that the exact solution is - $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc. + $sin(2 pi x) + sin(2 pi y)$ */ template class RHSPoly : public Function { @@ -112,18 +114,28 @@ class Solution : public Function { double RHSPoly<2>::operator () (const Point<2> &p) const { - return 2*(p(0)*(1-p(0)) + p(1)*(1-p(1))); + const double x = p(0), + y = p(1); + const double pi= 3.1415926536; + return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y)); }; double Solution<2>::operator () (const Point<2> &p) const { - return p(0)*(1-p(0))*p(1)*(1-p(1)); + const double x = p(0), + y = p(1); + const double pi= 3.1415926536; + return sin(2*pi*x)+sin(2*pi*y); }; Point<2> Solution<2>::gradient (const Point<2> &p) const { - return Point<2> ((1-2*p(0))*p(1)*(1-p(1)), (1-2*p(1))*p(0)*(1-p(0))); + const double x = p(0), + y = p(1); + const double pi= 3.1415926536; + return Point<2> (2*pi*cos(2*pi*x), + 2*pi*cos(2*pi*y)); }; @@ -175,8 +187,9 @@ void PoissonEquation::assemble (dVector &, template -PoissonProblem::PoissonProblem () : - tria(0), dof(0), rhs(0), boundary_values(0) {}; +PoissonProblem::PoissonProblem (unsigned int order) : + tria(0), dof(0), rhs(0), + boundary_values(0), order(order) {}; @@ -223,30 +236,59 @@ void PoissonProblem::create_new () { -#include - template -void PoissonProblem::run (const unsigned int level) { +int PoissonProblem::run (const unsigned int level) { create_new (); cout << "Refinement level = " << level + << ", using elements of order " << order << endl; cout << " Making grid... "; - tria->create_hypercube (); + tria->create_hyper_ball (); + HyperBallBoundary boundary_description; + tria->set_boundary (&boundary_description); + tria->begin_active()->set_refine_flag(); + (++(++(tria->begin_active())))->set_refine_flag(); + tria->execute_refinement (); tria->refine_global (level); cout << tria->n_active_cells() << " active cells." << endl; rhs = new RHSPoly(); - boundary_values = new ZeroFunction (); + boundary_values = new Solution (); - - FELinear fe; - PoissonEquation equation (*rhs); - QGauss3 quadrature; + + FiniteElement *fe; + PoissonEquation equation (*rhs); + Quadrature *quadrature; + Quadrature *boundary_quadrature; + switch (order) { + case 1: + fe = new FELinear(); + quadrature = new QGauss3(); + boundary_quadrature = new QGauss2(); + break; + case 2: + fe = new FEQuadraticSub(); + quadrature = new QGauss4(); + boundary_quadrature = new QGauss3(); + break; + case 3: + fe = new FECubicSub(); + quadrature = new QGauss5(); + boundary_quadrature = new QGauss4(); + break; + case 4: + fe = new FEQuarticSub(); + quadrature = new QGauss6(); + boundary_quadrature = new QGauss5(); + break; + default: + return 100000; + }; cout << " Distributing dofs... "; - dof->distribute_dofs (fe); + dof->distribute_dofs (*fe); cout << dof->n_dofs() << " degrees of freedom." << endl; n_dofs.push_back (dof->n_dofs()); @@ -256,7 +298,7 @@ void PoissonProblem::run (const unsigned int level) { ProblemBase::FunctionMap dirichlet_bc; dirichlet_bc[0] = boundary_values; - assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + assemble (equation, *quadrature, *fe, update_flags, dirichlet_bc); cout << " Solving..." << endl; solve (); @@ -264,13 +306,12 @@ void PoissonProblem::run (const unsigned int level) { 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... "; VectorTools::integrate_difference (*dof_handler, solution, sol, l1_error_per_cell, - q, fe, L1_norm); + *quadrature, *fe, L1_norm); cout << l1_error_per_cell.l1_norm() << endl; l1_error.push_back (l1_error_per_cell.l1_norm()); @@ -278,7 +319,7 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, l2_error_per_cell, - q, fe, L2_norm); + *quadrature, *fe, L2_norm); cout << l2_error_per_cell.l2_norm() << endl; l2_error.push_back (l2_error_per_cell.l2_norm()); @@ -286,7 +327,7 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, linfty_error_per_cell, - q, fe, Linfty_norm); + *quadrature, *fe, Linfty_norm); cout << linfty_error_per_cell.linfty_norm() << endl; linfty_error.push_back (linfty_error_per_cell.linfty_norm()); @@ -294,7 +335,7 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, h1_seminorm_error_per_cell, - q, fe, H1_seminorm); + *quadrature, *fe, H1_seminorm); cout << h1_seminorm_error_per_cell.l2_norm() << endl; h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); @@ -302,45 +343,72 @@ void PoissonProblem::run (const unsigned int level) { VectorTools::integrate_difference (*dof_handler, solution, sol, h1_error_per_cell, - q, fe, H1_norm); + *quadrature, *fe, H1_norm); cout << h1_error_per_cell.l2_norm() << endl; h1_error.push_back (h1_error_per_cell.l2_norm()); - if (level<=5) + if (dof->n_dofs()<=5000) { 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 (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."; + dVector projected_solution; + ConstraintMatrix constraints; + constraints.close (); + VectorTools::project (*dof, constraints, *fe, + StraightBoundary(), *quadrature, + sol, projected_solution, false, + *boundary_quadrature); + cout << " Calculating L2 error of projected solution... "; + VectorTools::integrate_difference (*dof_handler, + projected_solution, sol, + l2_error_per_cell, + *quadrature, *fe, L2_norm); + cout << l2_error_per_cell.l2_norm() << endl; + + + string filename; + filename = ('0'+order); + filename += "."; filename += ('0'+level); + filename += ".ucd"; cout << " Writing error plots to <" << filename << ">..." << endl; DataOut out; - ofstream gnuplot(filename.c_str()); + ofstream o(filename.c_str()); fill_data (out); + out.add_data_vector (projected_solution, "projected u"); 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 (); + out.write_ucd (o); + o.close (); } else cout << " Not writing error as grid." << endl; cout << endl; + + delete fe; + delete quadrature; + delete boundary_quadrature; + + return dof->n_dofs(); }; template -void PoissonProblem::print_history () const { - ofstream out("gnuplot.history"); +void PoissonProblem::print_history (string filename) const { + ofstream out(filename.c_str()); out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" << endl; for (unsigned int i=0; i::print_history () const { average_h1_semi /= (l1_error.size()-1); average_h1 /= (l1_error.size()-1); + cout << "==========================================================\n"; cout << "Average error reduction rates for h->h/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; + cout << "==========================================================\n"; + cout << "==========================================================\n"; }; int main () { - PoissonProblem<2> problem; - - for (unsigned int level=1; level<8; ++level) - problem.run (level); - - cout << endl << "Printing convergence history to ..." << endl; - problem.print_history (); - + for (unsigned int order=1; order<5; ++order) + { + PoissonProblem<2> problem (order); + + unsigned int level=0; + unsigned int n_dofs; + do + n_dofs = problem.run (level++); + while (n_dofs<25000); + + string filename; + switch (order) + { + case 1: + filename = "linear"; + break; + case 2: + filename = "quadratic"; + break; + case 3: + filename = "cubic"; + break; + case 4: + filename = "quartic"; + break; + }; + filename += ".history"; + + cout << endl << "Printing convergence history to <" + << filename << ">..." << endl; + problem.print_history (filename); + cout << endl << endl << endl; + }; + return 0; }; diff --git a/tests/big-tests/convergence/make_ps b/tests/big-tests/convergence/make_ps index ea538afd10..502062d53d 100644 --- a/tests/big-tests/convergence/make_ps +++ b/tests/big-tests/convergence/make_ps @@ -4,6 +4,45 @@ set data style linespoints set logscale xy set term postscript eps -set output "convergence.eps" +set output "linear.eps" + +plot "linear.history" using 1:2 title "L1 error","linear.history" using 1:3 title "L2 error","linear.history" using 1:4 title "Linfty error","linear.history" using 1:5 title "H1 seminorm error","linear.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "quadratic.eps" + +plot "quadratic.history" using 1:2 title "L1 error","quadratic.history" using 1:3 title "L2 error","quadratic.history" using 1:4 title "Linfty error","quadratic.history" using 1:5 title "H1 seminorm error","quadratic.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "cubic.eps" + +plot "cubic.history" using 1:2 title "L1 error","cubic.history" using 1:3 title "L2 error","cubic.history" using 1:4 title "Linfty error","cubic.history" using 1:5 title "H1 seminorm error","cubic.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "quartic.eps" + +plot "quartic.history" using 1:2 title "L1 error","quartic.history" using 1:3 title "L2 error","quartic.history" using 1:4 title "Linfty error","quartic.history" using 1:5 title "H1 seminorm error","quartic.history" using 1:6 title "H1 error" + + + +set term postscript eps +set output "l2error.eps" +set ylabel "L2-error" + +plot "linear.history" using 1:3 title "Linear elements", "quadratic.history" using 1:3 title "Quadratic elements", "cubic.history" using 1:3 title "Cubic elements", "quartic.history" using 1:3 title "Quartic elements" + + + +set term postscript eps +set output "h1error.eps" +set ylabel "H1-error" + +plot "linear.history" using 1:6 title "Linear elements", "quadratic.history" using 1:6 title "Quadratic elements", "cubic.history" using 1:6 title "Cubic elements", "quartic.history" using 1:6 title "Quartic elements" + -plot "gnuplot.history" using 1:2 title "L1 error","gnuplot.history" using 1:3 title "L2 error","gnuplot.history" using 1:4 title "Linfty error","gnuplot.history" using 1:5 title "H1 seminorm error","gnuplot.history" using 1:6 title "H1 error"