From 50125914480af7083c5b5c17ac03264bca3b439b Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 27 May 2020 11:50:35 -0700 Subject: [PATCH] Add scalar properties --- include/deal.II/particles/data_out.h | 8 +- source/particles/data_out.cc | 45 +++++++- tests/particles/data_out_02.cc | 148 +++++++++++++++++++++++++++ tests/particles/data_out_02.output | 57 +++++++++++ 4 files changed, 255 insertions(+), 3 deletions(-) create mode 100644 tests/particles/data_out_02.cc create mode 100644 tests/particles/data_out_02.output diff --git a/include/deal.II/particles/data_out.h b/include/deal.II/particles/data_out.h index 745a8aaf82..cce7a0718e 100644 --- a/include/deal.II/particles/data_out.h +++ b/include/deal.II/particles/data_out.h @@ -19,6 +19,8 @@ #include +#include + #include #include @@ -66,7 +68,11 @@ namespace Particles * @author Bruno Blais, Luca Heltai 2019 */ void - build_patches(const Particles::ParticleHandler &particles); + build_patches(const Particles::ParticleHandler &particles, + const std::vector & names = {}, + const std::vector< + DataComponentInterpretation::DataComponentInterpretation> + &data_component_interpretation = {}); protected: /** diff --git a/source/particles/data_out.cc b/source/particles/data_out.cc index b359db11b0..3ff8ab4bae 100644 --- a/source/particles/data_out.cc +++ b/source/particles/data_out.cc @@ -23,10 +23,39 @@ namespace Particles template void DataOut::build_patches( - const Particles::ParticleHandler &particles) + const Particles::ParticleHandler &particles, + const std::vector & data_component_names, + const std::vector + &data_component_interpretation) { + AssertThrow( + data_component_names.size() == data_component_interpretation.size(), + ExcMessage( + "When calling Particles::DataOut::build_patches with data component " + "names and interpretations you need to provide as many data component " + "names as interpretations. Provide the same name for components that " + "belong to a single vector or tensor.")); + + if (data_component_names.size() > 0) + { + AssertThrow( + data_component_names.size() == + particles.begin()->get_properties().size(), + ExcMessage( + "When calling Particles::DataOut::build_patches with data component " + "names and interpretations you need to provide as many data component " + "names as the particles have properties.")); + } + dataset_names.clear(); dataset_names.emplace_back("id"); + dataset_names.insert(dataset_names.end(), + data_component_names.begin(), + data_component_names.end()); + + const unsigned int n_property_components = data_component_names.size(); + const unsigned int n_data_components = dataset_names.size(); + patches.resize(particles.n_locally_owned_particles()); auto particle = particles.begin(); @@ -35,9 +64,21 @@ namespace Particles patches[i].vertices[0] = particle->get_location(); patches[i].patch_index = i; patches[i].n_subdivisions = 1; - patches[i].data.reinit(dataset_names.size(), 1); + + // We have one more data components than dataset_names (the particle id) + patches[i].data.reinit(n_data_components, 1); patches[i].data(0, 0) = particle->get_id(); + + if (n_data_components > 1) + { + const ArrayView properties = particle->get_properties(); + for (unsigned int property_index = 0; + property_index < n_property_components; + ++property_index) + patches[i].data(property_index + 1, 0) = + properties[property_index]; + } } } diff --git a/tests/particles/data_out_02.cc b/tests/particles/data_out_02.cc new file mode 100644 index 0000000000..6c0ada09ef --- /dev/null +++ b/tests/particles/data_out_02.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Create a particle handler then generate an output +// Tests that scalar particle properties are output correctly. + +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +create_regular_particle_distribution( + Particles::ParticleHandler &particle_handler, + const Triangulation & tr, + const unsigned int particles_per_direction = 3) +{ + for (unsigned int i = 0; i < particles_per_direction; ++i) + for (unsigned int j = 0; j < particles_per_direction; ++j) + { + Point position; + Point reference_position; + unsigned int id = i * particles_per_direction + j; + + position[0] = static_cast(i) / + static_cast(particles_per_direction - 1); + position[1] = static_cast(j) / + static_cast(particles_per_direction - 1); + + if (dim > 2) + for (unsigned int k = 0; k < particles_per_direction; ++k) + { + position[2] = static_cast(j) / + static_cast(particles_per_direction - 1); + id = i * particles_per_direction * particles_per_direction + + j * particles_per_direction + k; + Particles::Particle particle(position, + reference_position, + id); + + typename Triangulation::active_cell_iterator cell = + GridTools::find_active_cell_around_point( + tr, particle.get_location()); + + Particles::ParticleIterator pit = + particle_handler.insert_particle(particle, cell); + + for (unsigned int i = 0; i < spacedim; ++i) + pit->get_properties()[i] = pit->get_location()[i]; + } + else + { + Particles::Particle particle(position, + reference_position, + id); + + typename Triangulation::active_cell_iterator cell = + GridTools::find_active_cell_around_point(tr, + particle.get_location()); + + Particles::ParticleIterator pit = + particle_handler.insert_particle(particle, cell); + + for (unsigned int i = 0; i < spacedim; ++i) + pit->get_properties()[i] = pit->get_location()[i]; + } + } +} + +template +void +test() +{ + { + Triangulation tr; + + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + MappingQ mapping(1); + + // Create a particle handler using two manually created particles + Particles::ParticleHandler particle_handler(tr, + mapping, + spacedim); + + create_regular_particle_distribution(particle_handler, tr, 2); + + std::vector data_names; + std::vector + data_interpretations; + + for (unsigned int i = 0; i < spacedim; ++i) + { + data_names.emplace_back("property_coord_" + std::to_string(i)); + data_interpretations.emplace_back( + DataComponentInterpretation::component_is_scalar); + } + + Particles::DataOut particle_output; + particle_output.build_patches(particle_handler, + data_names, + data_interpretations); + particle_output.write_gnuplot(deallog.get_file_stream()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/data_out_02.output b/tests/particles/data_out_02.output new file mode 100644 index 0000000000..27e7f2288b --- /dev/null +++ b/tests/particles/data_out_02.output @@ -0,0 +1,57 @@ + +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 0.00000 0.00000 + +1.00000 0.00000 2.00000 1.00000 0.00000 + +0.00000 1.00000 1.00000 0.00000 1.00000 + +1.00000 1.00000 3.00000 1.00000 1.00000 + +DEAL:2d/2d::OK +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +1.00000 0.00000 0.00000 2.00000 1.00000 0.00000 0.00000 + +0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 + +1.00000 1.00000 0.00000 3.00000 1.00000 1.00000 0.00000 + +DEAL:2d/3d::OK +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 + +1.00000 0.00000 0.00000 4.00000 1.00000 0.00000 0.00000 + +1.00000 0.00000 0.00000 5.00000 1.00000 0.00000 0.00000 + +0.00000 1.00000 1.00000 2.00000 0.00000 1.00000 1.00000 + +0.00000 1.00000 1.00000 3.00000 0.00000 1.00000 1.00000 + +1.00000 1.00000 1.00000 6.00000 1.00000 1.00000 1.00000 + +1.00000 1.00000 1.00000 7.00000 1.00000 1.00000 1.00000 + +DEAL:3d/3d::OK -- 2.39.5