Continuity at interface between two subdomains

variational formulation
11

So, the way the existing plot function works is that most UFL objects you pass as the first argument are projected (in the sense of L^2) onto a CG1 (linear, continuous) finite element space for visualization. (If you pass a CG1 function, it is plotted directly, or, if you pass a DG0 function, that can also be plotted directly.) You can mimic this behavior for a piecewise function that is f1 on subdomain 1 and f2 on subdomain 2 by writing a custom projection function:

def split_project(f1,f2,V):
    u = TrialFunction(V)
    v = TestFunction(V)
    uh = Function(V)
    solve(u*v*dx == f1*v*dx(1) + f2*v*dx(2),uh)
    return uh

The return value of this will be a Function in the space V that is the L^2 projection of the piecewise function defined by f1 and f2. If V is a CG1 or DG0 space, that can be plotted with plot directly. (Otherwise, it will be projected onto CG1 to plot.) As @bhaveshshrimali pointed out, you can get a higher-fidelity visualization by choosing V to be a higher-degree DG space, but you would need to write it out to a file and plot with Paraview.

An extra note: L^2-projecting a discontinuous function onto a CG space will give oscillatory results; you see this in your original plot of grad(w)[1,1] (which was projected onto CG1, as explained above), where there is some over- and under-shoot near the discontinuity. One way around such oscillatory behavior is to use a lumped-mass L^2 projection, as demonstrated in my answer here; that is also more computationally-efficient than solving a linear system for the consistent L^2 projection.

2 Likes