From: Wolfgang Bangerth Date: Fri, 28 Aug 1998 12:08:42 +0000 (+0000) Subject: Rename target file. X-Git-Tag: v8.0.0~22718 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8f788994f41362dc93e82b8fa4ad21975a4169d6;p=dealii.git Rename target file. git-svn-id: https://svn.dealii.org/trunk@536 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/multigrid/Makefile b/deal.II/deal.II/Attic/examples/multigrid/Makefile index 022319dc4e..884d655e97 100644 --- a/deal.II/deal.II/Attic/examples/multigrid/Makefile +++ b/deal.II/deal.II/Attic/examples/multigrid/Makefile @@ -4,7 +4,7 @@ # Template for makefiles for the examples subdirectory. In principle, # everything should be done automatically if you set the target file # here correctly: -target = convergence +target = multigrid # All dependencies between files should be updated by the included # file Makefile.dep if necessary. Object files are compiled into diff --git a/tests/big-tests/multigrid/convergence.cc b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc similarity index 98% rename from tests/big-tests/multigrid/convergence.cc rename to deal.II/deal.II/Attic/examples/multigrid/multigrid.cc index 0aa8f561f7..b163635c41 100644 --- a/tests/big-tests/multigrid/convergence.cc +++ b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc @@ -3,9 +3,9 @@ #include -#include +#include #include -#include +#include #include #include #include @@ -14,7 +14,6 @@ #include #include #include -#include #include #include @@ -55,7 +54,7 @@ class PoissonEquation : public Equation { template -class PoissonProblem : public ProblemBase { +class PoissonProblem { public: PoissonProblem (unsigned int order); @@ -66,7 +65,7 @@ class PoissonProblem : public ProblemBase { protected: Triangulation *tria; - DoFHandler *dof; + MGDoFHandler *dof; Function *rhs; Function *boundary_values; @@ -231,7 +230,7 @@ void PoissonProblem::create_new () { clear (); tria = new Triangulation(); - dof = new DoFHandler (tria); + dof = new MGDoFHandler (tria); set_tria_and_dof (tria, dof); }; diff --git a/tests/big-tests/multigrid/Makefile b/tests/big-tests/multigrid/Makefile index 022319dc4e..884d655e97 100644 --- a/tests/big-tests/multigrid/Makefile +++ b/tests/big-tests/multigrid/Makefile @@ -4,7 +4,7 @@ # Template for makefiles for the examples subdirectory. In principle, # everything should be done automatically if you set the target file # here correctly: -target = convergence +target = multigrid # All dependencies between files should be updated by the included # file Makefile.dep if necessary. Object files are compiled into diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc new file mode 100644 index 0000000000..b163635c41 --- /dev/null +++ b/tests/big-tests/multigrid/multigrid.cc @@ -0,0 +1,510 @@ +/* $Id$ */ +/* Copyright W. Bangerth, University of Heidelberg, 1998 */ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + + + + +template +class PoissonEquation : public Equation { + public: + PoissonEquation (const Function &rhs) : + Equation(1), + right_hand_side (rhs) {}; + + virtual void assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + protected: + const Function &right_hand_side; +}; + + + + + + +template +class PoissonProblem { + public: + PoissonProblem (unsigned int order); + + void clear (); + void create_new (); + int run (unsigned int level); + void print_history (string filename) const; + + protected: + Triangulation *tria; + MGDoFHandler *dof; + + Function *rhs; + Function *boundary_values; + + vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; + vector n_dofs; + + unsigned int order; +}; + + + + + +/** + Right hand side constructed such that the exact solution is + $sin(2 pi x) + sin(2 pi y)$ + */ +template +class RHSPoly : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + +template +class Solution : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; + /** + * Return the gradient of the function + * at the given point. + */ + virtual Point gradient (const Point &p) const; +}; + + + + +double RHSPoly<2>::operator () (const Point<2> &p) const { + 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 { + 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 { + 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)); +}; + + + + + + +void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point +void PoissonEquation::assemble (dFMatrix &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + +template +void PoissonEquation::assemble (dVector &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + + + + + + + +template +PoissonProblem::PoissonProblem (unsigned int order) : + tria(0), dof(0), rhs(0), + boundary_values(0), order(order) {}; + + + + +template +void PoissonProblem::clear () { + if (tria != 0) { + delete tria; + tria = 0; + }; + + if (dof != 0) { + delete dof; + dof = 0; + }; + + if (rhs != 0) + { + delete rhs; + rhs = 0; + }; + + if (boundary_values != 0) + { + delete boundary_values; + boundary_values = 0; + }; + + ProblemBase::clear (); +}; + + + + +template +void PoissonProblem::create_new () { + clear (); + + tria = new Triangulation(); + dof = new MGDoFHandler (tria); + set_tria_and_dof (tria, dof); +}; + + + + +template +int PoissonProblem::run (const unsigned int level) { + create_new (); + + cout << "Refinement level = " << level + << ", using elements of type <"; + switch (order) + { + case 0: + cout << "criss-cross"; + break; + default: + cout << "Lagrange-" << order; + break; + }; + cout << ">" << endl; + + cout << " Making grid... "; + 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 Solution (); + + + FiniteElement *fe; + PoissonEquation equation (*rhs); + Quadrature *quadrature; + Quadrature *boundary_quadrature; + switch (order) { + case 0: + fe = new FECrissCross(); + quadrature = new QCrissCross1(); + boundary_quadrature = new QGauss2(); + break; + 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); + 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_JxW_values); + + ProblemBase::FunctionMap dirichlet_bc; + dirichlet_bc[0] = boundary_values; + assemble (equation, *quadrature, *fe, update_flags, dirichlet_bc); + + cout << " Solving..." << endl; + solve (); + + Solution sol; + dVector l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell; + dVector h1_seminorm_error_per_cell, h1_error_per_cell; + + cout << " Calculating L1 error... "; + VectorTools::integrate_difference (*dof_handler, + solution, sol, + l1_error_per_cell, + *quadrature, *fe, L1_norm); + cout << l1_error_per_cell.l1_norm() << endl; + l1_error.push_back (l1_error_per_cell.l1_norm()); + + cout << " Calculating L2 error... "; + VectorTools::integrate_difference (*dof_handler, + solution, sol, + l2_error_per_cell, + *quadrature, *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... "; + VectorTools::integrate_difference (*dof_handler, + solution, sol, + linfty_error_per_cell, + *quadrature, *fe, Linfty_norm); + cout << linfty_error_per_cell.linfty_norm() << endl; + linfty_error.push_back (linfty_error_per_cell.linfty_norm()); + + cout << " Calculating H1-seminorm error... "; + VectorTools::integrate_difference (*dof_handler, + solution, sol, + h1_seminorm_error_per_cell, + *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()); + + cout << " Calculating H1 error... "; + VectorTools::integrate_difference (*dof_handler, + solution, sol, + h1_error_per_cell, + *quadrature, *fe, H1_norm); + cout << h1_error_per_cell.l2_norm() << endl; + h1_error.push_back (h1_error_per_cell.l2_norm()); + + 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 (h1_error_per_cell, h1_error_per_dof); + +// 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 o(filename.c_str()); + fill_data (out); + out.add_data_vector (l1_error_per_dof, "L1-Error"); + out.add_data_vector (l2_error_per_dof, "L2-Error"); + out.add_data_vector (linfty_error_per_dof, "Linfty-Error"); + out.add_data_vector (h1_seminorm_error_per_dof, "H1-seminorm-Error"); + out.add_data_vector (h1_error_per_dof, "H1-Error"); + out.write_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 (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; ih/2:" << endl; + cout << " L1 error : " << 1./average_l1 << endl + << " L2 error : " << 1./average_l2 << endl + << " Linfty error : " << 1./average_linfty << endl + << " H1 seminorm error: " << 1./average_h1_semi << endl + << " H1 error : " << 1./average_h1 << endl; + cout << "==========================================================\n"; + cout << "==========================================================\n"; +}; + + + + +int main () { + for (unsigned int order=0; 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 0: + filename = "criss_cross"; + break; + 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; +};