Dear all,
I am quite confused by following this page: https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html.
I am trying to solve the heat equation in a region that consists of several mshr boxes put one on top of each other, so that there are surfaces, the interfaces between each boxes, where I do not want to apply any boundary conditions. And there are surfaces (the exposed ones) where I want to apply either Dirichlet, Robin or Neumann boundary conditions.
I already have the weak form and it would work if I could assign ds(n) to the nth surface of interest. I haven’t been able to do that, so far.
Here’s my current code
from dolfin import *
import mshr
block_length = 4
block_height = 3
block_sink_length = block_length
block_sink_height = 1
paralle_1 = mshr.Box(Point(0.0, 0.0, 0.0), Point(block_length, block_length, block_height))
paralle_2 = mshr.Box(Point(0.0, 0.0, block_height), Point(block_sink_length, block_sink_length, block_height + block_sink_height))
domain = paralle_1 + paralle_2
mesh = mshr.generate_mesh(domain, 15)
V = FunctionSpace(mesh, 'CG', 1)
t_f = 10
num_steps = 10
dt = t_f / num_steps
T_hot = 273.15+20 # in K.
T_room = 273.15 # in K.
tol = 1E-14
class Bottom(SubDomain):
def inside(self, x , on_boundary):
return on_boundary and near(x[2], 0.0, tol)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[2] >= 0.0 + tol and x[2] <= block_height - tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[2] >= block_height + tol and x[2] <= block_sink_height + 0.01 - tol
class top_block(SubDomain):
def inside(self, x , on_boundary):
return on_boundary and near(x[2], block_height, tol)
class top_blokus(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], block_height + block_sink_height, tol)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0) # mark all surface elements 0
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
subdomain_1 = top_block()
subdomain_1.mark(boundary_markers, 1)
subdomain_2 = top_blokus()
subdomain_2.mark(boundary_markers, 2)
materials = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomain_10 = Omega_0()
subdomain_11 = Omega_1()
subdomain_10.mark(materials, 10)
subdomain_11.mark(materials, 11)
k_0 = 30
k_1 = 22
class K(UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
kappa = K(materials, k_0, k_1)#, k_2, k_3, k_4)#, degree=1)
u_D0 = Constant(T_hot)
bc0 = DirichletBC(V, u_D0, Bottom())
u_0 = Expression('T_room', T_room = T_room, degree = 1)
u_n = interpolate(u_0, V)
u = TrialFunction(V)
v = TestFunction(V)
F = dt * dot(kappa*grad(u), grad(v)) * dx + v *(u-u_n)*dx
a, L = lhs(F), rhs(F)
u = Function(V)
t = 0
for n in range(num_steps):
u_D0.t = t
t += dt
solve(a == L, u, bcs = bc0)
u_n.assign(u)
So, as you can see I have defined several classes. Some are supposed to be volumes, and others are supposed to be surfaces (only the surfaces I care about and where I wish to apply ds(n) with them). However the top_block surface doesn’t seem to work well because when I apply Dirichlet boundary conditions on it by replacing the bc0 line by “bc0 = [DirichletBC(V, u_D0, Bottom()), DirichletBC(V, u_D0, top_block())]”, I get a warning that “*** Warning: Found no facets matching domain for boundary condition.”. So I am unable to do it properly… any idea what’s wrong? I am getting quite confused!
