How to correctly define a mixed finite element system in FreeFEM?

General Discussion
1

Now I need to solve a system of equations using FreeFEM, and I’m not sure if the format of the equations I’ve defined is correct. In fact, I only know how to define and solve individual equations using FreeFEM, so I tried to mimic the format used for a single equation and ended up with an obviously incorrect result.

find \left(\mathbf{u}_{h}, \gamma_{h}, d_{h}\right) \in V_{h} \times R_{h} \times R_{h}.

\begin{array}{ll} \left(\gamma_{h}, \mathbf{v}_{h}\right)-\left(\mathbf{u}_{h}, \operatorname{curl} \mathbf{v}_{h}\right)=0 & \forall \mathbf{v}_{h} \in R_{h}, \\ \left(d_{h}, \omega_{h}\right)+\left(\mathbf{u}_{h}, \nabla \omega_{h}\right)=0 & \forall \omega_{h} \in R_{h}, \\ \left(\operatorname{curl} \gamma_{h}, \mathbf{g}_{h}\right)-\left(\operatorname{grad} d_{h}, \mathbf{g}_{h}\right)=\left(\mathbf{f}, \mathbf{g}_{h}\right) & \forall \mathbf{g}_{h} \in V_{h} . \end{array}

R_h is the Lagrange elements with degree k ≥ 1 and V_h is the Raviart-Thomas elements with degree k.

Here is my code, and I’m not sure how to correctly define and solve a system of equations in FreeFEM. Can it only be solved through a matrix-based approach?

func f1=(2*pi*pi)*sin(pi*x)*sin(pi*y);
func f2=(2*pi*pi)*sin(pi*x)*sin(pi*y);
func re1=sin(pi*x)*sin(pi*y);
func re2=sin(pi*x)*sin(pi*y);
macro grad(u) [dx(u), dy(u)] //
macro curl(u) [dy(u),-dx(u)]//
macro div (u) [dx(u.x)+dy(u.y)]//
macro rot (u) [dx(u.y)-dy(u.x)]//
macro dot(u,v)[u.x*v.x+u.y*v.y]//
macro Grad(u) [dx(u), dy(u)]//
mesh Th = square(8,8); 
fespace Ph(Th,P1); Ph v,w,r,d,real1,real2;
fespace Rh(Th,RT0);Rh [u1,u2],[g1,g2];
problem laplaceMixte([u1,u2,r,d], [g1,g2,v,w]) = 
    int2d(Th)(
         r*v - u1*dy(v)+u2*dx(v)
        + d*w+u1*dx(w)+u2*dy(w)
        + dy(r)*g1-dx(r)*g2-dx(d)*g1-dy(d)*g2
    )-int2d(Th)(f1*g1+f2*g2)
	+ on(1,2,3,4, u1=0,u2=0); 
 laplaceMixte;
 real1=re1;
 real2=re2;
 plot([u1,u2],coef=0.1,value=true,wait=1); 
 plot([real1,real2],coef=0.1,value=true);
2

Again,I’m trying to solve the matrix equations by defining them in the form of a matrix equation, and I’d like to use the two components of the RT0 element separately, but I don’t know how to extract the components, so I’m stuck with the form a(Rh, Ph), only to be shown the error that the degrees of freedom don’t match when assembling the matrix:
– Square mesh : nb vertices =81 , nb triangles = 128 , nb boundary edges 32
current line = 21
Assertion fail : (Vh.N == Uh.N)
line :9476, in file problem.cpp
Assertion fail : (Vh.N == Uh.N)
line :9476, in file problem.cpp
err code 6 , mpirank 0
This is the code I’ve tried; does anyone have any suggestions on what I should do?

func f1=(2*pi*pi)*sin(pi*x)*sin(pi*y);
func f2=(2*pi*pi)*sin(pi*x)*sin(pi*y);
func re1=sin(pi*x)*sin(pi*y);
func re2=sin(pi*x)*sin(pi*y);
macro grad(u) [dx(u), dy(u)] //
macro curl(u) [dy(u),-dx(u)]//
macro div (u) [dx(u.x)+dy(u.y)]//
macro rot (u) [dx(u.y)-dy(u.x)]//
macro Grad(u) [dx(u), dy(u)]//
mesh Th = square(8,8); 
fespace Ph(Th,P1); Ph v,w,r,d,real1,real2;
fespace Ph1(Th,P2);
fespace Rh(Th,RT0);Rh [u1,u2],[g1,g2];
varf a(u,v)=int2d(Th)(u*v);
varf b([u1,u2],[v1,v2])=int2d(Th)(u1*v1+u2*v2);
varf a1([u1,u2],[v])=-int2d(Th)(u1*dy(v))+on(1,2,3,4,u1=0);
varf a2([u1,u2],[v])=int2d(Th)(u2*dx(v))+on(1,2,3,4,u2=0);
varf b1([u1,u2],[v])=int2d(Th)(u1*dx(v))+on(1,2,3,4,u1=0);
varf b2([u1,u2],[v])=int2d(Th)(u2*dy(v))+on(1,2,3,4,u2=0);
varf c1([u1,u2],[v])=int2d(Th)(u1*dy(v)-u2*dx(v));
varf c2([u1,u2],[v])=-int2d(Th)(u1*dx(v)+u2*dy(v));
matrix A1=a1(Rh,Ph);matrix A2=a2(Rh,Ph);matrix A3=a(Ph,Ph);// An error occurred.
matrix B1=b1(Rh,Ph);matrix B2=b2(Rh,Ph);matrix B4=a(Ph,Ph);
matrix C3=c1(Rh,Ph);matrix C4=c2(Rh,Ph);
cout << " A1 "<< A1.n << "x" <<A1.m << endl; 
cout << " B1 "<< B1.n << "x" <<B1.m << endl; 
cout <<"Ph"<<Ph.ndof<<endl;
matrix A=[[A1,A2,A3,0],[B1,B2,0,B4],[0,0,C3,C4]];
set(A,solver=UMFPACK);