Skip to content

Commit adfaad9

Browse files
bangerthbangerth
authored
Merge pull request #5857 from gassmoeller/simplify_particle_interpolators
Remove unnecessary check in particle interpolator. Remove deprecated function.
2 parents f7934a5 + 13e020f commit adfaad9

File tree

7 files changed

+35
-184
lines changed

7 files changed

+35
-184
lines changed

‎include/aspect/particle/interpolator/interface.h‎

Lines changed: 7 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -50,31 +50,6 @@ namespace aspect
5050
class Interface : public Plugins::InterfaceBase
5151
{
5252
public:
53-
/**
54-
* Perform an interpolation of the properties of the particles in
55-
* this cell onto a vector of positions in this cell.
56-
* Implementations of this function must return a vector of a vector
57-
* of doubles which contains a somehow computed
58-
* value of all particle properties at all given positions.
59-
*
60-
* @param [in] particle_handler Reference to the particle handler
61-
* that allows accessing the particles in the domain.
62-
* @param [in] positions The vector of positions where the properties
63-
* should be evaluated.
64-
* @param [in] cell An optional iterator to the cell containing the
65-
* particles. Not all callers will know the cell of the particles,
66-
* but providing the cell when known speeds up the interpolation
67-
* significantly.
68-
* @return A vector with as many entries as @p positions. Every entry
69-
* is a vector of interpolated particle properties at this position.
70-
*/
71-
DEAL_II_DEPRECATED
72-
virtual
73-
std::vector<std::vector<double>>
74-
properties_at_points(const ParticleHandler<dim> &particle_handler,
75-
const std::vector<Point<dim>> &positions,
76-
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell = typename parallel::distributed::Triangulation<dim>::active_cell_iterator()) const;
77-
7853
/**
7954
* Perform an interpolation of the properties of the particles in
8055
* this cell onto a vector of positions in this cell.
@@ -92,19 +67,20 @@ namespace aspect
9267
* should be evaluated.
9368
* @param [in] selected_properties A component mask that determines
9469
* which particle properties are interpolated in this function.
95-
* @param [in] cell An optional iterator to the cell containing the
96-
* particles. Not all callers will know the cell of the particles,
97-
* but providing the cell when known speeds up the interpolation
98-
* significantly.
99-
* @return A vector with as many entries as @p positions. Every entry
70+
* @param [in] cell An iterator to the cell containing the
71+
* positions.
72+
* @return A vector with as many entries as @p positions. Each entry
10073
* is a vector of interpolated particle properties at this position.
74+
* This property vector has as many entries as there are particle
75+
* properties, however entries that have not been selected in
76+
* @p selected_properties are filled with signalling NaNs.
10177
*/
10278
virtual
10379
std::vector<std::vector<double>>
10480
properties_at_points(const ParticleHandler<dim> &particle_handler,
10581
const std::vector<Point<dim>> &positions,
10682
const ComponentMask &selected_properties,
107-
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell = typename parallel::distributed::Triangulation<dim>::active_cell_iterator()) const = 0;
83+
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const = 0;
10884
};
10985

11086

‎source/particle/interpolator/bilinear_least_squares.cc‎

Lines changed: 11 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -48,32 +48,8 @@ namespace aspect
4848
ExcMessage("Internal error: the particle property interpolator was "
4949
"called without a specified component to interpolate."));
5050

51-
const Point<dim> approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>())
52-
/ static_cast<double> (positions.size());
53-
54-
typename parallel::distributed::Triangulation<dim>::active_cell_iterator found_cell;
55-
56-
if (cell == typename parallel::distributed::Triangulation<dim>::active_cell_iterator())
57-
{
58-
// We can not simply use one of the points as input for find_active_cell_around_point
59-
// because for vertices of mesh cells we might end up getting ghost_cells as return value
60-
// instead of the local active cell. So make sure we are well in the inside of a cell.
61-
Assert(positions.size() > 0,
62-
ExcMessage("The particle property interpolator was not given any "
63-
"positions to evaluate the particle cell_properties at."));
64-
65-
66-
found_cell =
67-
(GridTools::find_active_cell_around_point<> (this->get_mapping(),
68-
this->get_triangulation(),
69-
approximated_cell_midpoint)).first;
70-
}
71-
else
72-
found_cell = cell;
73-
7451
const typename ParticleHandler<dim>::particle_iterator_range particle_range =
75-
particle_handler.particles_in_cell(found_cell);
76-
52+
particle_handler.particles_in_cell(cell);
7753

7854
std::vector<std::vector<double>> cell_properties(positions.size(),
7955
std::vector<double>(n_particle_properties,
@@ -88,7 +64,7 @@ namespace aspect
8864
return fallback_interpolator.properties_at_points(particle_handler,
8965
positions,
9066
selected_properties,
91-
found_cell);
67+
cell);
9268

9369
// Noticed that the size of matrix A is n_particles x n_matrix_columns
9470
// which usually is not a square matrix. Therefore, we find the
@@ -138,7 +114,7 @@ namespace aspect
138114
if (use_linear_least_squares_limiter.n_selected_components() != 0)
139115
{
140116
std::vector<typename parallel::distributed::Triangulation<dim>::active_cell_iterator> active_neighbors;
141-
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(found_cell, active_neighbors);
117+
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(cell, active_neighbors);
142118
for (const auto &active_neighbor : active_neighbors)
143119
{
144120
if (active_neighbor->is_artificial())
@@ -153,23 +129,23 @@ namespace aspect
153129
}
154130
}
155131
}
156-
if (found_cell->at_boundary())
132+
if (cell->at_boundary())
157133
{
158-
const std::vector<double> cell_average_values = fallback_interpolator.properties_at_points(particle_handler, {positions[0]}, selected_properties, found_cell)[0];
159-
for (unsigned int face_id = 0; face_id < found_cell->reference_cell().n_faces(); ++face_id)
134+
const std::vector<double> cell_average_values = fallback_interpolator.properties_at_points(particle_handler, {positions[0]}, selected_properties, cell)[0];
135+
for (unsigned int face_id = 0; face_id < cell->reference_cell().n_faces(); ++face_id)
160136
{
161-
if (found_cell->at_boundary(face_id))
137+
if (cell->at_boundary(face_id))
162138
{
163139
const unsigned int opposing_face_id = GeometryInfo<dim>::opposite_face[face_id];
164-
const auto &opposing_cell = found_cell->neighbor(opposing_face_id);
140+
const auto &opposing_cell = cell->neighbor(opposing_face_id);
165141
if (opposing_cell.state() == IteratorState::IteratorStates::valid && opposing_cell->is_active() && opposing_cell->is_artificial() == false)
166142
{
167143
const auto neighbor_cell_average = fallback_interpolator.properties_at_points(particle_handler, {positions[0]}, selected_properties, opposing_cell)[0];
168144
for (unsigned int property_index = 0; property_index < n_particle_properties; ++property_index)
169145
{
170146
if (selected_properties[property_index] == true && use_boundary_extrapolation[property_index] == true)
171147
{
172-
Assert(found_cell->reference_cell().is_hyper_cube() == true, ExcNotImplemented());
148+
Assert(cell->reference_cell().is_hyper_cube() == true, ExcNotImplemented());
173149
const double expected_boundary_value = 1.5 * cell_average_values[property_index] - 0.5 * neighbor_cell_average[property_index];
174150
property_minimums[property_index] = std::min(property_minimums[property_index], expected_boundary_value);
175151
property_maximums[property_index] = std::max(property_maximums[property_index], expected_boundary_value);
@@ -192,7 +168,7 @@ namespace aspect
192168
return fallback_interpolator.properties_at_points(particle_handler,
193169
positions,
194170
selected_properties,
195-
found_cell);
171+
cell);
196172

197173
std::vector<Vector<double>> QTb(n_particle_properties, Vector<double>(n_matrix_columns));
198174
std::vector<Vector<double>> c(n_particle_properties, Vector<double>(n_matrix_columns));
@@ -233,7 +209,7 @@ namespace aspect
233209
std::size_t positions_index = 0;
234210
for (typename std::vector<Point<dim>>::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++positions_index)
235211
{
236-
Point<dim> relative_support_point_location = this->get_mapping().transform_real_to_unit_cell(found_cell, *itr);
212+
Point<dim> relative_support_point_location = this->get_mapping().transform_real_to_unit_cell(cell, *itr);
237213
double interpolated_value = c[property_index][0];
238214
for (unsigned int i = 1; i < n_matrix_columns; ++i)
239215
{

‎source/particle/interpolator/cell_average.cc‎

Lines changed: 2 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -38,30 +38,8 @@ namespace aspect
3838
const ComponentMask &selected_properties,
3939
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
4040
{
41-
typename parallel::distributed::Triangulation<dim>::active_cell_iterator found_cell;
42-
43-
if (cell == typename parallel::distributed::Triangulation<dim>::active_cell_iterator())
44-
{
45-
// We can not simply use one of the points as input for find_active_cell_around_point
46-
// because for vertices of mesh cells we might end up getting ghost_cells as return value
47-
// instead of the local active cell. So make sure we are well in the inside of a cell.
48-
Assert(positions.size() > 0,
49-
ExcMessage("The particle property interpolator was not given any "
50-
"positions to evaluate the particle properties at."));
51-
52-
const Point<dim> approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>())
53-
/ static_cast<double> (positions.size());
54-
55-
found_cell =
56-
(GridTools::find_active_cell_around_point<> (this->get_mapping(),
57-
this->get_triangulation(),
58-
approximated_cell_midpoint)).first;
59-
}
60-
else
61-
found_cell = cell;
62-
6341
const typename ParticleHandler<dim>::particle_iterator_range particle_range =
64-
particle_handler.particles_in_cell(found_cell);
42+
particle_handler.particles_in_cell(cell);
6543

6644
const unsigned int n_particles = std::distance(particle_range.begin(),particle_range.end());
6745
const unsigned int n_particle_properties = particle_handler.n_properties_per_particle();
@@ -92,7 +70,7 @@ namespace aspect
9270
else
9371
{
9472
std::vector<typename parallel::distributed::Triangulation<dim>::active_cell_iterator> neighbors;
95-
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(found_cell,neighbors);
73+
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(cell,neighbors);
9674

9775
unsigned int non_empty_neighbors = 0;
9876
for (unsigned int i=0; i<neighbors.size(); ++i)

‎source/particle/interpolator/harmonic_average.cc‎

Lines changed: 2 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -39,30 +39,8 @@ namespace aspect
3939
const ComponentMask &selected_properties,
4040
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
4141
{
42-
typename parallel::distributed::Triangulation<dim>::active_cell_iterator found_cell;
43-
44-
if (cell == typename parallel::distributed::Triangulation<dim>::active_cell_iterator())
45-
{
46-
// We can not simply use one of the points as input for find_active_cell_around_point
47-
// because for vertices of mesh cells we might end up getting ghost_cells as return value
48-
// instead of the local active cell. So make sure we are well in the inside of a cell.
49-
Assert(positions.size() > 0,
50-
ExcMessage("The particle property interpolator was not given any "
51-
"positions to evaluate the particle properties at."));
52-
53-
const Point<dim> approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>())
54-
/ static_cast<double> (positions.size());
55-
56-
found_cell =
57-
(GridTools::find_active_cell_around_point<> (this->get_mapping(),
58-
this->get_triangulation(),
59-
approximated_cell_midpoint)).first;
60-
}
61-
else
62-
found_cell = cell;
63-
6442
const typename ParticleHandler<dim>::particle_iterator_range particle_range =
65-
particle_handler.particles_in_cell(found_cell);
43+
particle_handler.particles_in_cell(cell);
6644

6745
const unsigned int n_particles = std::distance(particle_range.begin(),particle_range.end());
6846
const unsigned int n_particle_properties = particle_handler.n_properties_per_particle();
@@ -99,7 +77,7 @@ namespace aspect
9977
else
10078
{
10179
std::vector<typename parallel::distributed::Triangulation<dim>::active_cell_iterator> neighbors;
102-
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(found_cell,neighbors);
80+
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(cell,neighbors);
10381

10482
unsigned int non_empty_neighbors = 0;
10583
for (unsigned int i=0; i<neighbors.size(); ++i)

‎source/particle/interpolator/interface.cc‎

Lines changed: 0 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -27,17 +27,6 @@ namespace aspect
2727
{
2828
namespace Interpolator
2929
{
30-
template <int dim>
31-
std::vector<std::vector<double>>
32-
Interface<dim>::properties_at_points(const ParticleHandler<dim> &particle_handler,
33-
const std::vector<Point<dim>> &positions,
34-
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
35-
{
36-
return properties_at_points(particle_handler,positions,ComponentMask(), cell);
37-
}
38-
39-
40-
4130
// -------------------------------- Deal with registering models and automating
4231
// -------------------------------- their setup and selection at run time
4332

‎source/particle/interpolator/nearest_neighbor.cc‎

Lines changed: 2 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -37,29 +37,7 @@ namespace aspect
3737
const ComponentMask &selected_properties,
3838
const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell) const
3939
{
40-
typename parallel::distributed::Triangulation<dim>::active_cell_iterator found_cell;
41-
42-
if (cell->state() == IteratorState::invalid)
43-
{
44-
// We can not simply use one of the points as input for find_active_cell_around_point
45-
// because for vertices of mesh cells we might end up getting ghost_cells as return value
46-
// instead of the local active cell. So make sure we are well in the inside of a cell.
47-
Assert(positions.size() > 0,
48-
ExcMessage("The particle property interpolator was not given any "
49-
"positions to evaluate the particle properties at."));
50-
51-
const Point<dim> approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point<dim>())
52-
/ static_cast<double> (positions.size());
53-
54-
found_cell =
55-
(GridTools::find_active_cell_around_point<> (this->get_mapping(),
56-
this->get_triangulation(),
57-
approximated_cell_midpoint)).first;
58-
}
59-
else
60-
found_cell = cell;
61-
62-
const typename ParticleHandler<dim>::particle_iterator_range particle_range = particle_handler.particles_in_cell(found_cell);
40+
const typename ParticleHandler<dim>::particle_iterator_range particle_range = particle_handler.particles_in_cell(cell);
6341

6442
const unsigned int n_particles = std::distance(particle_range.begin(),particle_range.end());
6543
const unsigned int n_particle_properties = particle_handler.n_properties_per_particle();
@@ -91,7 +69,7 @@ namespace aspect
9169
else
9270
{
9371
std::vector<typename parallel::distributed::Triangulation<dim>::active_cell_iterator> neighbors;
94-
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(found_cell,neighbors);
72+
GridTools::get_active_neighbors<parallel::distributed::Triangulation<dim>>(cell,neighbors);
9573

9674
unsigned int nearest_neighbor_cell = numbers::invalid_unsigned_int;
9775
for (unsigned int i=0; i<neighbors.size(); ++i)

0 commit comments

Comments
 (0)