2
$\begingroup$

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: enter image description here

It has something to do with the If function. using a value instead of c works fine.

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

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}}];
eqn = Inactive[Div][ -c.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 1;
\[CapitalOmega] = Disk[];
d1 = DirichletCondition[u[x, y] == 0, x >= 0];
sol = NDSolveValue[{eqn, d1}, u, {x, y} \[Element] \[CapitalOmega]]

more Information can also be found in the FEM documentation - section "Partial Differential Equations with Variable Coefficients"

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.