From 7cdbee674b3b59e6089f29006d6c142ee84c113b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 17 Nov 2017 09:43:17 -0700 Subject: [PATCH] Add test. --- tests/numerics/error_estimator_02_complex.cc | 711 ++++++++++++++++++ .../error_estimator_02_complex.output | 62 ++ 2 files changed, 773 insertions(+) create mode 100644 tests/numerics/error_estimator_02_complex.cc create mode 100644 tests/numerics/error_estimator_02_complex.output diff --git a/tests/numerics/error_estimator_02_complex.cc b/tests/numerics/error_estimator_02_complex.cc new file mode 100644 index 0000000000..0539ed4ded --- /dev/null +++ b/tests/numerics/error_estimator_02_complex.cc @@ -0,0 +1,711 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2016 by the deal.II authors +// +// 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. +// +// --------------------------------------------------------------------- + + +// like the error_estimator_02 test but using a complex-valued vector + + +#include "../tests.h" + +// dealii +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// c++ +#include + +using namespace dealii; + +template +class MyFunction : public dealii::Function > +{ +public: + MyFunction(const double k); + + virtual std::complex value(const dealii::Point &point, + const unsigned int component = 0 ) const; + + double get_k() const; + +private: + const double k; +}; + +template +MyFunction::MyFunction(const double k) + : + dealii::Function >(1), + k(k) +{ + +} + +template +std::complex +MyFunction::value(const dealii::Point &point, + const unsigned int ) const +{ + const double x = point[0]-1.0; + + if (x < 0) + return 0.0; + else + return std::complex(0, k * x); +} + +template +double MyFunction::get_k() const +{ + return k; +} + +// neuman bc +template +class NeumanBC : public dealii::Function > +{ +public: + NeumanBC(const double c); + + virtual std::complex value(const dealii::Point &point, + const unsigned int component = 0 ) const; + + double get_c() const; + +private: + const double c; +}; + +template +NeumanBC::NeumanBC(const double c) + : + dealii::Function >(1), + c(c) +{ +} + +template +std::complex NeumanBC::value(const dealii::Point &point, + const unsigned int ) const +{ + return std::complex(0, c); +} + +template +double NeumanBC::get_c() const +{ + return c; +} + +// helper function to get diagonal and +// area of the squared element with lenght h +template +void get_h_area(double &h, double &a, const double L); + +template <> +void get_h_area<2>(double &h, double &a, const double L) +{ + h = L; + a = L; +} + +template <> +void get_h_area<3>(double &h, double &a, const double L) +{ + h = std::sqrt(2.0)*L; + a = L*L; +} + +// helper function to get diagonal and area of the +// h-refined face. +template +void get_h_area_sub(double &h, double &a, const double L); + +template <> +void get_h_area_sub<2>(double &h, double &a, const double L) +{ + h = L/2; + a = L/2; +} + +template <> +void get_h_area_sub<3>(double &h, double &a, const double L) +{ + h = std::sqrt(2.0)*L/2; + a = L*L/4.0; +} + +// output for inspection +template +void output(const std::string name, + const Triangulation &triangulation, + const hp::DoFHandler &dof_handler, + const Vector > &values, + const Vector &error) +{ + dealii::Vector fe_degrees(triangulation.n_active_cells()); + { + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + fe_degrees(index) = dof_handler.get_fe()[cell->active_fe_index()].degree; + } + + dealii::DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector(values, + std::string("function_interpolation")); + data_out.add_data_vector(fe_degrees, + std::string("fe_degree")); + data_out.add_data_vector(error, + std::string("error")); + data_out.build_patches (); + + std::ofstream output (name.c_str()); + data_out.write_vtu (output); +} + +// case 1) +template +void test_neumann(const NeumanBC &func) +{ + deallog << "NeumanBC case:"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + hp::QCollection face_quadrature_formula; + ConstraintMatrix constraints; + + const unsigned int p = 3; + + fe_collection.push_back(dealii::FE_Q(QIterated<1>(QTrapez<1>(),p))); + quadrature_formula.push_back(dealii::QGauss(p+5)); + face_quadrature_formula.push_back(dealii::QGauss(p+5)); + + const double L = 2.0; + + // set-up domain + { + GridGenerator::hyper_cube (triangulation,.0,L,/*colorize*/true); + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector > values(dof_handler.n_dofs()); + + dealii::deallog << "dof values:"< >::type function_map; + function_map[0] = &func; + + dealii::Vector error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + function_map, + values, + error, + dealii::ComponentMask(), + nullptr, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"<(h,A,L); + const double expected_value_squared = h*A*std::pow(func.get_c(),2)/p; + dealii::deallog << "expected:"<< std::endl <<" "<< std::sqrt(expected_value_squared) << std::endl; + + AssertThrow (std::fabs(std::sqrt(expected_value_squared) - error[0] ) < 1e-5, dealii::ExcInternalError()); + + dof_handler.clear(); +} + +// case 2) +template +void test_regular(const MyFunction &func) +{ + deallog << std::endl; + deallog << "Regular face:"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + hp::QCollection face_quadrature_formula; + ConstraintMatrix constraints; + + const unsigned int p1 = 1; + const unsigned int p2 = 2; + std::vector p_degree; + p_degree.push_back(p1); + p_degree.push_back(p2); + + for (unsigned int i=0; i(QIterated<1>(QTrapez<1>(),p))); + quadrature_formula.push_back(dealii::QGauss(p+5)); + face_quadrature_formula.push_back(dealii::QGauss(p+5)); + } + + const double L = 2.0; + + // set-up domain + { + std::vector repetitions(dim,1); + repetitions[0] = 2; + dealii::Point p1; + dealii::Point p2; + for (unsigned int d = 0; d < dim; d++) + { + p1[d] = 0.0; + p2[d] = L; + } + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + p1, + p2, + /*colorize*/ false); + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; cell++) + if (cell->center()[0] > 1.0) + { + cell->set_active_fe_index(1); + break; + } + + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector > values(dof_handler.n_dofs()); + dealii::VectorTools::interpolate (dof_handler, + func, + values); + + dealii::deallog << "dof values:"< error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + typename dealii::FunctionMap >::type(), + values, + error, + dealii::ComponentMask(), + nullptr, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"<(h,A,L); + const double expected_value_squared = h*A*std::pow(func.get_k(),2)/2.0/std::max(p1,p2); + dealii::deallog << "expected:"<< std::endl <<" "<< std::sqrt(expected_value_squared) << std::endl; + for (unsigned int i = 0; i < error.size(); i++) + AssertThrow (std::fabs(std::sqrt(expected_value_squared) - error[i] ) < 1e-6, dealii::ExcInternalError()); + + dof_handler.clear(); +} + +// case 3) +template +void test_irregular(const MyFunction &func) +{ + deallog << std::endl; + deallog << "Irregular face:"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + hp::QCollection face_quadrature_formula; + ConstraintMatrix constraints; + + const unsigned int p1 = 1; + const unsigned int p2 = 2; + const unsigned int p3 = 3; + std::vector p_degree; + p_degree.push_back(p1); + p_degree.push_back(p2); + p_degree.push_back(p3); + + for (unsigned int i=0; i(QIterated<1>(QTrapez<1>(),p))); + quadrature_formula.push_back(dealii::QGauss(p+5)); + face_quadrature_formula.push_back(dealii::QGauss(p+5)); + } + + const double L = 2.0; + + // set-up domain + { + std::vector repetitions(dim,1); + repetitions[0] = 2; + dealii::Point p1; + dealii::Point p2; + for (unsigned int d = 0; d < dim; d++) + { + p1[d] = 0.0; + p2[d] = L; + } + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + p1, + p2, + /*colorize*/ false); + // refine left side + { + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + } + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; cell++) + if (cell->center()[0] > 1.0) // right + { + cell->set_active_fe_index(0); + } + else if (cell->center()[1] > 1.0) // top + { + cell->set_active_fe_index(2); + } + else + { + cell->set_active_fe_index(1); + } + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector > values(dof_handler.n_dofs()); + dealii::VectorTools::interpolate (dof_handler, + func, + values); + + dealii::deallog << "dof values:"< error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + typename dealii::FunctionMap >::type(), + values, + error, + dealii::ComponentMask(), + nullptr, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"<(h,A,L); + // + const double expected_squared_1 = h*A*std::pow(func.get_k(),2)/2.0/std::max(p3,p1); + const double expected_squared_2 = h*A*std::pow(func.get_k(),2)/2.0/std::max(p2,p1); + const double expected_squared_3 = (dim==2) ? + expected_squared_1 + expected_squared_2: + 2*expected_squared_1 + 2*expected_squared_2; + + std::vector expected_error(error.size(),0.0); + + expected_error[0] = std::sqrt(expected_squared_3); + expected_error[2] = std::sqrt(expected_squared_2); + expected_error[4] = std::sqrt(expected_squared_1); + + if (dim ==3) + { + expected_error[6] = expected_error[2]; + expected_error[8] = expected_error[4]; + } + + dealii::deallog << "expected:"<< std::endl; + for (unsigned int i = 0; i < expected_error.size(); i++) + deallog<<" " << expected_error[i]; + deallog < +class MySecondFunction : public dealii::Function > +{ +public: + MySecondFunction(); + + virtual std::complex value(const dealii::Point &point, + const unsigned int component = 0 ) const; +}; + +template +MySecondFunction::MySecondFunction() + : + dealii::Function >(1) +{ + +} + +template +std::complex MySecondFunction::value(const dealii::Point &point, + const unsigned int ) const +{ + double f = 0.0; + const double &x = point[0]; + Assert (dim>1, dealii::ExcNotImplemented()); + const double &y = point[1]; + + return std::complex(0, (1.-x)*(1.-y)*(1.-y)+std::pow(1.0-y,4)*std::exp(-x)); +} + +template +void test(const MySecondFunction &func) +{ + deallog << std::endl; + deallog << "More complicated mesh:"< triangulation; + dealii::hp::DoFHandler dof_handler(triangulation); + dealii::hp::FECollection fe_collection; + dealii::hp::QCollection quadrature_formula; + dealii::hp::QCollection face_quadrature_formula; + dealii::ConstraintMatrix constraints; + for (unsigned int p = 1; p <=3; p++) + { + fe_collection.push_back(dealii::FE_Q(QIterated<1>(QTrapez<1>(),p))); + quadrature_formula.push_back(dealii::QGauss(p+1)); + face_quadrature_formula.push_back(dealii::QGauss(p+1)); + } + dealii::GridGenerator::hyper_cube (triangulation,0.0,1.0); // reference cell + + // refine + { + triangulation.refine_global(1); + + // otherwise the next set_active_fe_index + // will not carry to the child cells. + dof_handler.distribute_dofs (fe_collection); + + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; cell++) + { + bool in_top_left = true; + for (unsigned int d=0; d< dim; d++) + in_top_left = in_top_left && (cell->center()[d] < 0.5); + + if (in_top_left) + { + cell->set_active_fe_index(1); + cell->set_refine_flag(); + break; + } + } + + triangulation.prepare_coarsening_and_refinement(); + + triangulation.execute_coarsening_and_refinement(); + + cell = dof_handler.begin_active(); + + for (; cell != endc; cell++) + { + if (cell->center()[0] < 0.25) + { + cell->set_active_fe_index(2); + } + } + } + + dof_handler.distribute_dofs (fe_collection); + + // constraints + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + // interpolate some function + dealii::Vector > values(dof_handler.n_dofs()); + + dealii::VectorTools::interpolate (dof_handler, + func, + values); + + dealii::deallog << "dof values:"< error(dof_handler.n_dofs()); + dealii::KellyErrorEstimator::estimate(dof_handler, + face_quadrature_formula, + typename dealii::FunctionMap >::type (), + values, + error, + dealii::ComponentMask(), + nullptr, + dealii::numbers::invalid_unsigned_int, + dealii::numbers::invalid_subdomain_id, + dealii::numbers::invalid_material_id, + dealii::KellyErrorEstimator::face_diameter_over_twice_max_degree); + + dealii::deallog <<"error:"< func(1.25); + test_neumann(func); + } + + { + MyFunction<2> func(0.25); + test_regular(func); + } + + { + MyFunction<2> func(0.75); + test_irregular(func); + } + + deallog << "===3d==="< func(1.25); + test_neumann(func); + } + + { + MyFunction<3> func(0.25); + test_regular(func); + } + + { + MyFunction<3> func(0.75); + test_irregular(func); + } + + { + MySecondFunction<2> function; + test(function); + } + + + + dealii::deallog << "Ok"<