# The Edge03d element is used for numerical solution in freefem
**URL:** https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974
**Category:** General Discussion
**Created:** [June 18, 2025, 8:43am UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974 "2025-06-18T08:43:46Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@Dainy-Jia](https://community.freefem.org/u/Dainy-Jia)
#### Post date: [June 18, 2025, 8:43am UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/1 "2025-06-18T08:43:46Z")
I have a three-dimensional model：
E+u\times B-\curl B=g in \Omega, E \times n=0 on \partal \Omega
here, u represents the velocity field, E represents the electric field, and B represents the magnetic field.
The discrete format of this model is
(E\_h^n,F)+(u\_h^n \times B\_h^{n-1},F )-(\curl B\_h^{n},F)=(g,F),
In the code, I do the decoupling operation, which involves the exact solution with respect to the values of u and B. But unfortunately, I didn’t get the convergence order right.
load “Element\_Mixte3d”
load “msh3”
load “medit”
int MaxIter = 2;
real[int] errL2(MaxIter), orderL2(MaxIter);
macro C(t) 1./pi_cos(pi_t)-1.5_t //EOM
macro D(t) sin(pi_t)+1.5 //EOM
real cpu = clock();
real mu=1.;
real muc=2.;
real Rm=1.;
func real[int] uExact(real t){
real[int] uExactvec(3);
uExactvec(0) = (exp(-t)+1.)_sin(pi_x)_sin(pi_y)_sin(pi_z);
uExactvec(1) = 0.;
uExactvec(2) = (exp(-t)+1.)_(z^2-z)sin(pix)sin(piy);
return uExactvec;
};
func real[int] udxExact(real t){
real[int] udxExactvec(3);
udxExactvec(0) = (exp(-t)+1.)picos(pi_x)_sin(pi_y)_sin(pi_z);
udxExactvec(1) = 0.;
udxExactvec(2) = (exp(-t)+1.)_pi_cos(pi_x)_(z^2-z)_sin(pi_y);
return udxExactvec;
};
func real[int] udyExact(real t){
real[int] udyExactvec(3);
udyExactvec(0) = (exp(-t)+1.)_pi_cos(pi_y)sin(pix)sin(piz);
udyExactvec(1) = 0.;
udyExactvec(2) = (exp(-t)+1.)picos(pi_y)\*(z^2-z)_sin(pi_x);
return udyExactvec;
};
func real[int] udzExact(real t){
real[int] udzExactvec(3);
udzExactvec(0) = (exp(-t)+1.)_pi_cos(pi_z)sin(pix)sin(piy);
udzExactvec(1) = 0.;
udzExactvec(2) = (exp(-t)+1.)_(2.\*z-1.)_sin(pi_x)_sin(pi_y);
return udzExactvec;
};
func real[int] E1Exact(real t){
real[int] E1Exactvec(3);
E1Exactvec(0) = (sin(pi_t)+1.5)_(y^2-y)_sin(pi_z);
E1Exactvec(1) = (sin(pi_t)+1.5)_(x^2-x)_sin(pi_z);
E1Exactvec(2) = (sin(pi_t)+1.5)_(x^2-x)_sin(pi_y);
return E1Exactvec;
};
macro EExact(t) [(sin(pi_t)+1.5)_(y^2-y)_sin(pi_z), (sin(pi_t)+1.5)_(x^2-x)_sin(pi_z),(sin(pi_t)+1.5)_(x^2-x)_sin(pi_y)] //EOM
func real[int] EdxExact(real t){
real[int] EdxExactvec(3);
EdxExactvec(0) = 0.;
EdxExactvec(1) = (sin(pi_t)+1.5)_(2._x-1.)sin(piz);
EdxExactvec(2) = (sin(pi_t)+1.5)_(2.x-1.)sin(piy) ;
return EdxExactvec;
};
func real[int] EdyExact(real t){
real[int] EdyExactvec(3);
EdyExactvec(0) = (sin(pit)+1.5)_(2._y-1.)sin(piz);
EdyExactvec(1) = 0.;
EdyExactvec(2) = (sin(pi_t)+1.5)_(x^2-x)picos(pi_y);
return EdyExactvec;
};
func real[int] EdzExact(real t){
real[int] EdzExactvec(3);
EdzExactvec(0) = (sin(pi_t)+1.5)_(y^2-y)_pi_cos(pi_z);
EdzExactvec(1) = (sin(pi_t)+1.5)_(x^2-x)picos(pi_z);
EdzExactvec(2) = 0.;
return EdzExactvec;
};
func real[int] B1Exact(real t){
real[int] B1Exactvec(3);
B1Exactvec(0) = (1./pi_cos(pi_t)-1.5_t)pi(x^2-x)_(cos(pi_y)-cos(pi_z));
B1Exactvec(1) = (1./pi_cos(pi_t)-1.5_t)_(pi\*(y^2-y)_cos(pi_z)-(2._x-1.)sin(piy));
B1Exactvec(2) = (1./pi_cos(pi_t)-1.5_t)_(2.x-2.y)sin(piz);
return B1Exactvec;
};
macro BExact(t) [(1./picos(pit)-1.5_t)_pi_(x^2-x)_(cos(pi_y)-cos(pi_z)), (1./pi_cos(pi_t)-1.5_t)_(pi_(y^2-y)_cos(pi_z)-(2._x-1.)sin(piy)),(1./pi_cos(pi_t)-1.5_t)_(2.x-2.y)sin(piz)] // EOM
func real[int] BdxExact(real t){
real[int] BdxExactvec(3);
BdxExactvec(0) = (1./picos(pit)-1.5_t)_pi_(2._x-1.)_(cos(pi_y)-cos(pi_z));
BdxExactvec(1) = (1./pi_cos(pi_t)-1.5_t)_(-2.)_sin(pi_y);
BdxExactvec(2) = (1./pi_cos(pi_t)-1.5_t)_(2.)_sin(pi_z);
return BdxExactvec;
};
func real[int] BdyExact(real t){
real[int] BdyExactvec(3);
BdyExactvec(0) = (1./pi_cos(pi_t)-1.5_t)_(-pi^2)_(x^2-x)sin(piy);
BdyExactvec(1) = (1./pi_cos(pi_t)-1.5_t)_(pi_(2._y-1.)cos(piz)-pi_(2._x-1.)cos(piy));
BdyExactvec(2) = (1./pi_cos(pi_t)-1.5_t)\*(-2._sin(pi_z));
return BdyExactvec;
};
func real[int] BdzExact(real t){
real[int] BdzExactvec(3);
BdzExactvec(0) = (1./pi_cos(pi_t)-1.5_t)pi^2(x^2-x)sin(piz);
BdzExactvec(1) = (1./pi_cos(pi_t)-1.5_t)_(-pi^2)_(y^2-y)_sin(pi_z);
BdzExactvec(2) = (1./pi_cos(pi_t)-1.5\*t)_pi_(2.\*x-2.\*y)_cos(pi_z);
return BdzExactvec;
};
func real[int] BBcd(real t){
real[int] BBcdvec(3);
BBcdvec(0) = (1./pi_cos(pi_t)-1.5_t)pi(x^2-x)_(cos(pi_y)-cos(pi_z));
BBcdvec(1) = (1./pi_cos(pi_t)-1.5_t)_(pi\*(y^2-y)_cos(pi_z)-(2._x-1.)sin(piy));
BBcdvec(2) = (1./pi_cos(pi_t)-1.5_t)\*(2.\*x-2.\*y)_sin(pi_z);
return BBcdvec;
};
func real[int] EBcd(real t){
real[int] EBcdvec(3);
EBcdvec(0) = (sin(pi_t)+1.5)_(y^2-y)_sin(pi_z);
EBcdvec(1) = (sin(pi_t)+1.5)_(x^2-x)_sin(pi_z);
EBcdvec(2) = (sin(pi_t)+1.5)_(x^2-x)_sin(pi_y);
return EBcdvec;
};
func real[int] g(real t){
real[int] gvec(3);
gvec(0) = E1Exact(t)(0) + uExact(t)(1) \* B1Exact(t)(2) - uExact(t)(2) \* B1Exact(t)(1)
- 1. / Rm \* (BdyExact(t)(2) - BdzExact(t)(1));
gvec(1) = E1Exact(t)(1) + uExact(t)(2) \* B1Exact(t)(0) - uExact(t)(0) \* B1Exact(t)(2)
- 1. / Rm \* (BdzExact(t)(0) - BdxExact(t)(2));
gvec(2) = E1Exact(t)(2) + uExact(t)(0) \* B1Exact(t)(1) - uExact(t)(1) \* B1Exact(t)(0)
- 1. / Rm \* (BdxExact(t)(1) - BdyExact(t)(0));
return gvec;
};
//real cpu = clock();
// mesh iteration
for (int i = 0; i \< MaxIter; i++){
int n=pow(2, i + 2);
mesh3 Th = cube(n, n, n, [x, y, z]);
fespace Vh(Th, P1);
fespace Ph(Th, RT03d);
fespace Wh(Th, Edge03d);
Wh [E1, E2,E3],[F1,F2,F3];
Ph [B1,B2,B3], [B1old,B2old,B3old];
Vh u1old, u2old, u3old, u1,u2,u3;
```
int nt = pow(2*(i+2),2) / 4;
real tfinal = 1./4.;
real dt = tfinal / nt, alpha = 1. / dt;
real tnow = 0;
u1old = uExact(0.)(0);
u2old = uExact(0.)(1);
u3old = uExact(0.)(2);
[B1old,B2old,B3old] = BExact(0.);
problem solverE([E1,E2,E3],[F1,F2,F3], solver=CG) =
int3d(Th)(E1*F1+E2*F2+E3*F3)
+ int3d(Th)((u2*B3old-u3*B2old)*F1+(u3*B1old-u1*B3old)*F2+(u1*B2old-u2*B1old)*F3)
- int3d(Th)((1./Rm)*(F1*(dy(B3)-dz(B2))+F2*(dz(B1)-dx(B3))+F3*(dx(B2)-dy(B1))))
- int3d(Th)(g(tnow+dt)(0)*F1+g(tnow+dt)(1)*F2+g(tnow+dt)(2)*F3)
+ on(1, 2, 3, 4, 5, 6, E1=0.,E2=0.,E3=0.);
for (int i = 1; i <= nt; i++) {
real tnext = tnow + dt;
u1=uExact(tnext)(0);
u2=uExact(tnext)(1);
u3=uExact(tnext)(2);
[B1,B2,B3]=BExact(tnext);
solverE;
tnow = tnext;
[B1old,B2old,B3old]=BExact(tnow);
```
}
// E\_error
errL2[i] = sqrt(int3d(Th)((E1 - E1Exact(tfinal)(0)) ^ 2+ (E2 - E1Exact(tfinal)(1)) ^ 2+ (E3 - E1Exact(tfinal)(2)) ^ 2));
};
for (int i = 1; i \< MaxIter; i++){
orderL2[i] = log(errL2[i - 1] / errL2[i]) / log(2.);
};
cout \<\< “================OUTPUT===============” \<\< endl;
cout \<\< "CPU time = " \<\< clock() - cpu \<\< endl;
cout \<\< "errors and convergence rates: " \<\< endl;
cout \<\< "E = " \<\< endl;
cout \<\< “h\t L2 error\t order” \<\< endl;
for (int i = 0; i \< MaxIter; i++){
cout \<\< “1/” \<\< 2 ^ (i + 2) \<\< “:\t”
\<\< errL2[i] \<\< “\t” \<\< orderL2[i] \<\< “\t\t” \<\< endl;
};
---
### Author:  [@Monirul25](https://community.freefem.org/u/Monirul25)
#### Post date: [June 18, 2025, 12:09pm UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/2 "2025-06-18T12:09:05Z")
Hello, dear Dainy-Jia, Can you please provide your code as your\_code.edp format so that it can run directly.
---
### Author:  [@Dainy-Jia](https://community.freefem.org/u/Dainy-Jia)
#### Post date: [June 18, 2025, 12:41pm UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/3 "2025-06-18T12:41:51Z")
Sure. Thank you for your reply!
> Blockquote
> [solver\_E\_3D(1).edp](https://community.freefem.org/uploads/short-url/sBOezo6HtRJr9VokCGgtFS46IRY.edp) (6.8 KB)
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 20, 2025, 9:47pm UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/4 "2025-06-20T21:47:42Z")
Hello,
Your space Ph for the magnetic field chosen as RT03d seems inappropriate. I have taken P0 instead.
Then the curl B term needs to be integrated by parts, and the boundary condition E\times n=0 needs to be set correctly.
Another thing is that all macros need parenthesis () around the argument because they do literal replacement of the argument. When the argument is tnow+dt for example, without parenthesis it would not distribute the factor to each term tnow and dt.
[solver\_E\_3D.edp](https://community.freefem.org/uploads/short-url/tTTbiSkIj4OGqaTy65VZNH1DLca.edp) (7.1 KB)
---
### Author:  [@Dainy-Jia](https://community.freefem.org/u/Dainy-Jia)
#### Post date: [June 27, 2025, 4:00am UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/5 "2025-06-27T04:00:48Z")
Thank you very much for your reply！
Why does the magnetic field here use [P0,P0,P0] elements instead of RT03d elements?
In the model I presented above, these two equations are solved simultaneously. E belongs to the H(curl) space, and B belongs to the H(div) space.
The discrete format is as follows：
So naturally, the electric field (E\_h^n) is approximated using Edge03d(the lowest-order edge element), and the magnetic field (B\_h^n) is approximated using RT03d(the lowest-order face element).
By decoupling, I separately calculated the discrete format regarding the magnetic field as follows：
The code is
> Blockquote
> [slover\_B\_3D.edp](https://community.freefem.org/uploads/short-url/A1LuYoPqHsyjNsJqwXmmDPQPXXK.edp) (4.8 KB)
> The calculation result is:
> T=1/4; solver\_B\_3D
> ================OUTPUT===============
> CPU time = 645.957
> errors and convergence rates:
> B =
> h L2 error order divBh
> 1/4: 0.0448658 0 6.77796e-17
> 1/8: 0.0223948 1.00245 3.09493e-16
> 1/16: 0.0111921 1.00068 6.56293e-16
> This is a good result, but why does the order of convergence get wrong when we decouple the electric field? Looking forward to your reply！
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 27, 2025, 11:09am UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/6 "2025-06-27T11:09:44Z")
Dear Dainy-Jia,
I don’t know why RT03d does not work for B in the equation on E. Maybe one has to take finer meshes to see the convergence? Or maybe \nabla\times B, which is involved in the equation on E, is not well approximated when B in RT03d?
For the equation on B, RT03d works, and also P0, and Edge03d (I have tried with your code).
Thus I think that P0 for B will work for the coupled problem. Indeed with P0, the variational formulation on B is equivalent to setting directly
\frac{B^n\_h-B\_h^{n-1}}{\tau}+\nabla\times E\_h^n=0, because when E is in Edge03d, \nabla\times E is in P0. Thus you can replace B\_h^n in the equation on E by this expression and find a variational formulation for E\_h^n alone, involving (\nabla\times E\_h^n,\nabla\times F).
This formulation is good for E in Edge03d, and you have a single problem on E instead of a coupled problem in (E,B) to solve.
You could also do the converse: B in Edge03d and E in P0, and replace the value of E from the first equation into the second one (but then you have probably to make a correction with a “pressure type” Lagrange multiplier to get \nabla\cdot B=0 weakly, and a regularization term (\nabla\cdot B,\nabla\cdot C) ).
---
### Author:  [@Dainy-Jia](https://community.freefem.org/u/Dainy-Jia)
#### Post date: [June 27, 2025, 1:46pm UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/7 "2025-06-27T13:46:27Z")
Thank you for your reply and patient explanation. However, I still have one question. Why is the boundary condition chosen as:
+int2d(Th,1,3)(1.e30\*(E1_F1+E3_F3))
+int2d(Th,2,4)(1.e30\*(E2_F2+E3_F3))
+int2d(Th,5,6)(1.e30\*(E1_F1+E2_F2))
Looking forward to your reply！
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 27, 2025, 3:21pm UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/8 "2025-06-27T15:21:41Z")
This is for E\times n=0, knowing that
labels 1,3 correspond to n=(0,\pm 1,0)
labels 2,4 correspond to n=(\pm 1,0,0)
labels 5,6 correspond to n=(0,0,\pm1)
thus
for labels 1,3 we have to set E1=0 and E3=0,
for labels 2,4 we have to set E2=0 and E3=0,
for labels 5,6 we have to set E1=0 and E2=0
This is done by penalty.
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [June 30, 2025, 8:34am UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/9 "2025-06-30T08:34:04Z")
Please do not put 1e30 as a penalty term i integral because you have term off diagonal and it is true penalisation method,
1e6 is OK , I think. Remark on(1,3)(E1=0,E2=0,E3=0) impose E.tangent = 0 and tangent is (±1,0,0)
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 30, 2025, 4:37pm UTC](https://community.freefem.org/t/the-edge03d-element-is-used-for-numerical-solution-in-freefem/3974/10 "2025-06-30T16:37:58Z")
Dear Frédéric,
If I understand well, you say that because of the particular structure of the space `Edge03d`, setting `E1=0,E2=0,E3=0` indeed only imposes E.tangent=0.
This is good news, it simplifies the boundary statements!