From 253f73b3d7efe8feec678d7f8b0935864d0abcce Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 28 Feb 2020 11:15:02 -0700 Subject: [PATCH] Add tests. --- tests/numerics/nedelec_projection_mixed_01.cc | 67 +++ .../nedelec_projection_mixed_01.output | 5 + tests/numerics/project_bv_curl_conf_04.cc | 419 ++++++++++++++++++ tests/numerics/project_bv_curl_conf_04.output | 15 + 4 files changed, 506 insertions(+) create mode 100644 tests/numerics/nedelec_projection_mixed_01.cc create mode 100644 tests/numerics/nedelec_projection_mixed_01.output create mode 100644 tests/numerics/project_bv_curl_conf_04.cc create mode 100644 tests/numerics/project_bv_curl_conf_04.output diff --git a/tests/numerics/nedelec_projection_mixed_01.cc b/tests/numerics/nedelec_projection_mixed_01.cc new file mode 100644 index 0000000000..cf765609f7 --- /dev/null +++ b/tests/numerics/nedelec_projection_mixed_01.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// check VectorTools::project_boundary_values_curl_conforming_l2 for a +// combination element that has edge DoFs other than those just from +// the Nedelec base element + +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +int +main() +{ + initlog(); + + const unsigned int dim = 3; + + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, 0, 1, true); + MappingQ mapping(1); + AffineConstraints constraints; + + FESystem fe_system(FE_Q(2), 1, FE_Nedelec(0), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + + constraints.clear(); + VectorTools::project_boundary_values_curl_conforming_l2( + dof_handler, + 1, + Functions::ZeroFunction(dim + 1), + 0, + constraints, + mapping); + constraints.close(); + + constraints.print(deallog.get_file_stream()); +} diff --git a/tests/numerics/nedelec_projection_mixed_01.output b/tests/numerics/nedelec_projection_mixed_01.output new file mode 100644 index 0000000000..8d00f4cbc9 --- /dev/null +++ b/tests/numerics/nedelec_projection_mixed_01.output @@ -0,0 +1,5 @@ + + 9 = 0 + 17 = 0 + 25 = 0 + 29 = 0 diff --git a/tests/numerics/project_bv_curl_conf_04.cc b/tests/numerics/project_bv_curl_conf_04.cc new file mode 100644 index 0000000000..8f93314a55 --- /dev/null +++ b/tests/numerics/project_bv_curl_conf_04.cc @@ -0,0 +1,419 @@ +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include +using namespace dealii; + +const unsigned int dim = 3; + +void +apply_boundary_values(DoFHandler & dof_handler, + AffineConstraints &constraints, + unsigned int n_comps, + unsigned int start_comp, + Mapping & mapping, + Vector & dst) +{ + constraints.clear(); + for (unsigned int i = 0; i < 6; i++) + VectorTools::project_boundary_values_curl_conforming_l2( + dof_handler, + start_comp, + Functions::ConstantFunction(1, n_comps), + i, + constraints, + mapping); + constraints.close(); + + constraints.distribute(dst); +} + +bool +test_boundary_values(DoFHandler & dof_handler, + Mapping & mapping, + FiniteElement &fe, + unsigned int n_comps, + unsigned int start_comp, + Vector & vec) +{ + double ret = true; + // Initialize + QGaussLobatto quadrature(3); + FEFaceValues fe_values(mapping, + fe, + quadrature, + update_values | update_normal_vectors | + update_quadrature_points | update_JxW_values); + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + unsigned num_cells = 0; + std::vector> local_averages( + 6, Vector(dof_handler.n_dofs())); + const FEValuesExtractors::Vector nedelec(start_comp); + + // For testing we take only one cell! + assert(dof_handler.n_dofs() == cell->get_fe().dofs_per_cell); + + for (; cell != endc; ++cell) + { + // For testing we take only one cell! Otherwise, local to global is + // missing + assert(num_cells == 0); + for (unsigned int face = 0; + face < dealii::GeometryInfo::faces_per_cell; + ++face) + { + if (cell->face(face)->at_boundary()) + { + unsigned int boundary_number = cell->face(face)->boundary_id(); + Tensor<1, dim> expected_value; + expected_value[0] = 1; + expected_value[1] = 1; + expected_value[2] = 1; + assert(boundary_number < 6); + fe_values.reinit(cell, face); + + std::vector> local_values( + fe_values.n_quadrature_points, Vector(n_comps)); + fe_values.get_function_values(vec, local_values); + + for (unsigned int q_point = 0; + q_point < fe_values.n_quadrature_points; + q_point++) + { + // Compare n x v values + Tensor<1, dim> nedelec_values; + nedelec_values[0] = local_values[q_point](start_comp); + nedelec_values[1] = local_values[q_point](start_comp + 1); + nedelec_values[2] = local_values[q_point](start_comp + 2); + Tensor<1, dim> normal = fe_values.normal_vector(q_point); + + for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; + i++) + { + local_averages[boundary_number](i) += + ((cross_product_3d(normal, nedelec_values) - + cross_product_3d(normal, expected_value)) * + cross_product_3d( + normal, fe_values[nedelec].value(i, q_point))) * + fe_values.JxW(q_point); + } + } + } + } + num_cells++; + } + // Test if ok + for (unsigned int bc = 0; bc < 6; bc++) + { + for (unsigned int basis = 0; basis < dof_handler.n_dofs(); basis++) + { + if (fabs(local_averages[bc](basis)) > 1.e-10) + { + std::cout << "Wrong values on boundary " << bc + << " in basis function " << basis << " of size " + << fabs(local_averages[bc](basis)) << std::endl; + ret = false; + } + } + } + return ret; +} + +int +main(int /*argc*/, char ** /*argv*/) +{ + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, 0, 1, true); + MappingQ mapping(1); + AffineConstraints constraints; + Vector test_vec; + Vector test_vec2; + + { + // Test 1: Only FE_Nedelec + std::cout << "Testing FE_Nedelec(0): "; + FE_Nedelec fe(0); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + test_vec2.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec2); + if (test_boundary_values(dof_handler, mapping, fe, 3, 0, test_vec2)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 2: FESystem(FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Nedelec(0)): "; + FESystem fe_system(FE_Nedelec(0), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 3, 0, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + // Now, in the case of only one element both initializations should give + // identical vectors + std::cout << "Checking Consistency: "; + bool equal = (test_vec.size() == test_vec2.size()); + if (equal) + { + for (unsigned int i = 0; i < test_vec.size(); i++) + { + equal = equal && (test_vec(i) == test_vec2(i)); + } + } + if (equal) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 1b: Only FE_Nedelec + std::cout << "Testing FE_Nedelec(1): "; + FE_Nedelec fe(1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + test_vec2.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec2); + if (test_boundary_values(dof_handler, mapping, fe, 3, 0, test_vec2)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 2b: FESystem(FE_Nedelec(1)) + std::cout << "Testing FESystem(FE_Nedelec(1)): "; + FESystem fe_system(FE_Nedelec(1), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 3, 0, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 3, 0, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + // Now, in the case of only one element both initializations should give + // identical vectors + std::cout << "Checking Consistency: "; + bool equal = (test_vec.size() == test_vec2.size()); + if (equal) + { + for (unsigned int i = 0; i < test_vec.size(); i++) + { + equal = equal && (test_vec(i) == test_vec2(i)); + } + } + if (equal) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 3: FESystem(FE_Nedelec(0), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Nedelec(0)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Nedelec(0), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 6, 3, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 6, 3, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 4: FESystem(FE_Nedelec(1), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Nedelec(1), FE_Nedelec(0)): "; + FESystem fe_system(FE_Nedelec(1), 1, FE_Nedelec(0), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 6, 3, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 6, 3, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 5: FESystem(FE_Nedelec(0), FE_Nedelec(1)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Nedelec(1)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Nedelec(1), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 6, 3, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 6, 3, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 6: FESystem(FE_Nedelec(0), FE_Q(1)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Q(1)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Q(1), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 4, 0, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 4, 0, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 7: FESystem(FE_Q(1), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Q(1), FE_Nedelec(0)): "; + FESystem fe_system(FE_Q(1), 1, FE_Nedelec(0), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 4, 1, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 4, 1, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 8: FESystem(FE_Nedelec(0), FE_Q(2)) + std::cout << "Testing FESystem(FE_Nedelec(0), FE_Q(2)): "; + FESystem fe_system(FE_Nedelec(0), 1, FE_Q(2), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 4, 0, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 4, 0, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 9: FESystem(FE_Q(2), FE_Nedelec(0)) + std::cout << "Testing FESystem(FE_Q(2), FE_Nedelec(0)): "; + FESystem fe_system(FE_Q(2), 1, FE_Nedelec(0), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 4, 1, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 4, 1, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 10: FESystem(FE_Nedelec(1), FE_Q(1)) + std::cout << "Testing FESystem(FE_Nedelec(1), FE_Q(1)): "; + FESystem fe_system(FE_Nedelec(1), 1, FE_Q(1), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 4, 0, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 4, 0, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + { + // Test 11: FESystem(FE_Q(1), FE_Nedelec(1)) + std::cout << "Testing FESystem(FE_Q(1), FE_Nedelec(1)): "; + FESystem fe_system(FE_Q(1), 1, FE_Nedelec(1), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + test_vec.reinit(dof_handler.n_dofs()); + apply_boundary_values(dof_handler, constraints, 4, 1, mapping, test_vec); + if (test_boundary_values(dof_handler, mapping, fe_system, 4, 1, test_vec)) + { + std::cout << "OK" << std::endl; + } + else + { + std::cout << "Failed" << std::endl; + abort(); + } + } + + return 0; +} diff --git a/tests/numerics/project_bv_curl_conf_04.output b/tests/numerics/project_bv_curl_conf_04.output new file mode 100644 index 0000000000..b8d4974d03 --- /dev/null +++ b/tests/numerics/project_bv_curl_conf_04.output @@ -0,0 +1,15 @@ +Testing FE_Nedelec(0): OK +Testing FESystem(FE_Nedelec(0)): OK +Checking Consistency: OK +Testing FE_Nedelec(1): OK +Testing FESystem(FE_Nedelec(1)): OK +Checking Consistency: OK +Testing FESystem(FE_Nedelec(0), FE_Nedelec(0)): OK +Testing FESystem(FE_Nedelec(1), FE_Nedelec(0)): OK +Testing FESystem(FE_Nedelec(0), FE_Nedelec(1)): OK +Testing FESystem(FE_Nedelec(0), FE_Q(1)): OK +Testing FESystem(FE_Q(1), FE_Nedelec(0)): OK +Testing FESystem(FE_Nedelec(0), FE_Q(2)): OK +Testing FESystem(FE_Q(2), FE_Nedelec(0)): OK +Testing FESystem(FE_Nedelec(1), FE_Q(1)): OK +Testing FESystem(FE_Q(1), FE_Nedelec(1)): OK -- 2.39.5