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!
