# The rectangular region composed of two triangular regions contains a pde connecting the bc of the first and the second kind

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

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}]


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["NDSolveFEM"]
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}]


• great answer and thanks a lot!!! Jul 24, 2018 at 9:00
• 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. Jul 24, 2018 at 9:12
• hi,long time no see,can you explain some code "bmesh[Wireframe] mesh = ToElementMesh[bmesh]; mesh[Wireframe]" Aug 9, 2018 at 9:20
• @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". Aug 9, 2018 at 10:49
• @xzczd ok,thanks a lot! Aug 10, 2018 at 1:33