A FreeFEM++ code uses the Modified Weak Galerkin Finite Element Method to solve the Poisson Problem

General Discussion
1

Hello everyone,
I wrote a FreeFEM++ code to solve the Poisson Problem -∆u=f in Ω u=g on ∂Ω using the Modified Weak Galerkin Finite Element Method based on 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 3.

I modified a FreeFEM++ code for the Discontinuous Galerkin FEM to solve this problem to fit the Modified Weak Galerkin FEM, but the results are not the same as in the paper. I am wondering what’s wrong with my code.
MWGFEM.edp (1.4 KB)

1 Like
2

hello!! Can your code run with a crrect results now? I have the same problem with the MWG method on this paper >.<

4

The weak gradient was not defined correctly in this code (it must not be the usual gradient). Here is a code with the weak gradient
mwgfem.edp (1.4 KB)
However still it does not give the numerical values listed in the paper (Table 1). The results in Table 1 are doubtful because the error cannot be so small (even the L2 projection of uex gives an error to uex larger that stated!).

5

Respected sir, I have not so much idea about weak gradient but results coming if take a penaltity parameter 1000 in the intalledges term of the formulation using P1 element. I don’t know why it is happening. But i have not found any thing about penaltity parameter in the paper.
Here, is the Code>
My Code_wgm.edp (2.2 KB)

6

The weak gradient is defined in the paper mentioned above.
If we put a large penalty parameter I agree that the method becomes second order, indeed it approaches the conforming P1 framework, which is known.
Here the paper claims that with P1 discontinuous it is better than classical, with the penalty parameter 1/h. We don’t get it.

7

Yes Sir. I have also observed that. Thank you.

8

Respected sir, I wrote a freefem code with the weak gradient, and the error can be small enough like the paper mentioned.
MWG.edp (2.4 KB)
But When I try to change the exact solution’s boundary condition
u=0
to
u = u_{ex}, \partial_n u = (u_x, u_y)\cdot n ,
it doesn’t work. I don’t know why. TT
MWG_2.edp (2.6 KB)

9

Hello Dear orange, your code is giving very good results. But which formulation you used in the code??. Why you are defining two spaces??.
Do you used the formulation in the following:

image
image1207×382 60.8 KB

(equation: 3)
I have not understand your problem formulation: you have not taken (dxx(sigma1)+dyy(sigma2)) term still you have used others.

Blockquote
problem WGM([sigma1, sigma2, uh], [tau1, tau2, vh], solver=UMFPACK)
= -int2d(Th)( vh*(dx(sigma1)+dy(sigma2)) )
+intalledges(Th)( (1-nTonEdge)( mean(vh)(jump(sigma1)N.x + jump(sigma2)N.y) )/nTonEdge )
+intalledges(Th)((jump(uh)jump(vh))/lenEdge)
+int2d(Th)(sigma1tau1 + sigma2tau2)
+int2d(Th)( uh(dx(tau1)+dy(tau2)) )
-intalledges(Th)( (1-nTonEdge)( mean(uh)(jump(tau1)N.x + jump(tau2)N.y) )/nTonEdge )
//+int2d(Th)(ghvh)
-int2d(Th)(fvh)// linear form
+on(1,2,3,4,uh=0)
Blockquote

Does you taken the formulation of the following paper or anything else:
J Comp and app math 2014.pdf (377.6 KB)

Specifically, why you taken the following term:

Blockquote
+int2d(Th)(sigma1tau1 + sigma2tau2)
+int2d(Th)( uh*(dx(tau1)+dy(tau2)) )
-intalledges(Th)( (1-nTonEdge)( mean(uh)(jump(tau1)*N.x + jump(tau2)N.y) )/nTonEdge )
//+int2d(Th)(ghvh)
Blockquote

I have not found like any term in the above paper. Can you share your formulation pls??

Also, I have observed :

Blockquote
fespace Vh(Th,P2dc); // Discontinous P2 finite element
Vh uh,vh;

Blockquote
fespace Wh(Th,P1dc);
Wh sigma1,sigma2;
Wh tau1, tau2;

This two spaces Vh and Wh should be such that Wh has one less polynomial degree. Why you taken like this ??.
I need your clarification brother.
In the paper still they define with same degree>

image
image922×160 18.9 KB

Regarding boundary condition: +on(1,2,3,4,uh=0) . You have to take it as you used uex=sin(pi*x)sin(piy) which is zero on boundary of the domain D=[0,1]^2.

I have checked your code with P1, P2, and P3. All gives good results.

Thanks in advance.

10

Indeed my code works also if I replace +on(1,2,3,4,uh=0) by +int1d(Th)(1.e10*uh*vh). It seems that there is a bug in the implementation of on() on P1dc.

11

It is not a bug it is due to the logical support of the degree of freedom ( interior of the triangle not on boundary for discontinuity),.

12

Respected sir, Does the problem formulation in the following code is fine?. I am not getting the problem formulation in codes.
MWG.edp (2.4 KB)

If i understand well, i have to solve the following equation:

image
image1207×382 60.8 KB
(eqn. 3)
Does the code following it sir?.

Can please you explain it sir??. It will help to implement it on my works.

Thanks in advance sir.

13

The code is correct.
Apply Definition 1.1 to v\in V_h^0, taking q=\nabla_w u gives
(\nabla_w v,\nabla_w u)_T=-(v,\nabla\cdot(\nabla_w u))_T+(\{v\},\nabla_w u\cdot n)_{\partial T}
Replacing this in eq. (3) gives the variational formulation
-\sum_T\int_Tv\nabla\cdot(\nabla_w u) +\sum_T\int_{\partial T}\{v\}(\nabla_w u)\cdot n +\sum_e h^{-1}\int_e [[u]][[v]]-\int_\Omega fv=0
for all v\in V_h^0.
You write then Definition 1.1 for u, and call \nabla_wu=(\sigma_1,\sigma_2). Renaming the test function q as q=(\tau_1,\tau_2) you finally get the variational formulation of the code, with test functions v,\tau_1,\tau_2.
You may prefer this version
mwgfem.edp (1.5 KB)
which is exactly the above formulas ((\sigma_1,\sigma_2) is replaced by (Gx,Gy)).

14

Thanks for quick response sir. So, can i apply this code ideas in any time dependent problem in modified weak Galerkin method and weak galerkin method sir??.

15

A priori yes (but the proof of stability, convergence and error estimate may be missing).

16

So, why the proof of stability, convergence and error estimate may be missing sir??. Numerically, we are getting optimal results. So, theoritically should come.

17

It should, but as long as no proof is published, you are not sure.

19

As i am understand from your ideas sir,
i am getting two set of equations , but i am not getting from where the term - int2d(Th)(Gxqx+Gyqy) comes. Can tell sir??. If i am understand well, we do not add after getting two set of equations otherwise the terms + intalledges(Th)( jump(uh)*jump(vh)/h/nTonEdge ) and

  • int2d(Th)(f*vh) will comes twice. Can you tell me please sir??.
21

My understanding: (q, (Gx,Gy)) giving the term - int2d(Th)(Gxqx+Gyqy). Am i right sir??. But this should be with “+” sign why “-” is in formulation sir??. (\tau_1=qx, \tau_2=qy in the code )
Also, in this way i think we are adding many extra terms in formulation sir.