How to define right hand side L as integral over smaller circle?

variational formulation
1

Hello.

I want to solve a Poisson equation in 2D on the unit circle. So far, so clear. However, the right hand side L(v) in the weak form a(u,v)=L(v) is given by the curve integral over a smaller circle with radius 0.5. I.e., with I = \{(x,y)\in \mathbb{R^2}: x^2+y^2=0.5^2\}, the right hand side is given by
L(v) = \int_I v ds.

I read the Fenics tutorial book, but did not really find an answer on how to define such a right hand side L. How could this be done in FEniCS? Thanks in advance for any help and ideas.

2

If you want the right hand side to be a surface (curve integral), you need to resolve the circle in the mesh, i.e. create a mesh where this inner circle is aligning with the mesh facets. Then you can integrate over those facets using the ds measure.

3

Thanks for your quick answer. I am not sure, whether I understand you completely? Do you mean something like what I did in the code below? If yes, do you know what is still wrong? I sadly always get the constant 0 function as numerical solution.

from fenics import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
from mshr import * 

R_e = 1
R_0 = 0.5

domain = Circle(Point(0,0),R_e)
inner_circle = Circle(Point(0, 0), R_0)

domain.set_subdomain(1, inner_circle)

mesh = generate_mesh(domain, 10)
V = FunctionSpace(mesh, 'P', 1) 


boundary_markers = MeshFunction('size_t', mesh, dim=1)
boundary_markers.set_all(9999)

class inner_circ(SubDomain):                                                       
    def inside(self, x, on_boundary):
        return near(x[0]**2+x[1]**2, R_0**2, 1E-14)
bx0 = inner_circ()
bx0.mark(boundary_markers, 1)

class outer_circ(SubDomain):                                                       
    def inside(self, x, on_boundary):
        return near(x[0]**2+x[1]**2, R_e**2, 1E-14)                                   
bx1 = outer_circ()
bx1.mark(boundary_markers, 0)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

#Define VP
eps = 1
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u),grad(v))*dx +  pow(eps,-1)*(u*v*ds(0))
L = v*ds(1)

sol = Function(V)
solve(a==L, sol)
4

There are a few issues in your code.

  1. There is a difference between interior facets (a facet connected to two cells), and an exterior facet, (a facet connected to one cell). An interior facet integral is done by using the "dS"-measure, while an exterior facet is using the "ds"-measure.
  2. Your mesh is quite coarse, and this the circle is not exactly at the radius you are trying to use for marking. You can get the correct tolerance by reducing the tolerance.

In the following code I’ve fixed these things, and also added a print to verify the circumference.

from fenics import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
from mshr import * 
import numpy as np
R_e = 1
R_0 = 0.5

domain = Circle(Point(0,0),R_e)
inner_circle = Circle(Point(0, 0), R_0)

domain.set_subdomain(1, inner_circle)

mesh = generate_mesh(domain, 10)
V = FunctionSpace(mesh, 'P', 1) 


boundary_markers = MeshFunction('size_t', mesh, dim=1)
boundary_markers.set_all(9999)
tol = 5e-2

class inner_circ(SubDomain):                                                       
    def inside(self, x, on_boundary):
        return near(x[0]**2+x[1]**2, R_0**2, tol)
bx0 = inner_circ()
bx0.mark(boundary_markers, 1)

class outer_circ(SubDomain):                                                       
    def inside(self, x, on_boundary):
        return near(x[0]**2+x[1]**2, R_e**2, tol)                                   
bx1 = outer_circ()
bx1.mark(boundary_markers, 0)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dS = Measure("dS", domain=mesh, subdomain_data=boundary_markers)
#Define VP
eps = 1
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u),grad(v))*dx +  pow(eps,-1)*(u*v*ds(0))
L = v*ds(1)
print(assemble(1*dS(1))/(2*np.pi), assemble(1*ds(0))/(2*np.pi))
sol = Function(V)
solve(a==L, sol)
File("boundary.pvd")<<boundary_markers

Note that the best way of debugging boundary markers is by saving them to pvd and visualize them in Paraview.

A more robust way to get the markers between to sub domains of cells is presented in:

A final note is that mshr has been long deprecated, and it is recommended to use external meshing software, see: How i can make geometry and mesh in my program? - #2 by dokken

1 Like
5

Thanks a lot for this very instructive answer! One quick addendum: If I understood your explanations correctly, I would have to define L = v(‘+’)dS(1) instead of L = vds(1), correct?

6

Yes, you would have to select a side “+” or “-”. However, in your case, v is continuous, so you can choose either.

1 Like