From: Wolfgang Bangerth <bangerth@math.tamu.edu> Date: Mon, 16 Sep 2002 13:27:24 +0000 (+0000) Subject: New test. X-Git-Tag: v8.0.0~17552 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fb6d287bef602700463770b1f45908db9e4ede2c;p=dealii.git New test. git-svn-id: https://svn.dealii.org/trunk@6426 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/fe/Makefile b/tests/fe/Makefile index 03f3ae030d..5b272adc68 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -37,11 +37,12 @@ non_primitive_1.exe : non_primitive_1.go $(libraries) numbering.exe : numbering.go $(libraries) shapes.exe : shapes.go $(libraries) transfer.exe : transfer.go $(libraries) +up_and_down.exe : up_and_down.go $(libraries) tests = derivatives fe_data_test fe_traits fe_tools_test mapping \ mapping_c1 shapes numbering mapping_q1_eulerian \ transfer internals \ - non_primitive_1 nedelec nedelec_2 nedelec_3 + non_primitive_1 nedelec nedelec_2 nedelec_3 up_and_down ############################################################ diff --git a/tests/fe/up_and_down.cc b/tests/fe/up_and_down.cc new file mode 100644 index 0000000000..c4dc124d83 --- /dev/null +++ b/tests/fe/up_and_down.cc @@ -0,0 +1,233 @@ +// $Id$ +// (c) Wolfgang Bangerth +// +// Try something remarkably simple: take some arbitrary (discrete) +// function and transfer it to a coarser grid. transfer it back to the +// fine grid. store it and call it X. transfer down to the coarse grid +// and back up again. Should still be X, no? +// +// do this for a couple of possible finite elements, in particular to +// composed and vector valued ones, to stress test their embedding and +// restriction matrices. + +#include <base/quadrature_lib.h> +#include <base/logstream.h> +#include <lac/vector.h> +#include <grid/tria.h> +#include <grid/tria_iterator.h> +#include <dofs/dof_accessor.h> +#include <grid/grid_generator.h> +#include <grid/grid_tools.h> +#include <fe/fe_nedelec.h> +#include <fe/fe_q.h> +#include <fe/fe_dgq.h> +#include <fe/fe_dgp.h> +#include <fe/fe_system.h> +#include <fe/fe_values.h> + +#include <vector> +#include <fstream> +#include <string> + +#define PRECISION 2 + + +template <int dim> +Point<dim> transform (const Point<dim> p) +{ + switch (dim) + { + case 1: + return p; + case 2: + return Point<dim>(p(0)*(1+p(1)), p(1)*(1+p(0))); + case 3: + return Point<dim>(p(0)*(1+p(1))*(1+p(2)), + p(1)*(1+p(0))*(1+p(2)), + p(2)*(1+p(0))*(1+p(1))); + default: + Assert (false, ExcNotImplemented()); + return Point<dim>(); + }; +}; + + +template <int dim> +void check_element (const Triangulation<dim> &tr, + const FiniteElement<dim> &fe) +{ + DoFHandler<dim> dof_handler(tr); + dof_handler.distribute_dofs (fe); + + // create a mostly arbitrary + // function plus a trend on this + // grid + Vector<double> tmp(dof_handler.n_dofs()); + for (unsigned int i=0; i<tmp.size(); ++i) + tmp(i) = i;//(i + 13*i%17); + + // restrict this function to the + // next coarser level and + // distribute it again to the + // higher level + Vector<double> x(tmp.size()); + Vector<double> v(fe.dofs_per_cell); + for (typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin(); + cell!=dof_handler.end(); ++cell) + if (cell->has_children() && + cell->child(0)->active()) + { + // first make sure that what + // we do is reasonable. for + // this, _all_ children have + // to be active, not only + // some of them + for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c) + Assert (cell->child(c)->active(), ExcInternalError()); + + // then restrict and prolong + cell->get_interpolated_dof_values (tmp, v); + cell->set_dof_values_by_interpolation (v, x); + }; + + // now x is a function on the fine + // grid that is representable on + // the coarse grid. so another + // cycle should not alter it any + // more: + Vector<double> x2(x.size()); + for (typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin(); + cell!=dof_handler.end(); ++cell) + if (cell->has_children() && + cell->child(0)->active()) + { + cell->get_interpolated_dof_values (x, v); + cell->set_dof_values_by_interpolation (v, x2); + }; + + // then check that this is so: + x2 -= x; + const double relative_residual = (x2.l2_norm() / x.l2_norm()); + + const double threshold = 1e-6; + deallog << ", dofs_per_cell=" << fe.dofs_per_cell + << "; relative residual: " + << (relative_residual<threshold ? "ok" : "botched up!") + << std::endl; + + Assert (relative_residual < threshold*x.l2_norm(), ExcInternalError()); +}; + + +template <int dim> +void test () +{ + // make a coarse triangulation as a + // hypercube. if in more than 1d, + // distort it so that it is no more + // an affine image of the + // hypercube, to make things more + // difficult. then refine it twice + Triangulation<dim> tr; + GridGenerator::hyper_cube(tr, 0., 1.); + GridTools::transform (&transform<dim>, tr); + tr.refine_global (2); + + // now for a list of finite + // elements, for which we want to + // test. we happily waste tons of + // memory here, but who cares... + const FiniteElement<dim> *fe_list[] + = + { + // first for some scalar + // elements: + + // FE_Q + new FE_Q<dim>(1), + new FE_Q<dim>(2), + (dim<3 ? new FE_Q<dim>(3) : 0), + (dim<2 ? new FE_Q<dim>(4) : 0), + + // FE_DGQ + new FE_DGQ<dim>(0), + new FE_DGQ<dim>(1), + new FE_DGQ<dim>(2), + (dim<3 ? new FE_DGQ<dim>(3) : 0), + (dim<3 ? new FE_DGQ<dim>(4) : 0), + + // FE_DGP does not + // presently have the + // restriction matrices + // implemented, so comment + // out the following + // tests... +// new FE_DGP<dim>(0), +// new FE_DGP<dim>(1), +// new FE_DGP<dim>(2), +// new FE_DGP<dim>(3), + + // a non-primitive FE +// (dim != 1 ? new FE_Nedelec<dim>(1) : 0), + + // some composed elements + // of increasing + // complexity, to check the + // logics by which the + // matrices of the composed + // elements are assembled + // from those of the base + // elements. note that some + // of the base elements are + // additive, some not, so + // the result will be an + // element that is mixed in + // this respect + new FESystem<dim> (FE_Q<dim>(2), 2), + new FESystem<dim> (FE_Q<dim>(1), 2, + FE_DGQ<dim>(2), 2), + new FESystem<dim> (FE_Q<dim>(1), 2, + FE_DGQ<dim>(2), 2, + FE_DGQ<dim>(0), 1), + new FESystem<dim> (FE_Q<dim>(1), 2, + FESystem<dim> (FE_Q<dim>(1), 2, + FE_DGQ<dim>(2), 2, + FE_DGQ<dim>(2), 1), 2, + FE_DGQ<dim>(0), 1), + new FESystem<dim> (FE_Q<dim>(1), 2, + FESystem<dim> (FE_Q<dim>(1), 2, + FE_DGQ<dim>(2), 2, + FESystem<dim>(FE_DGQ<dim>(0), + 3), + 1), 2, + FE_DGQ<dim>(0), 1), + }; + + for (unsigned int i=0; i<sizeof(fe_list)/sizeof(fe_list[0]); ++i) + if (fe_list[i] != 0) + { + deallog << dim << "d, uniform grid, fe #" << i; + check_element (tr, *fe_list[i]); + }; +}; + + + +int +main() +{ + std::ofstream logfile ("up_and_down.output"); + logfile.precision (PRECISION); + logfile.setf(std::ios::fixed); + deallog.attach(logfile); + deallog.depth_console(0); + + test<1>(); + test<2>(); + test<3>(); + + return 0; +} + + + diff --git a/tests/fe/up_and_down.checked b/tests/fe/up_and_down.checked new file mode 100644 index 0000000000..f687a5662b --- /dev/null +++ b/tests/fe/up_and_down.checked @@ -0,0 +1,38 @@ + +DEAL::1d, uniform grid, fe #0, dofs_per_cell=2; relative residual: ok +DEAL::1d, uniform grid, fe #1, dofs_per_cell=3; relative residual: ok +DEAL::1d, uniform grid, fe #2, dofs_per_cell=4; relative residual: ok +DEAL::1d, uniform grid, fe #3, dofs_per_cell=5; relative residual: ok +DEAL::1d, uniform grid, fe #4, dofs_per_cell=1; relative residual: ok +DEAL::1d, uniform grid, fe #5, dofs_per_cell=2; relative residual: ok +DEAL::1d, uniform grid, fe #6, dofs_per_cell=3; relative residual: ok +DEAL::1d, uniform grid, fe #7, dofs_per_cell=4; relative residual: ok +DEAL::1d, uniform grid, fe #8, dofs_per_cell=5; relative residual: ok +DEAL::1d, uniform grid, fe #9, dofs_per_cell=6; relative residual: ok +DEAL::1d, uniform grid, fe #10, dofs_per_cell=10; relative residual: ok +DEAL::1d, uniform grid, fe #11, dofs_per_cell=11; relative residual: ok +DEAL::1d, uniform grid, fe #12, dofs_per_cell=31; relative residual: ok +DEAL::1d, uniform grid, fe #13, dofs_per_cell=31; relative residual: ok +DEAL::2d, uniform grid, fe #0, dofs_per_cell=4; relative residual: ok +DEAL::2d, uniform grid, fe #1, dofs_per_cell=9; relative residual: ok +DEAL::2d, uniform grid, fe #2, dofs_per_cell=16; relative residual: ok +DEAL::2d, uniform grid, fe #4, dofs_per_cell=1; relative residual: ok +DEAL::2d, uniform grid, fe #5, dofs_per_cell=4; relative residual: ok +DEAL::2d, uniform grid, fe #6, dofs_per_cell=9; relative residual: ok +DEAL::2d, uniform grid, fe #7, dofs_per_cell=16; relative residual: ok +DEAL::2d, uniform grid, fe #8, dofs_per_cell=25; relative residual: ok +DEAL::2d, uniform grid, fe #9, dofs_per_cell=18; relative residual: ok +DEAL::2d, uniform grid, fe #10, dofs_per_cell=26; relative residual: ok +DEAL::2d, uniform grid, fe #11, dofs_per_cell=27; relative residual: ok +DEAL::2d, uniform grid, fe #12, dofs_per_cell=79; relative residual: ok +DEAL::2d, uniform grid, fe #13, dofs_per_cell=67; relative residual: ok +DEAL::3d, uniform grid, fe #0, dofs_per_cell=8; relative residual: ok +DEAL::3d, uniform grid, fe #1, dofs_per_cell=27; relative residual: ok +DEAL::3d, uniform grid, fe #4, dofs_per_cell=1; relative residual: ok +DEAL::3d, uniform grid, fe #5, dofs_per_cell=8; relative residual: ok +DEAL::3d, uniform grid, fe #6, dofs_per_cell=27; relative residual: ok +DEAL::3d, uniform grid, fe #9, dofs_per_cell=54; relative residual: ok +DEAL::3d, uniform grid, fe #10, dofs_per_cell=70; relative residual: ok +DEAL::3d, uniform grid, fe #11, dofs_per_cell=71; relative residual: ok +DEAL::3d, uniform grid, fe #12, dofs_per_cell=211; relative residual: ok +DEAL::3d, uniform grid, fe #13, dofs_per_cell=163; relative residual: ok