From e65f90698bf5c14964bf76534f0ac48e0632819f Mon Sep 17 00:00:00 2001 From: Alexander Grayver Date: Thu, 26 Jul 2018 10:13:50 +0200 Subject: [PATCH] Add support for high-order VTK output --- include/deal.II/base/data_out_base.h | 13 +- source/base/data_out_base.cc | 230 +++++++++++++++++++++++++-- 2 files changed, 232 insertions(+), 11 deletions(-) diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h index 9100b73368..4177a1a1f4 100644 --- a/include/deal.II/base/data_out_base.h +++ b/include/deal.II/base/data_out_base.h @@ -205,6 +205,8 @@ class XDMFEntry; * *
  • Tecplot output by Benjamin Shelton Kirk * + *
  • Lagrange VTK output by Alexander Grayver + * * * * @ingroup output @@ -1146,6 +1148,14 @@ namespace DataOutBase */ ZlibCompressionLevel compression_level; + /** + * Flag determining whether to write patches as linear cells + * or as a high-order Lagrange cell. + * + * Default is false. + */ + bool write_lagrange_cells; + /** * Constructor. */ @@ -1153,7 +1163,8 @@ namespace DataOutBase const double time = std::numeric_limits::min(), const unsigned int cycle = std::numeric_limits::min(), const bool print_date_and_time = true, - const ZlibCompressionLevel compression_level = best_compression); + const ZlibCompressionLevel compression_level = best_compression, + const bool write_lagrange_cells = false); }; diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index a2717d421e..0ae8dc3554 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -692,12 +692,20 @@ namespace // DataOutBase in general Wolfgang // supposed that we don't need it in general, but however this choice avoids a // -Warray-bounds check warning - const unsigned int vtk_cell_type[5] = {1, - 3, - 9, - 12, + const unsigned int vtk_cell_type[5] = {1, // VTK_VERTEX + 3, // VTK_LINE + 9, // VTK_QUAD + 12, // VTK_HEXAHEDRON static_cast(-1)}; + // VTK cell ids defined in vtk_cell_type are used for linear cells, + // the ones defined below are used when Lagrange cells are written. + const unsigned int vtk_lagrange_cell_type[5] = {1, // VTK_VERTEX + 68, // VTK_LAGRANGE_CURVE + 70, // VTK_LAGRANGE_QUADRILATERAL + 72, // VTK_LAGRANGE_HEXAHEDRON + static_cast(-1)}; + //----------------------------------------------------------------------// // Auxiliary functions //----------------------------------------------------------------------// @@ -782,6 +790,116 @@ namespace } } + /** + * Given (i,j,k) coordinates within the Lagrange cell, return an + * offset into the local connectivity array + * + * Modified from + * https://github.com/Kitware/VTK/blob/265ca48a79a36538c95622c237da11133608bbe5/Common/DataModel/vtkLagrangeQuadrilateral.cxx#L558 + */ + int vtk_point_index_from_ijk(int i, int j, int, + const std::array &order) + { + bool ibdy = (i == 0 || i == order[0]); + bool jbdy = (j == 0 || j == order[1]); + // How many boundaries do we lie on at once? + int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0); + + if (nbdy == 2) // Vertex DOF + { // ijk is a corner node. Return the proper index (somewhere in [0,7]): + return (i ? (j ? 2 : 1) : (j ? 3 : 0)); + } + + int offset = 4; + if (nbdy == 1) // Edge DOF + { + if (!ibdy) + { // On i axis + return (i - 1) + (j ? order[0] - 1 + order[1] - 1 : 0) + offset; + } + + if (!jbdy) + { // On j axis + return (j - 1) + (i ? order[0] - 1 : 2 * (order[0] - 1) + order[1] - 1) + offset; + } + } + + offset += 2 * (order[0] - 1 + order[1] - 1); + // nbdy == 0: Face DOF + return offset + (i - 1) + (order[0] - 1) * ((j - 1)); + } + + /** + * Same as above, but for hexes. + * + * Modified from + * https://github.com/Kitware/VTK/blob/265ca48a79a36538c95622c237da11133608bbe5/Common/DataModel/vtkLagrangeHexahedron.cxx#L734 + */ + int vtk_point_index_from_ijk(int i, int j, int k, + const std::array &order) + { + bool ibdy = (i == 0 || i == order[0]); + bool jbdy = (j == 0 || j == order[1]); + bool kbdy = (k == 0 || k == order[2]); + // How many boundaries do we lie on at once? + int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0); + + if (nbdy == 3) // Vertex DOF + { // ijk is a corner node. Return the proper index (somewhere in [0,7]): + return (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0); + } + + int offset = 8; + if (nbdy == 2) // Edge DOF + { + if (!ibdy) + { // On i axis + return (i - 1) + (j ? order[0] - 1 + order[1] - 1 : 0) + (k ? 2 * (order[0] - 1 + order[1] - 1) : 0) + offset; + } + if (!jbdy) + { // On j axis + return (j - 1) + (i ? order[0] - 1 : 2 * (order[0] - 1) + order[1] - 1) + (k ? 2 * (order[0] - 1 + order[1] - 1) : 0) + offset; + } + // !kbdy, On k axis + offset += 4 * (order[0] - 1) + 4 * (order[1] - 1); + return (k - 1) + (order[2] - 1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) + offset; + } + + offset += 4 * (order[0] - 1 + order[1] - 1 + order[2] - 1); + if (nbdy == 1) // Face DOF + { + if (ibdy) // On i-normal face + { + return (j - 1) + ((order[1] - 1) * (k - 1)) + (i ? (order[1] - 1) * (order[2] - 1) : 0) + offset; + } + offset += 2 * (order[1] - 1) * (order[2] - 1); + if (jbdy) // On j-normal face + { + return (i - 1) + ((order[0] - 1) * (k - 1)) + (j ? (order[2] - 1) * (order[0] - 1) : 0) + offset; + } + offset += 2 * (order[2] - 1) * (order[0] - 1); + // kbdy, On k-normal face + return (i - 1) + ((order[0] - 1) * (j - 1)) + (k ? (order[0] - 1) * (order[1] - 1) : 0) + offset; + } + + // nbdy == 0: Body DOF + offset += 2 * ((order[1] - 1) * (order[2] - 1) + (order[2] - 1) * (order[0] - 1) + (order[0] - 1) * (order[1] - 1)); + return offset + (i - 1) + (order[0] - 1) * ((j - 1) + (order[1] - 1) * ((k - 1))); + } + + int vtk_point_index_from_ijk(int, int, int, + const std::array &) + { + ExcNotImplemented(); + return 0; + } + + int vtk_point_index_from_ijk(int, int, int, + const std::array &) + { + ExcNotImplemented(); + return 0; + } template @@ -1070,6 +1188,17 @@ namespace const unsigned int x_offset, const unsigned int y_offset, const unsigned int z_offset); + + /** + * Writes Lagrange type cell. The connectivity order of Lagrange + * points is given in the @p connectivity array, which are offset + * by the global index @p start. + */ + template + void + write_lagrange_cell(const unsigned int index, + const unsigned int start, + const std::vector &connectivity); }; @@ -1437,7 +1566,17 @@ namespace stream << '\n'; } - + template + void + VtkStream::write_lagrange_cell(const unsigned int, + const unsigned int start, + const std::vector &connectivity) + { + stream << connectivity.size(); + for(const auto& c: connectivity) + stream << '\t' << start + c; + stream << '\n'; + } VtuStream::VtuStream(std::ostream &out, const DataOutBase::VtkFlags &f) : StreamBase(out, f) @@ -2224,11 +2363,13 @@ namespace DataOutBase VtkFlags::VtkFlags(const double time, const unsigned int cycle, const bool print_date_and_time, - const VtkFlags::ZlibCompressionLevel compression_level) + const VtkFlags::ZlibCompressionLevel compression_level, + const bool write_lagrange_cells) : time(time) , cycle(cycle) , print_date_and_time(print_date_and_time) , compression_level(compression_level) + , write_lagrange_cells(write_lagrange_cells) {} @@ -2375,7 +2516,6 @@ namespace DataOutBase Assert(dim <= 3, ExcNotImplemented()); unsigned int count = 0; unsigned int first_vertex_of_patch = 0; - // Array to hold all the node numbers of a cell. 8 is sufficient for 3D for (typename std::vector>::const_iterator patch = patches.begin(); patch != patches.end(); @@ -2408,6 +2548,57 @@ namespace DataOutBase out.flush_cells(); } + template + void + write_lagrange_cells(const std::vector> &patches, StreamType &out) + { + Assert(dim <= 3 && dim > 1, ExcNotImplemented()); + unsigned int first_vertex_of_patch = 0; + unsigned int count = 0; + // Array to hold all the node numbers of a cell + std::vector connectivity; + // Array to hold cell order in each dimension + std::array cell_order; + + for (typename std::vector>::const_iterator patch = + patches.begin(); + patch != patches.end(); + ++patch) + { + const unsigned int n_subdivisions = patch->n_subdivisions; + const unsigned int n = n_subdivisions + 1; + + cell_order.fill(n_subdivisions); + connectivity.resize(Utilities::fixed_power(n)); + + // Length of loops in all dimensons + const unsigned int n1 = (dim > 0) ? n_subdivisions : 0; + const unsigned int n2 = (dim > 1) ? n_subdivisions : 0; + const unsigned int n3 = (dim > 2) ? n_subdivisions : 0; + // Offsets of outer loops + unsigned int d1 = 1; + unsigned int d2 = n; + unsigned int d3 = n * n; + for (unsigned int i3 = 0; i3 <= n3; ++i3) + for (unsigned int i2 = 0; i2 <= n2; ++i2) + for (unsigned int i1 = 0; i1 <= n1; ++i1) + { + const unsigned int local_offset = i3 * d3 + i2 * d2 + i1 * d1; + const unsigned int connectivity_idx = vtk_point_index_from_ijk(i1, i2, i3, cell_order); + connectivity[connectivity_idx] = local_offset; + } + + out.template write_lagrange_cell(count++, + first_vertex_of_patch, + connectivity); + + // finally update the number of the first vertex of this patch + first_vertex_of_patch += Utilities::fixed_power(n); + } + + out.flush_cells(); + } + template void @@ -5061,6 +5252,17 @@ namespace DataOutBase unsigned int n_nodes; unsigned int n_cells; compute_sizes(patches, n_nodes, n_cells); + + // If a user set to out Lagrange cells, we treat n_subdivions + // as cell order and adjust variables accordingly, otherwise + // each patch is written as a linear cell. + unsigned int n_points_per_cell = GeometryInfo::vertices_per_cell; + if(flags.write_lagrange_cells) + { + n_cells = patches.size(); + n_points_per_cell = n_nodes / n_cells; + } + // in gmv format the vertex coordinates and the data have an order that is a // bit unpleasant (first all x coordinates, then all y coordinate, ...; // first all data of variable 1, then variable 2, etc), so we have to copy @@ -5091,14 +5293,22 @@ namespace DataOutBase ///////////////////////////////// // now for the cells out << "CELLS " << n_cells << ' ' - << n_cells * (GeometryInfo::vertices_per_cell + 1) << '\n'; - write_cells(patches, vtk_out); + << n_cells * (n_points_per_cell + 1) << '\n'; + if(flags.write_lagrange_cells) + write_lagrange_cells(patches, vtk_out); + else + write_cells(patches, vtk_out); out << '\n'; // next output the types of the cells. since all cells are the same, this is // simple out << "CELL_TYPES " << n_cells << '\n'; + + // need to distinguish between linear and Lagrange cells + const unsigned int vtk_cell_id = flags.write_lagrange_cells + ? vtk_lagrange_cell_type[dim] + : vtk_cell_type[dim]; for (unsigned int i = 0; i < n_cells; ++i) - out << ' ' << vtk_cell_type[dim]; + out << ' ' << vtk_cell_id; out << '\n'; /////////////////////////////////////// // data output. -- 2.39.5