8

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
dcydhb
  • 615
  • 3
  • 9

1 Answers1

8

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • great answer and thanks a lot!!! – dcydhb Jul 24 '18 at 09: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. – Henrik Schumacher Jul 24 '18 at 09:12
  • hi,long time no see,can you explain some code "bmesh[Wireframe] mesh = ToElementMesh[bmesh]; mesh[Wireframe]" – dcydhb Aug 09 '18 at 09: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". – xzczd Aug 09 '18 at 10:49
  • @xzczd ok,thanks a lot! – dcydhb Aug 10 '18 at 01:33
  • and there is another question where you have used the boundary condition 8 in your solbetter – dcydhb Aug 10 '18 at 05:20
  • @dcydhb It's automatically satisfied, just as in sol, because we're using the same eq. – xzczd Aug 10 '18 at 07:28
  • @xzczd but have you ever define the bc 8 in sol? – dcydhb Aug 10 '18 at 08:01
  • @dcydhb As mentioned in my last comment, the 8th b.c. is satisfied automatically as long as eq is defined in the way above, this is a feature of "FiniteElement" method. Please check the link at the beginning of my answer for details. – xzczd Aug 10 '18 at 08:04
  • @xzczd but i have not found the 8th bc in your defined eq,so do you mean the 8th was in the fem method? – dcydhb Aug 11 '18 at 08:50
  • @dcydhb To some degree, yes. Once again, the 8th b.c. is automatically satisfied when 1. "FiniteElement" method is chosen, 2. the discontinuous coefficient is in proper position. Please check the linked post for details. – xzczd Aug 11 '18 at 09:15