Expand mesh deformation plugin diffusion to spherical model geometries

[anne-glerum]

At the moment, the mesh deformation plugin diffusion is only applicable to cartesian box geometries. However, we think it can be very useful for spherical models with a free surface as well.

For a 2D, cartesian ASPECT mesh, the plugin solves a 1D PDE of topography h(t,x). It does this by 1) taking the accumulated mesh deformation in the vertical direction as the current h, 2) projecting the shape function gradients onto the tangent space of the mesh surface, and 3) only assemblying the RHS and the matrix for the y-component of the DOFs.

After a discussion with @gassmoeller and @bangerth, I tried to generalize the cartesian case by 1) setting normal constraints on the deformation, and 2) assemblying all x and y terms of the equation. However, I realized this is not what we want: we do not want deformation in the normal direction of the surface - this normal changes over time and per surface element. We want deformation only in the vertical direction (i.e., the gravity direction). The image below shows the mesh for the hill_diffusion_field_boundary_conditions.prm test for 1) mainline ASPECT, 2) adding the compute_normal_flux_constraints, and 3) assembling for both the x and y-direction, resp.

I find it hard to wrap my head around another way of doing the spherical diffusion. I've tried projecting the shape function gradients to the plane normal to the vertical direction and vice versa to take the mesh deformation in the normal direction instead of the vertical direction. None of the different combinations give the mainline results though. I also can't think of a way to cancel out contributations to the horizontal direction without explicitely not assemblying them. Do you have an idea?

[gassmoeller]

Hey @anne-glerum, sorry for letting this slide for a while, it is a somewhat tricky problem. Do you have a version of the code where you apply normal flux constraints? If the problem is just to not apply normal flux constraints, but instead "vertical flux constraints" maybe we can create the necessary constraints manually? I.e. we can identify all the relevant constraints in x-, y-, and z-direction that make up the surface of the domain, and make sure that we constrain each of them in such a way that a certain deformation in x-direction has to create a certain direction in y-direction that happens to make sure the vertex can only move in "vertical" direction? To think more about that I would need to see how you integrated the normal flux constraints.