make computation of averages dimension agnostic

LibTrixi.jl/src/LibTrixi.jl

@@ -2,7 +2,7 @@ module LibTrixi
 
 using OrdinaryDiffEq: OrdinaryDiffEq, step!, check_error, DiscreteCallback
 using Trixi: Trixi, summary_callback, mesh_equations_solver_cache, nelements, nvariables,
-             wrap_array, eachelement, cons2prim, get_node_vars, eachnode
+             nnodes, wrap_array, eachelement, cons2prim, get_node_vars, eachnode
 using MPI: MPI, run_init_hooks, set_default_error_handler_return
 
 export trixi_initialize_simulation,

LibTrixi.jl/src/api_jl.jl

@@ -67,34 +67,43 @@ end
 
 function trixi_load_cell_averages_jl(data_, simstate)
     mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi)
-    nels = nelements(solver, cache)
-    nvars = nvariables(equations)
+    n_elements = nelements(solver, cache)
+    n_variables = nvariables(equations)
+    n_nodes = nnodes(solver)
+    n_dims = ndims(mesh)
 
     # convert C to julia
-    data = unsafe_wrap(Array, data_, nels*nvars)
+    data = unsafe_wrap(Array, data_, n_elements*n_variables)
 
     u_ode = simstate.integrator.u
     u = wrap_array(u_ode, mesh, equations, solver, cache)
 
+    # all permutations of nodes indices for arbitrary dimension
+    node_cis = CartesianIndices(ntuple(i -> n_nodes, n_dims))
+
     # temporary storage for mean value on current element for all variables
-    u_mean = zero(get_node_vars(u, equations, solver, 1, 1, 1))
+    u_mean = get_node_vars(u, equations, solver, node_cis[1], 1)
 
     for element in eachelement(solver, cache)
 
-        u_mean = zero(u_mean)
-
         # compute mean value using nodal dg values and quadrature
-        for j in eachnode(solver), i in eachnode(solver)
-            u_node_prim = cons2prim(get_node_vars(u, equations, solver, i, j, element), equations)
-            u_mean += u_node_prim * solver.basis.weights[i] * solver.basis.weights[j]
+        u_mean = zero(u_mean)
+        for node_ci in node_cis
+            u_node_prim = cons2prim(get_node_vars(u, equations, solver, node_ci, element), equations)
+            weight = 1.
+            for node_index in Tuple(node_ci)
+                weight *= solver.basis.weights[node_index]
+            end
+            u_mean += u_node_prim * weight
         end
 
         # normalize to unit element
-        u_mean = u_mean / 2^ndims(mesh)
+        u_mean = u_mean / 2^n_dims
 
         # copy to provided array
-        for ivar = 0:nvars-1
-            data[element + ivar * nels] = u_mean[ivar+1]
+        # all element values for first variable, then for second, ...
+        for ivar = 0:n_variables-1
+            data[element + ivar * n_elements] = u_mean[ivar+1]
         end
     end