From bb48d8dd5d65a5a1aef631c53e4c12a70ff70515 Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 25 Jul 2017 19:20:44 -0400 Subject: [PATCH] Describe the changes in a test and changelog entry. --- doc/news/changes/minor/20170725DiIlioEtAl | 4 + tests/mappings/mapping_q1_cartesian_grid.cc | 173 ++++++++++++++++++ .../mappings/mapping_q1_cartesian_grid.output | 9 + 3 files changed, 186 insertions(+) create mode 100644 doc/news/changes/minor/20170725DiIlioEtAl create mode 100644 tests/mappings/mapping_q1_cartesian_grid.cc create mode 100644 tests/mappings/mapping_q1_cartesian_grid.output diff --git a/doc/news/changes/minor/20170725DiIlioEtAl b/doc/news/changes/minor/20170725DiIlioEtAl new file mode 100644 index 0000000000..6a18dd8fac --- /dev/null +++ b/doc/news/changes/minor/20170725DiIlioEtAl @@ -0,0 +1,4 @@ +Fixed: MappingQGeneric::transform_real_to_unit_cell had a bug where, for certain +parallelograms, a problem with subtracting nearly equal floating point numbers +ruined most of the digits of accuracy in the coordinate transformation. +
(Giovanni Di Ilio, David Wells, Jean-Paul Pelteret, 2017/07/25) diff --git a/tests/mappings/mapping_q1_cartesian_grid.cc b/tests/mappings/mapping_q1_cartesian_grid.cc new file mode 100644 index 0000000000..f5ae71cc3e --- /dev/null +++ b/tests/mappings/mapping_q1_cartesian_grid.cc @@ -0,0 +1,173 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +// Test that MappingQ1 works on Cartesian grids. This is a regression test for +// a bug where, due to the quadratic equation solver in +// transform_unit_to_real_cell, we lost nearly all digits of precision due to +// roundoff. + +using namespace dealii; + +class RefineTest +{ + +public: + static constexpr int dim = 2 ; + +public: + + RefineTest(const unsigned int refine_case) + : + triangulation(), + dof_handler(triangulation), + refine_case(refine_case) + { + make_grid(); + test_mapping(); + } + +private: + void make_grid() + { + if (refine_case == 1) + { + std::vector< unsigned int > repetitions(dim); + repetitions[0] = repetitions[1] = 5; // Non-unity number of initial partitions + const std::array L = {{1.0, 1.0}}; + GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true); + triangulation.refine_global(1); // With refinement + } + else if (refine_case == 2) + { + std::vector< unsigned int > repetitions(dim); + repetitions[0] = repetitions[1] = 5; // Non-unity number of initial partitions + const std::array L = {{1.0, 1.0}}; + GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true); + triangulation.refine_global(0); // With no refinement + } + else if (refine_case == 3) + { + std::vector< unsigned int > repetitions(dim); + repetitions[0] = repetitions[1] = 10; // Finer initial mesh + const std::array L = {{1.0, 1.0}}; + GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true); + triangulation.refine_global(0); // With no refinement + } + else if (refine_case == 4) + { + std::vector< unsigned int > repetitions(dim); + repetitions[0] = repetitions[1] = 1; // Single cell + const std::array L = {{1.0, 1.0}}; + GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, Point(0.0,0.0), Point(L[0],L[1]), true); + triangulation.refine_global(4); // With pure global refinement + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + + void test_mapping() + { + const double tol = 1e-8; + const MappingQGeneric mapping(1); + + deallog << "Number of active cells: " << triangulation.n_active_cells() << std::endl; + + auto cell = dof_handler.begin_active(); + const auto endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + { + const auto P1 = cell->vertex(0); + const auto P2 = cell->vertex(1); + const auto P3 = cell->vertex(2); + + // Find the mid-point of the cell + Point test; + test[0] = (P1[0]+P2[0])/2.0; + test[1] = (P2[1]+P3[1])/2.0; + + // Map to and from the unit cell + const Point dp = mapping.transform_real_to_unit_cell(cell,test); + const Point dp_real = mapping.transform_unit_to_real_cell(cell,dp); + + // Now we start off at the mid-point of the unit cell + Point test_unit_mid; + test_unit_mid[0] = 0.5; + test_unit_mid[1] = 0.5; + + // And map to and from the real cell + const Point dp_real_mid = mapping.transform_unit_to_real_cell(cell,test_unit_mid); + const Point dp_test_mid = mapping.transform_real_to_unit_cell(cell,dp_real_mid); + + if ((test-dp_real).norm() > tol) + { + deallog << " " << std::endl; + deallog << "ERROR" << std::endl; + deallog << "cell = " << index << std::endl; + deallog << "cell vertex(0): " << cell->vertex(0) << std::endl; + deallog << "cell vertex(1): " << cell->vertex(1) << std::endl; + deallog << "cell vertex(2): " << cell->vertex(2) << std::endl; + deallog << "cell vertex(3): " << cell->vertex(3) << std::endl; + deallog << " " << std::endl; + deallog << "test point = " << test << " mapped point = " << dp << " back-mapping = " << dp_real << std::endl; + deallog << "test point (unit mid) = " << test_unit_mid << " mapped point = " << dp_real_mid << " forward-mapping = " << dp_test_mid << std::endl; + deallog << "=======================================================================================" << std::endl; + } + } + } + +private: + Triangulation triangulation; + DoFHandler dof_handler; + const unsigned int refine_case; +}; + +int main(/*int argc, char *argv[]*/) +{ + initlog(); + + { + deallog << "Refine case = 1" << std::endl; + RefineTest refineglobal(1); + } + { + deallog << "Refine case = 2" << std::endl; + RefineTest refineglobal(2); + } + { + deallog << "Refine case = 3" << std::endl; + RefineTest refineglobal(3); + } + { + deallog << "Refine case = 4" << std::endl; + RefineTest refineglobal(4); + } + return 0; +} diff --git a/tests/mappings/mapping_q1_cartesian_grid.output b/tests/mappings/mapping_q1_cartesian_grid.output new file mode 100644 index 0000000000..6265d14a74 --- /dev/null +++ b/tests/mappings/mapping_q1_cartesian_grid.output @@ -0,0 +1,9 @@ + +DEAL::Refine case = 1 +DEAL::Number of active cells: 100 +DEAL::Refine case = 2 +DEAL::Number of active cells: 25 +DEAL::Refine case = 3 +DEAL::Number of active cells: 100 +DEAL::Refine case = 4 +DEAL::Number of active cells: 256 -- 2.39.5