Skip to content

Commit 3afa9fa

Browse files
masterleinadmasterleinad
committed
Add ClosestSourcePointAdaptation
1 parent 842ee50 commit 3afa9fa

File tree

3 files changed

+248
-0
lines changed

3 files changed

+248
-0
lines changed

‎source/ClosestSourcePointAdaptation.hh‎

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

‎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: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
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+
{
65+
return closest_quad_point_adaptation.coarse_to_fine(parent,
66+
parent_values);
67+
},
68+
[&](const typename dealii::Triangulation<dim, spacedim>::cell_iterator
69+
&parent,
70+
const std::vector<std::vector<double>> child_values)
71+
{
72+
return closest_quad_point_adaptation.fine_to_coarse(parent,
73+
child_values);
74+
});
75+
76+
cell_data_transfer.prepare_for_coarsening_and_refinement(
77+
cell_quad_point_values);
78+
tria.execute_coarsening_and_refinement();
79+
80+
const unsigned int n_active_cells_new = tria.n_global_active_cells();
81+
82+
std::vector<std::vector<double>> cell_quad_point_values_new(
83+
n_active_cells_new,
84+
std::vector<double>(n_q_points,
85+
std::numeric_limits<double>::signaling_NaN()));
86+
cell_data_transfer.unpack(cell_quad_point_values_new);
87+
88+
for (auto &cell : tria.active_cell_iterators())
89+
if (cell->is_locally_owned())
90+
{
91+
cell->set_coarsen_flag();
92+
for (unsigned int i = 0; i < n_q_points; ++i)
93+
{
94+
const double quad_point_value =
95+
cell_quad_point_values_new[cell->active_cell_index()][i];
96+
BOOST_TEST(quad_point_value >= 1.);
97+
BOOST_TEST(quad_point_value <= 2.);
98+
}
99+
}
100+
101+
cell_data_transfer.prepare_for_coarsening_and_refinement(
102+
cell_quad_point_values_new);
103+
tria.execute_coarsening_and_refinement();
104+
105+
BOOST_TEST(tria.n_global_active_cells() == n_active_cells);
106+
std::vector<std::vector<double>> cell_quad_point_values_final(
107+
n_active_cells,
108+
std::vector<double>(n_q_points,
109+
std::numeric_limits<double>::signaling_NaN()));
110+
111+
cell_data_transfer.unpack(cell_quad_point_values_final);
112+
113+
for (auto &cell : tria.active_cell_iterators())
114+
if (cell->is_locally_owned())
115+
{
116+
cell->set_refine_flag();
117+
for (unsigned int i = 0; i < n_q_points; ++i)
118+
{
119+
BOOST_TEST(cell_quad_point_values_final[cell->active_cell_index()][i] ==
120+
cell_quad_point_values[cell->active_cell_index()][i]);
121+
}
122+
}
123+
}
124+
125+
BOOST_AUTO_TEST_CASE(closest_source_point_adaptation)
126+
{
127+
test<2, 2>();
128+
test<3, 3>();
129+
}

0 commit comments

Comments
 (0)