From: Rene Gassmoeller Date: Thu, 28 May 2020 00:28:09 +0000 (-0700) Subject: Add vector and tensor properties X-Git-Tag: v9.3.0-rc1~1523^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a6610fabe17dfe9b407265650f3511bd7464a7bc;p=dealii.git Add vector and tensor properties --- diff --git a/doc/news/changes/minor/20200527Gassmoeller b/doc/news/changes/minor/20200527Gassmoeller new file mode 100644 index 0000000000..8464738393 --- /dev/null +++ b/doc/news/changes/minor/20200527Gassmoeller @@ -0,0 +1,5 @@ +New: The class Particles::DataOut can now output particle properties as +scalars, vectors, or tensors, depending on the arguments handed over to the +Particles::DataOut::build_patches() function. +
+(Rene Gassmoeller, 2020/05/27) diff --git a/include/deal.II/particles/data_out.h b/include/deal.II/particles/data_out.h index cce7a0718e..cfbc789f56 100644 --- a/include/deal.II/particles/data_out.h +++ b/include/deal.II/particles/data_out.h @@ -64,15 +64,22 @@ namespace Particles * A dim=0 patch is built for each particle. The position of the particle is * used to build the node position and the ID of the particle is added as a * single data element. + * @param [in] data_component_names An optional vector of strings that + * describe the properties of each particle. Particle properties will only + * be written if this vector + * is provided. + * @param [in] data_component_interpretations An optional vector that + * controls if the particle properties are interpreted as scalars, vectors, + * or tensors. Has to be of the same length as @p data_component_names. * * @author Bruno Blais, Luca Heltai 2019 */ void build_patches(const Particles::ParticleHandler &particles, - const std::vector & names = {}, + const std::vector &data_component_names = {}, const std::vector< DataComponentInterpretation::DataComponentInterpretation> - &data_component_interpretation = {}); + &data_component_interpretations = {}); protected: /** @@ -89,18 +96,40 @@ namespace Particles virtual std::vector get_dataset_names() const override; + + /** + * Overload of the respective DataOutInterface::get_nonscalar_data_ranges() + * function. See there for a more extensive documentation. + * This function is a reimplementation of the function + * DataOut_DoFData::get_nonscalar_data_ranges(). + */ + virtual std::vector< + std::tuple> + get_nonscalar_data_ranges() const override; + private: /** - * This is a list of patches that is created each time build_patches() is + * This is a vector of patches that is created each time build_patches() is * called. These patches are used in the output routines of the base * classes. */ std::vector> patches; /** - * A list of field names for all data components stored in patches. + * A vector of field names for all data components stored in patches. */ std::vector dataset_names; + + /** + * A vector that for each of the data components of the + * current data set indicates whether they are scalar fields, parts of a + * vector-field, or any of the other supported kinds of data. + */ + std::vector + data_component_interpretations; }; } // namespace Particles diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index e35adfab56..a2642d1ef3 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -5421,7 +5421,7 @@ namespace DataOutBase for (unsigned int c = 0; c < size; ++c) { const auto ind = - Tensor<2, dim>::unrolled_to_component_indices(c); + Tensor<2, spacedim>::unrolled_to_component_indices(c); vtk_data[ind[0]][ind[1]] = data_vectors(first_component + c, n); } diff --git a/source/particles/data_out.cc b/source/particles/data_out.cc index 3ff8ab4bae..99d20548ed 100644 --- a/source/particles/data_out.cc +++ b/source/particles/data_out.cc @@ -16,6 +16,9 @@ #include #include +// We use some exceptions declared in this header +#include + DEAL_II_NAMESPACE_OPEN namespace Particles @@ -26,10 +29,10 @@ namespace Particles const Particles::ParticleHandler &particles, const std::vector & data_component_names, const std::vector - &data_component_interpretation) + &data_component_interpretations_) { - AssertThrow( - data_component_names.size() == data_component_interpretation.size(), + Assert( + data_component_names.size() == data_component_interpretations_.size(), ExcMessage( "When calling Particles::DataOut::build_patches with data component " "names and interpretations you need to provide as many data component " @@ -38,7 +41,7 @@ namespace Particles if (data_component_names.size() > 0) { - AssertThrow( + Assert( data_component_names.size() == particles.begin()->get_properties().size(), ExcMessage( @@ -53,6 +56,14 @@ namespace Particles data_component_names.begin(), data_component_names.end()); + data_component_interpretations.clear(); + data_component_interpretations.emplace_back( + DataComponentInterpretation::component_is_scalar); + data_component_interpretations.insert( + data_component_interpretations.end(), + data_component_interpretations_.begin(), + data_component_interpretations_.end()); + const unsigned int n_property_components = data_component_names.size(); const unsigned int n_data_components = dataset_names.size(); @@ -100,6 +111,140 @@ namespace Particles return dataset_names; } + + + template + std::vector< + std::tuple> + DataOut::get_nonscalar_data_ranges() const + { + std::vector< + std::tuple> + ranges; + + // Make sure the data structures were set up correctly. Since they + // can only be filled by build_patches() above, they should have + // been checked already. + Assert(dataset_names.size() == data_component_interpretations.size(), + ExcInternalError()); + + // collect the ranges of particle data + const unsigned int n_output_components = + data_component_interpretations.size(); + unsigned int output_component = 0; + for (unsigned int i = 0; i < n_output_components; /* i is updated below */) + // see what kind of data we have here. note that for the purpose of the + // current function all we care about is vector data + switch (data_component_interpretations[i]) + { + case DataComponentInterpretation::component_is_scalar: + { + // Just move component forward by one + ++i; + ++output_component; + + break; + } + case DataComponentInterpretation::component_is_part_of_vector: + { + // ensure that there is a continuous number of next space_dim + // components that all deal with vectors + Assert( + i + spacedim <= n_output_components, + Exceptions::DataOutImplementation::ExcInvalidVectorDeclaration( + i, dataset_names[i])); + for (unsigned int dd = 1; dd < spacedim; ++dd) + Assert( + data_component_interpretations[i + dd] == + DataComponentInterpretation::component_is_part_of_vector, + Exceptions::DataOutImplementation:: + ExcInvalidVectorDeclaration(i, dataset_names[i])); + + // all seems right, so figure out whether there is a common + // name to these components. if not, leave the name empty and + // let the output format writer decide what to do here + std::string name = dataset_names[i]; + for (unsigned int dd = 1; dd < spacedim; ++dd) + if (name != dataset_names[i + dd]) + { + name = ""; + break; + } + + // Finally add a corresponding range. + // + // This sort of logic is also explained in some detail in + // DataOut::build_one_patch(). + ranges.emplace_back(std::forward_as_tuple( + output_component, + output_component + spacedim - 1, + name, + DataComponentInterpretation::component_is_part_of_vector)); + + // increase the 'component' counter by the appropriate amount, + // same for 'i', since we have already dealt with all these + // components + output_component += spacedim; + i += spacedim; + + break; + } + + case DataComponentInterpretation::component_is_part_of_tensor: + { + const unsigned int size = spacedim * spacedim; + // ensure that there is a continuous number of next + // spacedim*spacedim components that all deal with tensors + Assert( + i + size <= n_output_components, + Exceptions::DataOutImplementation::ExcInvalidTensorDeclaration( + i, dataset_names[i])); + for (unsigned int dd = 1; dd < size; ++dd) + Assert( + data_component_interpretations[i + dd] == + DataComponentInterpretation::component_is_part_of_tensor, + Exceptions::DataOutImplementation:: + ExcInvalidTensorDeclaration(i, dataset_names[i])); + + // all seems right, so figure out whether there is a common + // name to these components. if not, leave the name empty and + // let the output format writer decide what to do here + std::string name = dataset_names[i]; + for (unsigned int dd = 1; dd < size; ++dd) + if (name != dataset_names[i + dd]) + { + name = ""; + break; + } + + // Finally add a corresponding range. + ranges.emplace_back(std::forward_as_tuple( + output_component, + output_component + size - 1, + name, + DataComponentInterpretation::component_is_part_of_tensor)); + + // increase the 'component' counter by the appropriate amount, + // same for 'i', since we have already dealt with all these + // components + output_component += size; + i += size; + break; + } + + default: + Assert(false, ExcNotImplemented()); + } + + return ranges; + } + } // namespace Particles #include "data_out.inst" diff --git a/tests/particles/data_out_03.cc b/tests/particles/data_out_03.cc new file mode 100644 index 0000000000..4644f346b2 --- /dev/null +++ b/tests/particles/data_out_03.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 vector 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_part_of_vector); + } + + Particles::DataOut particle_output; + particle_output.build_patches(particle_handler, + data_names, + data_interpretations); + particle_output.write_vtk(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_03.output b/tests/particles/data_out_03.output new file mode 100644 index 0000000000..25107a6db6 --- /dev/null +++ b/tests/particles/data_out_03.output @@ -0,0 +1,100 @@ + +# vtk DataFile Version 3.0 +#This file was generated +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 4 double +0.00000 0.00000 0 +1.00000 0.00000 0 +0.00000 1.00000 0 +1.00000 1.00000 0 + +CELLS 4 8 +1 0 +1 1 +1 2 +1 3 + +CELL_TYPES 4 + 1 1 1 1 +POINT_DATA 4 +VECTORS property_coord_0__property_coord_1 double +0.00000 0.00000 0 +1.00000 0.00000 0 +0.00000 1.00000 0 +1.00000 1.00000 0 +SCALARS id double 1 +LOOKUP_TABLE default +0.00000 2.00000 1.00000 3.00000 +DEAL:2d/2d::OK +# vtk DataFile Version 3.0 +#This file was generated +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 4 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +1.00000 1.00000 0.00000 + +CELLS 4 8 +1 0 +1 1 +1 2 +1 3 + +CELL_TYPES 4 + 1 1 1 1 +POINT_DATA 4 +VECTORS property_coord_0__property_coord_1__property_coord_2 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +1.00000 1.00000 0.00000 +SCALARS id double 1 +LOOKUP_TABLE default +0.00000 2.00000 1.00000 3.00000 +DEAL:2d/3d::OK +# vtk DataFile Version 3.0 +#This file was generated +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 8 double +0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 1.00000 +0.00000 1.00000 1.00000 +1.00000 1.00000 1.00000 +1.00000 1.00000 1.00000 + +CELLS 8 16 +1 0 +1 1 +1 2 +1 3 +1 4 +1 5 +1 6 +1 7 + +CELL_TYPES 8 + 1 1 1 1 1 1 1 1 +POINT_DATA 8 +VECTORS property_coord_0__property_coord_1__property_coord_2 double +0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 1.00000 +0.00000 1.00000 1.00000 +1.00000 1.00000 1.00000 +1.00000 1.00000 1.00000 +SCALARS id double 1 +LOOKUP_TABLE default +0.00000 1.00000 4.00000 5.00000 2.00000 3.00000 6.00000 7.00000 +DEAL:3d/3d::OK diff --git a/tests/particles/data_out_04.cc b/tests/particles/data_out_04.cc new file mode 100644 index 0000000000..42c30437bc --- /dev/null +++ b/tests/particles/data_out_04.cc @@ -0,0 +1,160 @@ +// --------------------------------------------------------------------- +// +// 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 tensor 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 * spacedim; ++i) + pit->get_properties()[i] = pit->get_location()[i % spacedim]; + } + } +} + +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 * spacedim); + + create_regular_particle_distribution(particle_handler, tr, 2); + + std::vector data_names; + std::vector + data_interpretations; + + for (unsigned int i = 0; i < spacedim * spacedim; ++i) + { + data_names.emplace_back("property_coord"); + data_interpretations.emplace_back( + DataComponentInterpretation::component_is_part_of_tensor); + } + + Particles::DataOut particle_output; + particle_output.build_patches(particle_handler, + data_names, + data_interpretations); + + // Write gnuplot to check the component output + particle_output.write_gnuplot(deallog.get_file_stream()); + + // Write to vtu to check that vtu can output tensors + std::stringstream test_stream; + particle_output.write_vtu(test_stream); + + // Write pvtu to check the component interpretation + // We can not write vtk like in data_out_03, because vtk + // only supports vectors, but not tensors (at the moment of writing this + // test) + particle_output.write_pvtu_record(deallog.get_file_stream(), + {"solution-01"}); + } + + 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_04.output b/tests/particles/data_out_04.output new file mode 100644 index 0000000000..8a69758b56 --- /dev/null +++ b/tests/particles/data_out_04.output @@ -0,0 +1,105 @@ + +# 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 2.00000 1.00000 0.00000 1.00000 0.00000 + +0.00000 1.00000 1.00000 0.00000 1.00000 0.00000 1.00000 + +1.00000 1.00000 3.00000 1.00000 1.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 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 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 + +0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 + +1.00000 1.00000 0.00000 3.00000 1.00000 1.00000 0.00000 1.00000 1.00000 0.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 0.00000 0.00000 0.00000 + +0.00000 0.00000 0.00000 1.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 4.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +1.00000 0.00000 0.00000 5.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +0.00000 1.00000 1.00000 2.00000 0.00000 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +0.00000 1.00000 1.00000 3.00000 0.00000 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +1.00000 1.00000 1.00000 6.00000 1.00000 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + +1.00000 1.00000 1.00000 7.00000 1.00000 1.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + + + + + + + + + + + + + + + +DEAL:3d/3d::OK