Hello FF!

I’m trying to distill from the example Laplace3d.edp the simplest example for applying mixed Neumann-Dirichlet boundary conditions to a general tetrahedral mesh, and to solve \Delta u = 0 in the interior (steady heat equation). There are regions 1, 5 for Dirichlet boundary conditions, and 2, 3, 4 for Neumann boundary conditions (it describes a beam whose temperatures at the ends are fixed). To make things simple, I set the function ue of the example Laplace3d.edp to 0 for the Neumann boundary part, basically dictating zero flow through the Neumann boundary. (I also checked with nonzero ue, it made the solution more regular but still wrong so I concluded that setting ue=0 does not seem to be the mistake). If it is of interest, the surface mesh used in the following script is here.

load “msh3”

load “gmsh”

load “tetgen”

meshS ThS = gmshloadS(“mainS.msh”); // a surface mesh file in msh file format

func newregion = 1 + (y > -500.0) + (y > -10.0) + (y > 10.0) + (y > 500.0); // define regions

ThS = change(ThS, fregion=newregion);

mesh3 Th=tetg(ThS,switch=“paAAQYY”);

fespace Vh(Th, P2);

Vh u, v;

macro Grad(u) [dx(u), dy(u), dz(u)] // EOM

problem Laplace(u, v, solver=CG)

= int3d(Th)(Grad(u)’ * Grad(v))

+int2d(Th, 2)(u * v)

+int2d(Th, 3)(u * v)

+int2d(Th, 4)(u * v)

+ on(1, u=20.0)

+ on(5, u=-13.1);

Laplace;

but the solution is incorrect, it seems the zero Neumann boundary conditions are not applied correctly:

Help appreciated!