Biharmonic equation using Argyris elements

[nilsfriess]

Hi, I'm currently going through the "Code generation for generally mapped finite elements" paper (https://arxiv.org/pdf/1808.05513). I was also looking at the code in the supplementary material and I then tried to solve the biharmonic equation using Argyris elements and the C0 IP method as described in the paper. Here's my code:

from firedrake import *

def delta(u):
    return div(grad(u))

def uexact(x, y):
    return Constant(1) * (x*(1-x)*y*(1-y))**2

mesh = UnitSquareMesh(16, 16)

if False:
    V = FunctionSpace(mesh, "Argyris", 5)

    u = TrialFunction(V)
    v = TestFunction(V)
    
    nu = Constant(0.25)    
    a = (inner(delta(u), delta(v))*dx -
         (1-nu)*(2*inner(u.dx(0).dx(0), v.dx(1).dx(1)) +
                 2*inner(u.dx(1).dx(1), v.dx(0).dx(0)) -
                 4*inner(u.dx(0).dx(1), v.dx(0).dx(1)))*dx)
else:
    V = FunctionSpace(mesh, "Lagrange", 3)
    u = TrialFunction(V)
    v = TestFunction(V)
    
    n = FacetNormal(mesh)
    h = CellSize(mesh)
    h_avg = (h('+') + h('-'))/2.0

    alpha = Constant(8.0)
    a = (inner(delta(u), delta(v))*dx -
         inner(avg(delta(u)), jump(grad(v), n))*dS -
         inner(jump(grad(u), n), avg(delta(v)))*dS +
         alpha/h_avg*inner(jump(grad(u), n), jump(grad(v), n))*dS)
    
x, y = SpatialCoordinate(mesh)
uex = uexact(x, y)
f = delta(delta(uex))
L = inner(f, v)*dx

uh = Function(V, name="Solution")

params = {"snes_type": "ksponly",
          "ksp_type": "preonly",
          "pc_type": "lu",
          "pc_factor_mat_solver_type": "pastix"}
solve(a == L, uh, solver_parameters=params)

error = sqrt(assemble((uex-uh)**2*dx))
print(f"Error: {error}")

outfile = VTKFile("biharmonic.pvd")

W = FunctionSpace(mesh, "Lagrange", 3)
uexf = Function(W, name="Exact solution")
outfile.write(uh, uexf.interpolate(uex))

Neither of the methods computes the correct solution. This is the first time I'm using firedrake and UFL code, so there might be some obvious error in my code that I'm just not seeing.

Edit: I actually forgot to specify boundary conditions for the interior penalty method. If I add them, it works fine. But for the Argyris element, I'm not sure what I'm missing.

[wence-]

In the argyris case, you need to apply strong boundary conditions using Nitsche's method (rather than DirichletBC). It looks like you are missing those Nitsche terms.

[nilsfriess]

That makes sense. I just copied some code from the supplementary material of the paper, and I thought that should probably just work. But I think the part I copied is not actually trying to solve the linear system, it's only measuring some setup cost or something like that. I'll try to incorporate the Dirichlet conditions weakly and see if that works.

[rckirby]

The code *should* be available (check the reproducilbity/zenodo bits at the end of the paper). But that was 5 years ago, so it would be interesting to see if the firedrake installer will *really* install the old branch on today's OS. It is also possible just to grab the biharmonic codes from the repo and see if they run against modern firedrake.

…

On Tue, Aug 27, 2024 at 2:25 PM Nils Friess ***@***.***> wrote: That makes sense. I just copied some code from the supplementary material of the paper, and I thought that should probably just work. But I think the part I copied is not actually trying to solve the linear system, it's only measuring some setup cost or something like that. I'll try to incorporate the Dirichlet conditions weakly and see if that works. — Reply to this email directly, view it on GitHub <#3746 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4XUL7XE6OMHMHRVZNKI3LZTTHEDAVCNFSM6AAAAABNF5HMK6VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTANBWGY4TIOA> . You are receiving this because you are subscribed to this thread.Message ID: <firedrakeproject/firedrake/repo-discussions/3746/comments/10466948@ github.com>

[nilsfriess]

I was just looking at the wrong example from the supplementary material. Using the correct bilinear form with the Nitsche terms gives the correct solution.

Different/Follow-up question: I want to use firedrake only for the matrix assembly, and then use a custom preconditioner to solve. Actually, it's not really a preconditioner and I'm also not really solving, I just need the matrix and in the end I'll have a bunch of PETSc Vecs that I'd like to visualise. What would be the right way to do this? I.e., what is the missing code below?

A = assemble(a).M.handle
x = MyExternalFunctionThatReturnsAPetscVec(A)
# ???
outfile = VTKFile("Result.pvd")
outfile.write(??)

And maybe another follow-up question because it feels like I'm doing something wrong: Solving the equation with solve(a == L, ...) takes ages (the first time I run it). For a 8x8 UnitSquareMesh (which end up being 694 dofs using the Argyris element) it takes about 6 minutes on my machine which seems too long. Is this just the expected overhead of the UFL compilation (is that the right term?)?