Skip to content

Commit

Permalink
Address comments
Browse files Browse the repository at this point in the history
  • Loading branch information
gassmoeller committed Aug 16, 2016
1 parent f790fda commit 54bea5d
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 20 deletions.
10 changes: 5 additions & 5 deletions include/aspect/particle/world.h
Original file line number Diff line number Diff line change
Expand Up @@ -376,7 +376,7 @@ namespace aspect
* process and in its current cell, or deleted (if it could not find
* its new process or cell).
*
* @param [in,out] particles_to_sort Vector containing all pairs of
* @param [in] particles_to_sort Vector containing all pairs of
* particles and their old cells that will be sorted into the
* 'particles' member variable in this function.
*/
Expand All @@ -399,9 +399,9 @@ namespace aspect
* that can not be found are discarded.
*/
void
move_particles_back_into_mesh(std::vector<std::pair<types::LevelInd, Particle<dim> > > &lost_particles,
std::vector<std::pair<types::LevelInd, Particle<dim> > > &moved_particles_cell,
std::multimap<types::subdomain_id, Particle<dim> > &moved_particles_domain);
move_particles_back_into_mesh(const std::vector<std::pair<types::LevelInd, Particle<dim> > > &lost_particles,
std::vector<std::pair<types::LevelInd, Particle<dim> > > &moved_particles_cell,
std::multimap<types::subdomain_id, Particle<dim> > &moved_particles_domain);

/**
* Transfer particles that have crossed subdomain boundaries to other
Expand All @@ -420,7 +420,7 @@ namespace aspect
* @param [in,out] received_particles List that stores all received
* particles. Note that it is not required nor checked that the list
* is empty, received particles are simply attached to the end of
* the list.
* the vector.
*/
void
send_recv_particles(const std::multimap<types::subdomain_id,Particle <dim> > &sent_particles,
Expand Down
31 changes: 16 additions & 15 deletions source/particle/world.cc
Original file line number Diff line number Diff line change
Expand Up @@ -537,8 +537,8 @@ namespace aspect
std::multimap<types::subdomain_id, Particle<dim> > moved_particles_domain;
std::vector<std::pair<types::LevelInd, Particle<dim> > > lost_particles;

sorted_particles.reserve(particles_to_sort.size()*2);
lost_particles.reserve(particles_to_sort.size());
sorted_particles.reserve(static_cast<unsigned int> (particles_to_sort.size()*1.25));
lost_particles.reserve(static_cast<unsigned int> (particles_to_sort.size()*0.25));

{
TimerOutput::Scope timer_section(this->get_computing_timer(), "Particles: Sort");
Expand Down Expand Up @@ -646,9 +646,9 @@ namespace aspect

template <int dim>
void
World<dim>::move_particles_back_into_mesh(std::vector<std::pair<types::LevelInd, Particle<dim> > > &lost_particles,
std::vector<std::pair<types::LevelInd, Particle<dim> > > &moved_particles_cell,
std::multimap<types::subdomain_id, Particle<dim> > &moved_particles_domain)
World<dim>::move_particles_back_into_mesh(const std::vector<std::pair<types::LevelInd, Particle<dim> > > &lost_particles,
std::vector<std::pair<types::LevelInd, Particle<dim> > > &moved_particles_cell,
std::multimap<types::subdomain_id, Particle<dim> > &moved_particles_domain)
{
// TODO: fix this to work with arbitrary meshes. Currently periodic boundaries only work for boxes.
// If the geometry is not a box, we simply discard particles that have left the
Expand All @@ -670,11 +670,12 @@ namespace aspect
for (; boundary != periodic_boundaries.end(); ++boundary)
periodic[boundary->second] = true;

typename std::vector<std::pair<types::LevelInd, Particle<dim> > >::iterator lost_particle = lost_particles.begin();
typename std::vector<std::pair<types::LevelInd, Particle<dim> > >::const_iterator lost_particle = lost_particles.begin();
for (; lost_particle != lost_particles.end(); ++lost_particle)
{
// modify the particle position if it crossed a periodic boundary
Point<dim> particle_position = lost_particle->second.get_location();
std::pair<types::LevelInd, Particle<dim> > cell_and_particle = *lost_particle;
Point<dim> particle_position = cell_and_particle.second.get_location();
for (unsigned int i = 0; i < dim; ++i)
{
if (periodic[i])
Expand All @@ -685,7 +686,7 @@ namespace aspect
particle_position[i] -= extent[i];
}
}
lost_particle->second.set_location(particle_position);
cell_and_particle.second.set_location(particle_position);

// Try again looking for the new cell with the updated position
typename parallel::distributed::Triangulation<dim>::active_cell_iterator cell;
Expand All @@ -695,9 +696,10 @@ namespace aspect
Point<dim> > current_cell_and_position =
GridTools::find_active_cell_around_point<> (this->get_mapping(),
this->get_triangulation(),
lost_particle->second.get_location());
cell_and_particle.second.get_location());
cell = current_cell_and_position.first;
lost_particle->second.set_reference_location(current_cell_and_position.second);
cell_and_particle.first = std::make_pair(cell->level(),cell->index());
cell_and_particle.second.set_reference_location(current_cell_and_position.second);
}
catch (GridTools::ExcPointNotFound<dim> &)
{
Expand All @@ -709,12 +711,9 @@ namespace aspect
// Reinsert the particle into our domain if we found its cell
// Mark it for MPI transfer otherwise
if (cell->is_locally_owned())
{
const types::LevelInd found_cell = std::make_pair(cell->level(),cell->index());
moved_particles_cell.push_back(std::make_pair(found_cell, lost_particle->second));
}
moved_particles_cell.push_back(cell_and_particle);
else
moved_particles_domain.insert(std::make_pair(cell->subdomain_id(),lost_particle->second));
moved_particles_domain.insert(std::make_pair(cell->subdomain_id(),cell_and_particle.second));
}
}
}
Expand Down Expand Up @@ -849,6 +848,8 @@ namespace aspect
recv_particle.set_reference_location(p_unit);
break;
}
// If the particle is not in this cell, do nothing and check
// the next neighbor cell
}
catch (typename Mapping<dim>::ExcTransformationFailed &)
{}
Expand Down

0 comments on commit 54bea5d

Please sign in to comment.