3
$\begingroup$

I have searched and read all previous questions but cannot get my head around this. I am new to mathematica. I have two regions in 2D where I want to solve PDE. Of the form:

Eqn=Laplacian[u[x,y], {x,y} - (alpha)u[x,y] == 0

Where alpha is equal to zero in one region (reduces to Laplace's eq) and equal to a constant (call it a so alpha=a). The regions are defined:

Omega = Rectangle[{100, 100}, {200, 101}];
Gamma = Rectangle[{0, 0}, {300, 200}];
Stigma = RegionDifference[Gamma, Omega];

So we have a narrow rectangle inside a much bigger rectangle. The ends of the inner rectange are equal to 1 and the outer borders are 0 as:

BCond1 = DirichletCondition[u[x, y] == 1., 
 x == 100 && 100 <= y <= 101 || 
 x == 200 && 100 <= y <= 101;
BCond2 = u[x, 0] == u[x, 200] == u[0, y] == u[300, y] == 0;

And my question, how do I solve for alpha=0 outside the inner box, i.e region Stigma and alpha=a=3.14(Example) inside the inner box, region Omega ? I have read loads of answers regarding this and it seems to stem on the Inactive func but can't for the life of me work out the notation or what exactly it is doing?

Any help greatly appreciated.

Improve this question
$\endgroup$
1
  • $\begingroup$ The ends of the inner rectangle are equal to 1 you only give B.C. for the left and right ends. What about the B.C. for top and bottom edges of the inner rectangle? $\endgroup$
    – Nasser
    Commented May 25, 2020 at 12:36

2 Answers 2

Reset to default
5
$\begingroup$

The following example under PeriodicBoundaryCondition should get you close. For example,

{length, height, xc, yc} = {300, 200, 0, 0};
{sx, sy, fx, fy} = {0, 0, length, height};
{srcxmin, srcxmax, srcymin, srcymax} = {100, 200, 100, 101};
{srcxmean, srcymean} = {(srcxmin + srcxmax)/2, (srcymin + srcymax)/2};
Ω = Rectangle[{sx, sy}, {fx, fy}];
source = If[srcxmin <= x <= srcxmax && srcymin <= y <= srcymax, 1., 
   0.];
op = ( Inactive[
      Div][({{-1, 0}, {0, -1}}.Inactive[Grad][u[x, y], {x, y}]), {x, 
      y}]) - source;
pde = op == 0;
dc = DirichletCondition[u[x, y] == 0, True];
ufun = NDSolveValue[{pde, dc}, u, {x, y} ∈ Ω];
ContourPlot[ufun[x, y], {x, y} ∈ Ω, 
 ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ Hello Tim, many thanks for your reply. Could you possibly explain in more detail the lines: op = ( Inactive[ Div][({{-1, 0}, {0, -1}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}]) - source;. More specifically what does the Inactive bit do? I've seen it many times and read the doc but I do not understand the syntax/what it is doing. $\endgroup$
    – BloodSexMagik
    Commented May 25, 2020 at 15:09
  • $\begingroup$ You don't actually use PeriodicBoundaryCondition here though?? $\endgroup$
    – chris
    Commented May 25, 2020 at 15:21
  • $\begingroup$ @chris Correct. I did not use the PeriodicBoundaryCondition here. I felt the non-BC parts (the operator and source terms) of the documented case were re-usable and that I would modify/experiment the boundary conditions. $\endgroup$
    – Tim Laska
    Commented May 25, 2020 at 16:19
  • $\begingroup$ @BloodSexMagik When possible, I try to cast my problem into coefficient form as described here. I find this format to be useful when comparing between other solvers. If you wanted to apply a Robin condition, you would need to convert to the Laplacian to coefficient form as described here. Inactive simplifies adjustment of c matrix Formal Differential Equations. $\endgroup$
    – Tim Laska
    Commented May 25, 2020 at 16:43
4
$\begingroup$

The documentation is a great place to get started, for example see here, here or here.

To get you started you can use:

Needs["NDSolve`FEM`"]
Omega = Rectangle[{100, 100}, {200, 101}];
Gam = Rectangle[{0, 0}, {300, 200}];
Stigma = RegionDifference[Gam, Omega];
(mesh = ToElementMesh[Stigma, "RegionHoles" -> None, 
    "RegionMarker" -> {{{150, 100.5}, 1, 0.25}, {{1, 1}, 2, 
       20}}])["Wireframe"]

enter image description here

mesh["Wireframe"[ElementMarker == 1]]

enter image description here

BCond1 = DirichletCondition[u[x, y] == 1., 
   x == 100 && 100 <= y <= 101 || x == 200 && 100 <= y <= 101];
BCond2 = u[x, 0] == u[x, 200] == u[0, y] == u[300, y] == 0;
if = NDSolveValue[{Laplacian[u[x, y], {x, y}] - 
     If[ElementMarker == 1, Pi, 0]*u[x, y] == 0, BCond1, BCond2}, 
  u, {x, y} \[Element] mesh]

![enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Hello, thank you for this. I am having trouble though. I have traced it back to an issue with RegionDifference. It seems that it is not computing the region correctly? If I plot the resulting region with RegionPlot[Stigma] then I get a 2 by 2 square, not at all what I expect! I have also tried in a seperate notebook with a variety of rectangles and get the same issue. Do you know of this? I am using 11.2. Thanks $\endgroup$
    – BloodSexMagik
    Commented May 25, 2020 at 16:19
  • $\begingroup$ @BloodSexMagik, you could try to make use of ToBoundaryMesh to get a region with the internal region. Have a look at the documentation. $\endgroup$
    – user21
    Commented May 25, 2020 at 20:02

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.