Skip to content

Commit fdc048b

Browse files
masterleinadmasterleinad
committed
Add ClosestSourcePointAdaptation
1 parent 842ee50 commit fdc048b

File tree

3 files changed

+241
-0
lines changed

3 files changed

+241
-0
lines changed

‎source/ClosestSourcePointAdaptation.hh

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
/* Copyright (c) 2016 - 2024, the adamantine authors.
2+
*
3+
* This file is subject to the Modified BSD License and may not be distributed
4+
* without copyright and license information. Please refer to the file LICENSE
5+
* for the text and further information on this license.
6+
*/
7+
8+
#include <deal.II/fe/fe_nothing.h>
9+
#include <deal.II/fe/fe_values.h>
10+
11+
#include <random>
12+
13+
namespace adamantine
14+
{
15+
16+
template <int dim, int spacedim, typename quad_point_value_type>
17+
class ClosestQuadPointAdaptation
18+
{
19+
public:
20+
template <class Quadrature>
21+
ClosestQuadPointAdaptation(const Quadrature &quadrature)
22+
: m_fe_values_coarse(m_fe_nothing, quadrature,
23+
dealii::update_quadrature_points),
24+
m_fe_values_fine(m_fe_nothing, quadrature,
25+
dealii::update_quadrature_points)
26+
{
27+
}
28+
29+
std::vector<std::vector<quad_point_value_type>> coarse_to_fine(
30+
const typename dealii::Triangulation<dim, spacedim>::cell_iterator
31+
&parent,
32+
const std::vector<quad_point_value_type> parent_values)
33+
{
34+
m_fe_values_coarse.reinit(parent);
35+
const std::vector<dealii::Point<spacedim>> &coarse_quad_points =
36+
m_fe_values_coarse.get_quadrature_points();
37+
38+
std::vector<std::vector<quad_point_value_type>> child_values(
39+
parent->n_children(),
40+
std::vector<quad_point_value_type>(coarse_quad_points.size()));
41+
for (unsigned int i = 0; i < parent->n_children(); ++i)
42+
{
43+
m_fe_values_fine.reinit(parent->child(i));
44+
const std::vector<dealii::Point<spacedim>> &fine_quad_points =
45+
m_fe_values_fine.get_quadrature_points();
46+
for (unsigned int j = 0; j < fine_quad_points.size(); ++j)
47+
{
48+
std::pair min_index_distance{-1, std::numeric_limits<double>::max()};
49+
for (unsigned int k = 0; k < coarse_quad_points.size(); ++k)
50+
{
51+
double candidate_distance =
52+
(fine_quad_points[j] - coarse_quad_points[k]).norm_square();
53+
if (candidate_distance < min_index_distance.second)
54+
{
55+
min_index_distance.first = k;
56+
min_index_distance.second = candidate_distance;
57+
}
58+
}
59+
child_values[i][j] = parent_values[min_index_distance.first];
60+
}
61+
}
62+
return child_values;
63+
}
64+
65+
std::vector<quad_point_value_type> fine_to_coarse(
66+
const typename dealii::Triangulation<dim, spacedim>::cell_iterator
67+
&parent,
68+
const std::vector<std::vector<quad_point_value_type>> &children_values)
69+
{
70+
m_fe_values_coarse.reinit(parent);
71+
const std::vector<dealii::Point<spacedim>> &coarse_quad_points =
72+
m_fe_values_coarse.get_quadrature_points();
73+
74+
std::vector<std::pair<std::pair<int, int>, double>> min_index_distances{
75+
coarse_quad_points.size(),
76+
std::pair{std::pair{-1, -1}, std::numeric_limits<double>::max()}};
77+
for (unsigned int i = 0; i < parent->n_children(); ++i)
78+
{
79+
m_fe_values_fine.reinit(parent->child(i));
80+
const std::vector<dealii::Point<spacedim>> &fine_quad_points =
81+
m_fe_values_fine.get_quadrature_points();
82+
for (unsigned int j = 0; j < fine_quad_points.size(); ++j)
83+
{
84+
for (unsigned int k = 0; k < coarse_quad_points.size(); ++k)
85+
{
86+
double candidate_distance =
87+
(fine_quad_points[j] - coarse_quad_points[k]).norm_square();
88+
if (candidate_distance < min_index_distances[k].second)
89+
{
90+
min_index_distances[k].first = std::pair{i, j};
91+
min_index_distances[k].second = candidate_distance;
92+
}
93+
}
94+
}
95+
}
96+
std::vector<quad_point_value_type> parent_values(coarse_quad_points.size());
97+
for (unsigned int k = 0; k < coarse_quad_points.size(); ++k)
98+
{
99+
parent_values[k] = children_values[min_index_distances[k].first.first]
100+
[min_index_distances[k].first.second];
101+
}
102+
103+
return parent_values;
104+
}
105+
106+
private:
107+
dealii::FE_Nothing<dim, spacedim> m_fe_nothing;
108+
dealii::FEValues<dim, spacedim> m_fe_values_coarse;
109+
dealii::FEValues<dim, spacedim> m_fe_values_fine;
110+
};
111+
112+
} // namespace adamantine

‎tests/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,8 @@ list(APPEND
4646
# test_integration_da is a special test that we need to test with more
4747
# processors than the others
4848
adamantine_ADD_BOOST_TEST(test_integration_da 1 2 3 6)
49+
# same goes for test_closest_source_point_adaptation
50+
adamantine_ADD_BOOST_TEST(test_closest_source_point_adaptation 1 2 3 6)
4951

5052
foreach(TEST_NAME ${UNIT_TESTS})
5153
adamantine_ADD_BOOST_TEST(${TEST_NAME})

‎tests/test_closest_source_point_adaptation.cc

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
/* Copyright (c) 2016 - 2024, the adamantine authors.
2+
*
3+
* This file is subject to the Modified BSD License and may not be distributed
4+
* without copyright and license information. Please refer to the file LICENSE
5+
* for the text and further information on this license.
6+
*/
7+
8+
#define BOOST_TEST_MODULE ClosestSourcePointAdaptation
9+
10+
#include "../source/ClosestSourcePointAdaptation.hh"
11+
12+
#include <deal.II/base/quadrature_lib.h>
13+
#include <deal.II/distributed/cell_data_transfer.templates.h>
14+
#include <deal.II/distributed/tria.h>
15+
#include <deal.II/fe/fe_nothing.h>
16+
#include <deal.II/fe/fe_values.h>
17+
#include <deal.II/grid/grid_generator.h>
18+
19+
#include <random>
20+
#include <string>
21+
22+
#include "main.cc"
23+
24+
template <int dim, int spacedim>
25+
void test()
26+
{
27+
dealii::parallel::distributed::Triangulation<dim, spacedim> tria(
28+
MPI_COMM_WORLD);
29+
30+
dealii::GridGenerator::subdivided_hyper_cube(tria, 2);
31+
tria.refine_global(3);
32+
33+
dealii::QGauss<dim> quadrature(2);
34+
const unsigned int n_q_points = quadrature.size();
35+
const unsigned int n_active_cells = tria.n_global_active_cells();
36+
37+
std::vector<std::vector<double>> cell_quad_point_values(
38+
n_active_cells,
39+
std::vector<double>(n_q_points,
40+
std::numeric_limits<double>::signaling_NaN()));
41+
42+
std::random_device rnd_device;
43+
std::mt19937 mersenne_engine{rnd_device()};
44+
std::uniform_real_distribution<double> dist{1., 2.};
45+
46+
for (auto &cell : tria.active_cell_iterators())
47+
if (cell->is_locally_owned())
48+
{
49+
cell->set_refine_flag();
50+
std::generate(cell_quad_point_values[cell->active_cell_index()].begin(),
51+
cell_quad_point_values[cell->active_cell_index()].end(),
52+
[&]() { return dist(mersenne_engine); });
53+
}
54+
55+
adamantine::ClosestQuadPointAdaptation<dim, spacedim, double>
56+
closest_quad_point_adaptation(quadrature);
57+
dealii::parallel::distributed::CellDataTransfer<
58+
dim, spacedim, std::vector<std::vector<double>>>
59+
cell_data_transfer(
60+
tria, false,
61+
[&](const typename dealii::Triangulation<dim, spacedim>::cell_iterator
62+
&parent,
63+
const std::vector<double> parent_values) {
64+
return closest_quad_point_adaptation.coarse_to_fine(parent,
65+
parent_values);
66+
},
67+
[&](const typename dealii::Triangulation<dim, spacedim>::cell_iterator
68+
&parent,
69+
const std::vector<std::vector<double>> child_values) {
70+
return closest_quad_point_adaptation.fine_to_coarse(parent,
71+
child_values);
72+
});
73+
74+
cell_data_transfer.prepare_for_coarsening_and_refinement(
75+
cell_quad_point_values);
76+
tria.execute_coarsening_and_refinement();
77+
78+
const unsigned int n_active_cells_new = tria.n_global_active_cells();
79+
80+
std::vector<std::vector<double>> cell_quad_point_values_new(
81+
n_active_cells_new,
82+
std::vector<double>(n_q_points,
83+
std::numeric_limits<double>::signaling_NaN()));
84+
cell_data_transfer.unpack(cell_quad_point_values_new);
85+
86+
for (auto &cell : tria.active_cell_iterators())
87+
if (cell->is_locally_owned())
88+
{
89+
cell->set_coarsen_flag();
90+
for (unsigned int i = 0; i < n_q_points; ++i)
91+
{
92+
const double quad_point_value =
93+
cell_quad_point_values_new[cell->active_cell_index()][i];
94+
BOOST_TEST(quad_point_value >= 1.);
95+
BOOST_TEST(quad_point_value <= 2.);
96+
}
97+
}
98+
99+
cell_data_transfer.prepare_for_coarsening_and_refinement(
100+
cell_quad_point_values_new);
101+
tria.execute_coarsening_and_refinement();
102+
103+
BOOST_TEST(tria.n_global_active_cells() == n_active_cells);
104+
std::vector<std::vector<double>> cell_quad_point_values_final(
105+
n_active_cells,
106+
std::vector<double>(n_q_points,
107+
std::numeric_limits<double>::signaling_NaN()));
108+
109+
cell_data_transfer.unpack(cell_quad_point_values_final);
110+
111+
for (auto &cell : tria.active_cell_iterators())
112+
if (cell->is_locally_owned())
113+
{
114+
cell->set_refine_flag();
115+
for (unsigned int i = 0; i < n_q_points; ++i)
116+
{
117+
BOOST_TEST(cell_quad_point_values_final[cell->active_cell_index()][i] ==
118+
cell_quad_point_values[cell->active_cell_index()][i]);
119+
}
120+
}
121+
}
122+
123+
BOOST_AUTO_TEST_CASE(closest_source_point_adaptation)
124+
{
125+
test<2, 2>();
126+
test<3, 3>();
127+
}

0 commit comments

Comments
 (0)