From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Wed, 27 May 2020 18:50:35 +0000 (-0700)
Subject: Add scalar properties
X-Git-Tag: v9.3.0-rc1~1523^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=50125914480af7083c5b5c17ac03264bca3b439b;p=dealii.git

Add scalar properties
---

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 <deal.II/base/data_out_base.h>
 
+#include <deal.II/numerics/data_component_interpretation.h>
+
 #include <string>
 #include <vector>
 
@@ -66,7 +68,11 @@ namespace Particles
      * @author Bruno Blais, Luca Heltai 2019
      */
     void
-    build_patches(const Particles::ParticleHandler<dim, spacedim> &particles);
+    build_patches(const Particles::ParticleHandler<dim, spacedim> &particles,
+                  const std::vector<std::string> &                 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 <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();
@@ -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<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];
+          }
       }
   }
 
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 <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();
+}
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.
+#
+# <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