Initial thought: could you be generating a NaN in your own code as the mesh is refined? E.g. if you're interpolating some expression that leads to a 0/0 at some points in the domain.
[Buesing, Henrik] I am dividing by Delta_h_facet (see [1]). So for the cell in the origin, the node should move closer to zero. But, since I am using DQ0 it should not be zero, right?
I am also dividing by dt. Should I reformulate s.t. I multiply by dt?
Thank you!
Henrik
[1]
x_func_expr = Expression("x[0]")
y_func_expr = Expression("x[1]")
z_func_expr = Expression("x[2]")
x_func = interpolate(x_func_expr, DG1)
y_func = interpolate(y_func_expr, DG1)
z_func = interpolate(z_func_expr, DG1)
Delta_h_facet = sqrt(jump(x_func)**2
+ jump(y_func)**2
+ jump(z_func)**2)