]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Restructure DataOutFilter
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Tue, 17 Oct 2017 21:27:40 +0000 (15:27 -0600)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Thu, 26 Oct 2017 16:35:45 +0000 (10:35 -0600)
include/deal.II/base/data_out_base.h
source/base/data_out_base.cc
tests/numerics/data_out_filter_01.cc [new file with mode: 0644]
tests/numerics/data_out_filter_01.with_hdf5=true.output [new file with mode: 0644]

index 125846ddbe8a07c4db7ee2e4bbfd25db90eef983..857c9e1c6e4243a79158473d7771ea8dc0e2e45e 100644 (file)
@@ -1274,76 +1274,17 @@ namespace DataOutBase
    */
   class DataOutFilter
   {
-  private:
+  public:
     /**
-     * Empty class to provide comparison function for Map3DPoint.
+     * Default constructor.
      */
-    struct Point3Comp
-    {
-      bool operator() (const Point<3> &lhs, const Point<3> &rhs) const
-      {
-        return (lhs(0) < rhs(0) || (!(rhs(0) < lhs(0)) && (lhs(1) < rhs(1) || (!(rhs(1) < lhs(1)) && lhs(2) < rhs(2)))));
-      }
-    };
-
-    typedef std::multimap<Point<3>, unsigned int, Point3Comp> Map3DPoint;
-
-    /// Flags used to specify filtering behavior
-    DataOutBase::DataOutFilterFlags   flags;
+    DataOutFilter();
 
     /**
-     * The number of space dimensions in which the vertices represented
-     * by the current object live. This corresponds to the usual
-     * @p dim argument, but since this class is not templated on the
-     * dimension, we need to store it here.
+     * Destructor with a given set of flags. See DataOutFilterFlags for
+     * possible flags.
      */
-    unsigned int      node_dim;
-
-    /**
-     * The number of vertices per cell. Equal to
-     * GeometryInfo<node_dim>::vertices_per_cell. We need to store
-     * it as a run-time variable here because the dimension
-     * node_dim is also a run-time variable.
-     */
-    unsigned int      vertices_per_cell;
-
-    /// Map of points to an internal index
-    Map3DPoint        existing_points;
-
-    /// Map of actual point index to internal point index
-    std::map<unsigned int, unsigned int>  filtered_points;
-
-    /// Map of cells to the filtered points
-    std::map<unsigned int, unsigned int>  filtered_cells;
-
-    /// Data set names
-    std::vector<std::string>    data_set_names;
-
-    /// Data set dimensions
-    std::vector<unsigned int>   data_set_dims;
-
-    /// Data set data
-    std::vector<std::vector<double> > data_sets;
-
-    /**
-     * Record a cell vertex index based on the internal reordering.
-     */
-    void internal_add_cell(const unsigned int &cell_index, const unsigned int &pt_index);
-
-  public:
-    DataOutFilter()
-      :
-      flags(false, true),
-      node_dim (numbers::invalid_unsigned_int),
-      vertices_per_cell (numbers::invalid_unsigned_int)
-    {}
-
-    DataOutFilter(const DataOutBase::DataOutFilterFlags &flags)
-      :
-      flags(flags),
-      node_dim (numbers::invalid_unsigned_int),
-      vertices_per_cell (numbers::invalid_unsigned_int)
-    {}
+    DataOutFilter(const DataOutBase::DataOutFilterFlags &flags);
 
     /**
      * Write a point with the specified index into the filtered data set. If
@@ -1351,13 +1292,18 @@ namespace DataOutBase
      * provided index will internally refer to another recorded point.
      */
     template <int dim>
-    void write_point(const unsigned int &index, const Point<dim> &p);
+    void write_point(const unsigned int index,
+                     const Point<dim> &p);
 
     /**
      * Record a deal.II cell in the internal reordered format.
      */
     template <int dim>
-    void write_cell(unsigned int index, unsigned int start, unsigned int d1, unsigned int d2, unsigned int d3);
+    void write_cell(const unsigned int index,
+                    const unsigned int start,
+                    const unsigned int d1,
+                    const unsigned int d2,
+                    const unsigned int d3);
 
     /**
      * Filter and record a data set. If there are multiple values at a given
@@ -1365,7 +1311,10 @@ namespace DataOutBase
      * chosen as the recorded value. In the future this can be expanded to
      * average/min/max multiple values at a given vertex.
      */
-    void write_data_set(const std::string &name, const unsigned int &dimension, const unsigned int &set_num, const Table<2,double> &data_vectors);
+    void write_data_set(const std::string &name,
+                        const unsigned int dimension,
+                        const unsigned int set_num,
+                        const Table<2,double> &data_vectors);
 
     /**
      * Resize and fill a vector with all the filtered node vertex points, for
@@ -1377,69 +1326,138 @@ namespace DataOutBase
      * Resize and fill a vector with all the filtered cell vertex indices, for
      * output to a file.
      */
-    void fill_cell_data(const unsigned int &local_node_offset, std::vector<unsigned int> &cell_data) const;
+    void fill_cell_data(const unsigned int local_node_offset,
+                        std::vector<unsigned int> &cell_data) const;
 
     /**
      * Get the name of the data set indicated by the set number.
      */
-    std::string get_data_set_name(const unsigned int &set_num) const
-    {
-      return data_set_names.at(set_num);
-    };
+    std::string get_data_set_name(const unsigned int set_num) const;
 
     /**
      * Get the dimensionality of the data set indicated by the set number.
      */
-    unsigned int get_data_set_dim(const unsigned int &set_num) const
-    {
-      return data_set_dims.at(set_num);
-    };
+    unsigned int get_data_set_dim(const unsigned int set_num) const;
 
     /**
      * Get the raw double valued data of the data set indicated by the set
      * number.
      */
-    const double *get_data_set(const unsigned int &set_num) const
-    {
-      return &data_sets[set_num][0];
-    };
+    const double *get_data_set(const unsigned int set_num) const;
 
     /**
      * Return the number of nodes in this DataOutFilter. This may be smaller
      * than the original number of nodes if filtering is enabled.
      */
-    unsigned int n_nodes() const
-    {
-      return existing_points.size();
-    };
+    unsigned int n_nodes() const;
 
     /**
      * Return the number of filtered cells in this DataOutFilter. Cells are
      * not filtered so this will be the original number of cells.
      */
-    unsigned int n_cells() const
-    {
-      return filtered_cells.size()/vertices_per_cell;
-    };
+    unsigned int n_cells() const;
 
     /**
      * Return the number of filtered data sets in this DataOutFilter. Data
      * sets are not filtered so this will be the original number of data sets.
      */
-    unsigned int n_data_sets() const
-    {
-      return data_set_names.size();
-    };
+    unsigned int n_data_sets() const;
 
     /**
      * Empty functions to do base class inheritance.
      */
-    void flush_points () {};
+    void flush_points ();
 
     /**
      * Empty functions to do base class inheritance.
      */
-    void flush_cells () {};
+    void flush_cells ();
+
+
+  private:
+    /**
+     * Empty class to provide comparison function for Map3DPoint.
+     */
+    struct Point3Comp
+    {
+      bool operator() (const Point<3> &one,
+                       const Point<3> &two) const
+      {
+        /*
+         * The return statement below is an optimized version of the following code:
+         *
+         * for (unsigned int d=0; d<3; ++d)
+         * {
+         *   if (one(d) < two(d))
+         *     return true;
+         *   else if (one(d) > two(d))
+         *     return false;
+         * }
+         * return false;
+         */
+
+        return (one(0) < two(0) || (!(two(0) < one(0)) && (one(1) < two(1) || (!(two(1) < one(1)) && one(2) < two(2)))));
+      }
+    };
+
+    typedef std::multimap<Point<3>, unsigned int, Point3Comp> Map3DPoint;
+
+    /**
+     * Flags used to specify filtering behavior.
+     */
+    DataOutBase::DataOutFilterFlags   flags;
+
+    /**
+     * The number of space dimensions in which the vertices represented
+     * by the current object live. This corresponds to the usual
+     * @p dim argument, but since this class is not templated on the
+     * dimension, we need to store it here.
+     */
+    unsigned int      node_dim;
+
+    /**
+     * The number of vertices per cell. Equal to
+     * GeometryInfo<node_dim>::vertices_per_cell. We need to store
+     * it as a run-time variable here because the dimension
+     * node_dim is also a run-time variable.
+     */
+    unsigned int      vertices_per_cell;
+
+    /**
+     * Map of points to an internal index.
+     */
+    Map3DPoint        existing_points;
+
+    /**
+     * Map of actual point index to internal point index.
+     */
+    std::map<unsigned int, unsigned int>  filtered_points;
+
+    /**
+     * Map of cells to the filtered points.
+     */
+    std::map<unsigned int, unsigned int>  filtered_cells;
+
+    /**
+     * Data set names.
+     */
+    std::vector<std::string>    data_set_names;
+
+    /**
+     * Data set dimensions.
+     */
+    std::vector<unsigned int>   data_set_dims;
+
+    /**
+     * Data set data.
+     */
+    std::vector<std::vector<double> > data_sets;
+
+    /**
+     * Record a cell vertex index based on the internal reordering.
+     */
+    void internal_add_cell(const unsigned int cell_index,
+                           const unsigned int pt_index);
   };
 
 
index 74af2846b930c9c7d6f382d89289137e56bf1695..4200ca4441edf74665cca4aad823e192a082a362 100644 (file)
@@ -453,124 +453,221 @@ namespace DataOutBase
                 ExcInternalError());
     }
   }
-}
 
-//----------------------------------------------------------------------//
-// DataOutFilter class member functions
-//----------------------------------------------------------------------//
 
-template <int dim>
-void DataOutBase::DataOutFilter::write_point(const unsigned int &index, const Point<dim> &p)
-{
-  Map3DPoint::const_iterator  it;
-  unsigned int      internal_ind;
-  Point<3>      int_pt;
 
-  for (int d=0; d<3; ++d) int_pt(d) = (d < dim ? p(d) : 0);
-  node_dim = dim;
-  it = existing_points.find(int_pt);
+  DataOutFilter::DataOutFilter()
+    :
+    flags(false, true),
+    node_dim (numbers::invalid_unsigned_int),
+    vertices_per_cell (numbers::invalid_unsigned_int)
+  {}
 
-  // If the point isn't in the set, or we're not filtering duplicate points, add it
-  if (it == existing_points.end() || !flags.filter_duplicate_vertices)
-    {
-      internal_ind = existing_points.size();
-      existing_points.insert(std::make_pair(int_pt, internal_ind));
-    }
-  else
-    {
-      internal_ind = it->second;
-    }
-  // Now add the index to the list of filtered points
-  filtered_points[index] = internal_ind;
-}
 
-void DataOutBase::DataOutFilter::internal_add_cell(const unsigned int &cell_index, const unsigned int &pt_index)
-{
-  filtered_cells[cell_index] = filtered_points[pt_index];
-}
 
-void DataOutBase::DataOutFilter::fill_node_data(std::vector<double> &node_data) const
-{
-  Map3DPoint::const_iterator  it;
+  DataOutFilter::DataOutFilter(const DataOutBase::DataOutFilterFlags &flags)
+    :
+    flags(flags),
+    node_dim (numbers::invalid_unsigned_int),
+    vertices_per_cell (numbers::invalid_unsigned_int)
+  {}
 
-  node_data.resize(existing_points.size()*node_dim);
 
-  for (it=existing_points.begin(); it!=existing_points.end(); ++it)
-    {
-      for (unsigned int d=0; d<node_dim; ++d)
-        node_data[node_dim*it->second+d] = it->first(d);
-    }
-}
 
-void DataOutBase::DataOutFilter::fill_cell_data(const unsigned int &local_node_offset, std::vector<unsigned int> &cell_data) const
-{
-  std::map<unsigned int, unsigned int>::const_iterator  it;
+  template <int dim>
+  void
+  DataOutFilter::write_point(const unsigned int index,
+                             const Point<dim> &p)
+  {
+    node_dim = dim;
 
-  cell_data.resize(filtered_cells.size());
+    Point<3> int_pt;
+    for (unsigned int d=0; d<dim; ++d)
+      int_pt(d) = p(d);
 
-  for (it=filtered_cells.begin(); it!=filtered_cells.end(); ++it)
-    {
-      cell_data[it->first] = it->second+local_node_offset;
-    }
-}
+    const Map3DPoint::const_iterator it = existing_points.find(int_pt);
+    unsigned int internal_ind;
 
-template <int dim>
-void
-DataOutBase::DataOutFilter::write_cell(
-  unsigned int index,
-  unsigned int start,
-  unsigned int d1,
-  unsigned int d2,
-  unsigned int d3)
-{
-  unsigned int base_entry = index * GeometryInfo<dim>::vertices_per_cell;
-  vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
-  internal_add_cell(base_entry+0, start);
-  if (dim>=1)
-    {
-      internal_add_cell(base_entry+1, start+d1);
-      if (dim>=2)
-        {
-          internal_add_cell(base_entry+2, start+d2+d1);
-          internal_add_cell(base_entry+3, start+d2);
-          if (dim>=3)
-            {
-              internal_add_cell(base_entry+4, start+d3);
-              internal_add_cell(base_entry+5, start+d3+d1);
-              internal_add_cell(base_entry+6, start+d3+d2+d1);
-              internal_add_cell(base_entry+7, start+d3+d2);
-            }
-        }
-    }
-}
+    // If the point isn't in the set, or we're not filtering duplicate points, add it
+    if (it == existing_points.end() || !flags.filter_duplicate_vertices)
+      {
+        internal_ind = existing_points.size();
+        existing_points.insert(std::make_pair(int_pt, internal_ind));
+      }
+    else
+      {
+        internal_ind = it->second;
+      }
+    // Now add the index to the list of filtered points
+    filtered_points[index] = internal_ind;
+  }
 
-void DataOutBase::DataOutFilter::write_data_set(const std::string &name, const unsigned int &dimension, const unsigned int &set_num, const Table<2,double> &data_vectors)
-{
-  unsigned int    num_verts = existing_points.size();
-  unsigned int    i, r, d, new_dim;
 
-  // HDF5/XDMF output only supports 1D or 3D output, so force rearrangement if needed
-  if (flags.xdmf_hdf5_output && dimension != 1) new_dim = 3;
-  else new_dim = dimension;
 
-  // Record the data set name, dimension, and allocate space for it
-  data_set_names.push_back(name);
-  data_set_dims.push_back(new_dim);
-  data_sets.emplace_back(new_dim*num_verts);
+  void
+  DataOutFilter::internal_add_cell(const unsigned int cell_index,
+                                   const unsigned int pt_index)
+  {
+    filtered_cells[cell_index] = filtered_points[pt_index];
+  }
 
-  // TODO: averaging, min/max, etc for merged vertices
-  for (i=0; i<filtered_points.size(); ++i)
-    {
-      for (d=0; d<new_dim; ++d)
-        {
-          r = filtered_points[i];
-          if (d < dimension) data_sets.back()[r*new_dim+d] = data_vectors(set_num+d, i);
-          else data_sets.back()[r*new_dim+d] = 0;
-        }
-    }
+
+
+  void
+  DataOutFilter::fill_node_data(std::vector<double> &node_data) const
+  {
+    node_data.resize(existing_points.size()*node_dim);
+
+    for (Map3DPoint::const_iterator it=existing_points.begin(); it!=existing_points.end(); ++it)
+      {
+        for (unsigned int d=0; d<node_dim; ++d)
+          node_data[node_dim*it->second+d] = it->first(d);
+      }
+  }
+
+
+
+  void
+  DataOutFilter::fill_cell_data(const unsigned int local_node_offset,
+                                std::vector<unsigned int> &cell_data) const
+  {
+    cell_data.resize(filtered_cells.size());
+
+    for (std::map<unsigned int, unsigned int>::const_iterator it=filtered_cells.begin(); it!=filtered_cells.end(); ++it)
+      {
+        cell_data[it->first] = it->second+local_node_offset;
+      }
+  }
+
+
+
+  std::string
+  DataOutFilter::get_data_set_name(const unsigned int set_num) const
+  {
+    return data_set_names.at(set_num);
+  }
+
+
+
+  unsigned int
+  DataOutFilter::get_data_set_dim(const unsigned int set_num) const
+  {
+    return data_set_dims.at(set_num);
+  }
+
+
+
+  const double *
+  DataOutFilter::get_data_set(const unsigned int set_num) const
+  {
+    return &data_sets[set_num][0];
+  }
+
+
+
+  unsigned int
+  DataOutFilter::n_nodes() const
+  {
+    return existing_points.size();
+  }
+
+
+
+  unsigned int
+  DataOutFilter::n_cells() const
+  {
+    return filtered_cells.size()/vertices_per_cell;
+  }
+
+
+
+  unsigned int
+  DataOutFilter::n_data_sets() const
+  {
+    return data_set_names.size();
+  }
+
+
+
+  void
+  DataOutFilter::flush_points ()
+  {}
+
+
+
+  void
+  DataOutFilter::flush_cells ()
+  {}
+
+
+
+  template <int dim>
+  void
+  DataOutFilter::write_cell(const unsigned int index,
+                            const unsigned int start,
+                            const unsigned int d1,
+                            const unsigned int d2,
+                            const unsigned int d3)
+  {
+    const unsigned int base_entry = index * GeometryInfo<dim>::vertices_per_cell;
+    vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
+    internal_add_cell(base_entry+0, start);
+    if (dim>=1)
+      {
+        internal_add_cell(base_entry+1, start+d1);
+        if (dim>=2)
+          {
+            internal_add_cell(base_entry+2, start+d2+d1);
+            internal_add_cell(base_entry+3, start+d2);
+            if (dim>=3)
+              {
+                internal_add_cell(base_entry+4, start+d3);
+                internal_add_cell(base_entry+5, start+d3+d1);
+                internal_add_cell(base_entry+6, start+d3+d2+d1);
+                internal_add_cell(base_entry+7, start+d3+d2);
+              }
+          }
+      }
+  }
+
+
+
+  void DataOutFilter::write_data_set(const std::string &name,
+                                     const unsigned int dimension,
+                                     const unsigned int set_num,
+                                     const Table<2,double> &data_vectors)
+  {
+    unsigned int new_dim;
+
+    // HDF5/XDMF output only supports 1D or 3D output, so force rearrangement if needed
+    if (flags.xdmf_hdf5_output && dimension != 1)
+      new_dim = 3;
+    else
+      new_dim = dimension;
+
+    // Record the data set name, dimension, and allocate space for it
+    data_set_names.push_back(name);
+    data_set_dims.push_back(new_dim);
+    data_sets.emplace_back(new_dim * existing_points.size());
+
+    // TODO: averaging, min/max, etc for merged vertices
+    for (unsigned int i=0; i<filtered_points.size(); ++i)
+      {
+        const unsigned int r = filtered_points[i];
+
+        for (unsigned int d=0; d<new_dim; ++d)
+          {
+            if (d < dimension)
+              data_sets.back()[r*new_dim+d] = data_vectors(set_num+d, i);
+            else
+              data_sets.back()[r*new_dim+d] = 0;
+          }
+      }
+  }
 }
 
 
+
 //----------------------------------------------------------------------//
 //Auxiliary data
 //----------------------------------------------------------------------//
diff --git a/tests/numerics/data_out_filter_01.cc b/tests/numerics/data_out_filter_01.cc
new file mode 100644 (file)
index 0000000..f6e90b2
--- /dev/null
@@ -0,0 +1,94 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// This test shows the ability of DataOutFilter to merge vertex duplicates taking
+// into account floating point inprecision.
+// The test creates a distorted triangulation with vertices very close to each other
+// that are subsequently merged by the DataOutFilter class.
+// Note that while the test case will likely never happen, there are cases
+// in which vertices during mesh output have slightly different positions
+// that would not be merged if the points would simply be compared for equality.
+
+#include "../tests.h"
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/numerics/data_out.h>
+
+
+
+
+template <int dim>
+void
+test ()
+{
+  Triangulation<dim> tria(Triangulation<dim>::MeshSmoothing::none,true);
+  GridGenerator::hyper_cube(tria, 1., 1. + 1000.0 * std::numeric_limits<double>::epsilon());
+
+  FE_Q<dim> fe1(1);
+
+  DoFHandler<dim> dof1(tria);
+  dof1.distribute_dofs(fe1);
+
+  Vector<double> v1(dof1.n_dofs());
+  for (unsigned int i=0; i<v1.size(); ++i) v1(i) = i;
+
+  DataOut<dim> data_out;
+  data_out.add_data_vector (dof1, v1, "linear");
+  data_out.build_patches ();
+
+  DataOutBase::DataOutFilter data_filter
+  (DataOutBase::DataOutFilterFlags (true, false));
+
+  data_out.write_filtered_data (data_filter);
+
+  deallog << "Number of filtered nodes: " << data_filter.n_nodes() << std::endl;
+
+  Triangulation<dim> tria2(Triangulation<dim>::MeshSmoothing::none,true);
+  GridGenerator::hyper_cube(tria2, 1., 1 + 10.0 * std::numeric_limits<double>::epsilon());
+
+  FE_Q<dim> fe2(1);
+
+  DoFHandler<dim> dof2(tria2);
+  dof2.distribute_dofs(fe2);
+
+  Vector<double> v2(dof2.n_dofs());
+  for (unsigned int i=0; i<v2.size(); ++i) v2(i) = i;
+
+  DataOut<dim> data_out2;
+  data_out2.add_data_vector (dof2, v2, "linear");
+  data_out2.build_patches ();
+
+  DataOutBase::DataOutFilter data_filter2
+  (DataOutBase::DataOutFilterFlags (true, false));
+
+  data_out2.write_filtered_data (data_filter2);
+
+  deallog << "Number of filtered nodes: " << data_filter2.n_nodes() << std::endl;
+
+  deallog << "ok" << std::endl;
+}
+
+
+int main(int argc, char *argv[])
+{
+  initlog();
+  test<2>();
+  test<3>();
+
+  return 0;
+}
diff --git a/tests/numerics/data_out_filter_01.with_hdf5=true.output b/tests/numerics/data_out_filter_01.with_hdf5=true.output
new file mode 100644 (file)
index 0000000..69858ee
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::Number of filtered nodes: 4
+DEAL::Number of filtered nodes: 1
+DEAL::ok
+DEAL::Number of filtered nodes: 8
+DEAL::Number of filtered nodes: 1
+DEAL::ok

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.