Prescribe discontinuity on internal boundary

1

Here is a model problem for what I want to do. Consider a unit square \Omega with a smaller square \Omega_i fully contained in \Omega. Let \Omega_e = \Omega \setminus \Omega_i be the region inside the unit square but outside the smaller square. We want to solve Poisson’s equation in \Omega_e and \Omega_i:

- \nabla ^2 \Theta_i = 0 \text{ in } \Omega_i \\ - \nabla ^2 \Theta_e = 0 \text{ in } \Omega_e

I want the normal derivative to be smooth across \partial \Omega_i:

\mathbf{n} \cdot \nabla \Theta_i = - \mathbf{n} \cdot \nabla \Theta_e \text{ on } \partial \Omega_i

I also want to specify the discontinuity in \Theta across the boundary:

\Theta_e - \Theta_i = V_m \text{ on } \partial \Omega_i

for a given function V_m defined on \partial \Omega_i.

I believe the weak form of this problem (not including that last condition on the discontinuity across \partial \Omega_i) is:

\int_{\Omega_i} \nabla \Theta_i \cdot \nabla v_i \, dx + \int_{\Omega_e} \nabla \Theta_e \cdot \nabla v_e \, dx + \int_{\partial \Omega_i} (\nabla \Theta_e \cdot v_e - \nabla \Theta_i \cdot v_i) ds = 0

So far, I think I can implement this in fenics with the rudimentary knowledge I’ve pieced together from tutorials (in this example code I impose Dirichlet BC on the left and bottom of the outer square, just to make the solution nontrivial):

from fenics import *
from mshr import *

# Define constants for marking subdomains/boundaries                                                                                                                                                      
OMEGA_E = 1
OMEGA_I = 2
INTERFACE = 3
LEFTEDGE = 4
BOTTOMEDGE = 5

# Create domain and mark subdomains                                                                                                                                                                       
Omega = Rectangle(Point(0, 0), Point(1, 1))
Omega_i = Rectangle(Point(0.2, 0.2), Point(0.4, 0.4))
Omega.set_subdomain(OMEGA_E, Omega - Omega_i)
Omega.set_subdomain(OMEGA_I, Omega_i)
mesh = generate_mesh(Omega, 256)
subdomains = MeshFunction("size_t", mesh, 2, mesh.domains())

# Mark boundaries                                                                                                                                                                                         
boundaries = MeshFunction("size_t", mesh, 1)
for f in facets(mesh):
    p0 = Vertex(mesh, f.entities(0)[0])
    p1 = Vertex(mesh, f.entities(0)[1])
    x0, y0 = p0.x(0), p0.x(1)
    x1, y1 = p1.x(0), p1.x(1)
    on_vert_edge = lambda x, y: (near(x, .2) or near(x, .4)) and y > .2 and y < .4
    on_horiz_edge = lambda x, y: (near(y, .2) or near(y, .4)) and x > .2 and x < .4
    on_edge = lambda x, y: on_vert_edge(x, y) or on_horiz_edge(x, y)
    if on_edge(x0, y0) and on_edge(x1, y1):
        boundaries[f] = INTERFACE
    elif x0 < DOLFIN_EPS and x1 < DOLFIN_EPS:
        boundaries[f] = LEFTEDGE
    elif y0 < DOLFIN_EPS and y1 < DOLFIN_EPS:
        boundaries[f] = BOTTOMEDGE

# TEST                                                                                                                                                                                                    
# Define a function of ones on the mesh:                                                                                                                                                                  
# V = FunctionSpace(mesh, "CG", 1)                                                                                                                                                                        
# f = Function(V)                                                                                                                                                                                         
# f.vector()[:] = 1.                                                                                                                                                                                      

# dStest = dS(subdomain_data=boundaries)                                                                                                                                                                  
# perim = assemble(f*dStest(INTERFACE))                                                                                                                                                                   
# print("Perimeter from integrating 1 on the interface: {} (expected 0.8)".format(perim))                                                                                                                 
# END TEST                                                                                                                                                                                                

V = FunctionSpace(mesh, "Lagrange", 1)

dx = Measure("dx")(subdomain_data=subdomains)
dS = Measure("dS")(subdomain_data=boundaries)
ds = Measure("ds")(subdomain_data=boundaries)

Theta = Function(V)
v = TestFunction(V)

n = FacetNormal(mesh)

LHS = inner(grad(Theta), grad(v)) * dx(OMEGA_E)
LHS += inner(grad(Theta), grad(v)) * dx(OMEGA_I)
LHS += inner(n('-'), grad(Theta('-'))) * v('-') * dS(INTERFACE)
LHS -= inner(n('+'), grad(Theta('+'))) * v('+') * dS(INTERFACE)
                                                                                                              
solve(LHS == 0, Theta, [DirichletBC(V, 1.0, boundaries, LEFTEDGE),
                        DirichletBC(V, 2.0, boundaries, BOTTOMEDGE)])

Now I need to implement that condition \Theta_e - \Theta_i = V_m \text{ on } \partial \Omega_i, and I have no idea where to start. One hint is that the paper I got this problem from specifies the weak form of that condition as

\int_{\partial\Omega_i} (\Theta_i - \Theta_e) u \, ds= \int_{\partial\Omega_i} V_m u \, ds

But I’m not sure how to apply that condition in fenics. I’m grateful for any advice that might help point me in the right direction!

1 Like
2

As you want your solution to be discontinuous, you need to use the discontinuous galerkin methods. See for instance http://www.maths.nuigalway.ie/~eoghan/GroupPresentations/Roisin-27Oct2017.pdf

Alternatively, I believe you can use @cdaversin framework for mixed dimensional finite element methods, see: https://arxiv.org/pdf/1911.01166.pdf

1 Like
3

Hi Jorgen, thank you for your reply, and apologies for my delay in writing back.

I need to learn more about discontinuous Galerkin methods, as simply changing V = FunctionSpace(mesh, "Lagrange", 1) to V = FunctionSpace(mesh, "DG", 1) resulted in a null solution. Presumably, a new weak form of the problem is needed. When using DG method, how does one penalize some jumps (away from the boundaries) but prescribe others (on the boundaries)? Essentially, It’s not clear to me how to obtain the weak form on slide 8 of the presentation you linked.

Also, since the discontinuity in this problem is at a fixed location, and the solution should be continuous elsewhere, I’m wondering whether it makes sense for this problem to be solved using continuous methods, but separate meshes for \Omega_e and \Omega_i. The condition which specifies the discontinuity between the two meshes would be difficult to write, and I’m not sure if it’s even possible in fenics. I may be barking up the wrong tree here, but let me know how this sounds to you.

4

To do DG, nitsches method for enforcing Dirichlet conditions can be used. There is a vast amount of litterature on this (search for nitsches method and weak enforcement of Dirichlet conditions on Google scholar) and the only difference from the other internal interfaces will be that jump(u)=g. This can be restricted in dolfin using Measure and MeshFunction’s. You can for instance see this undocumented dolfin demo: https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/undocumented/dg-poisson/demo_dg-poisson.py

You could also read up on MultiMesh fem, Which imposes g=0 at an interface between two meshes: https://arxiv.org/pdf/1804.06455.pdf

However, as i also mentioned in the previous post, you could use the Mixed dimensional Finite element framework (paper linked to above). Here is a link to examples for that paper: https://zenodo.org/record/3525001#.XuR0feRHklQ

5

This is all very helpful, thanks for the resources and key terms. I’m going to take some time to digest this and do some research.

I had completely misunderstood that your suggestion about mixed dim FEM was not an addition to the DG methods, but rather an alternative, so I’d like to go back and see if I can get something like that to work.

6

The mixed-dim framework looks promising right now, but I’m curious whether the MultiMesh* classes in fenics allow users to define forms that integrate variables from distinct meshes of the same topological dimension, as I’m trying to do.

This may be a moot point since it can be done by @cdaversin 's mixed-dimensional framework (with examples no less), but I’m just curious and trying to get acquainted with fenics a bit. What are the use cases for the MultiMesh classes?

7

There is plenty of literature on the MultiMesh FEm available. It is most commonly used to enforce continuity between non-matching meshes, with application to shape optimization and FSI: https://scholar.google.no/scholar?start=0&q=MultiMesh+finite+element+method+fenics&hl=en&as_sdt=0,5 (the first two pages contains plenty of relevant papers, those by Logg, Johansson, Dokken etc)

8

Ah I see - the MultiMesh* classes are the fenics implementation of the MultiMesh FEM you mentioned earlier. Thanks again for the quick reply. Sounds like I have a lot of options and will not be able to escape learning the Nitsche method, try as I might… :wink:

9

Hi Vyassa,

I am trying to solve a similar problem and I am stuck.
Were you able to solve yours?

Thanks

10

Unfortunately I did not pursue this much further, but I would be interested to hear if you succeed. Good luck