I want to solve the Laplace equation in 2 different medium in 2D, each with its own dielectric constant. This is the extremely simplified geometry. The problem is the same error (which you will find out in the end!) for my more complex geometry in 3D.

I start with defining the region:

Needs["NDSolve`FEM`"]  bmesh = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {.5, 0.}, {1., 0.},  {1., 1.}, {.5, 1.}, {0., 1.}},  "BoundaryElements" -> {LineElement[{{1,2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 1}, {2, 5}}]}];  bmesh["Wireframe"]  mesh = ToElementMesh[bmesh]; mesh["Wireframe"] 

The result is:

You can see there is a boundary at x==0.5 and the mesh of the region is made considering that there is a boundary at x==0.5.

My attempt to solve Laplace equation with only Dirichlet boundary condition is as follows, note that I just modified the wolfram’s FEM tutorial example for my problem:

\[Epsilon]r = If[x <= .5, {{1., 0.}, {0., 1.}}, {{10., 0.}, {0., 10.}}] op = Inactive[Div][-\[Epsilon]r.Inactive[Grad][u[x, y], {x, y}], {x,y}] \[CapitalGamma]D = {DirichletCondition[u[x, y] == 1., x == 0.], DirichletCondition[u[x, y] == 2., x == 1.],  DirichletCondition[u[x, y] == .1, 0. <= x <= .5 && y == 0.], DirichletCondition[u[x, y] == 0., 0.5 <= x <= 1. && y == 0.], DirichletCondition[u[x, y] == 0., 0. <= x <= .5 && y == 1.], DirichletCondition[u[x, y] == .1, .5 <= x <= 1 && y == 1.]}  \[CapitalGamma]N = NeumannValue[2, .45 < x < 0.55 ]; (*This is  simplified condition. I will explain the B.C. I have in mind at the  end.*)  ufun = NDSolveValue[{op == 0, \[CapitalGamma]D}, u, {x, y} \[Element] mesh] 

The result without using the Neumann condition is:

Show[ContourPlot[ufun[x, y], {x, y} \[Element] mesh,  ColorFunction -> "TemperatureMap"], bmesh["Wireframe"]] 

However, when I use Neumann boundary condition on the boundary at x==0.5 this happens:

ufun = NDSolveValue[{op == \[CapitalGamma]N, \[CapitalGamma]D}, u,  {x,y} \[Element] mesh]  (*NDSolveValue::bcnop: No places were found on the boundary where  x==0.5 was True, so NeumannValue[2,x==0.5] will effectively be  ignored.*) 

Please note that the Neumann condition that I have in my mind for the boundary at x==0.5 is

\[Epsilon]r(*for x < 0.5*)*D[u[x,y],x] == \[Epsilon]r(*for x > 0.5*)*D[u[x,y],x] 

and

D[u[x,y],y](*u[x,y] for left side of the region i.e. x < 0.5*) ==  D[u[x,y],y](*u[x,y] for right side of the region i.e. x > 0.5*) 

I appreciate your help to resolve this issue and correctly implement the Neumann boundary condition at x==0.5.