From: wolf Date: Fri, 2 Apr 2004 03:58:58 +0000 (+0000) Subject: New test. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=65e9a6f24b916cdb999e7d73fd4fb2217aaeda89;p=dealii-svn.git New test. git-svn-id: https://svn.dealii.org/trunk@8957 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/bits/error_estimator_01.cc b/tests/bits/error_estimator_01.cc new file mode 100644 index 0000000000..0dca7448f7 --- /dev/null +++ b/tests/bits/error_estimator_01.cc @@ -0,0 +1,214 @@ +//---------------------------- error_estimator_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- error_estimator_01.cc --------------------------- + +// check that the kelly error estimator returns the same result whether we +// compute all indicators at once, or on different subdomains separately + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +class MySquareFunction : public Function +{ + public: + MySquareFunction () : Function(2) {}; + + virtual double value (const Point &p, + const unsigned int component) const + { return (component+1)*p.square(); }; + + virtual void vector_value (const Point &p, + Vector &values) const + { values(0) = value(p,0); + values(1) = value(p,1); }; +}; + + + +template +Quadrature & +get_q_face (Function&) +{ + static QGauss4 q; + return q; +} + + +Quadrature<0> & +get_q_face (Function<1>&) +{ + Quadrature<0> *q = 0; + return *q; +} + + + +template +void make_mesh (Triangulation &tria) +{ + + GridGenerator::hyper_cube(tria, -1, 1); + + // refine the mesh in a random way so as to + // generate as many cells with + // hanging nodes as possible + tria.refine_global (4-dim); + const double steps[3] = { 9, 7, 3 }; + for (unsigned int i=0; i::active_cell_iterator + cell = tria.begin_active(); + for (unsigned int index=0; cell != tria.end(); ++cell, ++index) + if (index % (3*dim) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement (); + } + + // we now have a number of cells, + // flag them with some subdomain + // ids based on their position, in + // particular we take the quadrant + // (octant) + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end (); + for (; cell!=endc; ++cell) + { + unsigned int subdomain = 0; + for (unsigned int d=0; dcenter()(d) > 0) + subdomain |= (1<set_subdomain_id (subdomain); + } +} + + + + +template +void +check () +{ + Functions::CosineFunction function; + + Triangulation tria; + make_mesh (tria); + + FE_Q element(3); + DoFHandler dof(tria); + dof.distribute_dofs(element); + + MappingQ mapping(3); + Quadrature &q_face = get_q_face(function); + + std::map*> neumann_bc; + neumann_bc[0] = &function; + + Vector v (dof.n_dofs()); + VectorTools::interpolate (mapping, dof, function, v); + + Vector error1 (tria.n_active_cells()); + Vector error2 (tria.n_active_cells()); + + // compute error by looking at all cells at + // once and output this as the base line + // results. scale results so that they show + // up in the output file as a reasonable + // number + KellyErrorEstimator::estimate (mapping, dof, q_face, neumann_bc, + v, error1); + const double scaling_factor = 10000.*error1.size()/error1.l1_norm(); + error1 *= scaling_factor; + + deallog << "Estimated error indicators:" << std::endl; + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +}