# Magnetostatic issue

**URL:** <https://community.freefem.org/t/magnetostatic-issue/2102>\
**Category:** General Discussion\
**Created:** [November 4, 2022, 11:16am UTC](https://community.freefem.org/t/magnetostatic-issue/2102 "2022-11-04T11:16:30Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![dega](https://avatars.discourse-cdn.com/v4/letter/d/a183cd/32.png) [@dega](https://community.freefem.org/u/dega)\
**Post date:** [November 4, 2022, 11:16am UTC](https://community.freefem.org/t/magnetostatic-issue/2102/1 "2022-11-04T11:16:30Z")

</div>

Hello,  
I am trying to make a simple magnetostatic simulation of a coil, I started from the example in:  
[https://modules.freefem.org/modules/magnetostatic/](https://modules.freefem.org/modules/magnetostatic/)

I created a simple gemometry with 4 coils and the results seem ok, when I add one more, the simulation seems to break

This is the code with the problem:

```auto
//Parameters
real Mu0 = 4.*pi*1.e-7;	//Vacuum magnetic permeability
real MuC = 1.25e-6; //Copper magnetic permeability
real J0 = 500.e4; //Current density

//Mesh
int nn = 25; //Mesh quality
real BoxRadius = 2.; //"Inifite" boundary

real copperRadius = 8e-3;
real CoilYSpacing = 0.018;
real CoilRadius = 0.015;

real pipeHeight =0.3;
real pipeThickness =0.01;
real pipeID =0.04;

int BoxWall = 10;

int nnL = max(2., 5*BoxRadius*nn/2.);
border b1(t=0., 2.*pi){x=BoxRadius*cos(t); y=BoxRadius*sin(t); label=BoxWall;};

border b2(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=-(CoilYSpacing*2.)+copperRadius*sin(t);};
border b3(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=-CoilYSpacing+copperRadius*sin(t);};
border b4(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=copperRadius*sin(t);};
border b5(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=CoilYSpacing+copperRadius*sin(t);};
border b6(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=(CoilYSpacing*2.)+copperRadius*sin(t);};

border b7(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=(-CoilYSpacing*2.)+copperRadius*sin(t);};
border b8(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=-CoilYSpacing+copperRadius*sin(t);};
border b9(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=copperRadius*sin(t);};
border b10(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=CoilYSpacing+copperRadius*sin(t);};
border b11(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=(CoilYSpacing*2.)+copperRadius*sin(t);};

int nnC = max(2., 900*copperRadius*nn/2.);

mesh Th = buildmesh(
	  b1(nnL) + b2(nnC) + b3(nnC) + b4(nnC) + b5(nnC) + b6(nnC) + b7(nnC) + b8(nnC) + b9(nnC) + b10(nnC) + b11(nnC)
);

plot(Th,wait=1);

//cout << b13;

int Coil1 = Th(CoilRadius, -2.*CoilYSpacing).region;
int Coil2 = Th(CoilRadius, -CoilYSpacing).region;
int Coil3 = Th(CoilRadius, 0.).region;
int Coil4 = Th(CoilRadius, CoilYSpacing).region;
int Coil5 = Th(CoilRadius, 2.*CoilYSpacing).region;

int Coil6 = Th(-CoilRadius, -2.*CoilYSpacing).region;
int Coil7 = Th(-CoilRadius, -CoilYSpacing).region;
int Coil8 = Th(-CoilRadius, 0.).region;
int Coil9 = Th(-CoilRadius, CoilYSpacing).region;
int Coil10 = Th(-CoilRadius, 2.*CoilYSpacing).region;

int Box = Th(0., 0.).region;

//Fespace
func Pk = P2;
fespace Ah(Th, Pk);
Ah Az;

//Macro
macro Curl(A) [dy(A), -dx(A)] //

//Current
Ah Jz = J0*(-1.*(Th(x,y).region == Coil1) - 1.*(Th(x,y).region == Coil2) - 1.*(Th(x,y).region == Coil3) - 1.*(Th(x,y).region == Coil4) - 1.*(Th(x,y).region == Coil5)
            +1.*(Th(x,y).region == Coil6) + 1.*(Th(x,y).region == Coil7) + 1.*(Th(x,y).region == Coil8) + 1.*(Th(x,y).region == Coil9) + 1.*(Th(x,y).region == Coil10));
plot(Jz, nbiso=30, fill=true, value=true, cmm="J - Current Density");

//Problem
func Mu = Mu0 + (MuC-Mu0)*((region==Coil1)+(region==Coil2)+(region==Coil3)+(region==Coil4)+(region==Coil5)+(region==Coil6)+(region==Coil7)+(region==Coil8)+(region==Coil9)+(region==Coil10));
						  
func Nu = 1./Mu;
varf vLaplacian (Az, AAz)
	= int2d(Th)(
		  Nu * Curl(Az)' * Curl(AAz)
	)
	- int2d(Th)(
		  Jz * AAz
	)
	+ on(BoxWall, Az=0)
	;
matrix<real> Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver);
real[int] LaplacianBoundary = vLaplacian(0, Ah);
Az[] = Laplacian^-1 * LaplacianBoundary;

//Magnetic induction
Ah Bx, By;
Bx = dy(Az);
By = -dx(Az);
Ah B = sqrt(Bx^2 + By^2);

//Magnetic field
Ah Hx, Hy;
Hx = Nu * Bx;
Hy = Nu * By;
Ah H = sqrt(Hx^2 + Hy^2);

//Current
Ah J = dx(Hy) - dy(Hx);

//Lorentz force
Ah Lx, Ly;
Lx = -J * By;
Ly = J * Bx;
Ah L = sqrt(Lx^2 + Ly^2);

//Plot
plot(Az, nbiso=30, fill=true, value=true, cmm="A");
plot(B, nbiso=30, fill=true, value=true, cmm="B - Magnetic induction");
plot(H, nbiso=30, fill=true, value=true, cmm="H - Magnetic field");
plot(J, nbiso=30, fill=true, value=true, cmm="J - Current");
plot(L, nbiso=30, fill=true, value=true, cmm="Lorentz force");
Ah MuVal = Mu;
plot(MuVal, nbiso=30, fill=true, value=true, cmm="Permeability");

```

---

<div class="post-metadata">

**Author:** ![dega](https://avatars.discourse-cdn.com/v4/letter/d/a183cd/32.png) [@dega](https://community.freefem.org/u/dega)\
**Post date:** [November 4, 2022, 11:19am UTC](https://community.freefem.org/t/magnetostatic-issue/2102/3 "2022-11-04T11:19:48Z")

</div>

For reference,  
this is the code with just 4 coils, that doesn’t seem to have the issue:

```auto
//Parameters
real Mu0 = 4.*pi*1.e-7;	//Vacuum magnetic permeability
real MuC = 1.25e-6; //Copper magnetic permeability
real J0 = 500.e4; //Current density

//Mesh
int nn = 25; //Mesh quality
real BoxRadius = 2.; //"Inifite" boundary

real copperRadius = 8e-3;
real CoilYSpacing = 0.018;
real CoilRadius = 0.015;

int BoxWall = 10;

int nnL = max(2., 5*BoxRadius*nn/2.);
border b1(t=0., 2.*pi){x=BoxRadius*cos(t); y=BoxRadius*sin(t); label=BoxWall;};

border b2(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=-CoilYSpacing+copperRadius*sin(t);};
border b3(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=copperRadius*sin(t);};
border b4(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=CoilYSpacing+copperRadius*sin(t);};
border b5(t=0., 2.*pi){x=CoilRadius+copperRadius*cos(t); y=(CoilYSpacing*2.)+copperRadius*sin(t);};

border b6(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=-CoilYSpacing+copperRadius*sin(t);};
border b7(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=copperRadius*sin(t);};
border b8(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=CoilYSpacing+copperRadius*sin(t);};
border b9(t=0., 2.*pi){x=-CoilRadius+copperRadius*cos(t); y=(CoilYSpacing*2.)+copperRadius*sin(t);};

int nnC = max(2., 900*copperRadius*nn/2.);

mesh Th = buildmesh(
	  b1(nnL) + b2(nnC) + b3(nnC) + b4(nnC) + b5(nnC) + b6(nnC) + b7(nnC) + b8(nnC) + b9(nnC)
);

plot(Th,wait=1);

//cout << b13;

int Coil1 = Th(CoilRadius, -CoilYSpacing).region;
int Coil2 = Th(CoilRadius, 0.).region;
int Coil3 = Th(CoilRadius, CoilYSpacing).region;
int Coil4 = Th(CoilRadius, 2.*CoilYSpacing).region;

int Coil5 = Th(-CoilRadius, -CoilYSpacing).region;
int Coil6 = Th(-CoilRadius, 0.).region;
int Coil7 = Th(-CoilRadius, CoilYSpacing).region;
int Coil8 = Th(-CoilRadius, 2.*CoilYSpacing).region;

int Box = Th(0., 0.).region;

//Fespace
func Pk = P2;
fespace Ah(Th, Pk);
Ah Az;

//Macro
macro Curl(A) [dy(A), -dx(A)] //

//Current
Ah Jz = J0*(-1.*(Th(x,y).region == Coil1) - 1.*(Th(x,y).region == Coil2) - 1.*(Th(x,y).region == Coil3) - 1.*(Th(x,y).region == Coil4)
            +1.*(Th(x,y).region == Coil5) + 1.*(Th(x,y).region == Coil6) + 1.*(Th(x,y).region == Coil7) + 1.*(Th(x,y).region == Coil8));
plot(Jz, nbiso=30, fill=true, value=true, cmm="J - Current Density");

//Problem
func Mu = Mu0 + (MuC-Mu0)*((region==Coil1)+(region==Coil2)+(region==Coil3)+(region==Coil4)+(region==Coil5)+(region==Coil6)+(region==Coil7)+(region==Coil8));
						  
func Nu = 1./Mu;
varf vLaplacian (Az, AAz)
	= int2d(Th)(
		  Nu * Curl(Az)' * Curl(AAz)
	)
	- int2d(Th)(
		  Jz * AAz
	)
	+ on(BoxWall, Az=0)
	;
matrix<real> Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver);
real[int] LaplacianBoundary = vLaplacian(0, Ah);
Az[] = Laplacian^-1 * LaplacianBoundary;

//Magnetic induction
Ah Bx, By;
Bx = dy(Az);
By = -dx(Az);
Ah B = sqrt(Bx^2 + By^2);

//Magnetic field
Ah Hx, Hy;
Hx = Nu * Bx;
Hy = Nu * By;
Ah H = sqrt(Hx^2 + Hy^2);

//Current
Ah J = dx(Hy) - dy(Hx);

//Lorentz force
Ah Lx, Ly;
Lx = -J * By;
Ly = J * Bx;
Ah L = sqrt(Lx^2 + Ly^2);

//Plot
plot(Az, nbiso=30, fill=true, value=true, cmm="A");
plot(B, nbiso=30, fill=true, value=true, cmm="B - Magnetic induction");
plot(H, nbiso=30, fill=true, value=true, cmm="H - Magnetic field");
plot(J, nbiso=30, fill=true, value=true, cmm="J - Current");
plot(L, nbiso=30, fill=true, value=true, cmm="Lorentz force");
Ah MuVal = Mu;
plot(MuVal, nbiso=30, fill=true, value=true, cmm="Permeability");

```

---

<div class="post-metadata">

**Author:** ![dega](https://avatars.discourse-cdn.com/v4/letter/d/a183cd/32.png) [@dega](https://community.freefem.org/u/dega)\
**Post date:** [November 4, 2022, 1:11pm UTC](https://community.freefem.org/t/magnetostatic-issue/2102/4 "2022-11-04T13:11:32Z")

</div>

A screenshot to show the odd solution for the 5 coil:

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

---

<div class="post-metadata">

**Author:** ![dega](https://avatars.discourse-cdn.com/v4/letter/d/a183cd/32.png) [@dega](https://community.freefem.org/u/dega)\
**Post date:** [November 4, 2022, 1:12pm UTC](https://community.freefem.org/t/magnetostatic-issue/2102/5 "2022-11-04T13:12:10Z")

</div>

The 4 coil one seems to have computed ok:

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

---

<div class="post-metadata">

**Author:** ![gero](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@gero](https://community.freefem.org/u/gero)\
**Post date:** [November 5, 2022, 9:36pm UTC](https://community.freefem.org/t/magnetostatic-issue/2102/6 "2022-11-05T21:36:08Z")

</div>

I don’t know if this helps:

I could not figure out what is wrong with your code… so I reworked the script to automatically generate the mesh for multiple windings (via nWinding; using array of border {see documentation}) and to reassign region numbers (as Frederic Hecht suggested in a recent answer to a different issue).

It seems to work for nWinding 3…7. At least the results look symmetric.

[revN5.edp](https://community.freefem.org/uploads/short-url/hpli9fI12fM7slqVTRveW2qk4KI.edp) (2.9 KB)

Perhaps it helps you or somebody else to figure out what goes wrong with your scripts.

Kind regards  
Gerald

---

<div class="post-metadata">

**Author:** ![dega](https://avatars.discourse-cdn.com/v4/letter/d/a183cd/32.png) [@dega](https://community.freefem.org/u/dega)\
**Post date:** [November 6, 2022, 10:33am UTC](https://community.freefem.org/t/magnetostatic-issue/2102/7 "2022-11-06T10:33:49Z")

</div>

Hi Gero,  
I like your code, way neater than my raw one.

Maybe it has somethign to do with region names, a small bug lurking somewhere…

Thank you very much!  
Marco

---

<div class="post-metadata">

**Author:** ![gero](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@gero](https://community.freefem.org/u/gero)\
**Post date:** [November 8, 2022, 12:44pm UTC](https://community.freefem.org/t/magnetostatic-issue/2102/8 "2022-11-08T12:44:02Z")

</div>

Just to conclude the matter:

The original problem is caused by a poor choice of the value of BoxWall  
int BoxWall = 10;  
as for higher numbers of windings 10 is also used for an automatically assigned label.

Thus the boundary condition is applied incorrectly.

Lesson learned:  
Use reasonably high numbers for manual assigned labels…  
so e.g.  
int BoxWall = 99;  
fixes the problem…

All the best  
Gerald

---

<div class="post-metadata">

**Author:** ![dega](https://avatars.discourse-cdn.com/v4/letter/d/a183cd/32.png) [@dega](https://community.freefem.org/u/dega)\
**Post date:** [November 8, 2022, 1:31pm UTC](https://community.freefem.org/t/magnetostatic-issue/2102/9 "2022-11-08T13:31:41Z")

</div>

Ah I see…  
Defnitely my mistake and no bug then!

Thanks again,  
Marco