From 4a6b7b632327f2c08e454804d76a098b7d8767b1 Mon Sep 17 00:00:00 2001 From: Alexander Grayver Date: Fri, 18 Sep 2015 15:46:08 +0200 Subject: [PATCH] Add new tests --- tests/fe/fe_project_2d.cc | 265 ++++++++++++++++++++++++++++++++++ tests/fe/fe_project_2d.output | 79 ++++++++++ 2 files changed, 344 insertions(+) create mode 100644 tests/fe/fe_project_2d.cc create mode 100644 tests/fe/fe_project_2d.output diff --git a/tests/fe/fe_project_2d.cc b/tests/fe/fe_project_2d.cc new file mode 100644 index 0000000000..861b613525 --- /dev/null +++ b/tests/fe/fe_project_2d.cc @@ -0,0 +1,265 @@ +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +/* + * This program projects a function into FE spaces defined on meshes + * consisting of: rectangular cells, affine cells, non-affine cells + * + * The error, curl and divergence are then numerically calculated on a + * series of globally refined meshes and get output. + * + * Among FE spaces tested are: FE_ABF, FE_Nedelec, FE_RaviartThomas, + * FE_Q^dim (via FESystem) + * + * Alexander Grayver + */ + +using namespace dealii; + +static const Point<2> vertices_nonaffine[] = +{ + Point<2> (-1., -1.), + Point<2> (0., -1.), + Point<2> (1., -1.), + + Point<2> (-1., 0.), + Point<2> (0.3, 0.3), + Point<2> (1., 0.), + + Point<2> (-1., 1.), + Point<2> (0., 1.), + Point<2> (1., 1.), +}; + +static const Point<2> vertices_affine[] = +{ + Point<2> (-1.4, -1.), + Point<2> (-0.4, -1.), + Point<2> (0.6, -1.), + + Point<2> (-1.2, 0.), + Point<2> (-0.2, 0.), + Point<2> (0.8, 0.), + + Point<2> (-1., 1.), + Point<2> (0., 1.), + Point<2> (1., 1.), +}; + +static const Point<2> vertices_rectangular[] = +{ + Point<2> (-1., -1.), + Point<2> (0., -1.), + Point<2> (1., -1.), + + Point<2> (-1., 0.), + Point<2> (0., 0.), + Point<2> (1., 0.), + + Point<2> (-1., 1.), + Point<2> (0., 1.), + Point<2> (1., 1.), +}; + +static const unsigned n_vertices = sizeof(vertices_rectangular) / sizeof(vertices_rectangular[0]); +static const unsigned n_cells = 4; + +template +class VectorFunction : public Function +{ +public: + VectorFunction() : Function(dim) {} + virtual double value (const Point &p, const unsigned int component) const; + virtual void vector_value(const Point &p, Vector &values) const; +}; + +template<> +double VectorFunction<2>::value(const Point<2> &p, const unsigned int component) const +{ + Assert (component < 2, ExcIndexRange (component, 0, 1)); + + const double PI = numbers::PI; + double val = 0.0; + switch (component) + { + case 0: + val = cos(PI*p(0))*sin(PI*p(1)); + break; + case 1: + val = -sin(PI*p(0))*cos(PI*p(1)); + break; + } + return val; +} + +template +void VectorFunction::vector_value(const Point &p, Vector &values) const +{ + for (int i = 0; i < dim; ++i) + values(i) = value(p, i); +} + +void create_tria(Triangulation<2> &triangulation, const Point<2> *vertices_parallelograms) +{ + const std::vector > vertices (&vertices_parallelograms[0], + &vertices_parallelograms[n_vertices]); + + static const int cell_vertices[][GeometryInfo<2>::vertices_per_cell] = + { + {0, 1, 3, 4}, + {1, 2, 4, 5}, + {3, 4, 6, 7}, + {4, 5, 7, 8} + }; + + std::vector > cells (n_cells, CellData<2>()); + for (unsigned i = 0; i::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, cells, SubCellData()); + triangulation.refine_global(1); +} + +template +void test(const FiniteElement &fe, unsigned n_cycles, bool global, const Point *vertices_parallelograms) +{ + deallog << "dim: " << dim << "\t" << fe.get_name() << std::endl; + deallog << "DoFs\t\t||u-u_h||\tcurl(u_h)\tdiv(u_h)" << std::endl; + + Triangulation triangulation; + create_tria(triangulation, vertices_parallelograms); + + DoFHandler dof_handler(triangulation); + + VectorFunction fe_function; + const FEValuesExtractors::Vector vec (0); + const QGauss quadrature (fe.degree+1); + const unsigned int n_q_points = quadrature.size (); + MappingQ mapping(1); + //MappingQ1 mapping; + std::vector div_v(n_q_points); + std::vector::curl_type> curl_v(n_q_points); + + for (unsigned cycle = 0; cycle < n_cycles; ++cycle) + { + dof_handler.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + Vector v(dof_handler.n_dofs()); + VectorTools::project(mapping, dof_handler, constraints, quadrature, fe_function, v); + + Vector diff(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, dof_handler, v, fe_function, diff, + QGauss(fe.degree + 2), VectorTools::L2_norm); + + typename FEValuesViews::Vector::curl_type total_curl; + double total_div = 0; + total_curl *= 0.; + + FEValues fe_values (mapping, fe, quadrature, update_JxW_values | update_quadrature_points | update_values | update_gradients); + unsigned int cell_index = 0; + + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); + cell != dof_handler.end (); ++cell, ++cell_index) + { + fe_values.reinit (cell); + const std::vector &JxW_values = fe_values.get_JxW_values (); + fe_values[vec].get_function_divergences (v, div_v); + fe_values[vec].get_function_curls (v, curl_v); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + total_div += JxW_values[q_point] * div_v[q_point]; + total_curl += JxW_values[q_point] * curl_v[q_point]; + } + } + + deallog << dof_handler.n_dofs() << "\t\t" + << diff.l2_norm() << "\t" + << total_curl << "\t" + << total_div << std::endl; + + if (global) + triangulation.refine_global(); + else + { + GridRefinement::refine_and_coarsen_fixed_number(triangulation, diff, 0.3, 0.0); + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement(); + } + } +} + +int main () +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(6); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double (1e-8); + + const static unsigned dim = 2; + unsigned order = 1; + unsigned n_cycles = 4; + + deallog << "2d\nRectangular grid:\n"; + + const Point *vertices = &vertices_rectangular[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); + test(FE_ABF (order), n_cycles, true, vertices); + + deallog << "\nAffine grid:\n"; + + vertices = &vertices_affine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); + test(FE_ABF (order), n_cycles, true, vertices); + + deallog << "\nNon-affine grid:\n"; + + vertices = &vertices_nonaffine[0]; + test(FE_Nedelec (order), n_cycles, true, vertices); + test(FE_RaviartThomas (order), n_cycles, true, vertices); + test(FESystem (FE_Q(order), dim), n_cycles, true, vertices); + test(FE_ABF (order), n_cycles, true, vertices); + + deallog << std::endl; +} diff --git a/tests/fe/fe_project_2d.output b/tests/fe/fe_project_2d.output new file mode 100644 index 0000000000..b30affbef5 --- /dev/null +++ b/tests/fe/fe_project_2d.output @@ -0,0 +1,79 @@ +DEAL::2d +Rectangular grid: +dim: 2 FE_Nedelec<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::144 0.127197 0.000000 0 +DEAL::544 0.032385 -0.000000 0 +DEAL::2112 0.008122 0.000000 0 +DEAL::8320 0.002032 -0.000000 0 +DEAL::dim: 2 FE_RaviartThomas<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::144 0.127676 0.000000 0 +DEAL::544 0.032425 0.000000 0 +DEAL::2112 0.008124 0.000000 0 +DEAL::8320 0.002032 0.000000 0 +DEAL::dim: 2 FESystem<2>[FE_Q<2>(1)^2] +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::50 0.229086 -0.000000 0 +DEAL::162 0.049159 -0.000000 0 +DEAL::578 0.011701 0.000000 0 +DEAL::2178 0.002887 0.000000 0 +DEAL::dim: 2 FE_ABF<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::208 0.126166 0.000000 0 +DEAL::800 0.032829 -0.000000 0 +DEAL::3136 0.008676 -0.000000 0 +DEAL::12416 0.002544 -0.000000 0 +DEAL:: +Affine grid: +dim: 2 FE_Nedelec<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::144 0.126899 -0.010605 0 +DEAL::544 0.032013 -0.001829 0 +DEAL::2112 0.007983 -0.000255 0 +DEAL::8320 0.001993 -0.000033 0 +DEAL::dim: 2 FE_RaviartThomas<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::144 0.142775 -0.000000 0.018305 +DEAL::544 0.036564 0.000000 0.003444 +DEAL::2112 0.009186 0.000000 0.000490 +DEAL::8320 0.002299 0.000000 0.000063 +DEAL::dim: 2 FESystem<2>[FE_Q<2>(1)^2] +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::50 0.250113 -0.021467 -0.051342 +DEAL::162 0.052634 -0.002128 -0.005106 +DEAL::578 0.012425 -0.000250 -0.000600 +DEAL::2178 0.003058 -0.000031 -0.000074 +DEAL::dim: 2 FE_ABF<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::208 0.150632 0.000000 0 +DEAL::800 0.038793 -0.000000 0 +DEAL::3136 0.010523 -0.000000 0 +DEAL::12416 0.003304 0.000000 0 +DEAL:: +Non-affine grid: +dim: 2 FE_Nedelec<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::144 0.153434 0.020872 0 +DEAL::544 0.037846 0.005462 0 +DEAL::2112 0.009393 0.000834 0 +DEAL::8320 0.002340 0.000110 0 +DEAL::dim: 2 FE_RaviartThomas<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::144 0.163372 -0.033531 0 +DEAL::544 0.042792 -0.006665 0 +DEAL::2112 0.010774 -0.001430 0 +DEAL::8320 0.002695 -0.000330 0 +DEAL::dim: 2 FESystem<2>[FE_Q<2>(1)^2] +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::50 0.303441 -0.149674 0 +DEAL::162 0.064869 -0.009370 0 +DEAL::578 0.014893 -0.000911 0 +DEAL::2178 0.003620 -0.000108 0 +DEAL::dim: 2 FE_ABF<2>(1) +DEAL::DoFs ||u-u_h|| curl(u_h) div(u_h) +DEAL::208 0.191291 -0.048789 0 +DEAL::800 0.047830 -0.008916 0 +DEAL::3136 0.012716 -0.002205 0 +DEAL::12416 0.003947 -0.000586 0 +DEAL:: -- 2.39.5