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$
load “msh3”
load “tetgen”
load “medit”
//
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;
macro grad(u) [dx(u),dy(u),dz(u)]// def of grad operator
// Solve
solve potential(psi, wx)
= int3d(Th)(
grad(psi)'grad(wx)) // scalar product
//+int2d(Thsrc) (2wx)
+int2d(Th,1) (.200wx)
+int2d(Th,2) (-.200wx)
//+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 ,
df<-read.table(“psi.txt”,header=F)
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”)