The region is splited into Omega1 and Omega2 by a interface Gamma. When the solution is discontinuous on Gamma, we need to consider the jump conditions in weak formulation, assume [\beta*u \dot n]=g, how I calculate <g,v>, where <> denotes L2 inner product on Gamma, v is a test function. My code is as follows, I think the part the confused me is wrong.
from fenics import *
from mshr import *
import numpy as np
import time
# Define variables
ax, ay, bx, by = -1.0, -1.0, 1.0, 1.0
Nx = Ny = 20
r0 = 1 / 2
# Define geometry for background and subdomains
domain = Rectangle(Point(ax, ay), Point(bx, by))
out = Rectangle(Point(ax, ay), Point(bx, by)) - Circle(Point(0, 0), r0)
circle = Circle(Point(0, 0), r0)
# Set subdomains
domain.set_subdomain(2, out)
domain.set_subdomain(1, circle)
start_time = time.time()
# Generate mesh
mesh = generate_mesh(domain, Nx)
meshfile = File("mesh.pvd")
meshfile << mesh
# Display max mesh size
hmax = mesh.hmax()
print('max h', hmax, 'for mesh N', Nx)
# Function space and subdomain markers
V = FunctionSpace(mesh, 'P', 1)
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
# Define material properties
alfa = 2
beta1 = 1
beta2 = 1
class Beta(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = beta1
else:
values[0] = beta2
beta = Beta(markers, degree=0)
# Define w
class W(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = 10*(x[0]+x[1])-(5*np.exp(pow(x[0],2)+pow(x[1],2))+20)
else:
values[0] = 0
w = W(markers, degree=0)
# Define boundary condition
p_D = Expression('10*(x[0]+x[1])', degree=2)
def boundary(x):
tol = 1e-15
return abs(x[0] - ax) < tol or abs(x[0] - bx) < tol or abs(x[1] - ay) < tol or abs(x[1] - by) < tol
bc = DirichletBC(V, p_D, boundary)
# Define measures with subdomain data
dx = Measure('dx', domain=mesh, subdomain_data=markers)
ds = Measure("ds", domain=mesh, subdomain_data=None)
# Define source term
class F(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = (20*pow(x[0],2)+10)*np.exp(pow(x[0],2)+pow(x[1],2))
else:
values[0] = 0
f = F(markers, degree=0)
u2 = Expression("10*(x[0]+x[1])", degree=2)
u1 = Expression("5*np.exp(pow(x[0],2)+pow(x[1],2))+20", degree=2)
# Define variational problem
n = FacetNormal(mesh)
p = TrialFunction(V)
v = TestFunction(V)
a = beta * dot(grad(p), grad(v)) * dx
L = f * v * dx(1) + f * v * dx(2) + beta * dot(grad(w), grad(v)) * dx - dot(u2-u1, n) * v *ds
# Compute solution
p = Function(V)
solve(a == L, p, bc)
u = p-w
# Save solution to file
vtkfile = File('u.pvd')
vtkfile << u
end_time = time.time()
elapsed_time = end_time - start_time
print('time:', elapsed_time)
# Define exact solution
class u_exact(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = 5*np.exp(pow(x[0],2)+pow(x[1],2))+20
else:
values[0] = 10*(x[0]+x[1])
u_e = u_exact(markers, degree=2)
# Interpolate exact solution for error calculation
u_e_ = interpolate(u_e, V)
file = File("u_e.pvd")
file << u_e_
# Compute error in L2 norm
error_L2 = errornorm(u_e, u, 'L2')
u_L2 = norm(u, 'L2')
# Compute maximum error at vertices
vertex_values_u_e = u_e_.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u_e - vertex_values_u))
r_L2 = error_L2/u_L2
num_dofs = V.dim()
print("自由度的个数:", num_dofs)
# Print errors
print('error_L2:', error_L2, 'for N =', Nx)
print('relative_error_L2:', r_L2, 'for N =', Nx)
print('error_max:', error_max, 'for N =', Nx)
