DataOutFilter::DataOutFilter()
: flags(false, true)
, node_dim(numbers::invalid_unsigned_int)
- , vertices_per_cell(numbers::invalid_unsigned_int)
+ , num_cells(numbers::invalid_unsigned_int)
{}
DataOutFilter::DataOutFilter(const DataOutBase::DataOutFilterFlags &flags)
: flags(flags)
, node_dim(numbers::invalid_unsigned_int)
- , vertices_per_cell(numbers::invalid_unsigned_int)
+ , num_cells(numbers::invalid_unsigned_int)
{}
const unsigned int pt_index)
{
filtered_cells[cell_index] = filtered_points[pt_index];
+
+ // (Re)-initialize counter at any first call to this method.
+ if (cell_index == 0)
+ num_cells = 1;
}
unsigned int
DataOutFilter::n_cells() const
{
- return filtered_cells.size() / vertices_per_cell;
+ return num_cells;
}
const unsigned int d2,
const unsigned int d3)
{
+ ++num_cells;
+
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)
{
const unsigned int start,
const unsigned int n_points)
{
- (void)index;
- (void)start;
- (void)n_points;
+ ++num_cells;
- Assert(false, ExcNotImplemented());
+ const unsigned int base_entry = index * n_points;
+
+ for (unsigned int i = 0; i < n_points; ++i)
+ {
+ internal_add_cell(base_entry + i, start + i);
+ }
}
n_cells = 0;
for (const auto &patch : patches)
{
- // The following formula doesn't work for non-tensor products
+ // The following formula doesn't hold for non-tensor products.
if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim))
{
n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
unsigned int d3)
{
stream << GeometryInfo<dim>::vertices_per_cell << '\t' << start;
+
if (dim >= 1)
stream << '\t' << start + d1;
{
{
const unsigned int n_subdivisions = patch.n_subdivisions;
const unsigned int n = n_subdivisions + 1;
- // Length of loops in all dimensons
+ // Length of loops in all dimensions
const unsigned int n1 = (dim > 0) ? n_subdivisions : 1;
const unsigned int n2 = (dim > 1) ? n_subdivisions : 1;
const unsigned int n3 = (dim > 2) ? n_subdivisions : 1;
}
// max. and min. height of solution
- Assert(patches.size() > 0, ExcInternalError());
+ Assert(patches.size() > 0, ExcNoPatches());
double hmin = patches[0].data(0, 0);
double hmax = patches[0].data(0, 0);
for (const auto &patch : patches)
{
- // special treatment of simplices since they are not subdivided
- if (patch.reference_cell_type != ReferenceCell::get_hypercube(dim))
- {
- n_nodes += patch.data.n_cols();
- n_cells += 1;
- n_points_an_n_cell += patch.data.n_cols() + 1;
- }
- else
+ // The following formulas don't hold for non-tensor products.
+ if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim))
{
n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
(1 + GeometryInfo<dim>::vertices_per_cell);
}
}
+ else
+ {
+ n_nodes += patch.data.n_cols();
+ n_cells += 1;
+ n_points_an_n_cell += patch.data.n_cols() + 1;
+ }
}
// in gmv format the vertex coordinates and the data have an order that is a
<< " <Grid Name=\"CellTime\" GridType=\"Collection\" CollectionType=\"Temporal\">\n";
// Write out all the entries indented
+ const auto &patches = get_patches();
+ Assert(patches.size() > 0, DataOutBase::ExcNoPatches());
+
for (it = entries.begin(); it != entries.end(); ++it)
- xdmf_file << it->get_xdmf_content(3);
+ {
+ xdmf_file << it->get_xdmf_content(3, patches[0].reference_cell_type);
+ }
xdmf_file << " </Grid>\n";
xdmf_file << " </Domain>\n";
template <int dim, int spacedim>
void
DataOutBase::write_hdf5_parallel(
- const std::vector<Patch<dim, spacedim>> & /*patches*/,
- const DataOutBase::DataOutFilter &data_filter,
- const bool write_mesh_file,
- const std::string & mesh_filename,
- const std::string & solution_filename,
- MPI_Comm comm)
+ const std::vector<Patch<dim, spacedim>> &patches,
+ const DataOutBase::DataOutFilter & data_filter,
+ const bool write_mesh_file,
+ const std::string & mesh_filename,
+ const std::string & solution_filename,
+ MPI_Comm comm)
{
AssertThrow(
spacedim >= 2,
#ifndef DEAL_II_WITH_HDF5
// throw an exception, but first make sure the compiler does not warn about
// the now unused function arguments
+ (void)patches;
(void)data_filter;
(void)write_mesh_file;
(void)mesh_filename;
AssertThrow(false, ExcMessage("HDF5 support is disabled."));
#else
# ifndef DEAL_II_WITH_MPI
+ (void)comm;
+# endif
+
// verify that there are indeed patches to be written out. most of the times,
// people just forget to call build_patches when there are no patches, so a
// warning is in order. that said, the assertion is disabled if we support MPI
// since then it can happen that on the coarsest mesh, a processor simply has
// no cells it actually owns, and in that case it is legit if there are no
// patches
- Assert(data_filter.n_nodes() > 0, ExcNoPatches());
- (void)comm;
-# endif
+ Assert(patches.size() > 0, ExcNoPatches());
+
+ const auto &cell_info =
+ ReferenceCell::internal::Info::get_cell(patches[0].reference_cell_type);
hid_t h5_mesh_file_id = -1, h5_solution_file_id, file_plist_id, plist_id;
hid_t node_dataspace, node_dataset, node_file_dataspace,
AssertThrow(node_dataspace >= 0, ExcIO());
cell_ds_dim[0] = global_node_cell_count[1];
- cell_ds_dim[1] = GeometryInfo<dim>::vertices_per_cell;
+ cell_ds_dim[1] = cell_info.n_vertices();
cell_dataspace = H5Screate_simple(2, cell_ds_dim, nullptr);
AssertThrow(cell_dataspace >= 0, ExcIO());
// And repeat for cells
count[0] = local_node_cell_count[1];
- count[1] = GeometryInfo<dim>::vertices_per_cell;
+ count[1] = cell_info.n_vertices();
offset[0] = global_node_cell_offsets[1];
offset[1] = 0;
cell_memory_dataspace = H5Screate_simple(2, count, nullptr);
std::string
XDMFEntry::get_xdmf_content(const unsigned int indent_level) const
+{
+ return get_xdmf_content(indent_level,
+ ReferenceCell::get_hypercube(dimension));
+}
+
+
+
+std::string
+XDMFEntry::get_xdmf_content(
+ const unsigned int indent_level,
+ const ReferenceCell::Type &reference_cell_type) const
{
if (!valid)
return "";
<< "\" NumberOfElements=\"" << num_cells
<< "\" NodesPerElement=\"2\">\n";
else if (dimension == 2)
- ss << indent(indent_level + 1) << "<Topology TopologyType=\""
- << "Quadrilateral"
- << "\" NumberOfElements=\"" << num_cells << "\">\n";
+ {
+ Assert(reference_cell_type == ReferenceCell::Type::Quad ||
+ reference_cell_type == ReferenceCell::Type::Tri,
+ ExcNotImplemented());
+
+ ss << indent(indent_level + 1) << "<Topology TopologyType=\"";
+ if (reference_cell_type == ReferenceCell::Type::Quad)
+ {
+ ss << "Quadrilateral"
+ << "\" NumberOfElements=\"" << num_cells << "\">\n"
+ << indent(indent_level + 2) << "<DataItem Dimensions=\""
+ << num_cells << " " << (1 << dimension);
+ }
+ else // if (reference_cell_type == ReferenceCell::Type::Tri)
+ {
+ ss << "Triangle"
+ << "\" NumberOfElements=\"" << num_cells << "\">\n"
+ << indent(indent_level + 2) << "<DataItem Dimensions=\""
+ << num_cells << " " << 3;
+ }
+ }
else if (dimension == 3)
- ss << indent(indent_level + 1) << "<Topology TopologyType=\""
- << "Hexahedron"
- << "\" NumberOfElements=\"" << num_cells << "\">\n";
+ {
+ Assert(reference_cell_type == ReferenceCell::Type::Hex ||
+ reference_cell_type == ReferenceCell::Type::Tet,
+ ExcNotImplemented());
+
+ ss << indent(indent_level + 1) << "<Topology TopologyType=\"";
+ if (reference_cell_type == ReferenceCell::Type::Hex)
+ {
+ ss << "Hexahedron"
+ << "\" NumberOfElements=\"" << num_cells << "\">\n"
+ << indent(indent_level + 2) << "<DataItem Dimensions=\""
+ << num_cells << " " << (1 << dimension);
+ }
+ else // if (reference_cell_type == ReferenceCell::Type::Tet)
+ {
+ ss << "Tetrahedron"
+ << "\" NumberOfElements=\"" << num_cells << "\">\n"
+ << indent(indent_level + 2) << "<DataItem Dimensions=\""
+ << num_cells << " " << 4;
+ }
+ }
- ss << indent(indent_level + 2) << "<DataItem Dimensions=\"" << num_cells
- << " " << (1 << dimension)
- << "\" NumberType=\"UInt\" Format=\"HDF\">\n";
+ ss << "\" NumberType=\"UInt\" Format=\"HDF\">\n";
ss << indent(indent_level + 3) << h5_mesh_filename << ":/cells\n";
ss << indent(indent_level + 2) << "</DataItem>\n";
ss << indent(indent_level + 1) << "</Topology>\n";
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test DataOut with HDF5 for simplex meshes.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+ RightHandSideFunction(const unsigned int n_components)
+ : Function<dim>(n_components)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const
+ {
+ return p[component % dim] * p[component % dim];
+ }
+};
+
+template <int dim, int spacedim = dim>
+void
+test(const FiniteElement<dim, spacedim> &fe, const unsigned int n_components)
+{
+ Triangulation<dim, spacedim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2);
+
+ DoFHandler<dim> dof_handler(tria);
+
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> solution(dof_handler.n_dofs());
+
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+ VectorTools::interpolate(mapping,
+ dof_handler,
+ RightHandSideFunction<dim>(n_components),
+ solution);
+
+ static unsigned int counter = 0;
+
+ for (unsigned int n_subdivisions = 1; n_subdivisions <= 2; ++n_subdivisions)
+ {
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+
+ data_out.build_patches(mapping, n_subdivisions);
+
+ const std::string output_basename("test." + std::to_string(dim) + "." +
+ std::to_string(counter++));
+
+ DataOutBase::DataOutFilter data_filter(
+ DataOutBase::DataOutFilterFlags(true, true));
+ data_out.write_filtered_data(data_filter);
+ data_out.write_hdf5_parallel(data_filter,
+ output_basename + ".h5",
+ MPI_COMM_SELF);
+
+ std::vector<XDMFEntry> xdmf_entries({data_out.create_xdmf_entry(
+ data_filter, output_basename + ".h5", 0, MPI_COMM_SELF)});
+
+ data_out.write_xdmf_file(xdmf_entries,
+ output_basename + ".xdmf",
+ MPI_COMM_SELF);
+
+ data_out.clear();
+
+ // Sadly hdf5 is binary and we can not use hd5dump because it might
+ // not be in the path. At least we can make sure that both the xdmf and
+ // the h5 file are created.
+ std::ifstream h5((output_basename + ".h5").c_str());
+ AssertThrow(h5.good(), ExcIO());
+
+ std::ifstream xdmf((output_basename + ".xdmf").c_str());
+ AssertThrow(h5.good(), ExcIO());
+
+ deallog << "Files " << output_basename + ".h5"
+ << " and " << output_basename + ".xdmf"
+ << " created succesfully!" << std::endl;
+ }
+}
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ {
+ const unsigned int dim = 2;
+ test<dim>(Simplex::FE_P<dim>(2) /*=degree*/, 1);
+ test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/), dim), dim);
+ test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/),
+ dim,
+ Simplex::FE_P<dim>(1 /*=degree*/),
+ 1),
+ dim + 1);
+ }
+ {
+ const unsigned int dim = 3;
+ test<dim>(Simplex::FE_P<dim>(2) /*=degree*/, 1);
+ test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/), dim), dim);
+ test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/),
+ dim,
+ Simplex::FE_P<dim>(1 /*=degree*/),
+ 1),
+ dim + 1);
+ }
+}
projects / dealii.git / commitdiff