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 ) )
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
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.
[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
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.
There are some examples of calculations of the modified Galerkin FEM for the Poisson Problem in d = 2, 3.
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(pix) sin (piy); // exact solution
Vh f = 2sin (pix)pi^2 sin (piy); // corresponding RHS
varf a(uh ,vh) = int2d (Th)( Grad (uh)’ * Grad (vh) ) // bilinearform
- int2d (Th)(fvh) // 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)(fvh); // 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;}
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 (3.7 KB)
// 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;
}
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 , 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.xdx(u)+N.ydy(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(pix) sin (piy); // exact solution
Vh f = 2.0 sin (pix)(pi^2) * sin (piy); // corresponding RHS
// func fh = 2.0 sin (pix)(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))
+ penajump(uh)jump(vh) ) / nTonEdge);
matrix A = a(Vh, Vh); // build the matrix
varf L(unused, vh) = int2d (Th)(fvh); // 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; }
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.
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 (1.2 KB)
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.
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.
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 (3.1 KB)
Does this FF example look a lot like their eqn 2? You just need to get
from 2 to 3 it seems
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.
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.