Calculation of discrete weak gradient

We are trying to use freefem++ to address computational aspects of weak Galerkin finite element method.

To do so, in the first step, we need to solve the following problem in the finite dimensional settings which reads as:
Screenshot 2023-05-17 163156

We have used the following finite element spaces: For u0 we are using P0 elements, ub , we are using P1 elements and for H(Div ,Ω), we are using [P1]^2 elements.

In order to compute the discrete weak gradient we are proceeding as follows:

Note that to calculate the discrete weak Gradient, we need to solve the system (1) in finite dimensional settings.

We are solving the problem in the unit square.

int P=10;
mesh Th=square(P,P);
fespace Vh(Th, [P1,P1]);
fespace Wh1(Th, P1);
fespace Wh2(Th, P0);
Wh1 ub,vb; Wh2 u0,v0;
Vh [v1,v2],[w1,w2];
Wh2 div;
Wh1 vn;
vn = v1*N.x+v2*N.y;
solve dgrad(w1,w2,u0,ub,v1,v2,div,vn)= int2d(Th)(w1*v1+w2*v2) +int2d(Th)(u0*div)- int1d(Th)(ub*vn)+on(1,2,3,4,ub=0) ;

But, after compilation, we are observing the following error:

Error umpfack umfpack_di_numeric status 1
Error umfpack_di_solve status 1
– Solve :
min nan max nan
min 0 max 0
min nan max nan

I would be happy if anybody can help me with this.
Thanks in advance.

I’m not sure if you can mix and match vector fespaces like that.
Personally I would start with all scalar fespaces and see
if that works although there may be better ways to user the vector
my FF is down right now with fortran issues or I would run ia
modified version…

I’m new to FF and don’t have a lot of experience with FEM but there are a number
of things here that are not real common.
I guess in terms of the FF code, generally you pass 2n parameters to problem solve
or varf where the first n are the unknowns and the latter are the placeholders for the
test functions determined by the fespace elements. In any case, div and vn are AFAICT
identically zero. But do you have more context around (1) and how it
realtes to the FF code? This appears to have a vector test function with a scalar
unknown The first two volume integrals can be related by integration by parts.