Plot 2D surfaces in 3D space using P4estMesh{NDIMS, NDIMS_AMBIENT}

[tristanmontoya]

This PR adapts #82 to use the NDIMS_AMBIENT type parameter, which I added to P4estMesh in trixi-framework/Trixi.jl#2068. This was released in Trixi.jl v0.8.9, so I have removed compat for versions of Trixi before 0.9 (@sloede or @ranocha please let me know if this is not necessary, or if you'd like to maintain compatibility with earlier versions of Trixi.jl). The main change is in the function below:

Trixi2Vtk.jl/src/vtktools.jl

Lines 234 to 246 in 6fffc3e

	# Version of calc_node_coordinates for a P4estMesh representing a manifold of dimension
	# NDIMS embedded within an ambient space of dimension NDIMS_AMBIENT. This provides support
	# for standard 2D and 3D meshes (i.e. NDIMS = NDIMS_AMBIENT) as well as 2D surfaces in 3D
	# space (i.e. NDIMS = 2 and NDIMS_AMBIENT = 3).
	function calc_node_coordinates(mesh::P4estMesh{NDIMS, NDIMS_AMBIENT}, nodes,
	n_visnodes) where {NDIMS, NDIMS_AMBIENT}
	node_coordinates = Array{Float64, NDIMS+2}(undef, NDIMS_AMBIENT,
	ntuple(_ -> n_visnodes, NDIMS)...,
	Trixi.ncells(mesh))
	return Trixi.calc_node_coordinates!(node_coordinates, mesh, nodes)
	end

I have kept the call to size(node_coordinates, 1) in calc_vtk_points_cells from #82 and added a comment explaining why. This is needed because the argument node_coordinates is an AbstractArray{<:Any,4}, so we can't access any further information about the dimensions from type parameters. I don't expect this to be particularly problematic, because we already have to do this for the other dimensions of node_coordinates to determine the sizes of the new arrays vtk_points and vtk_cells. See below:

Trixi2Vtk.jl/src/vtktools.jl

Lines 589 to 621 in 6fffc3e

	function calc_vtk_points_cells(node_coordinates::AbstractArray{<:Any,4})
	n_elements = size(node_coordinates, 4)
	size_ = size(node_coordinates)
	n_points = prod(size_[2:end])
	# Linear indices to access points by node indices and element id
	linear_indices = LinearIndices(size_[2:end])
	# Use Lagrange nodes as VTK points. Note that we call size(node_coordinates, 1) in order
	# to provide support for standard two-dimensional meshes as well as meshes representing
	# 2D surfaces in 3D space, which are implemented using P4estMesh{2, 3}.
	vtk_points = reshape(node_coordinates, (size(node_coordinates, 1), n_points))
	vtk_cells = Vector{MeshCell}(undef, n_elements)
	# Create cell for each element
	for element in 1:n_elements
	vertices = [linear_indices[1, 1, element],
	linear_indices[end, 1, element],
	linear_indices[end, end, element],
	linear_indices[1, end, element]]
	edges = vcat(linear_indices[2:end-1, 1, element],
	linear_indices[end, 2:end-1, element],
	linear_indices[2:end-1, end, element],
	linear_indices[1, 2:end-1, element])
	faces = vec(linear_indices[2:end-1, 2:end-1, element])
	point_ids = vcat(vertices, edges, faces)
	vtk_cells[element] = MeshCell(VTKCellTypes.VTK_LAGRANGE_QUADRILATERAL, point_ids)
	end
	return vtk_points, vtk_cells
	end

Note that even when using the covariant solver in TrixiAtmo.jl, which can use an exact representation of the curved surface within the solver internals, the reinterpolation functionality in Trixi2Vtk will always use a polynomial representation of the mesh, and there will therefore be some (usually small) deviation from the true solution on the surface. This is due to the fact that Trixi2Vtk.jl calls Trixi.jl's calc_node_coordinates, which is based on a polynomial mesh representation.

[ranocha]

Do you have a test for the new functionality (2D surface in 3D)?

[tristanmontoya]

Do you have a test for the new functionality (2D surface in 3D)?

I will add one! At first, I didn't understand how this could be done without adding TrixiAtmo.jl as a dependency, but now I see the https://github.com/trixi-framework/Trixi2Vtk_reference_files repository, which I could add some TrixiAtmo.jl output to.

Edit: Would it make more sense to generate the .h5 files for Trixi2Vtk.jl by running TrixiAtmo.jl within the test, or could those be stored in https://github.com/trixi-framework/Trixi2Vtk_reference_files alongside the reference.vtu files? There currently is no straightforward way to run simulations on 2D surfaces in 3D outside of TrixiAtmo.jl.

[tristanmontoya]

Running TrixiAtmo.jl within the test is probably not feasible right now, as we haven't had our first release yet. Therefore, I think the only way to properly test this for now is to store both the .h5 and .vtu files in https://github.com/trixi-framework/Trixi2Vtk_reference_files.

[ranocha]

In the long run, we should run TrixiAtmo.jl for the tests. If you do not want to make a public release, you could store both files in the reference repo as you suggested.

[tristanmontoya]

In the long run, we should run TrixiAtmo.jl for the tests. If you do not want to make a public release, you could store both files in the reference repo as you suggested.

I agree; once there's a public release of TrixiAtmo.jl (which will be quite soon), it will be run similarly to Trixi.jl within the tests. In the meantime, I've added tests which use .h5 files stored in the reference repo.