Merge pull request #6073 from gassmoeller/fix_multiple_particle_worlds

Fix multiple particle worlds

source/material_model/rheology/strain_dependent.cc

@@ -412,8 +412,11 @@ namespace aspect
  // Currently this functionality only works in field composition
  if (healing_mechanism != no_healing && this->n_particle_worlds() > 0)
  {
- const Particle::Property::Manager<dim> &particle_property_manager = this->get_particle_world(0).get_property_manager();
- AssertThrow(particle_property_manager.plugin_name_exists("viscoplastic strain invariants") == false, ExcMessage("This healing mechanism currently does not work if the strain is tracked on particles."));
+ for (unsigned int i=0; i<this->n_particle_worlds(); ++i)
+ {
+ const Particle::Property::Manager<dim> &particle_property_manager = this->get_particle_world(i).get_property_manager();
+ AssertThrow(particle_property_manager.plugin_name_exists("viscoplastic strain invariants") == false, ExcMessage("This healing mechanism currently does not work if the strain is tracked on particles."));
+ }
  }
 
  // Temperature dependent strain healing requires that adiabatic surface temperature is non zero

source/mesh_refinement/particle_density.cc

@@ -36,16 +36,19 @@ namespace aspect
  "postprocessor plugin `particles' to be selected. Please activate the "
  "particles or deactivate this mesh refinement plugin."));
 
- const Particle::ParticleHandler<dim> &particle_handler = this->get_particle_world(0).get_particle_handler();
-
  for (const auto &cell : this->get_dof_handler().active_cell_iterators())
  if (cell->is_locally_owned())
  {
  // Note that this refinement indicator will level out the number
  // of particles per cell, therefore creating fine cells in regions
  // of high particle density and coarse cells in low particle
  // density regions.
- indicators(cell->active_cell_index()) = static_cast<float>(particle_handler.n_particles_in_cell(cell));
+ indicators(cell->active_cell_index()) = 0.0;
+ for (unsigned int i=0; i<this->n_particle_worlds(); ++i)
+ {
+ const Particle::ParticleHandler<dim> &particle_handler = this->get_particle_world(i).get_particle_handler();
+ indicators(cell->active_cell_index()) += static_cast<float>(particle_handler.n_particles_in_cell(cell));
+ }
  }
  }
  }

source/postprocess/crystal_preferred_orientation.cc

@@ -172,16 +172,30 @@ namespace aspect
  std::pair<std::string,std::string>
  CrystalPreferredOrientation<dim>::execute (TableHandler &statistics)
  {
+ unsigned int particle_world_index = numbers::invalid_unsigned_int;
+ for (unsigned int i = 0; i < this->n_particle_worlds(); ++i)
+ {
+ if (this->get_particle_world(i).get_property_manager().template has_matching_active_plugin<Particle::Property::CrystalPreferredOrientation<dim>>())
+ {
+ Assert(particle_world_index == numbers::invalid_unsigned_int,
+ ExcMessage("Multiple particle worlds with CPO data found. This is not supported."));
 
- const Particle::Property::Manager<dim> &manager = this->get_particle_world(0).get_property_manager();
- const Particle::Property::ParticleHandler<dim> &particle_handler = this->get_particle_world(0).get_particle_handler();
+ particle_world_index = i;
+ }
+ }
 
- const bool cpo_elastic_decomposition_plugin_exists = manager.plugin_name_exists("elastic tensor decomposition");
+ AssertThrow(particle_world_index != numbers::invalid_unsigned_int,
+ ExcMessage("No CPO particle properties found. Make sure that the CPO particle property plugin is selected."));
 
- // Get a reference to the CPO particle property.
+ // get particle data and a reference to the CPO particle property
+ const Particle::Property::Manager<dim> &manager = this->get_particle_world(particle_world_index).get_property_manager();
+ const Particle::Property::ParticleHandler<dim> &particle_handler = this->get_particle_world(particle_world_index).get_particle_handler();
+ const Particle::Property::ParticlePropertyInformation &property_information = manager.get_data_info();
  const Particle::Property::CrystalPreferredOrientation<dim> &cpo_particle_property =
  manager.template get_matching_active_plugin<Particle::Property::CrystalPreferredOrientation<dim>>();
 
+ const bool cpo_elastic_decomposition_plugin_exists = manager.plugin_name_exists("elastic tensor decomposition");
+
  const unsigned int n_grains = cpo_particle_property.get_number_of_grains();
  const unsigned int n_minerals = cpo_particle_property.get_number_of_minerals();
 
@@ -219,14 +233,6 @@ namespace aspect
  + "hexagonal_norm_square_p1 hexagonal_norm_square_p2 hexagonal_norm_square_p3 "
  + "isotropic_norm_square") : "") << std::endl;
 
- // get particle data
- const Particle::Property::ParticlePropertyInformation &property_information = this->get_particle_world(0).get_property_manager().get_data_info();
-
- AssertThrow(property_information.fieldname_exists("cpo mineral 0 type") ,
- ExcMessage("No CPO particle properties found. Make sure that the CPO particle property plugin is selected."));
-
-
-
  const unsigned int cpo_data_position = property_information.n_fields() == 0
  ?
  0