From: Wolfgang Bangerth Date: Thu, 28 Oct 2021 03:15:50 +0000 (-0600) Subject: Add the ability to output triangles with gnuplot. X-Git-Tag: v9.4.0-rc1~874^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bd009472d4e4caf5e3612a3e77b24a1c65117ca9;p=dealii.git Add the ability to output triangles with gnuplot. --- diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index 72aff84ea5..cadbe6b389 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -710,9 +710,9 @@ namespace const unsigned int zstep, const unsigned int n_subdivisions) { - Point node; if (patch.points_are_available) { + Assert(n_subdivisions == patch.n_subdivisions, ExcNotImplemented()); unsigned int point_no = 0; switch (dim) { @@ -735,20 +735,22 @@ namespace default: Assert(false, ExcNotImplemented()); } + Point node; for (unsigned int d = 0; d < spacedim; ++d) node[d] = patch.data(patch.data.size(0) - spacedim + d, point_no); + return node; } else { if (dim == 0) - node = patch.vertices[0]; + return patch.vertices[0]; else { // perform a dim-linear interpolation const double stepsize = 1. / n_subdivisions, xfrac = xstep * stepsize; - node = + Point node = (patch.vertices[1] * xfrac) + (patch.vertices[0] * (1 - xfrac)); if (dim > 1) { @@ -770,9 +772,9 @@ namespace zfrac; } } + return node; } } - return node; } // For a given patch, compute the nodes for arbitrary (non-hypercube) cells. @@ -3685,10 +3687,6 @@ namespace DataOutBase (n_data_sets + spacedim) : n_data_sets, patch.data.n_rows())); - Assert(patch.data.n_cols() == - Utilities::fixed_power(n_points_per_direction), - ExcInvalidDatasetSize(patch.data.n_cols(), n_subdivisions + 1)); - auto output_point_data = [&out, &patch, n_data_sets](const unsigned int point_index) mutable { @@ -3702,6 +3700,10 @@ namespace DataOutBase { Assert(patch.reference_cell == ReferenceCells::Vertex, ExcInternalError()); + Assert(patch.data.n_cols() == 1, + ExcInvalidDatasetSize(patch.data.n_cols(), + n_subdivisions + 1)); + // compute coordinates for this patch point out << compute_hypercube_node(patch, 0, 0, 0, n_subdivisions) @@ -3716,6 +3718,10 @@ namespace DataOutBase { Assert(patch.reference_cell == ReferenceCells::Line, ExcInternalError()); + Assert(patch.data.n_cols() == + Utilities::fixed_power(n_points_per_direction), + ExcInvalidDatasetSize(patch.data.n_cols(), + n_subdivisions + 1)); for (unsigned int i1 = 0; i1 < n_points_per_direction; ++i1) { @@ -3735,24 +3741,72 @@ namespace DataOutBase case 2: { - Assert(patch.reference_cell == ReferenceCells::Quadrilateral, - ExcNotImplemented()); - - for (unsigned int i2 = 0; i2 < n_points_per_direction; ++i2) + if (patch.reference_cell == ReferenceCells::Quadrilateral) { - for (unsigned int i1 = 0; i1 < n_points_per_direction; ++i1) + Assert(patch.data.n_cols() == Utilities::fixed_power( + n_points_per_direction), + ExcInvalidDatasetSize(patch.data.n_cols(), + n_subdivisions + 1)); + + for (unsigned int i2 = 0; i2 < n_points_per_direction; ++i2) { - // compute coordinates for this patch point - out << compute_hypercube_node( - patch, i1, i2, 0, n_subdivisions) - << ' '; + for (unsigned int i1 = 0; i1 < n_points_per_direction; + ++i1) + { + // compute coordinates for this patch point + out << compute_hypercube_node( + patch, i1, i2, 0, n_subdivisions) + << ' '; - output_point_data(i1 + i2 * n_points_per_direction); + output_point_data(i1 + i2 * n_points_per_direction); + out << '\n'; + } + // end of row in patch out << '\n'; } - // end of row in patch + } + else if (patch.reference_cell == ReferenceCells::Triangle) + { + Assert(n_subdivisions == 1, ExcNotImplemented()); + + Assert(patch.data.n_cols() == 3, ExcInternalError()); + + // Gnuplot can only plot surfaces if each facet of the + // surface is a bilinear patch, or a subdivided bilinear + // patch with equally many points along each row of the + // subdivision. This is what the code above for + // quadrilaterals does. We emulate this by repeating the + // third point of a triangle twice so that there are two + // points for that row as well -- i.e., we write a 2x2 + // bilinear patch where two of the points are collapsed onto + // one vertex. + // + // This also matches the example here: + // https://stackoverflow.com/questions/42784369/drawing-triangular-mesh-using-gnuplot + out << compute_arbitrary_node(patch, 0) << ' '; + output_point_data(0); + out << '\n'; + + out << compute_arbitrary_node(patch, 1) << ' '; + output_point_data(1); + out << '\n'; + out << '\n'; // end of one row of points + + out << compute_arbitrary_node(patch, 2) << ' '; + output_point_data(2); out << '\n'; + + out << compute_arbitrary_node(patch, 2) << ' '; + output_point_data(2); + out << '\n'; + out << '\n'; // end of the second row of points + out << '\n'; // end of the entire patch } + else + // There aren't any other reference cells in 2d than the + // quadrilateral and the triangle. So whatever we got here + // can't be any good + Assert(false, ExcInternalError()); // end of patch out << '\n'; @@ -3763,6 +3817,10 @@ namespace DataOutBase { Assert(patch.reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); + Assert(patch.data.n_cols() == + Utilities::fixed_power(n_points_per_direction), + ExcInvalidDatasetSize(patch.data.n_cols(), + n_subdivisions + 1)); // for all grid points: draw lines into all positive coordinate // directions if there is another grid point there