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"

// 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);

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
// Saving Mesh
// Saving Finite Element Space

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

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.