Utility for symmetric metric

[ddundo]

I think that it would be useful to add a utility to make the metric symmetric, as we'd want to sometimes have a symmetric mesh (when the problem is symmetric). I guess we can do this easily by copying the metric, rotating it and then taking the mean of the two?

[jwallwork23]

Is it not sufficient to just apply the enforce_spd method?

[jwallwork23]

Oh hang on, do you mean symmetric in space rather than within the entries at a given point?

[ddundo]

Ah sorry for not being clear! I did mean in space, yes - so that the generated mesh would be symmetric.

[ddundo]

And I think my suggestion would only work if the original mesh is symmetric so that the new and old nodes would still align after the transformation

Also saw this repo now https://github.com/MmgTools/MirrorMesh

[ddundo]

Hey @jwallwork23, could you please tell me what you think about the workflow below? I modified the simple_metric.py animate demo to get this. In particular, I was wondering if there is a way to make a deepcopy of a mesh rather than redefining it (labelled below with TODO)?

from firedrake import *
from firedrake.pyplot import *
import matplotlib.pyplot as plt
from animate import *

mesh = UnitSquareMesh(40, 40)
P1_ten = TensorFunctionSpace(mesh, "CG", 1)
metric = RiemannianMetric(P1_ten)

x, y = SpatialCoordinate(mesh)
r = 0.4
h1 = 0.02
h2 = 0.05
high = as_matrix([[1 / h1**2, 0], [0, 1 / h1**2]])
low = as_matrix([[1 / h2**2, 0], [0, 1 / h2**2]])
metric.interpolate(conditional(sqrt((x - 0.3) ** 2 + (y - 0.3) ** 2) < r, high, low))
density, *_ = metric.density_and_quotients()

# TEMPORARY METRIC ON A ROTATED DOMAIN
mesh_rotated = UnitSquareMesh(40, 40)  # TODO: Can we make a deepcopy of mesh?
P1_ten_rotated = TensorFunctionSpace(mesh_rotated, "CG", 1)
metric_rotated_temp = RiemannianMetric(P1_ten_rotated)
metric_rotated_temp.interpolate(metric)
mesh_rotated.coordinates.dat.data[:] = 1 - mesh_rotated.coordinates.dat.data

# INTERPOLATE TEMPORARY ROTATED METRIC ON THE ORIGINAL MESH
metric_rotated = RiemannianMetric(P1_ten)
metric_rotated.interpolate(metric_rotated_temp)
density_rotated, *_ = metric_rotated.density_and_quotients()

# COMPUTE SYMMETRIC METRIC ON THE ORIGINAL MESH
metric_symmetric = RiemannianMetric(P1_ten)
metric_symmetric.assign(0.5 * (metric + metric_rotated))
density_symmetric,*_ = metric_symmetric.density_and_quotients()

amesh = adapt(mesh, metric)
amesh_rotated = adapt(mesh, metric_rotated)
amesh_symmetric = adapt(mesh, metric_symmetric)

fig, axes = plt.subplots(2, 3, figsize=(10, 6.5), sharex=True, sharey=True)
axes[0, 0].set_title("Metric density")
tripcolor(density, axes=axes[0, 0])
axes[0, 1].set_title("Rotated m. density")
tripcolor(density_rotated, axes=axes[0, 1])
axes[0, 2].set_title("(Density + Rotated density)/2")
tripcolor(density_symmetric, axes=axes[0, 2])

triplot(amesh, axes=axes[1, 0])
triplot(amesh_rotated, axes=axes[1, 1])
triplot(amesh_symmetric, axes=axes[1, 2])
for ax in axes.flatten():
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.set_aspect("equal")
plt.show()

[jwallwork23]

Would it not be sufficient to just interpolate the rotated metric rather than rotating the mesh and interpolating between meshes?

To copy a mesh I normally do Mesh(old_mesh.coordinates.copy(deepcopy=True)).

[ddundo]

Would it not be sufficient to just interpolate the rotated metric rather than rotating the mesh and interpolating between meshes?

That won't work if old_mesh isn't symmetric to begin with, since rotated vertices wouldn't match up with existing (not rotated) vertices. So in that case we need to interpolate between meshes. Or am I misunderstanding? If so, could you code/pseudocode what you mean please? :)

To copy a mesh I normally do Mesh(old_mesh.coordinates.copy(deepcopy=True)).

Thanks, but that doesn't seem to create a deepcopy unfortunately (the above example doesn't work). I'll dig more :)

[ddundo]

And just a sidenote... having a symmetric metric doesn't guarantee at all that the mesh will be symmetric. But it should be close enough.

[jwallwork23]

Thanks, but that doesn't seem to create a deepcopy unfortunately (the above example doesn't work). I'll dig more :)

How about Mesh(Function(old_mesh.coordinates))?

[jwallwork23]

That won't work if old_mesh isn't symmetric to begin with, since rotated vertices wouldn't match up with existing (not rotated) vertices. So in that case we need to interpolate between meshes. Or am I misunderstanding? If so, could you code/pseudocode what you mean please? :)

But even with rotating the mesh you still interpolate back into $\mathbb{P}1$ space on the old mesh, so it isn't going to be fully symmetric. I can see that a similar approach would work if you took a supermesh of the old mesh and rotated mesh.

[jwallwork23]

Please could you clarify in what sense you want the metric to be symmetric? I was under the impression you wanted it to be symmetric (e.g.) in the $x$-axis, i.e., invariant to reflection top <-> bottom. But the discussion above refers to rotation, which is a different kind of transformation.

[ddundo]

Will start by saying that if you think this isn't something worth spending time on at the moment, please don't waste time thinking about it. I am just working on this for my lock exchange case and thought I'd implement it here while it's fresh, but can get back to it later :)

How about Mesh(Function(old_mesh.coordinates))?

Sadly not. When I do mesh_rotated = Mesh(Function(mesh.coordinates.copy(deepcopy=True))), the mesh_rotated still seems to use mesh.coordinates.

But even with rotating the mesh you still interpolate back into P 1 space on the old mesh, so it isn't going to be fully symmetric. I can see that a similar approach would work if you took a supermesh of the old mesh and rotated mesh.

Ah great point with supermeshing! And yes, sorry - I am assuming that the original mesh isn't horribly asymmetric.

Please could you clarify in what sense you want the metric to be symmetric? I was under the impression you wanted it to be symmetric (e.g.) in the x -axis, i.e., invariant to reflection top <-> bottom. But the discussion above refers to rotation, which is a different kind of transformation.

In my specific case (lock exchange) I indeed want to have 180deg rotational symmetry. But I thought that this approach is quite general if we wanted to add this utility to animate: transform the metric in whatever way (reflect, rotate...) and take the mean.