pip install -U festim

my_model = F.HTransportProblemDiscontinuous()

interface_1 = 0.2
interface_2 = 0.8

vertices = np.concatenate(
    [
        np.linspace(0, interface_1, num=100),
        np.linspace(interface_1, interface_2, num=100),
        np.linspace(interface_2, 1, num=100),
    ]
)

my_model.mesh = F.Mesh1D(vertices)

material_left = F.Material(D_0=2.0, E_D=0.1)
material_mid = F.Material(D_0=2.0, E_D=0.1)
material_right = F.Material(D_0=2.0, E_D=0.1)

material_left.K_S_0 = 2.0
material_left.E_K_S = 0
material_mid.K_S_0 = 4.0
material_mid.E_K_S = 0
material_right.K_S_0 = 6.0
material_right.E_K_S = 0

left_domain = F.VolumeSubdomain1D(3, borders=[0, interface_1], material=material_left)
middle_domain = F.VolumeSubdomain1D(4, borders=[interface_1, interface_2], material=material_mid)
right_domain = F.VolumeSubdomain1D(5, borders=[interface_2, 1], material=material_right)

left_surface = F.SurfaceSubdomain1D(id=1, x=vertices[0])
right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])

# the ids here are arbitrary in 1D, you can put anything as long as it's not the same as the surfaces
my_model.interfaces = {6: [left_domain, middle_domain], 7: [middle_domain, right_domain]}
my_model.subdomains = [left_domain, middle_domain, right_domain, left_surface, right_surface]
my_model.surface_to_volume = {right_surface: right_domain, left_surface: left_domain}

H = F.Species("H", mobile=True)

my_model.species = [H]

H.subdomains = [left_domain, middle_domain, right_domain]

my_model.boundary_conditions = [
    F.DirichletBC(left_surface, value=0.05, species=H),
    F.DirichletBC(right_surface, value=0.2, species=H),
]


my_model.temperature = lambda x: 300 + 100 * x[0]

my_model.settings = F.Settings(atol=None, rtol=None, transient=False)


my_model.initialise()
my_model.run()