template <int dim, int spacedim>
void
DataOut<dim, spacedim>::build_patches(
- const Particles::ParticleHandler<dim, spacedim> &particles)
+ const Particles::ParticleHandler<dim, spacedim> &particles,
+ const std::vector<std::string> & data_component_names,
+ const std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ &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();
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<double> 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];
+ }
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/data_out_base.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/particles/data_out.h>
+#include <deal.II/particles/particle_handler.h>
+
+#include "../tests.h"
+
+template <int dim, int spacedim>
+void
+create_regular_particle_distribution(
+ Particles::ParticleHandler<dim, spacedim> &particle_handler,
+ const Triangulation<dim, spacedim> & 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<spacedim> position;
+ Point<dim> reference_position;
+ unsigned int id = i * particles_per_direction + j;
+
+ position[0] = static_cast<double>(i) /
+ static_cast<double>(particles_per_direction - 1);
+ position[1] = static_cast<double>(j) /
+ static_cast<double>(particles_per_direction - 1);
+
+ if (dim > 2)
+ for (unsigned int k = 0; k < particles_per_direction; ++k)
+ {
+ position[2] = static_cast<double>(j) /
+ static_cast<double>(particles_per_direction - 1);
+ id = i * particles_per_direction * particles_per_direction +
+ j * particles_per_direction + k;
+ Particles::Particle<dim, spacedim> particle(position,
+ reference_position,
+ id);
+
+ typename Triangulation<dim, spacedim>::active_cell_iterator cell =
+ GridTools::find_active_cell_around_point(
+ tr, particle.get_location());
+
+ Particles::ParticleIterator<dim, spacedim> 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<dim, spacedim> particle(position,
+ reference_position,
+ id);
+
+ typename Triangulation<dim, spacedim>::active_cell_iterator cell =
+ GridTools::find_active_cell_around_point(tr,
+ particle.get_location());
+
+ Particles::ParticleIterator<dim, spacedim> pit =
+ particle_handler.insert_particle(particle, cell);
+
+ for (unsigned int i = 0; i < spacedim; ++i)
+ pit->get_properties()[i] = pit->get_location()[i];
+ }
+ }
+}
+
+template <int dim, int spacedim>
+void
+test()
+{
+ {
+ Triangulation<dim, spacedim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ tr.refine_global(2);
+ MappingQ<dim, spacedim> mapping(1);
+
+ // Create a particle handler using two manually created particles
+ Particles::ParticleHandler<dim, spacedim> particle_handler(tr,
+ mapping,
+ spacedim);
+
+ create_regular_particle_distribution(particle_handler, tr, 2);
+
+ std::vector<std::string> data_names;
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ 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<dim, spacedim> 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();
+}
--- /dev/null
+
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <id> <property_coord_0> <property_coord_1>
+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.
+#
+# <x> <y> <z> <id> <property_coord_0> <property_coord_1> <property_coord_2>
+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.
+#
+# <x> <y> <z> <id> <property_coord_0> <property_coord_1> <property_coord_2>
+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