8
$\begingroup$

enter image description here

enter image description here

I'm going to solve the Laplacian equation of the electrostatic field, which consists of two triangular regions, a rectangular region, a square, and on the intersection of the two regions of $$y=x$$, there are the first and second boundary conditions. How to set the correct Boundary condition and solve the problem

tried

Ω1 = DiscretizeRegion@Triangle[{{0, 0}, {0, 1}, {1, 1}}];
(*RegionPlot[Ω1]*)
Ω2 = DiscretizeRegion@Triangle[{{0, 0}, {1, 0}, {1, 1}}];
(*RegionPlot[Ω2]*)

nv1 = NeumannValue[0, x == 0];
nv2 = NeumannValue[0, x == 1];
nv3 = NeumannValue[0, x == y];
nv4 = NeumannValue[0, x == y];

sol1 = NDSolveValue[{D[u1[x, y], x, x] + D[u1[x, y], y, y] == 
    nv1 + nv3, 
   DirichletCondition[u1[x, y] == 10, y == 1 && 0 <= x <= 1]}, 
  u1, {x, y} ∈ Ω1]
sol2 = NDSolveValue[{D[u2[x, y], x, x] + D[u2[x, y], y, y] == 
    nv2 + nv4, 
   DirichletCondition[u2[x, y] == 0, y == 0 && 0 <= x <= 1]}, 
  u2, {x, y} ∈ Ω2]

DensityPlot[sol1[x, y], {x, y} ∈ Ω1, 
 Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]
DensityPlot[sol2[x, y], {x, y} ∈ Ω2, 
 Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]

but the right answer should be this enter image description here

Improve this question
$\endgroup$

1 Answer 1

Reset to default
8
$\begingroup$

Your translation for the 8th equation i.e. continuity condition is wrong. This issue has been discussed in detail here so I'd like to omit corresonding explanation and simply give the solution. Values of $\gamma_1$ and $\gamma_2$ aren't given in the question so I've picked them casually. 

gamma1 = 1; gamma2 = 2;
gamma = Piecewise[{{gamma1, y > x}}, gamma2];

With[{phi = phi[x, y]},
  eq = gamma Laplacian[phi, {x, y}] == 0;
  (* Alternatively: *)
  (* eq= Inactive[Div][{{gamma, 0}, {0, gamma}}.Inactive[Grad][phi,{x,y}],{x,y}] == 0; *)
  bc = {phi == 10 /. y -> 1, phi == 0 /. y -> 0};]


sol = NDSolveValue[{eq, bc}, phi, {x, 0, 1}, {y, 0, 1}];

ContourPlot[sol[x, y], {x, 0, 1}, {y, 0, 1}]~Show~Plot[x, {x, 0, 1}]

Mathematica graphics

The quality of solution above isn't that good actually, because NDSolve won't take the internal boundary at $y=x$ into consideration when discretizing the region automatically. To improve the quality, we can:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh["Coordinates" -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}], 
     LineElement[{{3, 1}}]}];
bmesh[Wireframe]
mesh = ToElementMesh[bmesh];
mesh[Wireframe]
solbetter = NDSolveValue[{eq, bc}, phi, {x, y} ∈ mesh];

Plot[{sol[x, 1/2], solbetter[x, 1/2]}, {x, 0, 1}]

Mathematica graphics

Improve this answer
$\endgroup$
11
  • $\begingroup$ great answer and thanks a lot!!! $\endgroup$
    – dcydhb
    Commented Jul 24, 2018 at 9:00
  • $\begingroup$ Very neat. As often, the problem lies in the fact that people want (i) to stick to strong formulation and (ii) to write the diffusion constants in front. $\endgroup$
    – Henrik Schumacher
    Commented Jul 24, 2018 at 9:12
  • $\begingroup$ hi,long time no see,can you explain some code "bmesh[Wireframe] mesh = ToElementMesh[bmesh]; mesh[Wireframe]" $\endgroup$
    – dcydhb
    Commented Aug 9, 2018 at 9:20
  • $\begingroup$ @dcydhb Please check the document of ToElementMesh and the tutorial FEMDocumentation/tutorial/ElementMeshCreation, especially the part starting from "A boundary element mesh may have internal structure; for example, to represent two material regions". $\endgroup$
    – xzczd
    Commented Aug 9, 2018 at 10:49
  • $\begingroup$ @xzczd ok,thanks a lot! $\endgroup$
    – dcydhb
    Commented Aug 10, 2018 at 1:33

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