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

variational formulation
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