Hi everyone.
Just curious about a few topics. The summary of my problem is a domain (rectangle) is divided in 3 subdomains where there are jumps across each interface that divides the domains. I have a few questions. My code is below. My main questions are really with respect to this normal. I believe I should have a normal for each subdomain instead of just saying n = FacetNormal(mesh), my current error is ufl.log.UFLException: Shape mismatch. - anyone know how I would go about this? Secondly I'm breaking my problem into 3 problems - one for each subdomain. Is there a neat way that I can just put the entire problem together - or is the way I'm doing this completely wrong?
from dolfin import *
from mshr import *
a=0
b=5.0
c=0.0
d=10.0
tol = 1E-14
v0=100.0
e_1=1.0
e_2=1.0
e_3=1.0
def set_mesh(ref,w,h):
h2=h/2.0
w2=w/2.0
p1 = dolfin.Point(0.,0.)
p2 = dolfin.Point(w,h)
o21=p1
o22 = dolfin.Point(w2,h2)
domain = Rectangle(p1,p2)
o11 = dolfin.Point(0.,h2)
o12 = dolfin.Point(w,h)
o31 = dolfin.Point(w2,0.0)
o32 = dolfin.Point(w,h2)
domain.set_subdomain(1, Rectangle(o11, o12)) #Omega 1
domain.set_subdomain(2, Rectangle(p1, o22)) #Omega 2
domain.set_subdomain(3, Rectangle(o31, o32)) #Omega 2
mesh2d = generate_mesh(domain, ref)
return (mesh2d)
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= ( (c+d)/2.0 - tol)
class Omega_2(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= (c+d)/2.0 - tol and x[0] <= (a+b)/2.0 - tol
class Omega_3(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= (c+d)/2.0 - tol and x[0] >= (a+b)/2.0 +tol
def bound_top(x, on_boundary):
return on_boundary and near(x[1], 10, tol)
def bound_bottom(x, on_boundary):
return on_boundary and near(x[1], 0, tol)
def bound_left(x, on_boundary):
return on_boundary and near(x[0], 0, tol)
def bound_right(x, on_boundary):
return on_boundary and near(x[0], 5, tol)
u_top = Expression("v", degree=2,v=v0)
u_left = Expression("0", degree=2)
u_right = Expression("0", degree=2)
u_bottom = Expression("0", degree=2)
mesh= set_mesh(128,5.,10.)
facet_marker = FacetFunction("size_t", mesh)
facet_marker.set_all(10)
Omega_1().mark(facet_marker, 1)
Omega_2().mark(facet_marker, 2)
Omega_3().mark(facet_marker, 3)
n = FacetNormal(mesh)
ds = ds(subdomain_data=facet_marker)
V = FunctionSpace(mesh, "Lagrange", 2)
bc_t = DirichletBC(V, u_top, bound_top)
bc_b = DirichletBC(V, u_bottom , bound_bottom)
bc_l = DirichletBC(V, u_left , bound_left)
bc_r = DirichletBC(V, u_right, bound_right)
bcs=[bc_t,bc_b,bc_l,bc_r]
u1 = TrialFunction(V)
u2 = TrialFunction(V)
u3 = TrialFunction(V)
n = FacetNormal(mesh)
v = TestFunction(V)
f = Expression("0.0", degree=2)
b1= (e_2/e_1)*inner(dot(grad(u2), n),grad(v) )*ds(2) + (e_3/e_1)*inner(dot(grad(u3), n),grad(v) )*ds(3)
b2= (e_1/e_2)*inner(dot(grad(u1), n),grad(v) )*ds(1) + (e_3/e_2)*inner(dot(grad(u3), n),grad(v) )*ds(3)
b3= (e_1/e_3)*inner(dot(grad(u1), n),grad(v) )*ds(1) + (e_2/e_3)*inner(dot(grad(u2), n),grad(v) )*ds(2)
a1=-inner(grad(u1), grad(v))*dx + b1
a2 = -inner(grad(u2), grad(v))*dx + b2
a3 = -inner(grad(u3), grad(v))*dx + b3
L = f*v*dx
# Compute solution
u1 = Function(V)
solve(a1 == L, u1, bcs)
u2 = Function(V)
solve(a2 == L, u2, bcs)
u3 = Function(V)
solve(a3 == L, u3, bcs)
