diff --git a/source/ClosestSourcePointAdaptation.hh b/source/ClosestSourcePointAdaptation.hh new file mode 100644 index 00000000..fad45050 --- /dev/null +++ b/source/ClosestSourcePointAdaptation.hh @@ -0,0 +1,117 @@ +/* Copyright (c) 2016 - 2024, the adamantine authors. + * + * This file is subject to the Modified BSD License and may not be distributed + * without copyright and license information. Please refer to the file LICENSE + * for the text and further information on this license. + */ + +#ifndef CLOSEST_SOURCE_POINT_ADAPTATION_HH +#define CLOSEST_SOURCE_POINT_ADAPTATION_HH + +#include +#include + +#include + +namespace adamantine +{ + +template +class ClosestQuadPointAdaptation +{ +public: + template + ClosestQuadPointAdaptation(const Quadrature &quadrature) + : m_fe_values_coarse(m_fe_nothing, quadrature, + dealii::update_quadrature_points), + m_fe_values_fine(m_fe_nothing, quadrature, + dealii::update_quadrature_points) + { + } + + std::vector> coarse_to_fine( + const typename dealii::Triangulation::cell_iterator + &parent, + const std::vector parent_values) + { + m_fe_values_coarse.reinit(parent); + const std::vector> &coarse_quad_points = + m_fe_values_coarse.get_quadrature_points(); + + std::vector> child_values( + parent->n_children(), + std::vector(coarse_quad_points.size())); + for (unsigned int i = 0; i < parent->n_children(); ++i) + { + m_fe_values_fine.reinit(parent->child(i)); + const std::vector> &fine_quad_points = + m_fe_values_fine.get_quadrature_points(); + for (unsigned int j = 0; j < fine_quad_points.size(); ++j) + { + std::pair min_index_distance{-1, std::numeric_limits::max()}; + for (unsigned int k = 0; k < coarse_quad_points.size(); ++k) + { + double candidate_distance = + (fine_quad_points[j] - coarse_quad_points[k]).norm_square(); + if (candidate_distance < min_index_distance.second) + { + min_index_distance.first = k; + min_index_distance.second = candidate_distance; + } + } + child_values[i][j] = parent_values[min_index_distance.first]; + } + } + return child_values; + } + + std::vector fine_to_coarse( + const typename dealii::Triangulation::cell_iterator + &parent, + const std::vector> &children_values) + { + m_fe_values_coarse.reinit(parent); + const std::vector> &coarse_quad_points = + m_fe_values_coarse.get_quadrature_points(); + + std::vector, double>> min_index_distances{ + coarse_quad_points.size(), + std::pair{std::pair{-1, -1}, std::numeric_limits::max()}}; + for (unsigned int i = 0; i < parent->n_children(); ++i) + { + m_fe_values_fine.reinit(parent->child(i)); + const std::vector> &fine_quad_points = + m_fe_values_fine.get_quadrature_points(); + for (unsigned int j = 0; j < fine_quad_points.size(); ++j) + { + for (unsigned int k = 0; k < coarse_quad_points.size(); ++k) + { + double candidate_distance = + (fine_quad_points[j] - coarse_quad_points[k]).norm_square(); + if (candidate_distance < min_index_distances[k].second) + { + min_index_distances[k].first = std::pair{i, j}; + min_index_distances[k].second = candidate_distance; + } + } + } + } + std::vector parent_values(coarse_quad_points.size()); + for (unsigned int k = 0; k < coarse_quad_points.size(); ++k) + { + parent_values[k] = children_values[min_index_distances[k].first.first] + [min_index_distances[k].first.second]; + } + + return parent_values; + } + +private: + dealii::FE_Nothing m_fe_nothing; + dealii::FEValues m_fe_values_coarse; + dealii::FEValues m_fe_values_fine; +}; + +} // namespace adamantine + +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8486c32e..c493ec7c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -46,6 +46,8 @@ list(APPEND # test_integration_da is a special test that we need to test with more # processors than the others adamantine_ADD_BOOST_TEST(test_integration_da 1 2 3 6) +# same goes for test_closest_source_point_adaptation +adamantine_ADD_BOOST_TEST(test_closest_source_point_adaptation 1 2 3 6) foreach(TEST_NAME ${UNIT_TESTS}) adamantine_ADD_BOOST_TEST(${TEST_NAME}) diff --git a/tests/test_closest_source_point_adaptation.cc b/tests/test_closest_source_point_adaptation.cc new file mode 100644 index 00000000..3e16183e --- /dev/null +++ b/tests/test_closest_source_point_adaptation.cc @@ -0,0 +1,129 @@ +/* Copyright (c) 2016 - 2024, the adamantine authors. + * + * This file is subject to the Modified BSD License and may not be distributed + * without copyright and license information. Please refer to the file LICENSE + * for the text and further information on this license. + */ + +#define BOOST_TEST_MODULE ClosestSourcePointAdaptation + +#include "../source/ClosestSourcePointAdaptation.hh" + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "main.cc" + +template +void test() +{ + dealii::parallel::distributed::Triangulation tria( + MPI_COMM_WORLD); + + dealii::GridGenerator::subdivided_hyper_cube(tria, 2); + tria.refine_global(3); + + dealii::QGauss quadrature(2); + const unsigned int n_q_points = quadrature.size(); + const unsigned int n_active_cells = tria.n_global_active_cells(); + + std::vector> cell_quad_point_values( + n_active_cells, + std::vector(n_q_points, + std::numeric_limits::signaling_NaN())); + + std::random_device rnd_device; + std::mt19937 mersenne_engine{rnd_device()}; + std::uniform_real_distribution dist{1., 2.}; + + for (auto &cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->set_refine_flag(); + std::generate(cell_quad_point_values[cell->active_cell_index()].begin(), + cell_quad_point_values[cell->active_cell_index()].end(), + [&]() { return dist(mersenne_engine); }); + } + + adamantine::ClosestQuadPointAdaptation + closest_quad_point_adaptation(quadrature); + dealii::parallel::distributed::CellDataTransfer< + dim, spacedim, std::vector>> + cell_data_transfer( + tria, false, + [&](const typename dealii::Triangulation::cell_iterator + &parent, + const std::vector parent_values) + { + return closest_quad_point_adaptation.coarse_to_fine(parent, + parent_values); + }, + [&](const typename dealii::Triangulation::cell_iterator + &parent, + const std::vector> child_values) + { + return closest_quad_point_adaptation.fine_to_coarse(parent, + child_values); + }); + + cell_data_transfer.prepare_for_coarsening_and_refinement( + cell_quad_point_values); + tria.execute_coarsening_and_refinement(); + + const unsigned int n_active_cells_new = tria.n_global_active_cells(); + + std::vector> cell_quad_point_values_new( + n_active_cells_new, + std::vector(n_q_points, + std::numeric_limits::signaling_NaN())); + cell_data_transfer.unpack(cell_quad_point_values_new); + + for (auto &cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->set_coarsen_flag(); + for (unsigned int i = 0; i < n_q_points; ++i) + { + const double quad_point_value = + cell_quad_point_values_new[cell->active_cell_index()][i]; + BOOST_TEST(quad_point_value >= 1.); + BOOST_TEST(quad_point_value <= 2.); + } + } + + cell_data_transfer.prepare_for_coarsening_and_refinement( + cell_quad_point_values_new); + tria.execute_coarsening_and_refinement(); + + BOOST_TEST(tria.n_global_active_cells() == n_active_cells); + std::vector> cell_quad_point_values_final( + n_active_cells, + std::vector(n_q_points, + std::numeric_limits::signaling_NaN())); + + cell_data_transfer.unpack(cell_quad_point_values_final); + + for (auto &cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->set_refine_flag(); + for (unsigned int i = 0; i < n_q_points; ++i) + { + BOOST_TEST(cell_quad_point_values_final[cell->active_cell_index()][i] == + cell_quad_point_values[cell->active_cell_index()][i]); + } + } +} + +BOOST_AUTO_TEST_CASE(closest_source_point_adaptation) +{ + test<2, 2>(); + test<3, 3>(); +}