# Modified Weak Galerkin Finite Element Method

**URL:** <https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656>\
**Category:** General Discussion\
**Created:** [August 15, 2023, 8:59am UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656 "2023-08-15T08:59:58Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![ahbabil1983](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ahbabil1983/32/1860_2.png) [@ahbabil1983](https://community.freefem.org/u/ahbabil1983)\
**Post date:** [August 15, 2023, 8:59am UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/1 "2023-08-15T08:59:58Z")

</div>

Hello everyone, I want to write a FreeFEM++ program about the Poisson Problem  
-∆u=f in Ω  
u=g on ∂Ω  
using the Modified Weak Galerkin Finite Element Method  
∑\_T〈∇\_w u\_h,∇\_w v\_h 〉_T + ∑\_e h^(-1) 〈[[u\_h]], [[v\_h]]〉e = 〈f, v\_h 〉, ∀ v\_h ∈ V\_h^0  
where [[u]] = u\_0 |(T\_1 ) n\_1+u\_0 |_(T\_2 ) n\_2 is the jump function across the edge e, and the discrete gradient of u on T is given by  
〈∇\_w u\_h,q〉\_T = -〈u\_h,∇.q〉_T + 〈{u\_h },q.n〉∂T,  
{u\_h } = (1/2) (u\_0 |(T\_1 ) +u\_0 |_(T\_2 ) )

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 15, 2023, 12:26pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/2 "2023-08-15T12:26:56Z")

</div>

Something like this[1]? There are discontinuous elements  
This is a first google hit,  
[https://community.freefem.org/t/problem-in-3d-discontinuous-galerkin-computation/2015/5](https://community.freefem.org/t/problem-in-3d-discontinuous-galerkin-computation/2015/5)

I’ve never used these but curious, If you can post a worked example with  
solution I can try to do it in FF if I get a chance.

```auto
[1]A modified weak Galerkin finite element method for the Stokes equations, Mu , Lin[...] Ye , Xiu; Journal of Computational and Applied Mathematics. 2015 
https://www.sciencedirect.com/science/article/pii/S037704271400363X

```

---

<div class="post-metadata">

**Author:** ![ahbabil1983](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ahbabil1983/32/1860_2.png) [@ahbabil1983](https://community.freefem.org/u/ahbabil1983)\
**Post date:** [August 15, 2023, 1:40pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/3 "2023-08-15T13:40:23Z")

</div>

> [@marchywka](#):
>
> I’ve never used these but curious, If you can post a worked example with  
> solution I can try to do it in FF if I get a chance

Thanks for your reply. I found this problem in this paper:  
[1] A modified weak Galerkin finite element method. X. Wang, N.S. Malluwawadu, F. Gao, T.C. McMillan. [https://www.sciencedirect.com/science/article/pii/S0377042714002131](https://www.sciencedirect.com/science/article/pii/S0377042714002131).  
There are some examples of calculations of the modified Galerkin FEM for the Poisson Problem in d = 2, 3.

---

<div class="post-metadata">

**Author:** ![ahbabil1983](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ahbabil1983/32/1860_2.png) [@ahbabil1983](https://community.freefem.org/u/ahbabil1983)\
**Post date:** [August 15, 2023, 2:25pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/4 "2023-08-15T14:25:03Z")

</div>

I used the following FreeFEM++ program to claculate the poisson problem in the classical Galerkin method, but I do not know how to use it to claculate the poisson problem in the modified weak Galerkin method:  
int nref = 4;  
real [int] L2error ( nref ); // initialize the L2 error array  
for ( int n=0;n\< nref ;n++) {  
int N =2^( n +4) ; // space discretization  
mesh Th= square (N,N); // mesh generation of a square  
fespace Vh(Th ,P1); // space of P1 Finite Elements  
Vh uh ,vh; // uh and vh belongs to Vh  
macro Grad (u)[dx(u),dy(u)] //  
Vh uex = sin(pi_x)_ sin (pi_y); // exact solution  
Vh f = 2_sin (pi_x)pi^2 sin (pi_y); // corresponding RHS  
varf a(uh ,vh) = int2d (Th)( Grad (uh)’ \* Grad (vh) ) // bilinearform  
- int2d (Th)(f_vh) // linear form  
+ on (1 ,2 ,3 ,4 , uh = 0); // Dirichlet B.  
matrix A = a(Vh ,Vh); // build the matrix  
varf l( unused ,vh) = int2d (Th)(f_vh); // linear form  
Vh F; F = l(0, Vh); // build the right hand side vector  
set (A, solver = sparsesolver );  
uh = A^ -1\*F ;  
L2error [n] = sqrt ( int2d (Th)((uh - uex )^2));  
}  
for ( int n = 0; n \< nref; n ++){  
int N = 2^(n + 1);  
cout \<\< " N " \<\< " = " \<\< N \<\< endl;  
cout \<\< " L2error " \<\< n \<\< " = " \<\< L2error [n] \<\< endl;}  
for ( int n = 1; n \< nref; n ++){  
cout \<\<" convergence rate = " \<\< log ( L2error [n -1]/ L2error [n])/log (2.) \<\< endl;}

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [August 17, 2023, 12:29am UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/5 "2023-08-17T00:29:04Z")

</div>

See if this helps. Its a copy of an example using the FF discontinuous  
element. It gives a decent error vs the other approach but there are a couple of problems.  
First I think there is a sign mistake I can’t find meaning the result may  
be spurious but I don’t have time to find it right now. If you have not used  
this before the link in the code may be a good starting point. I think in the paper  
you cited, also link in the code, they say discrete L2 and I’m not sure  
if they mean integral or dot product although I guess with a uniform grid it may not  
matter much. It does look like they use base 2 logs for convergence but  
I was playing with that too…  
Also you may want to get the syntax highlighter for FF code. Some of these  
variables like “N” are used elsewhere and if you define your own some  
code will give confusing error messages. If it is highlighted don’t use it  
for your own 🙂

[wg.edp](https://community.freefem.org/uploads/short-url/vgwCTgRSLrjkm2IWiaTdqNAerIw.edp) (3.7 KB)

```auto

// https://www.um.es/freefem/ff++/pmwiki.php?n=Main.DiscontinuousGalerkin
// Discontinous Galerlin Method
// based on paper from 
// Riviere, Beatrice; Wheeler, Mary F.; Girault, Vivette
// title: 
// A priori error estimates for finite element 
// methods based on discontinuous approximation spaces
// for elliptic problems.
// SIAM J. Numer. Anal. 39 (2001), no. 3, 902--931 (electronic).
// ---------------------------------
// Formulation given by Vivette Girault
// ------ 
// Author: F. Hecht , december 2003
// -------------------------------
// nonsymetric bilinear form

int nref = 3;
real [int] L2error ( nref ); // initialize the L2 error array
real [int] L2errord ( nref ); // initialize the L2 error array
for ( int n=0;n< nref ;n++) {
int NN =2^( n +4) ; // space discretization
mesh Th= square (NN,NN); // mesh generation of a square
fespace Vh(Th ,P1); // space of P1 Finite Elements
Vh uh ,vh; // uh and vh belongs to Vh
macro Grad (u)[dx(u),dy(u)] //
Vh uex = - sin(pi*x)* sin (pi*y); // exact solution
Vh f =2.0* sin (pi*x)*(pi^2) * sin (pi*y); // corresponding RHS
varf a(uh ,vh,init=1) = int2d (Th)( Grad (uh)' * Grad (vh) ) // bilinearform
- int2d (Th)(f*vh) // linear form
+ on (1 ,2 ,3 ,4 , uh = 0); // Dirichlet B.
{
/////////////////////////
macro dn(u) (N.x*dx(u)+N.y*dy(u) ) // def the normal derivative 

//fespace Xh(Th,P2);
fespace Vhh(Th,P2dc); // Discontinous P2 finite element
// if param = 0 => Vh must be P2 otherwise we need some penalisation  
real pena=0; // a paramater to add penalisation 
Vhh u,v;
Vhh uexh = - sin(pi*x)* sin (pi*y); // exact solution
//Vhh fh =2.0* sin (pi*x)*(pi^2) * sin (pi*y); // corresponding RHS
func fh =2.0* sin (pi*x)*(pi^2) * sin (pi*y); // corresponding RHS
solve Ans(u,v,solver=sparsesolver)=
//solve Ans(u,v)=
   int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v) )
 + intalledges(Th)(// loop on all edge of all triangle 
       // the edge are see nTonEdge times so we / nTonEdge
       // remark: nTonEdge =1 on border edge and =2 on internal 
       // we are in a triange th normal is the exterior normal
       // def: jump = external - internal value; on border exter value =0
       // mean = (external + internal value)/2, on border just internal value
            ( jump(v)*mean(dn(u)) - jump(u)*mean(dn(v))
          + pena*jump(u)*jump(v) ) / nTonEdge
)
// FIXME this sign is wrong AFAICT but gives good error :) 
+ int2d (Th)(fh*v) // linear form
//+ on (1 ,2 ,3 ,4 , u = 0); // Dirichlet B.
;
//solve Ans(u,v);
//plot(u,wait=1,fill=0,value=1);
L2errord [n] = sqrt ( int2d (Th)((u - uexh )^2));
// I'm not sure if they mean this by "discreate L2-norm in 
// https://www.sciencedirect.com/science/article/pii/S0377042714002131
real tot=0;
for(int i=0; i<u[].n; ++i){ tot+=(u[][i]-uexh[][i])* (u[][i]-uexh[][i]); } 
//L2errord [n] = sqrt ( tot);

}

///////////////////////////////////////////////

matrix A = a(Vh ,Vh); // build the matrix
//varf l( unused ,vh) = int2d (Th)(f*vh); // linear form
//Vh F; 
real[int] F = a(0, Vh); // build the right hand side vector
set (A, solver = sparsesolver );
uh[] = A^ -1*F ;
L2error [n] = sqrt ( int2d (Th)((uh - uex )^2));
}
for ( int n = 0; n < nref; n ++){
int NN = 2^(n + 1);
cout << " NN " << " = " << NN << endl;
cout << " L2error " << n << " = " << L2error [n] << endl;
cout << " L2errord " << n << " = " << L2errord [n] << endl;}
for ( int n = 1; n < nref; n ++)
{
cout <<" convergence rate = " <<( log ( L2error [n -1]/ L2error [n])/log(exp(1.0)) )<< endl;
cout <<" convergence rate d = " << log ( L2errord [n -1]/ L2errord [n])/log (exp(1.0)) << endl;
cout <<" convergence rate e vs d base 2 = " 
	<<( log ( L2error [n -1]/ L2error [n])/log (2.0)) // << endl;
<<" "	<<( log ( L2errord [n -1]/ L2errord [n])/log (2.0)) << endl;

}

```

---

<div class="post-metadata">

**Author:** ![ahbabil1983](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ahbabil1983/32/1860_2.png) [@ahbabil1983](https://community.freefem.org/u/ahbabil1983)\
**Post date:** [August 18, 2023, 8:04pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/6 "2023-08-18T20:04:17Z")

</div>

Thanks for your reply, I just modified your code a little to avoid the sign issue. But I am still confused about how to calculate the algorithm for the Modified Weak Galerkin FEM mentioned in the paper  
[1] A modified weak Galerkin finite element method. X. Wang, N.S. Malluwawadu, F. Gao, T.C. McMillan. [https://www.sciencedirect.com/science/article/pii/S0377042714002131](https://www.sciencedirect.com/science/article/pii/S0377042714002131), where the definitions of the modified discrete gradient and the jump function in the Modified Weak Galerkin FEM are different than the ones given in the discontinuous Galerkin FEM. Moreover, they are not defined or used by anyone before in FreeFEM++.

int nref = 3;  
real [int] L2error ( nref ); // initialize the L2 error array  
real [int] L2errord ( nref ); // initialize the L2 error d array  
for ( int n = 0; n \< nref; n++) {  
int Nn = 2^( n +4) ; // space discretization  
mesh Th = square (Nn,Nn); // mesh generation of a square  
macro dn(u) (N.x_dx(u)+N.y_dy(u) ) // def the normal derivative  
fespace Vh(Th,P2dc); // Discontinous P2 finite element  
// if param = 0 =\> Vh must be P2 otherwise we need some penalisation  
real pena=0; // a paramater to add penalisation  
Vh uh,vh;  
Vh uex = sin(pi_x)_ sin (pi_y); // exact solution  
Vh f = 2.0_ sin (pi_x)_(pi^2) \* sin (pi_y); // corresponding RHS  
// func fh = 2.0_ sin (pi_x)_(pi^2) \* sin (pi\*y); // corresponding RHS  
varf a(uh,vh, solver=sparsesolver)=  
//solve Ans(u,v)=  
int2d(Th)(dx(uh)\*dx(vh)+dy(uh)\*dy(vh) )

- intalledges(Th)(// loop on all edge of all triangle  
// the edge are see nTonEdge times so we / nTonEdge  
// remark: nTonEdge =1 on border edge and =2 on internal  
// we are in a triange th normal is the exterior normal  
// def: jump = external - internal value; on border exter value =0  
// mean = (external + internal value)/2, on border just internal value  
( jump(vh)\*mean(dn(uh)) - jump(uh)_mean(dn(vh))  
+ pena_jump(uh)_jump(vh) ) / nTonEdge);  
matrix A = a(Vh, Vh); // build the matrix  
varf L(unused, vh) = int2d (Th)(f_vh); // linear form  
Vh F; F = L(0, Vh); // build the right hand side vector  
// set(A, solver=sparsesolver);  
uh = A^-1 \* F ;  
L2errord [n] = sqrt ( int2d (Th)((uh - uex )^2));  
}  
for ( int n = 0; n \< nref; n ++){  
int Nn = 2^(n + 4);  
cout \<\< " N " \<\< " = " \<\< Nn \<\< endl;  
cout \<\< " L2errord " \<\< n \<\< " = " \<\< L2errord [n] \<\< endl;}  
for ( int n = 1; n \< nref; n ++) {  
int Nn = 2^(n + 4);  
cout \<\< " N " \<\< " = " \<\< Nn \<\< endl;  
cout \<\<" convergence rate d = " \<\< log ( L2errord [n - 1]/ L2errord [n])/log (2.) \<\< endl; }

---

<div class="post-metadata">

**Author:** ![ahbabil1983](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ahbabil1983/32/1860_2.png) [@ahbabil1983](https://community.freefem.org/u/ahbabil1983)\
**Post date:** [September 1, 2023, 5:42pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/7 "2023-09-01T17:42:55Z")

</div>

I did the following code, but the results are confusing. They are not the same as in the paper:  
[1] A modified weak Galerkin finite element method. X. Wang, N.S. Malluwawadu, F. Gao, T.C. McMillan. [https://www.sciencedirect.com/science/article/pii/S0377042714002131](https://www.sciencedirect.com/science/article/pii/S0377042714002131).  
I am not sure about the jump function and the weak discrete gradient is the same in both the discontinuous Galerkin FEM and the modified Weak Galerkin FEM which is not defined at all in FreeFEM++.

[MWGFEM\_1.edp](https://community.freefem.org/uploads/short-url/xx09ELtGHE2xYY4foGgnYn8NARz.edp) (1.2 KB)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 2, 2023, 1:34am UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/8 "2023-09-02T01:34:31Z")

</div>

I’m trying to make a more complete response but can’t entirely get  
it to work. I did note you have 1/Nn where you need 1.0/Nn in both h  
and the expression for h in the jump term. However, that itself doesn’t  
fix everything. For the continuous case you need to boundary conditions  
as the value of u does not appear anywhere. I always set tgv when  
doing this to keep the matrix clean for the solvers The P2 continuous beats  
their L2 errors. The best I could do with the P2dc with a factor of two fudge  
was a bit worse and IIRC in that case I used “lenEdge” instead of h.  
And I wasn’t sure if h should be hmax which would intridouce a sqrt(2)  
but that did not help AFAICT.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 2, 2023, 4:55pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/9 "2023-09-02T16:55:42Z")

</div>

I don’t know if you checked the L2 from the example DC but it is pretty close to  
the paper although their method appears slightly better its not much.  
I can’t get their method to work yet but you may want to look at this for comparison.  
See if this helps. I think though the error along x=y may be informative. I’m doing this  
to test mjmdatascope but there may be a similar way with existing FF code.  
Obvoously you can just delete that or I can put everything up on github if there is interest.

 ![mwgferr](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/ec81338f313c52610c167fa78da87746a53d7e68.jpeg)

```auto
load "msh3"
load "medit"
 load "mjm_dscope_freefem"

// mjm - I added the formulation from this, 
// and it is very close but slightly worse than the paper
// https://www.sciencedirect.com/science/article/pii/S0377042714002131#fd000050
// although I still can't get that to work. ZZ
// https://www.um.es/freefem/ff++/pmwiki.php?n=Main.DiscontinuousGalerkin
// Discontinous Galerlin Method
// based on paper from 
// Riviere, Beatrice; Wheeler, Mary F.; Girault, Vivette
// title: 
// A priori error estimates for finite element 
// methods based on discontinuous approximation spaces
// for elliptic problems.
// SIAM J. Numer. Anal. 39 (2001), no. 3, 902--931 (electronic).
// ---------------------------------
// Formulation given by Vivette Girault
// ------ 
// Author: F. Hecht , december 2003
// -------------------------------

macro dn(u) (N.x*dx(u)+N.y*dy(u) ) // def the normal derivative 

int nref = 6;
real [int] L2errorm ( nref ); // initialize the L2 error d array
for ( int nn = 0; nn < nref; nn++) {
int Nn = 2^(nn + 1) ; // space discretization
// "2" means 2 squares or a side of h but hmax has sqrt(2)
real h = 1.0/Nn;
mesh Th = square (Nn,Nn); // mesh generation of a square
//medit(""+Nn,Th);
//fespace Vh(Th,P2); // piecewise continuous quadratic finite element.
//fespace Vh(Th,P1); // Discontinous P2 finite element
fespace Vh(Th,P2dc); // Discontinous P2 finite element
Vh uh,vh;
Vh uex = sin(pi*x) * sin (pi*y); // exact solution
Vh f = 2.0 * sin(pi*x) *(pi^2) * sin (pi*y); // corresponding RHS
//
varf a(uh,vh, solver=sparsesolver,tgv=-1) = int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh) )  
// this seems to be important? for the continuous case anyway 
+on(0,1,2,3,4,uh=0)

 +intalledges(Th)( ( jump(vh)*mean(dn(uh)) - jump(uh)*mean(dn(vh))
// + pena*jump(u)*jump(v)
 ) / nTonEdge)

// + intalledges(Th)(sqrt(1)*(1.0/Nn)^-1 * jump(uh)*jump(vh) ); // bilinear form                            
// + intalledges(Th)((lenEdge)^-1 * jump(uh)*jump(vh) ); // bilinear form                            
;
matrix A = a(Vh, Vh,tgv=-1); // build the matrix
//set(A,solver=CG);
set(A,solver=sparsesolver);
varf L(uh, vh,tgv=-1) = int2d (Th)(f*vh)+ on(0,1,2,3,4,uh=0) ; // linear form
Vh F; F[] = L(0, Vh,tgv=-1); // build the right hand side vector
uh[] = A^-1 * F [];
int na=uh[].n;
//real s=0;
//for(int i=0; i<na; ++i) s=s+(uh[][i]-uex[][i])^2; 
//L2errorm [nn] = sqrt(s)*h*h; // sqrt ( int2d (Th)((uh - uex )^2));
L2errorm [nn] = sqrt ( int2d (Th)((uh - uex )^2));

meshL Th1=Sline(1024,[x]);
fespace Mh(Th1,P1);
Mh asol=uh(x,x);
Mh anex=uex(x,x);
Mh diff=anex-asol; // =uex(x,x);
mjmdscopeL(Th1,"names asol anex diff "+nn,asol,anex,diff); // warning do not forget () 

}
//
for ( int n = 0; n < nref; n ++){
int Nn = 2^(n + 1);
cout << " N " << " = " << Nn << endl;
cout << " L2errorm " << n << " = " << L2errorm [n] << endl;}
for ( int n = 1; n < nref; n ++) {
int Nn = 2^(n + 1);
cout << " N " << " = " << Nn << endl;
cout <<" convergence rate d = " << log ( L2errorm [n - 1]/ L2errorm [n])/log (2.) << endl; }

```

[MWGFEM\_1.edp](https://community.freefem.org/uploads/short-url/ghWoaytskYgqg6L2d0yv6PqVWWR.edp) (3.1 KB)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 2, 2023, 7:25pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/10 "2023-09-02T19:25:37Z")

</div>

Does this FF example look a lot like their eqn 2? You just need to get  
from 2 to 3 it seems 🙂

---

<div class="post-metadata">

**Author:** ![ahbabil1983](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ahbabil1983/32/1860_2.png) [@ahbabil1983](https://community.freefem.org/u/ahbabil1983)\
**Post date:** [September 6, 2023, 2:28pm UTC](https://community.freefem.org/t/modified-weak-galerkin-finite-element-method/2656/11 "2023-09-06T14:28:55Z")

</div>

I think we can move form eq. 2 to eq. 3 if u is continuous. But if u is discontinuous, this will be another problem.  
Moreover, there is a difference between the weak formulations for both the DGFEM and the MWGFEM. And we can consider MWGFEM as a special case of DGFEM if we consider the penalty parameter = 1 and neglect the second and third terms of the weak formulation for the DGFEM. But the jump and average functions used in both DGFEM and MWGFEM are not the same.  
Please, have a look to the code made by F. Hecht :  
[FreeFem++ users | Main \>\> DiscontinuousGalerkin](https://www.um.es/freefem/ff++/pmwiki.php?n=Main.DiscontinuousGalerkin).  
The jump function and the average function used there are different from the ones defined in the paper.  
Even though the solution obtained by the above code is close to the one in the paper, it is still not the same, because of the previous differences.