Dual mortar

[dewenyushu]

Reason

The dual basis is often utilized in the dual Lagrange approach. Here, the shape functions for the Lagrange multipliers (LMs) satisfy the biorthogonal condition across the interface. The biorthogonal condition results in a diagonalized local mass matrix, which allows to condense out the LMs and transform the saddle point problem in to a positive definite one.

Design

• Add dual shape functions in libmesh (see issue 2550)

• Proof of concept using diffusion with EqualValueConstraint

• Enable dual mortar for mechanical contact

• Interface for condensing out the LMs

Impact

Improve the performance for mortar-based contact problems by producing a positive definite system and allows a broader range of preconditioners/solvers

[dewenyushu]

Mortar formulation results in a singular Jacobian for conforming mesh while using Dirichlet BCs. Two (nearly) zero singular values are associated with the Lagrange multipliers DOFs at the boundary. Refs #15216 (comment)

[dewenyushu]

Here I am including some initial verifications of dual mortar implementation using MMS. The result is shown for a 2D diffusion problem in a unit square with analytical solution u(x,y) = sin(2pix)sin(2piy). The domain is composed of two subdomains (left and right), where solution continuity is enforced via using Lagrange multipliers (LMs) along the vertical interface. In this case the LM, \lambda=-2pisin(2pi*y).

The analytical solution is shown below:

The LM also converges with refined mesh:

The upper-right sub-block of Jacobian (which couples slave and LM) is simplified/diagonalized:

[dewenyushu]

@lindsayad @jwpeterson here are some initial verification results. I will be working on different element types and convergence studies as @lindsayad suggested.

[lindsayad]

@lindsayad @jwpeterson here are some initial verification results. I will be working on different element types and convergence studies as @lindsayad suggested.

Yea some L2/H1 error plots would be great. Very neat figures!

[jwpeterson]

@dewenyushu these results look very promising indeed, the sparsity patterns are especially interesting to see. I think the relative sizes of the u1 and u2 blocks indicate that you have chosen the LM/slave variable to be on the coarser mesh? From what I recall, this approach gives better results than the converse, but I'm also curious whether it makes any difference to your results. When the method is used in practice, sometimes it is not possible to pick a "coarse" side and obviously users can always set up problems incorrectly.

Also I was asking @lindsayad about it previously, but for the dual LM method, the dual basis functions are piecewise-discontinuous, so I would expect the LM solution to also be piecewise discontinuous. I can't tell from your figure if that's the case or whether the LM variable is actually continuous? Can you comment on this/help me understand this aspect better?

[lindsayad]

I think the relative sizes of the u1 and u2 blocks indicate that you have chosen the LM/slave variable to be on the coarser mesh? From what I recall, this approach gives better results than the converse, but I'm also curious whether it makes any difference to your results. When the method is used in practice, sometimes it is not possible to pick a "coarse" side and obviously users can always set up problems incorrectly.

I could have sworn that you advised picking the finer side as the slave side, but going back through your technical report, I actually don't see any statement regarding about a choice of side. In your test case 9.2 the slave side is about twice as fine as the master side.

[dewenyushu]

think the relative sizes of the u1 and u2 blocks indicate that you have chosen the LM/slave variable to be on the coarser mesh? From what I recall, this approach gives better results than the converse, but I'm also curious whether it makes any difference to your results. When the method is used in practice, sometimes it is not possible to pick a "coarse" side and obviously users can always set up problems incorrectly.

Yes, I am using a coarser mesh for the slave. Actually I did not choose it intentionally or know that would give me a better performance. So I agree it would be beneficial to try a finer mesh on the slave side. I will remember to report back what I see!

for the dual LM method, the dual basis functions are piecewise-discontinuous, so I would expect the LM solution to also be piecewise discontinuous. I can't tell from your figure if that's the case or whether the LM variable is actually continuous? Can you comment on this/help me understand this aspect better?

Great question! The dual basis I am using is actually piecewise continuous for the first order Lagrange. According to my understanding, the dual basis must satisfy the biorthogonal condition:

where \phi is the primal basis and \psi is the dual basis. Therefore, theoretically there are many choices of the dual basis. Here, I am referring to Alexander Popp's paper and assuming the dual basis to be linear combinations of the primal basis, i.e.,

This is why my dual basis is actually continuous, and so is the LM field. I think this choice should have some benefit in terms of convergence, but I did not find it explicitly stated in literature. But yea, if you choose the primal basis to be discontinuous, then you will get discontinuous dual basis as well.

[dewenyushu]

Here is the 1st order and 2nd order LAGRANGE primal-dual basis comparison:

[jwpeterson]

In your biorthogonality equation, is the domain "\Gamma_e" a single edge in the discretization (Popp's Eqn (18)) or is it the entire contact interface (Popp's Eqn (17))?

I think you are referring to Popp's Eqn (17), in which case the dual basis will be continuous, but you have to solve a "global" (i.e. on the entire contact interface) L2-projection problem to get the continuous dual basis function coefficients. From what I can tell, you have actually implemented Popp's (18) in libmesh, i.e. local L2-projections on an element-by-element basis, which should in general lead to a discontinuous dual basis.

The dual Lagrange basis functions you posted above are discontinuous because, unlike the primal Lagrange basis functions, they do not go to zero at the element boundary on either side of the element... This is difficult to explain in words, but you can also look at e.g. Fig. 3.2 in (Wohmuth 2000) to try and get some idea of what I am talking about.

All that being said, your method seems to work so there is probably more that I am just not understanding about the approach...