resolve merge

src/equations/reference_data.jl

@@ -71,6 +71,7 @@ the test suite described in the following paper:
     return SVector(h, vx, vy, vz, 0.0f0)
 end
 
+# Version for spherical coordinates
 @inline function initial_condition_gaussian(x, t,
                                             ::AbstractCovariantEquations{2, 3,
                                                                          GlobalSphericalCoordinates})

src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl

@@ -131,7 +131,7 @@ function init_elements_2d_manifold_in_3d!(elements,
     (; node_coordinates, jacobian_matrix,
     contravariant_vectors, inverse_jacobian) = elements
 
-    calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis)
+    Trixi.calc_node_coordinates!(node_coordinates, mesh, basis)
 
     for element in 1:Trixi.ncells(mesh)
 
@@ -172,58 +172,6 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis.nodes)
 end
 
-# Interpolate tree_node_coordinates to each quadrant at the specified nodes
-function calc_node_coordinates_2d_shell!(node_coordinates,
-                                         mesh::P4estMesh{2},
-                                         nodes::AbstractVector)
-    # We use `StrideArray`s here since these buffers are used in performance-critical
-    # places and the additional information passed to the compiler makes them faster
-    # than native `Array`s.
-    tmp1 = Trixi.StrideArray(undef, real(mesh),
-                             Trixi.StaticInt(3),
-                             Trixi.static_length(nodes),
-                             Trixi.static_length(mesh.nodes))
-    matrix1 = Trixi.StrideArray(undef, real(mesh),
-                                Trixi.static_length(nodes),
-                                Trixi.static_length(mesh.nodes))
-    matrix2 = similar(matrix1)
-    baryweights_in = Trixi.barycentric_weights(mesh.nodes)
-
-    # Macros from `p4est`
-    p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL
-    p4est_quadrant_len(l) = 1 << (Trixi.P4EST_MAXLEVEL - l)
-
-    trees = Trixi.unsafe_wrap_sc(Trixi.p4est_tree_t, mesh.p4est.trees)
-
-    for tree in eachindex(trees)
-        offset = trees[tree].quadrants_offset
-        quadrants = Trixi.unsafe_wrap_sc(Trixi.p4est_quadrant_t, trees[tree].quadrants)
-
-        for i in eachindex(quadrants)
-            element = offset + i
-            quad = quadrants[i]
-
-            quad_length = p4est_quadrant_len(quad.level) / p4est_root_len
-
-            nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+
-                           quad.x / p4est_root_len) .- 1
-            nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+
-                           quad.y / p4est_root_len) .- 1
-            Trixi.polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x,
-                                                   baryweights_in)
-            Trixi.polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y,
-                                                   baryweights_in)
-
-            Trixi.multiply_dimensionwise!(view(node_coordinates, :, :, :, element),
-                                          matrix1, matrix2,
-                                          view(mesh.tree_node_coordinates, :, :, :,
-                                               tree), tmp1)
-        end
-    end
-
-    return node_coordinates
-end
-
 # Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping,
 # using eq (12) of :
 #   Giraldo, F. X. (2001). A spectral element shallow water model on spherical geodesic grids.