How to solve coupled PDEs with discontinuous boundary conditions?

So I have a system of coupled PDEs that represent the movement of an interface of fluid. The interfaces are functions of time and space, ie. h_1(x,t), h_2(x,t), h_3(x,t). The PDEs are horrendous (so I won’t write in full here) but they are of the form:

$h_{1_t} = f(h_1,h_2,h_3,h_{1_x},h_{2_x},h_{3_x},h_{1_{xx}},h_{2_{xx}},h_{3_{xx}}) \
h_{2_t} = g(h_1,h_2,h_3,h_{1_x},h_{2_x},h_{3_x},h_{1_{xx}},h_{2_{xx}},h_{3_{xx}})
\
h_{3_t} = p(h_1,h_2,h_3,h_{1_x},h_{2_x},h_{3_x},h_{1_{xx}},h_{2_{xx}},h_{3_{xx}})$

I wish to solve in a domain |x| \leq 0.5 and for time say t \in (0,1). The initial conditions are as follows:

$h_1(x,0) = a(x) \quad \text{for} , -0.5 < x \leq 0.1 \
h_2(x,0) = b(x) \quad \text{for} , 0 \leq x \leq 0.1 \
h_3(x,0) = c(x) \quad \text{for} , 0 \leq x < 0.5$

a(x), b(x), c(x) are known functions and join up continuously, such that a(0.1) = b(0.1) and b(0) = c(0). Also a(-1/2) = c(1/2) = 0 to conform with the following boundary conditions:

$h_1(-0.5,t) = 0 \
h_3(0.5,t) = 0 \
h_2(\max(x),t) = h_3(\max(x),t) \
h_2(\min(x),t) = h_1(\min(x),t)$

such that the ends of h_1 and h_3 are pinned at the end of my domain, and h_2 is always smoothly connecting h_1 and h_3.

I wish to solve this system numerically with NDSolve or NDSolveValue. However, I am having some issues with implementation. I am unsure as to how to input the discontinuous boundary conditions?

My current go is something like this:

NDSolve[{
   D[h1[x, t], t] == stuff1,
   D[h2[x, t], t] == stuff,
   D[h3[x, t], t] == stuff,
   h1[x,0] == a[x], h2[x,0] == b[x], h3[x,0] == c[x], h1[-0.5,t] == 0, h3[0.5,t] == 0
  },
 {h1,h2,h3}, {x,-0.5,0.5}, {t,0,1}
]

Clearly this isn’t working, or handing the h_2 boundary condition. The functions a(x), c(x), b(x) are all InterpolatingFunction running for their respective ranges (ie. a(x) is defined from -0.5 < x < 0.1 etc).

Any ideas on how to deal with this type of boundary condition?