Ordering of vertices for wedges and pyramids

[bangerth]

@peterrum -- primarily a question for you:

• In ReferenceCell::vertex(), we have this code for the ordering of vertices of wedges:

  else if ((dim == 3) && (*this == ReferenceCells::Wedge))
    {
      static const Point<dim> vertices[6] = {Point<dim>{1.0, 0.0, 0.0},
                                             Point<dim>{0.0, 1.0, 0.0},
                                             Point<dim>{0.0, 0.0, 0.0},
                                             Point<dim>{1.0, 0.0, 1.0},
                                             Point<dim>{0.0, 1.0, 1.0},
                                             Point<dim>{0.0, 0.0, 1.0}};
      return vertices[v];

It is a little bit surprising that the origin and its translate in the z-direction are the third and sixth vertices, rather than the first and fourth. It would have been nice if the wedge were simply the tensor product of the vertices of a triangle with the z-direction, in the same way as the cube is simply the tensor product of the square with the z-direction. Is there a reason for doing it that way? What would happen if we still changed that?

• For pyramids, we have this in ReferenceCell::vertex():

      static const Point<dim> vertices[5] = {Point<dim>{-1.0, -1.0, 0.0},
                                             Point<dim>{+1.0, -1.0, 0.0},
                                             Point<dim>{-1.0, +1.0, 0.0},
                                             Point<dim>{+1.0, +1.0, 0.0},
                                             Point<dim>{+0.0, +0.0, 1.0}};

which follows the tensor-product principle for the construction of the base square, but then we have the following in data_out_dof_data.templates.h:

              std::vector<Point<dim>> points;
              points.emplace_back(-1.0, -1.0, 0.0);
              points.emplace_back(+1.0, -1.0, 0.0);
              points.emplace_back(+1.0, +1.0, 0.0);
              points.emplace_back(-1.0, +1.0, 0.0);
              points.emplace_back(+0.0, +0.0, 1.0);

which uses a different order (counterclockwise) for the points on the pyramid's base. Do you remember what the reason for the discrepancy is? Can that still be changed? If the difference is between our own numbering and what VTK wants, my preference would be to use our own numbering and make sure that we translate where necessary, as we do in other places where specific output formats want specific orderings.

[peterrum]

Regarding, the wedges: they should have been tensor products of triangles and lines. See also:

dealii/source/base/polynomials_wedge.cc

Lines 60 to 61 in 6e182af

	static const constexpr ndarray<unsigned int, 6, 2> wedge_table_1{
	{{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};

dealii/source/fe/fe_wedge_p.cc

Lines 106 to 114 in 6e182af

	if (degree == 1)
	{
	this->unit_support_points.emplace_back(0.0, 0.0, 0.0);
	this->unit_support_points.emplace_back(1.0, 0.0, 0.0);
	this->unit_support_points.emplace_back(0.0, 1.0, 0.0);
	this->unit_support_points.emplace_back(0.0, 0.0, 1.0);
	this->unit_support_points.emplace_back(1.0, 0.0, 1.0);
	this->unit_support_points.emplace_back(0.0, 1.0, 1.0);
	}

I have created PR #12845 to make the code consistent. No test has failed after this change. I guess the actual positions of the points are not tested but implicitly testes via L2-norms ...

Regarding pyramids: DataOut contains the vertices in order how VTK needs it. See https://www.paraview.org/Wiki/images/5/51/VTK-File-Formats.pdf. My guess is that the reordering should happen in the base class:

dealii/source/base/data_out_base.cc

Lines 774 to 802 in 6dce9bf

	// For a given patch, compute the nodes for arbitrary (non-hypercube) cells.
	// If the points are saved in the patch.data member, return the saved point
	// instead.
	template <int dim, int spacedim>
	inline Point<spacedim>
	compute_arbitrary_node(const DataOutBase::Patch<dim, spacedim> &patch,
	const unsigned int point_no)
	{
	Point<spacedim> node;
	if (patch.points_are_available)
	{
	for (unsigned int d = 0; d < spacedim; ++d)
	node[d] = patch.data(patch.data.size(0) - spacedim + d, point_no);
	return node;
	}
	else
	{
	AssertDimension(patch.n_subdivisions, 1);
	Assert(
	patch.reference_cell != ReferenceCells::Pyramid,
	ExcMessage(
	"Pyramids need different ordering of the vertices, which is not implemented yet here."));
	node = patch.vertices[point_no];
	}
	return node;
	}

[drwells]

This is similar to what I wrote in #11124.

[bangerth]

I think we fixed this!