# Why FEM not giving output at the boundary

I am trying to solve simple 3D Poisson's equation of thermal conduction due to heat generation (q) in a spherical composite, where from 0 to a (inner radius) conductivity is k1 and from a to b (outer radius) conductivity is k2. And BC is temperature becomes zero at r=b. I am using FEM to solve this problem, but the result of NDSolveValue is coming only till b=4.99 not b=5. Please help. I guess the output is also not correct.

  Needs["NDSolveFEM"]
k1=10;k2=1;a=1;q=1;b=5;
kk[x_, y_, z_] :=
Boole[x^2 + y^2 + z^2 < a^2] k1 + Boole[x^2 + y^2 + z^2 >= a^2] k2;
area = ImplicitRegion[x^2 + y^2 + z^2 - b^2 <= 0, {x, y, z}];
bc1 = DirichletCondition[t[x, y, z] == 0, x^2 + y^2 + z^2 - b^2 == 0];
op = kk[x, y, z] Laplacian[t[x, y, z], {x, y, z}] + q;
tfun = NDSolveValue[{op == 0, bc1}, t, {x, y, z} ∈ area];
Plot[tfun[x, 0, 0], {x, 0, 5}]

• tfun[x,y,z] is evaluable only, if the point x,y,z lies inside the grid(mesh) of your FEM-solution. Obviously point 5,0,0lies outside! Jun 21, 2020 at 21:13
• Clarification: why do you think output is not correct? Also variable "area" is misleading. Domain is a volume. Like the use of Boole. Jun 21, 2020 at 23:10
• Thank you for replying. @Ulrich Neumann Sorry I didn't understand why 5,0,0 lies outside the mesh. Jun 22, 2020 at 4:16
• Thank you for replying. @PaulCommentary I think output is not correct because when I am trying to solve by DSolve by assuming temperature and heat flux continuity at the interface, it is giving me different result. Jun 22, 2020 at 4:19

When a symbolic region is meshed there can be discrepancies between what the symbolic region represents and what the numerical approximation represent; that is especially the case for curved regions. "OpenCascade" typically does a good job at preserving the boundary. So if you can, it's advisable to represent the region with graphics primitives such that "OpenCascade" can work like so:

Needs["NDSolveFEM"]
k1 = 10; k2 = 1; a = 1; q = 1; b = 5;
kk[x_, y_, z_] :=
Boole[x^2 + y^2 + z^2 < a^2] k1 + Boole[x^2 + y^2 + z^2 >= a^2] k2;
area = Ball[{0, 0, 0}, b];
bc1 = DirichletCondition[t[x, y, z] == 0, True];
op = kk[x, y, z] Laplacian[t[x, y, z], {x, y, z}] + q;
mesh = ToElementMesh[area, "BoundaryMeshGenerator" -> "OpenCascade"];
tfun = NDSolveValue[{op == 0, bc1}, t, {x, y, z} \[Element] mesh];
(*Plot[tfun[x,0,0],{x,0,5}]*)
tfun[b, 0, 0]

(* 1.13485*10^-16 *)


Other options would be to refine the mesh. Such that the point {b,0,0} is actually part of the mesh or change the interpolating function property such that extrapolation will be performed.

• Thank you so much. Jun 22, 2020 at 14:36
• @Nilabh, you are welcome. Jun 22, 2020 at 19:02