]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add adaptive refinement.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 4 May 1998 12:04:03 +0000 (12:04 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 4 May 1998 12:04:03 +0000 (12:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@242 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/error-estimation/Makefile
deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc
deal.II/deal.II/Attic/examples/error-estimation/make_ps
tests/big-tests/error-estimation/Makefile
tests/big-tests/error-estimation/error-estimation.cc
tests/big-tests/error-estimation/make_ps

index ee6a9b67138c205b2f769d5cf4cc2f9887d89ae5..772ffb9c4367e8ac92d5f89f33fcee69232419b3 100644 (file)
@@ -28,7 +28,7 @@ additional-run-action = gnuplot make_ps
 # 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
+delete-files = gnuplot* *inp *.ps *.eps
 
 
 
index 346fc734ccc0c6c4d1264fbc185d7c51e94b0c85..c530fc2001f7383da219eb11486a9edb4415303d 100644 (file)
@@ -56,12 +56,17 @@ class PoissonEquation :  public Equation<dim> {
 template <int dim>
 class PoissonProblem : public ProblemBase<dim> {
   public:
+    enum RefineMode {
+         global, true_error, error_estimator
+    };
+    
     PoissonProblem ();
 
     void clear ();
     void create_new ();
-    void run (unsigned int level);
-    void print_history () const;
+    void run (const unsigned int start_level,
+             const RefineMode refine_mode);
+    void print_history (const RefineMode refine_mode) const;
     
   protected:
     Triangulation<dim> *tria;
@@ -211,6 +216,12 @@ void PoissonProblem<dim>::clear () {
       boundary_values = 0;
     };
 
+  l2_error.clear ();
+  linfty_error.clear ();
+  h1_error.clear ();
+  estimated_error.clear();
+  n_dofs.clear ();
+  
   ProblemBase<dim>::clear ();
 };
 
@@ -232,17 +243,32 @@ void PoissonProblem<dim>::create_new () {
 
 
 template <int dim>
-void PoissonProblem<dim>::run (const unsigned int start_level) {
+void PoissonProblem<dim>::run (const unsigned int start_level,
+                              const RefineMode refine_mode) {
   create_new ();
-  
-  cout << "Refinement level = " << start_level
-       << endl;
-  
-  cout << "    Making grid... ";
+
+  cout << "======================================="
+       << "=======================================" << endl
+       << "===== Doing computation with refinement criterion: ";
+  switch (refine_mode) 
+    {
+      case global:
+           cout << "global";
+           break;
+      case true_error:
+           cout << "true error";
+           break;
+      case error_estimator:
+           cout << "error estimator";
+           break;
+    };
+  cout << endl
+       << "======================================="
+       << "=======================================" << endl;
+  cout << "Making initial grid... " << endl;
   tria->set_boundary (&boundary);
   tria->create_hyper_ball ();
   tria->refine_global (start_level);
-  cout << tria->n_active_cells() << " active cells." << endl;
 
   rhs             = new RHSPoly<dim>();
   boundary_values = new Solution<dim> ();
@@ -251,96 +277,144 @@ void PoissonProblem<dim>::run (const unsigned int start_level) {
   FELinear<dim>                   fe;
   PoissonEquation<dim>            equation (*rhs);
   QGauss3<dim>                    quadrature;
+
+  unsigned int refine_step = 0;
+  while (tria->n_active_cells() < 5000)
+    {
+      cout << "Refinement step " << refine_step
+          << ", using " << tria->n_active_cells() << " active cells."
+          << endl;
+      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;
+      UpdateFlags update_flags = UpdateFlags(update_q_points  | update_gradients |
+                                            update_jacobians | update_JxW_values);
   
-  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;
-  UpdateFlags update_flags = UpdateFlags(update_q_points  | update_gradients |
-                                        update_jacobians | update_JxW_values);
-  
-  ProblemBase<dim>::FunctionMap dirichlet_bc;
-  dirichlet_bc[0] = boundary_values;
-  assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
+      ProblemBase<dim>::FunctionMap dirichlet_bc;
+      dirichlet_bc[0] = boundary_values;
+      assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
 
-  cout << "    Solving..." << endl;
-  solve ();
+      cout << "    Solving..." << endl;
+      solve ();
 
-  Solution<dim> sol;
-  dVector       l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell;
-  dVector       estimated_error_per_cell;
-  QGauss3<dim>  q;
+      Solution<dim> sol;
+      dVector       l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell;
+      dVector       estimated_error_per_cell;
+      QGauss3<dim>  q;
   
-  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 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());
-
-  cout << "    Estimating H1 error... ";
-  KellyErrorEstimator<dim> ee;
-  QSimpson<dim-1> eq;
-  ee.estimate_error (*dof, eq, fe, boundary,
-                    KellyErrorEstimator<dim>::FunctionMap(),
-                    solution,
-                    estimated_error_per_cell);
-  cout << estimated_error_per_cell.l2_norm() << endl;
-  estimated_error.push_back (estimated_error_per_cell.l2_norm());
-
-  dVector l2_error_per_dof, linfty_error_per_dof;
-  dVector h1_error_per_dof, estimated_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_error_per_cell, h1_error_per_dof);
-  dof->distribute_cell_to_dof_vector (estimated_error_per_cell,
-                                     estimated_error_per_dof);
+      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 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());
+
+      cout << "    Estimating H1 error... ";
+      KellyErrorEstimator<dim> ee;
+      QSimpson<dim-1> eq;
+//      ee.estimate_error (*dof, eq, fe, boundary,
+//                      KellyErrorEstimator<dim>::FunctionMap(),
+//                      solution,
+//                      estimated_error_per_cell);
+      cout << estimated_error_per_cell.l2_norm() << endl;
+      estimated_error.push_back (estimated_error_per_cell.l2_norm());
+
+      dVector l2_error_per_dof, linfty_error_per_dof;
+      dVector h1_error_per_dof, estimated_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_error_per_cell, h1_error_per_dof);
+//      dof->distribute_cell_to_dof_vector (estimated_error_per_cell,
+//                                       estimated_error_per_dof);
   
   
-  DataOut<dim> out;
-  fill_data (out);
-  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_error_per_dof, "H1-Error");
-  out.add_data_vector (estimated_error_per_dof, "Estimated Error");
-  string filename = "gnuplot.";
+      DataOut<dim> out;
+      fill_data (out);
+      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_error_per_dof, "H1-Error");
+//      out.add_data_vector (estimated_error_per_dof, "Estimated Error");
+      string filename = "gnuplot.";
 //  string filename = "ee.";
-  filename += ('0'+start_level);
-//  filename += ".inp";
-  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();
+      switch (refine_mode) 
+       {
+         case global:
+               filename += "global.";
+               break;
+         case true_error:
+               filename += "true_error.";
+               break;
+         case error_estimator:
+               filename += "estimated_error.";
+               break;
+       };
+      filename += ('0'+start_level+refine_step);
+      filename += ".inp";                            //*
+      cout << "    Writing error plots to <" << filename << ">..." << endl;
   
-  tria->refine_fixed_number (estimated_error_per_cell, 0.3);
-  tria->execute_refinement ();
-  string ref_filename = "gnuplot.ref.";
-  ref_filename += ('0'+start_level);
-  ofstream oxx(ref_filename.c_str());
-  tria->print_gnuplot (oxx);
+//      ofstream gnuplot(filename.c_str());
+//      out.write_gnuplot (gnuplot);
+//      gnuplot.close ();
+      ofstream ucd(filename.c_str());
+      out.write_ucd (ucd);
+      ucd.close();
+
+      cout << "    Refining triangulation...";
+      switch (refine_mode) 
+       {
+         case global:
+               tria->refine_global (1);
+               break;
+         case true_error:
+               tria->refine_fixed_number (h1_error_per_cell, 0.3);
+               tria->execute_refinement ();
+               break;
+         case error_estimator:
+               tria->execute_refinement ();
+               tria->refine_fixed_number (estimated_error_per_cell, 0.3);
+               break;
+       };
+      cout << endl << endl;
+      cout << endl;
+      ++refine_step;
+    };
   
-  cout << endl;
+  print_history (refine_mode);
 };
 
 
 template <int dim>
-void PoissonProblem<dim>::print_history () const {
-  ofstream out("gnuplot.history");
+void PoissonProblem<dim>::print_history (const RefineMode refine_mode) const {
+  string filename("history.");
+  switch (refine_mode) 
+    {
+      case global:
+           filename += "global.";
+           break;
+      case true_error:
+           filename += "true_error.";
+           break;
+      case error_estimator:
+           filename += "estimated_error.";
+           break;
+    };
+  filename += "gnuplot";
+  
+  cout << endl << "Printing convergence history to" << filename << "..."
+       << endl;
+  ofstream out(filename.c_str());
   out << "# n_dofs    l2_error linfty_error "
       << "h1_error estimated_error"
       << endl;
@@ -384,11 +458,8 @@ void PoissonProblem<dim>::print_history () const {
 int main () {
   PoissonProblem<2> problem;
 
-  for (unsigned int level=2; level<6; ++level)
-    problem.run (level);
-
-  cout << endl << "Printing convergence history to <gnuplot.history>..." << endl;
-  problem.print_history ();
+  problem.run (2, PoissonProblem<2>::global);
+  problem.run (2, PoissonProblem<2>::true_error);
 
   return 0;
 };
index f633e712352eaaedf4ada1b860c7f82bf0d108e3..8fa3f5b4eee12553daefbfef625ee9b59fdc6f1d 100644 (file)
@@ -4,7 +4,20 @@ set data style linespoints
 set logscale xy
 
 set term postscript eps
-set output "error-estimation.eps"
 
-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","gnuplot.history" using 1:7 title "Estimated H1 error"
+
+set output "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"
+
+
+set output "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"
+
+
+
+set output "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"
 
index ee6a9b67138c205b2f769d5cf4cc2f9887d89ae5..772ffb9c4367e8ac92d5f89f33fcee69232419b3 100644 (file)
@@ -28,7 +28,7 @@ additional-run-action = gnuplot make_ps
 # 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
+delete-files = gnuplot* *inp *.ps *.eps
 
 
 
index 346fc734ccc0c6c4d1264fbc185d7c51e94b0c85..c530fc2001f7383da219eb11486a9edb4415303d 100644 (file)
@@ -56,12 +56,17 @@ class PoissonEquation :  public Equation<dim> {
 template <int dim>
 class PoissonProblem : public ProblemBase<dim> {
   public:
+    enum RefineMode {
+         global, true_error, error_estimator
+    };
+    
     PoissonProblem ();
 
     void clear ();
     void create_new ();
-    void run (unsigned int level);
-    void print_history () const;
+    void run (const unsigned int start_level,
+             const RefineMode refine_mode);
+    void print_history (const RefineMode refine_mode) const;
     
   protected:
     Triangulation<dim> *tria;
@@ -211,6 +216,12 @@ void PoissonProblem<dim>::clear () {
       boundary_values = 0;
     };
 
+  l2_error.clear ();
+  linfty_error.clear ();
+  h1_error.clear ();
+  estimated_error.clear();
+  n_dofs.clear ();
+  
   ProblemBase<dim>::clear ();
 };
 
@@ -232,17 +243,32 @@ void PoissonProblem<dim>::create_new () {
 
 
 template <int dim>
-void PoissonProblem<dim>::run (const unsigned int start_level) {
+void PoissonProblem<dim>::run (const unsigned int start_level,
+                              const RefineMode refine_mode) {
   create_new ();
-  
-  cout << "Refinement level = " << start_level
-       << endl;
-  
-  cout << "    Making grid... ";
+
+  cout << "======================================="
+       << "=======================================" << endl
+       << "===== Doing computation with refinement criterion: ";
+  switch (refine_mode) 
+    {
+      case global:
+           cout << "global";
+           break;
+      case true_error:
+           cout << "true error";
+           break;
+      case error_estimator:
+           cout << "error estimator";
+           break;
+    };
+  cout << endl
+       << "======================================="
+       << "=======================================" << endl;
+  cout << "Making initial grid... " << endl;
   tria->set_boundary (&boundary);
   tria->create_hyper_ball ();
   tria->refine_global (start_level);
-  cout << tria->n_active_cells() << " active cells." << endl;
 
   rhs             = new RHSPoly<dim>();
   boundary_values = new Solution<dim> ();
@@ -251,96 +277,144 @@ void PoissonProblem<dim>::run (const unsigned int start_level) {
   FELinear<dim>                   fe;
   PoissonEquation<dim>            equation (*rhs);
   QGauss3<dim>                    quadrature;
+
+  unsigned int refine_step = 0;
+  while (tria->n_active_cells() < 5000)
+    {
+      cout << "Refinement step " << refine_step
+          << ", using " << tria->n_active_cells() << " active cells."
+          << endl;
+      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;
+      UpdateFlags update_flags = UpdateFlags(update_q_points  | update_gradients |
+                                            update_jacobians | update_JxW_values);
   
-  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;
-  UpdateFlags update_flags = UpdateFlags(update_q_points  | update_gradients |
-                                        update_jacobians | update_JxW_values);
-  
-  ProblemBase<dim>::FunctionMap dirichlet_bc;
-  dirichlet_bc[0] = boundary_values;
-  assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
+      ProblemBase<dim>::FunctionMap dirichlet_bc;
+      dirichlet_bc[0] = boundary_values;
+      assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
 
-  cout << "    Solving..." << endl;
-  solve ();
+      cout << "    Solving..." << endl;
+      solve ();
 
-  Solution<dim> sol;
-  dVector       l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell;
-  dVector       estimated_error_per_cell;
-  QGauss3<dim>  q;
+      Solution<dim> sol;
+      dVector       l2_error_per_cell, linfty_error_per_cell, h1_error_per_cell;
+      dVector       estimated_error_per_cell;
+      QGauss3<dim>  q;
   
-  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 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());
-
-  cout << "    Estimating H1 error... ";
-  KellyErrorEstimator<dim> ee;
-  QSimpson<dim-1> eq;
-  ee.estimate_error (*dof, eq, fe, boundary,
-                    KellyErrorEstimator<dim>::FunctionMap(),
-                    solution,
-                    estimated_error_per_cell);
-  cout << estimated_error_per_cell.l2_norm() << endl;
-  estimated_error.push_back (estimated_error_per_cell.l2_norm());
-
-  dVector l2_error_per_dof, linfty_error_per_dof;
-  dVector h1_error_per_dof, estimated_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_error_per_cell, h1_error_per_dof);
-  dof->distribute_cell_to_dof_vector (estimated_error_per_cell,
-                                     estimated_error_per_dof);
+      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 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());
+
+      cout << "    Estimating H1 error... ";
+      KellyErrorEstimator<dim> ee;
+      QSimpson<dim-1> eq;
+//      ee.estimate_error (*dof, eq, fe, boundary,
+//                      KellyErrorEstimator<dim>::FunctionMap(),
+//                      solution,
+//                      estimated_error_per_cell);
+      cout << estimated_error_per_cell.l2_norm() << endl;
+      estimated_error.push_back (estimated_error_per_cell.l2_norm());
+
+      dVector l2_error_per_dof, linfty_error_per_dof;
+      dVector h1_error_per_dof, estimated_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_error_per_cell, h1_error_per_dof);
+//      dof->distribute_cell_to_dof_vector (estimated_error_per_cell,
+//                                       estimated_error_per_dof);
   
   
-  DataOut<dim> out;
-  fill_data (out);
-  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_error_per_dof, "H1-Error");
-  out.add_data_vector (estimated_error_per_dof, "Estimated Error");
-  string filename = "gnuplot.";
+      DataOut<dim> out;
+      fill_data (out);
+      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_error_per_dof, "H1-Error");
+//      out.add_data_vector (estimated_error_per_dof, "Estimated Error");
+      string filename = "gnuplot.";
 //  string filename = "ee.";
-  filename += ('0'+start_level);
-//  filename += ".inp";
-  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();
+      switch (refine_mode) 
+       {
+         case global:
+               filename += "global.";
+               break;
+         case true_error:
+               filename += "true_error.";
+               break;
+         case error_estimator:
+               filename += "estimated_error.";
+               break;
+       };
+      filename += ('0'+start_level+refine_step);
+      filename += ".inp";                            //*
+      cout << "    Writing error plots to <" << filename << ">..." << endl;
   
-  tria->refine_fixed_number (estimated_error_per_cell, 0.3);
-  tria->execute_refinement ();
-  string ref_filename = "gnuplot.ref.";
-  ref_filename += ('0'+start_level);
-  ofstream oxx(ref_filename.c_str());
-  tria->print_gnuplot (oxx);
+//      ofstream gnuplot(filename.c_str());
+//      out.write_gnuplot (gnuplot);
+//      gnuplot.close ();
+      ofstream ucd(filename.c_str());
+      out.write_ucd (ucd);
+      ucd.close();
+
+      cout << "    Refining triangulation...";
+      switch (refine_mode) 
+       {
+         case global:
+               tria->refine_global (1);
+               break;
+         case true_error:
+               tria->refine_fixed_number (h1_error_per_cell, 0.3);
+               tria->execute_refinement ();
+               break;
+         case error_estimator:
+               tria->execute_refinement ();
+               tria->refine_fixed_number (estimated_error_per_cell, 0.3);
+               break;
+       };
+      cout << endl << endl;
+      cout << endl;
+      ++refine_step;
+    };
   
-  cout << endl;
+  print_history (refine_mode);
 };
 
 
 template <int dim>
-void PoissonProblem<dim>::print_history () const {
-  ofstream out("gnuplot.history");
+void PoissonProblem<dim>::print_history (const RefineMode refine_mode) const {
+  string filename("history.");
+  switch (refine_mode) 
+    {
+      case global:
+           filename += "global.";
+           break;
+      case true_error:
+           filename += "true_error.";
+           break;
+      case error_estimator:
+           filename += "estimated_error.";
+           break;
+    };
+  filename += "gnuplot";
+  
+  cout << endl << "Printing convergence history to" << filename << "..."
+       << endl;
+  ofstream out(filename.c_str());
   out << "# n_dofs    l2_error linfty_error "
       << "h1_error estimated_error"
       << endl;
@@ -384,11 +458,8 @@ void PoissonProblem<dim>::print_history () const {
 int main () {
   PoissonProblem<2> problem;
 
-  for (unsigned int level=2; level<6; ++level)
-    problem.run (level);
-
-  cout << endl << "Printing convergence history to <gnuplot.history>..." << endl;
-  problem.print_history ();
+  problem.run (2, PoissonProblem<2>::global);
+  problem.run (2, PoissonProblem<2>::true_error);
 
   return 0;
 };
index f633e712352eaaedf4ada1b860c7f82bf0d108e3..8fa3f5b4eee12553daefbfef625ee9b59fdc6f1d 100644 (file)
@@ -4,7 +4,20 @@ set data style linespoints
 set logscale xy
 
 set term postscript eps
-set output "error-estimation.eps"
 
-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","gnuplot.history" using 1:7 title "Estimated H1 error"
+
+set output "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"
+
+
+set output "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"
+
+
+
+set output "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"
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.