# Toy Maxwell equation code

**URL:** https://community.freefem.org/t/toy-maxwell-equation-code/2937
**Category:** General Discussion
**Created:** [February 7, 2024, 11:12am UTC](https://community.freefem.org/t/toy-maxwell-equation-code/2937 "2024-02-07T11:12:58Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![Alexei](https://avatars.discourse-cdn.com/v4/letter/a/c77e96/32.png) [@Alexei](https://community.freefem.org/u/Alexei)
#### Post date: [February 7, 2024, 11:12am UTC](https://community.freefem.org/t/toy-maxwell-equation-code/2937/1 "2024-02-07T11:12:58Z")

</div>

Hi, I made a FreeFem code to test a toy 2D MaxWell equation model which is proposed as a benchmark in the presentation by Monique Dauge & Martin Costabel  
( [https://perso.univ-rennes1.fr/monique.dauge/publis/MaxwellMD.pdf](https://perso.univ-rennes1.fr/monique.dauge/publis/MaxwellMD.pdf) )

Correct eigenvalues seem independent of the div . div regularization term s while spurious eigenvalues disperse linearly - which is the expected behavior from the presentation.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c75e261b361010133296bda938bb511212500d73.png)

The FreeFem++ script is listed below.

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

int Left=99;  
int Right=98;  
int Top=96;  
int Bottom=97;

border D1(t=0., pi){x=t; y=0; label=Bottom;}  
border D2(t=0, pi){x=pi; y=t; label=Right;}  
border D3(t=pi, 0){x=t; y=pi; label=Top;}  
border D4(t=pi, 0){x=0; y=t; label=Left;}

int n = 5;  
mesh Th=buildmesh( D1(10_n) + D2(10_n) + D3(10_n) + D4(10_n) );

fespace Vh(Th, [P2, P2]); // change to P2 if convergence problems

Vh [Ex, Ey];  
Vh [Ux, Uy];  
real sigma = 1.;  
real s;

ofstream fout(“evalout.dat”);

for (s = 0; s \<= 4.01; s += 0.05) {

varf TE([Ex, Ey], [Ux, Uy])= int2d(Th)( (dx(Uy) - dy(Ux)) \* (dx(Ey) - dy(Ex) ) )

- int2d(Th)( s \* (dx(Ux) + dy(Uy)) \* (dx(Ex) + dy(Ey) ) )
  - int2d(Th)( sigma \* ( Ex \* Ux + Ey \* Uy ) )  
+on(Left,Ey=0) + on(Right, Ey=0) + on(Top, Ex=0) + on(Bottom, Ex = 0)  
;

varf TEid([Ex, Ey], [Ux, Uy])=int2d(Th)( Ex \* Ux + Ey \* Uy);

int nev = 30;  
real[int] ev(nev);

matrix A = TE(Vh, Vh,solver=UMFPACK);  
matrix B = TEid(Vh, Vh, solver=UMFPACK);

// Vh[int] eVx, eVy;  
// int k = EigenValue(A, B, sym=true, sigma=sigma, value=ev, vector=eVx, tol=1e-10, maxit=3, ncv=0);

int k = EigenValue(A, B, sym=true, sigma=sigma, value=ev, tol=1e-10);

for (int i = 0; i \< k; i++){  
fout \<\< s \<\< " " \<\< i \<\< " " \<\< ev[i] \<\< endl;  
}  
fout \<\< endl;  
fout.flush;

//plot( eVx[25], fill=true, value=true, wait= 1);  
//plot( eVy[25], fill=true, value=true, wait= 1);

}