# Region and label

In movemesh23, the code I copied had a “region” parameter but I also added
a “label” parameter that appears to work with the boundary conditions. Do these
conflict or are they different? The resulting code seems to work as I could include
neumann conditions on the X planes and verify in R the solution was mostly planar
although the slope at the one end was closer to .1 than .2 . Did I setup the
weak formulation Neumann conditions right? It looks like you get zero by default.

Also, when dumping the psi.txt, I used code that I copied from somewhere else
and it looks a bit odd and is probably dumping the same point multiple times.
Is the mesh indexing discussed somewhere?

This code is supposed to be two cubes ( cars ) in an airstream ( copied the code
from the airfoil example in the manual ). I was debating about trying a full Navier
Stokes with this example or moving on the Dirac or Landau-Lifschitz- Gilbert.
Has anyone encountered issues with weak formulations that may be in the examples?
Thanks, so far so good.

///home/ubuntu/dev/freefem/FreeFem-sources-master/examples/3d$vi Poisson-cube-ballon.edp //marchywka@happy:/home/ubuntu/dev/freefem/FreeFem-sources-master/examples/3d$

//
real volumetet; // use in tetg.
meshS Thx0,Thx1, Thy0, Thy1, Thz0, Thz1;
int nx=15,ny=15,nz=15;
func meshS mycube ( real x0, real x1, real y0,real y1, real z0, real z1, int orient, int lb)
{
meshS ThHex;
// first build the 6 faces of the hex.
//real x0=-1,x1=1;
//real y0=-1.1,y1=1.1;
//real z0=-1.2,z1=1.2;

// a volume of on tet.
volumetet= (x1-x0)(y1-y0)(z1-z0)/ (nxnyny) /6.;

mesh Thx = square(ny,nz,[y0+(y1-y0)*x,z0+(z1-z0)*y]);
mesh Thy = square(nx,nz,[x0+(x1-x0)*x,z0+(z1-z0)*y]);
mesh Thz = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

int[int] refz=[0,5]; // bas
int[int] refZ=[0,6]; // haut
int[int] refy=[0,3]; // devant
int[int] refY=[0,4]; // derriere
int[int] refx=[0,1]; // gauche
int[int] refX=[0,2]; // droite
int f1=-1;
int f2=1;
if (orient==1) { f1=1; f2=-1; }
int[int] l0=[0,lb+0];
//meshS
Thx0 = movemesh23(Thx,transfo=[x0,x,y],label=l0,orientation=f1,region=refx,removeduplicate=false);
++l0[1];
//meshS
Thx1 = movemesh23(Thx,transfo=[x1,x,y],label=l0,orientation=f2,region=refX,removeduplicate=false);
++l0[1];
//meshS
Thy0 = movemesh23(Thy,transfo=[x,y0,y],label=l0,orientation=f2,region=refy,removeduplicate=false);
++l0[1];
//meshS
Thy1 = movemesh23(Thy,transfo=[x,y1,y],label=l0,orientation=f1,region=refY,removeduplicate=false);
++l0[1];
//meshS
Thz0 = movemesh23(Thz,transfo=[x,y,z0],label=l0,orientation=f1,region=refz,removeduplicate=false);
++l0[1];
//meshS
Thz1 = movemesh23(Thz,transfo=[x,y,z1],label=l0,orientation=f2,region=refZ,removeduplicate=false);

cout<<“l0=”<<l0<<endl;

//medit(" — ", Thx0,Thx1,Thy0,Thy1,Thz0,Thz1);
ThHex = Thx0+Thx1+Thy0+Thy1+Thz0+Thz1;
return ThHex;
}

real x0=-1,x1=1; real yy0=-1.1,yy1=.1; real z0=-1.2,z1=.2;
real xa0=-.5; real ya0=-.8; real za0=-1;
real xb0=-.1; real yb0=-.8; real zb0=-1;
real sz=.25;

meshS ThHex = mycube(x0,x1,yy0,yy1,z0,z1,1,0);
real vol=volumetet;
meshS Thsrc=Thx0;
nx=10; ny=10; nz=10;
meshS Tha = mycube(xa0,xa0+sz,ya0,ya0+sz,za0,za0+sz,-1,10);
meshS Thb = mycube(xb0,xb0+sz,yb0,yb0+sz,zb0,zb0+sz,-1,20);
meshS Thc= Tha+Thb+ThHex;
//medit(" asdf ",Thc);

real[int] domaine = [-.9,-1,-1.1,.5,vol];
mesh3 Th = tetg(Thc,switch=“pqaAAYYQ”,nbofregions=1,regionlist=domaine);
// Tetrahelize the interior of the cube with tetgen
//medit(“tetg”,Th,wait=1);

fespace Xh(Th,P23d);
Xh psi, wx,dxpsi,dypsi;
// Solve
solve potential(psi, wx)
= int3d(Th)(
//+int2d(Thsrc) (2
wx)
+int2d(Th,1) (.200wx)
+int2d(Th,2) (-.200
wx)
//+on(Thl,psi=1)
//+on(102,psi=1)
+on(1,psi=0)
//+on(1,psi=0)
//+on(1,psi=1)
//+on(3,dx(psi)=0)
//+on(4,psi=1)
//+on(5,psi=1)
//+on(6,psi=1)
;

dxpsi=dx(psi);
dypsi=dy(psi);
{
ofstream ff(“psi.txt”);
for (int ix = 0; ix < Th.nv; ix++)
{
// this does produce dups, maybe confusing everything…
for (int j = 0; j < 3; j++) // foo[ix]+=Thx[ix][j]*Thx[ix][j];
//ff << Th[ix][j].x << " “<< Th[ix][j].y << " " << psi[Xh(ix,j)] << endl;
ff << Th[ix][j].x << " “<< Th[ix][j].y << " "
<<Th[ix][j].z<<” "
<< dxpsi[Xh(ix,j)]<<” “<<dypsi[Xh(ix,j)]<<” "<<psi[Xh(ix,j)] << endl;
//ff << Th[ix][0].x << " " << Th[ix][0].y << " " << psi[Xh(ix,0)] << “\n\n\n”;
}
} // scoping

//real[int] viso=[1.2,1.15,1.1,1.02,1.01];
//real[int] viso=[2,1.1,1.05,1.01,1,.999,.99,.9];
//real[int] viso=[1.7,1.1,1,.999,.99,.9];
real[int] viso=[0,.01,.05,.1,.2,.3];
plot(psi, viso=viso, value=true, wait=1);
//plot(psi, value=true, wait=1);
//real[int] xxx=[dx(psi),dy(psi)];
//plot(xxx, value=true, wait=1);

R-code ,

library(“rgl”)
x=which((df$V3> -.9)&(df$V3< -.8))
plot3d(df$V1,df$V2,df$V4) plot3d(df$V1,df$V2,df$V5,add=T, col=“red”)
plot3d(df$V1,df$V2,df$V6) plot3d(df$V1,df$V2,df$V6)
history()
?history
savehistory(file=“xxx.hist”)

This is probably behaving as described in the manual
I can probably figure out the variable dump format …

No big problem with the weak forms. But for LLG you’ll have hard times to figure out how to work with a vector field with constant norm. At least I did not find a convenient way to deal wit this…

angles? Thanks. There does seem to be a lot of literature on numerical
analysis of LLG. Since you seem familiar with it, curious what manifestations
there may be of spin inertia on linger time scales.

One other thing I was considering was modelling coagulation or maybe some
classical electrochemistry ( drift and diffusion and reactions etc ). Just offhand
I did not expect any issues, just a learning curve.

Well, I think angles is a complicated option because angle are compact variables. So this means you can have discontinuity of you angular variable due to the identification between zero and 2\pi

What people usually do is to work in Cartesian coordinates for the vector field. More precisely, the magnetization say \boldsymbol{m} that has value on S^2, that is \boldsymbol{m}\cdot\boldsymbol{m}=1 everywhere. In general, algorithms will not preserve the norm of \boldsymbol{m}, so typically after each update, you have to project back \boldsymbol{m} on S^2. Computation-wise, this is quite time consuming I think…

That said, whenever you drop the constraint of normalization, everything is fine

Thanks. I’ll have to look at it some more but does theta actually appear anywhere
or just cos and sin? Overall tracking rotations will likely come up somewhere
Is norm conservation fundamentally different from energy conservation?

It will probably be a while before I get there,. Right now still playing with navier
stokes, thinking about coagulation or some eletrochemistry or a simple
configuration of metal ( like LC tank or antenna ) and solving for vector potential
A and rho. Previously I had been playing with libmesh but this
seems simpler for starting out and I was also looking at
some front end code for FEM and this looks like a good example
of that if nothing else.
Thanks.