# Problem with Internal Boundaries in NDSolve / FEM using the If Function

I'm still fairly new to Mathematica and am learning how to use NDSolve and FEM.

On youtube I found an interesting presentation done by Paritosh Mokhasi about finite elements that I'm trying to follow. When trying to implement an internal boundary shown in the presentation (~6:30 in the video), I run into an error. The code is working fine in Paritosh's presentation, is there a specific reason to why it isn't working for me:

c = If[x^2 + y^2 <= 1/4, {{10, 0}, {0, 10}}, {{1, 0}, {0, 1}}];
eqn = Div[-c.Grad[u[x, y], {x, y}], {x, y}] == 1;
Ω = Disk[];
d1 = DirichletCondition[u[x, y] == 0, x >= 0];
sol = NDSolveValue[{eqn, d1}, u, {x, y} ∈ Ω]


As the copied code is fairly messy, here's a screenshot: It has something to do with the If function. using a value instead of c works fine.

• The basic problem is that eqn evaluates before c evaluates. May 4 '18 at 10:50
• You might want to try using Inactive[Grad] and Inactive[Div] in your eqn. May 4 '18 at 14:44
• Thanks @Silvia it works fine!
– p4di
May 7 '18 at 8:11

Using Inactive[Grad] and Inactive[Div] in eqn fixes the Problem:
c = If[x^2 + y^2 <= 1/4, {{10, 0}, {0, 10}}, {{1, 0}, {0, 1}}];