From 77247368473f432a0558ec0762b3c0d06d735c0b Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 24 Oct 2007 20:11:24 +0000 Subject: [PATCH] Boundary projection for RT git-svn-id: https://svn.dealii.org/trunk@15377 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/deal.II/Makefile | 1 + tests/deal.II/project_boundary_rt_01.cc | 192 ++++++++++++++++++ .../project_boundary_rt_01/cmp/generic | 7 + 3 files changed, 200 insertions(+) create mode 100644 tests/deal.II/project_boundary_rt_01.cc create mode 100644 tests/deal.II/project_boundary_rt_01/cmp/generic diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile index 071d849107..ba01ac2875 100644 --- a/tests/deal.II/Makefile +++ b/tests/deal.II/Makefile @@ -62,6 +62,7 @@ tests_x = block_matrices \ have_same_coarse_mesh_* \ get_finest_common_cells_* \ project_*[0-9] \ + project_boundary_* \ minimal_cell_diameter \ maximal_cell_diameter \ union_triangulation \ diff --git a/tests/deal.II/project_boundary_rt_01.cc b/tests/deal.II/project_boundary_rt_01.cc new file mode 100644 index 0000000000..44bda3e0db --- /dev/null +++ b/tests/deal.II/project_boundary_rt_01.cc @@ -0,0 +1,192 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2007 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. +// +//---------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +/** + * A vector-valued polynomial for testing RT elements. + */ +template +class TestFunction : public Function +{ + public: + + /// Construct a polynomial of degree p + TestFunction(unsigned int degree); + + virtual void vector_value_list (const std::vector > &points, + std::vector > &values) const; + private: + unsigned int degree; +}; + + +template +TestFunction::TestFunction (unsigned int p) + : + Function(dim), degree(p) +{} + + +template +void +TestFunction::vector_value_list (const std::vector > &points, + std::vector > &values) const +{ + if (degree < 2) + { + Assert(false, ExcNotImplemented()); + } + else + { + for (unsigned int k=0;k p = points[k]; + for (unsigned int dd=0;dd +double integrate_error(const DoFHandler& dof, + FEFaceValues& fe, + const Vector& u, + const Function& f) +{ + double result = 0.; + std::vector > f_values(fe.n_quadrature_points, Vector(dim)); + std::vector > fe_values(fe.n_quadrature_points, Vector(dim)); + + for (typename DoFHandler::active_cell_iterator cell = dof.begin_active(); + cell != dof.end(); ++cell) + { + for (unsigned int face=0 ; face != GeometryInfo::faces_per_cell; ++face) + { + if (!cell->at_boundary(face)) continue; + + fe.reinit(cell, face); + f.vector_value_list(fe.get_quadrature_points(), f_values); + fe.get_function_values(u, fe_values); + for (unsigned int k=0;k +void test_projection (const Triangulation& tr, + const FiniteElement& fe) +{ + deallog << fe.get_name() << std::endl; + + const unsigned int degree = fe.tensor_degree(); + + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + QGauss quadrature(degree+2); + MappingQ1 mapping; + + TestFunction f(degree-1); + std::map boundary_constraints; + typename FunctionMap::type boundary_map; + for (unsigned char i=0;i<255;++i) + boundary_map[i] = &f; + VectorTools::project_boundary_values(mapping, dof, boundary_map, quadrature, + boundary_constraints); + + // Fill a vector with the projected + // boundary values + Vector u(dof.n_dofs()); + for (typename std::map::const_iterator + i = boundary_constraints.begin(); i != boundary_constraints.end(); ++i) + u(i->first) = i->second; + + FEFaceValues feval(mapping, fe, quadrature, + update_quadrature_points + | update_normal_vectors + | update_JxW_values + | update_values); + double err = integrate_error(dof, feval, u, f); + deallog << err << std::endl; + + DataOut dout; + dout.attach_dof_handler(dof); + std::ofstream of("T.gnuplot"); + dout.add_data_vector(u, "u"); + dout.build_patches(3); + dout.write_gnuplot(of); +} + + +template +void test_hyper_cube(const FiniteElement& fe) +{ + Triangulation tr; + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + test_projection(tr, fe); +} + + +int main() +{ + std::ofstream logfile ("project_boundary_rt_01/output"); + logfile.precision (2); + logfile.setf(std::ios::fixed); + deallog.attach(logfile); + deallog.depth_console (0); + deallog.threshold_double(1.e-12); + + FE_RaviartThomasNodal<2> rt22(2); + test_hyper_cube(rt22); + FE_RaviartThomasNodal<2> rt23(3); + test_hyper_cube(rt23); + FE_RaviartThomasNodal<2> rt24(4); + test_hyper_cube(rt24); + + +} diff --git a/tests/deal.II/project_boundary_rt_01/cmp/generic b/tests/deal.II/project_boundary_rt_01/cmp/generic new file mode 100644 index 0000000000..5bef325afb --- /dev/null +++ b/tests/deal.II/project_boundary_rt_01/cmp/generic @@ -0,0 +1,7 @@ + +DEAL::FE_RaviartThomasNodal<2>(2) +DEAL::0 +DEAL::FE_RaviartThomasNodal<2>(3) +DEAL::0 +DEAL::FE_RaviartThomasNodal<2>(4) +DEAL::0 -- 2.39.5