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, 2018 at 10:50
  • 3
    $\begingroup$ You might want to try using Inactive[Grad] and Inactive[Div] in your eqn. $\endgroup$
    – Silvia
    May 4, 2018 at 14:44
  • $\begingroup$ Thanks @Silvia it works fine! $\endgroup$
    – p4di
    May 7, 2018 at 8:11

1 Answer 1

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 and acknowledge you have read our privacy policy.

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