Finite element computations using tmop meshes

[tsokar]

Dear all,
I am currently playing with tmop (both for mesh optimization and mesh fitting).
Everything seemed to be OK until I started to do some finite element computations based on the optimized meshes, as I could not retrieve optimal convergence for higher orders.

As an example, I considered a Laplace problem whose (regular) solution is $u = \exp(-(r - r_0) ^2) \sin(\theta)$ in polar coordinates (centre at (0.5,0.5)). The problem is solved (using mfem) on a square mesh of unit length containing a circular hole at the centre.

Exact Dirichlet BCs are enforced on the outer boundary (the hole is flux-free in this case, so nothing is prescribed). The approximation order is 2.

Two series of meshes are compared:

1. Using second order meshes produced by gmsh.
2. Using meshes obtained through tmop (v4.7): the gmsh meshes are taken, then optimized using mesh-optimizer. The boundary is frozen, and the default settings are used in terms of metric and target ( -mid 2 -tid 2 -fix-bnd). In practice, changing the metric or the target has no obvious influence on the optimized mesh.

The results of the convergence study are depicted below:

• Using plain gmsh meshes, we recover the optimal $O(dofs^{p/2})$ convergence rate;
• Using the tmop optimized meshes, the optimal convergence is lost

The only difference that I can spot is that tmop improves the mesh mainly by moving the mid-side nodes, so that the jacobian is somehow continuous over the mesh, leading to oscillations in its evolution.

Maybe there is a bug in my finite element code, but I do not see what kind of bug could be triggered by only moving the mid-side nodes of the mesh.

I would be interested in your thoughts ?
(And btw, thanks a lot for your work on mfem, the library is very nice !)

[vladotomov]

Hi @tsokar,

We'll take a look at this. Could you please send me the mesh files of the gmsh and tmop meshes?

Do you optimize the starting mesh and then refine, or do you optimize separately every refinement?

In theory, improving the mesh quality wrt certain metric does not guarantee better convergence. And this may be aggravated with curved meshes. But we must take a closer look at the meshes before being certain.

[tsokar]

Thanks a lot for taking time to look at it !
Here are the meshes (testMeshX.msh and optMeshX.mesh for gmsh resp. optimized meshes).
meshes.zip

In practice, I optimize each refinement separately. To give you a bigger picture, I am actually working on mesh fitting, but I was able to trace my problem down to mesh optimization only.

I would be interested in your experience regarding to the choice of the right metric for a given computation objective. I think I understand how the metric can "filter" some desired mesh features, but your experience would be valuable !
I tried to see if this behaviour could be linked to some specific metric / target choice, but I wasn't able to find anything interesting.

[tsokar]

I also tested solving the very same problem with another finite element solver, and I get the same convergence mismatch, so it seems to discard the solver bug...

T.

[termi-official]

Maybe a stupid guess, but have you tried a different quadrature rule or increasing the quadrature order?

[tsokar]

Hi, @termi-official, not a stupid guess, but it does not work (I tried 😉 ).
The problem is that the basis in the physical space which is no more polynomial, depending on the mapping on the element (see below)

T.

[vladotomov]

Hi @tsokar,

I tried the meshes with simple interpolation of a smooth function, and I observed a similar degradation of the rate.

Theoretically this is not so surprising, see details in this paper. The relevant takeaway from the paper is that to control the maximum interpolation error, one must limit the high-order derivatives of the geometric mapping, i.e., keep the mesh less curved. The TMOP optimization probably slightly increases the curvature (as the metric does not focus on curvature), which increases the errors.

Note that the interpolation bounds in this paper are worst-case theoretical estimates. That is, decreasing curvature would not always lead to better errors. There are many examples where curved meshes lead to better accuracy (e.g. fitting a curved interface). But one must use them appropriately, based on the PDE properties.

To guarantee that a mesh optimization procedure improves the solution of a given PDE (that will be solved on the optimized mesh), one should involve the PDE solve in the optimization, i.e., move information both ways (the mesh depends on the PDE solution, which depends on the mesh). We have actually been working on such methods and will publish them soon.

BTW, starting with the coarse meshes and refining them gives better interpolation errors, for both gmsh and tmop, than generating / optimizing each resolution. That's because refinement tends to preserve the starting curvature.

[tsokar]

Dear @vladotomov, thank you for the answer.
In fact, the paper you linked reminded a very nice paper from Sevilla et al.. The curvature can have an influence, but more generally, the regularity of the mapping between the reference and physical element can kill the polynomial reproducibility of the basis functions written in the physical space.
This behaviour can prevent the satisfaction of higher-order (>2) patch tests, although the effect can be unnoticeable upon refinement for real problems.

What I do not understand here is why, even for straight-sided elements, tmop tends to shift the mid-side nodes towards one of the vertices of the triangle. It seems to me that this creates "strongly" non-polynomial bases in the physical space (left gmsh, right tmop):

(position of the nodes (left gmsh, right tmop))

(determinant of the jacobian (Sorry, left is tmop, right is gmsh for this figure !))

Did I mess up with the choice of the metric / target ? Or would it be possible to somehow penalize the movement of the nodes along the edges ?

I also coded a small "higher order patch-test" where I project a quadratic solution on the approximation defined on these meshes. As expected, the error is around machine precision for the initial mesh, and non-zero for the optimized one.

T.

[vladotomov]

Hi @tsokar,

It moves them because this results in a better total integral of the metric.

You can try with other shape metrics, or shape+size metrics like TMOP_Metric_007 / TMOP_Metric_009 / TMOP_Metric_094 with IDEAL_SHAPE_EQUAL_SIZE target. These will produce different meshes.

You can also check the displacement-limiting functionality in TMOP_Integrator::EnableLimiting.

But as I said before, none of these have any notion of optimal polynomial representation or convergence.

We will have methods that address this situation soon.

[tsokar]

Thanks a lot for your answer @vladotomov.
I will take a look at it, and I am looking forward to reading your next paper describing the new method to fix (mitigate ?) this issue !

T.