I have been trying to learn FreeFem++. I have a code that shows no error. However, the code is not working.
The domain is $$ \Omega=[-5,5]\times[0,10]. $$
And the problem I am trying to solve is
$$ \Delta p=0\,\text{on } \Omega, $$ $$ p=y\,\text{ if $y=10$}, $$ $$ \frac{\partial p}{\partial n}=0\,\text{ if }y=0, x=-5 \text{ or }x=5. $$
Notice that the solution is $p\equiv10$.
My code reads
///////////////Code///////////////
border Gb(t=-5,5) {x=t;y=0;label=2;};//flat bottom
border Hr(t=0,10) {x=5;y=t;label=3;}; //boundary horizontal variable
border Gr(t=5,pi) {x=t;y=10;label=4;};//interface
border G(t=-pi,pi) {x=-t;y=10;label=5;};//interface
border Gl(t=-pi,-5) {x=t;y=10;label=4;};//interface
border Hl(t=0,10) {x=-5;y=10-t;label=1;};//boundary horizontal variable
func z=y;
mesh Th=buildmesh(Gl(25)+Gr(25)+G(50)+Gb(100)+Hr(200)+Hl(200));
fespace Ph(Th,P1); //impermeable BC on x
Ph p,phi;
solve Laplace(p,phi)=int2d(Th)(dx(phi)*dx(p)+dy(phi)*dy(p))-int1d(Th,Gb)(0*phi)-int1d(Th,Hr)(0*phi)-int1d(Th,Hl)(0*phi)+on(G,p=z)+on(Gr,p=z)+on(Gl,p=z); //weak formulation
plot(p,wait=1,value=true,fill=true);
/////////////End of Code/////////////
Then FreeFem++ returns the following message
-- mesh: Nb of Triangles = 37728, Nb of Vertices 19165 -- Solve : min 1.49772e-31 max 10 times: compile 0.186s, execution 6.051s, mpirank:0 CodeAlloc : nb ptr 2727, size :165264 mpirank: 0 Ok: Normal End
together with a plot. The problem is that the plot shows a nonconstant solution!! Any comment is very welcome.
The problem lies in the way you write the boundary conditions. First of all what do the pieces of the boundary with parameter $\pi$ stand for? Note that you put some boundary conditions on them also, which have nothing to do with the problem...
Since you want Neumann boundary conditions on all the sides except the upper one, you do not need to put any boundary conditions on them. The default boundary condition is the Neumann one (this is a direct consequence of the weak formulation). You only need to put the Dirichlet boundary condition on the upper side of the square which is
on(G,p=z)(remove all other boundary conditions...). With this modification I do obtain a constant solution.