# Magnetostatic Simulation of Dipole

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 and some examples. I have also read many other threads in this forum but I believe I am treading in uncharted waters here.

My problem can be viewed in detail here. 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:

``````    + 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.

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

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

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 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 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
+ 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.