Adress subdomain of external mesh for diffiusion coeficient

Hi,
Apologies for the beginner-level questions:
there is an external mesh which is imported to the Fenics. the elements of the Mesh are tagged with physical groups.

i wrote a simple diffusion code. But I want to have different diffusion coefficient in in different materials which is showed in the picture with mark 1,2 and 3.
I know that I have to use and define subdomain but unfortunately i could not do it with external mesh.

from fenics import *

t_end = 3
dt = 0.1
k = 3 #diffusion Coefficient
u_in = 20
u_out = -5

xml_file = "01Mesh.xml"
mesh = Mesh(xml_file)
fd = MeshFunction('size_t', mesh, "01Mesh_facet_region.xml");

V = FunctionSpace(mesh, 'P', 7)

bc1 = DirichletBC(V, Constant(u_in), fd, 97)
bc2 = DirichletBC(V, Constant(u_out), fd, 2)
bc = [bc1, bc2]

u = TrialFunction(V)
v = TestFunction(V)
u_n = Function(V)

F = u*v*dx + dt*k*dot(grad(u), grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)

u = Function(V)
t = 0
vtkfile = File('output/output.pvd')

num_steps = int(t_end/dt)
for n in range(num_steps):
    t += dt
    solve(a == L, u, bc)
    u_n.assign(u)
    vtkfile << (u, t)

Many thanks.

See: How to define different materials / importet 3D geometry (gmsh) - #2 by dokken
Here, you can replace the SubDomains that generate the MeshFunction by reading in the MeshFunction from file.

1 Like