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)/ (nx*ny*ny) /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) (2*wx)

+int2d(Th,1) (.200

*wx)*

+int2d(Th,2) (-.200wx)

+int2d(Th,2) (-.200

//+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”)