Hi

i have to plot the Following function:

q(x/eps,y/eps,z/eps)=1 if (x/eps,y/eps,z/eps) \in B_{\delta}+ Z^3, else 0.

where B_{delta} is the ball with radius

\delta= ((3/4)*(1/(w*pi)))^{1/3}

we choise eps=0.1 and w=2000.

The problem is that in .vtu we Don’t ditingue the ball when q=1 to the rest when q=0.

I search in the doc but i Don’t find how we save .vtu in 3d (i choice 3d in paraview). So there is an solution?

And the image after copilation Don’t show the ball.

Thank you in advance to the help. There is the code:

’ ’ ’ freefem

load “UMFPACK64”

defaulttoUMFPACK64();

load “iovtk”

load “msh3”

verbosity=1;

func real q (real t, real s, real c,real w) {

real t1 = t - rint(t);

real s1 = s - rint(s);

real c1 = c -rint©;

real delta= ((3./4.)*(1./(pi*w)))^(1/3);

if ( square(s1)+square(t1)+square(c1) <= square(delta))

return 1;

else

return 0.0;

};

real eps=0.2, w=2000.;

mesh3 Th=cube(30,30,30);

fespace Wh(Th,P2,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);

Wh uplot;

int[int] Order= [1];

uplot=q(x/eps,y/eps,z/eps,w);

plot(uplot,wait=1,dim=3);

savevtk(“q.vtu”,Th,uplot,dataname=“qdelta”, order=Order);

Kind regards