*/
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
* 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
* 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
* 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);
};
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
//----------------------------------------------------------------------//
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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