For simplicity, consider

-\Delta u_s = f \quad \text{in }\Omega_s,
\\
-\Delta u_f = 0 \quad \text{in }\Omega_f,
\\
u_s = u_f \quad \text{on }\Gamma,
\\
u_s = 0 \quad \text{on }\partial\Omega_s\setminus\Gamma,
\\
u_f = 0 \quad \text{on }\partial\Omega_s\setminus\Gamma,

where \Gamma is the interface between the two domains \Omega_s (solid) and \Omega_f (fluid).

Let V_s=\{ u \mid u=0 \text{ on } \partial\Omega_s\setminus\Gamma \}, V_f=\{ u \mid u=0 \text{ on } \partial\Omega_f\setminus\Gamma \}, and V_\Gamma = \{ \text{functions on }\Gamma \}. Then, the weak form is: find (u_s,u_f,\lambda)\in V_s\times V_f\times V_\Gamma such that

\int_{\Omega_s} \nabla u_s\cdot \nabla v_s dx + \int_{\Omega_s} \nabla u_f\cdot \nabla v_f dx + \int_\Gamma \lambda (v_s-v_f) ds + \int_\Gamma \varphi(u_s-u_f)ds = \int_{\Omega_s} fv_s dx

for all test functions (v_s,v_f,\varphi)\in V_s\times V_f\times V_\Gamma. To solve this, we need multiple FE spaces corresponding to V_s, V_f, and V_\Gamma. A possible way is to use composite finite element spaces. Here is an example:

```
load "msh3"
// solid
mesh Ths = square(50,50, [-1+x, -0.5+y]);
// fluid
mesh Thf = square(50,50, [x, -0.5+y]);
// interface
int[int] labs = [2];
meshL ThL = extract(Ths,label=labs);
fespace Vhs(Ths, P1);
Vhs us;
fespace Vhf(Thf, P1);
Vhf uf;
fespace Lh(ThL, P1);
Lh lambda;
{
// Test functions
Vhs vs;
Vhf vf;
Lh phi;
func f = 1; // source term
solve a(<[us],[uf],[lambda]>,<[vs],[vf],[phi]>) =
int2d(Ths)(dx(us)*dx(vs)+dy(us)*dy(vs))
+ int2d(Thf)(dx(uf)*dx(vf)+dy(uf)*dy(vf))
+ int1d(ThL)(lambda*(vs-vf))
+ int1d(ThL)( phi*(us-uf))
- int2d(Ths)(f*vs)
+ on(1,3,4,us=0)
+ on(1,2,3,uf=0);
}
```

I hope this answer fits your question.