# Magnetostatic Simulation of Dipole

**URL:** <https://community.freefem.org/t/magnetostatic-simulation-of-dipole/903>\
**Category:** General Discussion\
**Created:** [April 9, 2021, 6:22pm UTC](https://community.freefem.org/t/magnetostatic-simulation-of-dipole/903 "2021-04-09T18:22:44Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![bazzil](https://avatars.discourse-cdn.com/v4/letter/b/f1d935/32.png) [@bazzil](https://community.freefem.org/u/bazzil)\
**Post date:** [April 9, 2021, 6:22pm UTC](https://community.freefem.org/t/magnetostatic-simulation-of-dipole/903/1 "2021-04-09T18:22:44Z")

</div>

Hello Everyone. I have been learning how to use FreeFEM++ for the last few weeks to solve a research problem but first I am attempting to validate a model to gain confidence in my implementation. I have done lots of reading in the [documentation](https://doc.freefem.org/documentation/index.html) and [some examples](https://modules.freefem.org/modules/magnetostatic/). I have also read [many other threads](https://community.freefem.org/t/curl-in-cylindrical-coordinate/305) in this forum but I believe I am treading in uncharted waters here.

My problem can be viewed in detail [here](https://physics.stackexchange.com/questions/623234/solving-laplaces-equation-for-magnetostatic-problem-in-the-absence-of-currents?noredirect=1&lq=1). I did eventually solve this problem by identifying the correct boundary conditions, but I do not understand the FreeFEM implementation.

The shorthand is:

I am simulating a paramagnetic sphere in a static magnetic field which should distort the field in the shape of a dipole. I am first testing a 2D case before moving on to more complicated 3D configurations. In 2D, the magnetic field points in y, but I am working in z and x coordinates. I export the data to Matlab for further processing and analysis.

My main issues are as follows:  
-Why am I not able to take the partial derivative with respect to z? If I change all instances of dy to dz, the compiler throws an error. I need this because the definition of the magnetic field H in absence of currents is H=-grad(potential\_M)  
Despite this I was still able to obtain the expected dipole shape. Still, I find this odd.

-The boundary condition describing the continuity of the potential at the boundary seems to be defined using these two lines:

```auto
    + on(1, phi = z) // Dirichlet
    + on(2, phi = x) // Dirichlet

```

Why is this so? I found this totally by accident, but not sure I understand why it worked.

```auto
// 2D solution magnetic dipole
// Circular object in magnetic field (sphere)

include "D:\\FreeFem++\\examples\\ffpmatlib\\demos\\ffmatlib.idp"

// Reading in 2D mesh
// meshS Th=readmeshS("circle_in_square.mesh");

border a(t=0, 2*pi){x=4*cos(t); y=4*sin(t); label=1;};
border b(t=0, 2*pi){x=1*cos(t); y=1*sin(t); label=2;};
// plot(a(50) + b(30)); //to see a plot of the border mesh
mesh Th = buildmesh(a(50) + b(30));

// plot(Th);

//Parameters
real Mu0 = 4.*pi*1.e-7;	//Vacuum magnetic permeability
real MuC = 1.25e-6; //Iron magnetic permeability
real MuW = 1e-8; // Water permeability

real BoxWidth = 4.; //Box (square) size
real Radius = 1.41; //Radius of circular object

real B0 = 1.; // Tesla, magnetic field strength

//Fespace (Finite Element Space) in 2D
func Pk = P2;
fespace Ah(Th, Pk);
Ah phi, dphi; // Magnetic Potential scalar

func g = 1;

// Macro
macro grad(r) [dx(r),dy(r)] // Gradient function
macro div(r) dx(r) + dy(r) // Divergence function

// Defining equation parameters
// func f = B0; // RHS of Poisson's equation

int[int] lab = labels(Th);

func Mu1 = Mu0 + (MuW-Mu0)*(region==1);
func Nu1 = 1./Mu1;

func Mu2 = Mu0 + (MuC-Mu0)*(region==2); // circle
func Nu2 = 1./Mu2;

  cout << "labels:" << lab << " End labels" << endl;
  // 1 = outer, 2 = inner circle

// Problem -- Solve Laplace's equation for magnetic potential inside circle
solve magnetostatics(phi,dphi,solver=CG) // GMRES, CG, sparsesolver
	= int2d(Th)( grad(phi)' * grad(dphi) )-int1d(Th,1)( Mu1*Nu2*dphi ) // equation
    + on(1, phi = z) // Dirichlet
    + on(2, phi = x) // Dirichlet
	;

//Magnetic field H
Ah Hx, Hz;
Hx = -dx(phi);
Hz = -dy(phi);
Ah Htot = sqrt(Hx^2 + Hz^2);

plot(Hz, nbarrow=30, fill=true, value=true, cmm="Bz",wait=2);
plot(Hx, nbarrow=30, fill=true, value=true, cmm="Bx",wait=2);
plot(Htot, nbarrow=30, fill=true, value=true, cmm="Btot");

//Save a 2D vector field
ffSaveData2(Hx,Hz,"dipole_2D.txt");
// Saving Mesh
savemesh(Th,"circle_dipole_2D.msh");
// Saving Finite Element Space
ffSaveVh(Th,Ah,"dipole_2D_vh.txt");

```

Any insight would be appreciated. Thank you very much for your time and patience.

---

<div class="post-metadata">

**Author:** ![SamuelJosephs](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/samueljosephs/32/1046_2.png) [@SamuelJosephs](https://community.freefem.org/u/SamuelJosephs)\
**Post date:** [March 7, 2022, 7:02pm UTC](https://community.freefem.org/t/magnetostatic-simulation-of-dipole/903/2 "2022-03-07T19:02:43Z")

</div>

Where can I download the mesh?  
As for not being able to use z I would think that is due to the fact FreeFem assumes int2d is in x and y only, so the function dz() is not implemented for int2d.

To get round this I would either make my system fully 3D (As that is what maxwell’s equations assume) or I would just use y in place of z and pretend that it is instead z.