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