# ToElementMesh Error/General meshing strategies

I am trying to create a mesh using ToBoundaryMesh as opposed to DiscretizeRegion. I was hoping this would be a better strategy since most of the action in the problem I am trying to solve happens in a very small region along the right boundary; however I am getting the following error from ToElementMesh

ToElementMesh: The mesh elements are not valid. A set of valid mesh element incidents needs to be positive integers and be able to form a complete sequence starting from 1 to the largest incident present. There are missing incidents; a complete sequence cannot be formed.

Why am I getting this error? As far as I can tell, the mesh elements incidents do form a complete sequence.

Additionally, there are several different ways to create a mesh in Mathematica and I am not familiar with all of them and what types of problems they are best suited for. Could anyone suggest the best strategy given the problem I am trying to solve? The aspect ratio of the domain is extremely high and the segments of the boundary on the right side where things happen are on the nano scale. So this problem requires extremely fine mesh elements in that region.

I've included the full code, including the problem I am trying to solve once the mesh is working, below. I would appreciate any help or insight I can get!

ClearAll["Global*"]
Needs["NDSolveFEM"]
(* 1) Define Constants*)
e = 1.60217662*10^-19;
F = 96485;
kb = 1.381*10^-23;
sigi = 18;
sigini = 0;
sigeni = 2*10^6;
T = 1000;
n = -0.02;
c = 1;
pH2 = 0.2;
pH2O = 1 - pH2;
pO2 = 1.52*^-19;
l = 10*10^-6;
a = 100*10^-9;
b = 50*10^-9;
d = 300*10^-9;
y1 = 0.01;
y2 = 0.5*y1;
y3 = y2 + a;
y4 = y3 + d;
y5 = y4 + b;
mu1 = 0;
mu2 = -5.98392*^-19;
phi1 = 0;

(* 2) Create mesh*)
m = 0.1*l;
size1 = 10^-17;
size2 = 10^-15;
size3 = 10^-7;
pts = {{0, 0}, {l, 0}, {l, y2}, {l, y3}, {l, y4}, {l, y5}, {l,
y1}, {0, y1}};
incidents = Partition[FindShortestTour[pts][[2]], 2, 1];
markers = {1, 2, 3, 4, 5, 6, 1, 7};
bcEle = {LineElement[incidents, markers]};
bmesh = ToBoundaryMesh["Coordinates" -> pts,
"BoundaryElements" -> bcEle];
mrf = With[{rmf =
RegionMember[
Region@RegionUnion[Disk[{l, y2}, m], Disk[{l, y3}, m],
Disk[{l, y4}, m], Disk[{l, y5}, m]]]},
Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
Which[rmf[{x, y}],
area > size1, (0 <= x <= l && y2 - l <= y <= y2 + l),
area > size2, (0 <= x <= l && y3 - l <= y <= y3 + l),
area > size2, (0 <= x <= l && y4 - l <= y <= y4 + l),
area > size2, (0 <= x <= l && y5 - l <= y <= y5 + l),
area > size2, True, area > size3]]]];
mesh = ToElementMesh[bmesh, MeshRefinementFunction -> mrf];

(* 3) Solve for mu*)
bcmu = {DirichletCondition[mu[x, y] == mu1, (x == 0 && 0 < y < y1)],
DirichletCondition[
mu[x, y] ==
mu2, (x == l && y2 <=  y <=  y3) || (x == l && y4 <= y <= y5)]};
solmu = NDSolve[{Laplacian[mu[x, y], {x, y}] ==
0 + NeumannValue[0, y == 0 || y == y1 ||
(x == l && 0 <= y < y2) || (x == l &&
y3 < y < y4) || (x == l && y5 < y < y1)], bcmu},
mu, {x, y} \[Element] mesh, WorkingPrecision -> 50];

(* 4) Solve for electronic conductivity everywhere*)
pO2data = Exp[(mu[x, y] /. solmu)/kb/T];
sige0 = 2.77*10^-7;
sigedata = Piecewise[{{sige0*pO2data^(-1/4), 0 <= x <= l - m},
{sige0*pO2data^(-1/4), (l - m < x <= l && 0 <= y < y2)},
{(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/m*(x - (l - m)) +
sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l &&
y2 <=  y <= y3)},
{sige0*pO2data^(-1/4), (l - m < x <= l && y3 < y < y4)},
{(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/m*(x - (l - m)) +
sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l &&
y4 <= y <= y5)},
{sige0*pO2data^(-1/4), (l - m < x <= l && y5 < y <= y1)}}];

(* 5) Solve for phi*)
Irxn = -(2*F)*(c*pO2^n );
A = (Irxn - sigi/(4*e)*(D[mu[x, y] /. solmu, x] /. x -> l))/(-sigi);
B = sigi/(4*e)*(D[mu[x, y] /. solmu, x] /.
x -> l)/(sigi + sigedata /. x -> l - m);
bcphi = DirichletCondition[phi[x, y] == phi1, (x == 0 && 0 < y < y1)];
solphi = NDSolve[{Laplacian[phi[x, y], {x, y}] ==
0 + NeumannValue[0,
y == 0 ||
y == y1 || (x == l && 0 <= y < y2) || (x == l &&
y3 < y < y4) || (x == l && y5 < y < y1)] +
NeumannValue[-A[[1]], (x == l && y2 <= y <= y3)] +
NeumannValue[-B[[1]], (x == l && y4 <= y <= y5)], bcphi},
phi, {x, y} \[Element] mesh, WorkingPrecision -> 50];

• Related: (226486). Jul 25, 2020 at 14:57
• Using ToBoundaryMesh won't resolve your problem. NDSolve will still have to mesh the whole domain. Once again, the boundary element method that I mentioned in my comment on the other post is not built into Mathematica. Jul 25, 2020 at 14:59
• @kjcole, what problem do you try to solve? Jul 25, 2020 at 22:27

We can solve this problem using

ClearAll["Global*"]
Needs["NDSolveFEM"]

(*1) Define Constants*)
e = 1.60217662*10^-19;
F = 96485;
kb = 1.381*10^-23;
sigi = 18;
sigini = 0;
sigeni = 2*10^6;
T = 1000;
n = -0.02;
c = 1;
pH2 = 0.2;
pH2O = 1 - pH2;
pO2 = 1.52*^-19;
l = 10*10^-6;
a = 100*10^-7;
b = 50*10^-7;
d = 300*10^-7;
y1 = 0.01;
y2 = 0.5*y1;
y3 = y2 + a;
y4 = y3 + d;
y5 = y4 + b;
mu1 = 0;
mu2 = -5.98392*^-19;
phi1 = 0;

(*2) Create mesh*)
m = 0.1*l;
size1 = 10^-16;
size2 = 10^-15;
size3 = 10^-7;
mesh = ToElementMesh[Rectangle[{0, 0}, {l, y1}],
MaxCellMeasure -> 5 10^-12]

(*3) Solve for mu*)
bcmu = {DirichletCondition[mu[x, y] == mu1, (x == 0 && 0 < y < y1)],
DirichletCondition[
mu[x, y] ==
mu2, (x == l && y2 <= y <= y3) || (x == l && y4 <= y <= y5)]};
solmu = NDSolve[{Laplacian[mu[x, y], {x, y}] ==
0 + NeumannValue[0,
y == 0 ||
y == y1 || (x == l && 0 <= y < y2) || (x == l &&
y3 < y < y4) || (x == l && y5 < y < y1)], bcmu},
mu, {x, y} \[Element] mesh];

(*4) Solve for electronic conductivity everywhere*)
pO2data = Exp[(mu[x, y] /. solmu)/kb/T];
sige0 = 2.77*10^-7;
sigedata =
Piecewise[{{sige0*pO2data^(-1/4),
0 <= x <= l - m}, {sige0*
pO2data^(-1/4), (l - m < x <= l &&
0 <= y <
y2)}, {(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/
m*(x - (l - m)) +
sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l &&
y2 <= y <= y3)}, {sige0*
pO2data^(-1/4), (l - m < x <= l &&
y3 < y <
y4)}, {(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/
m*(x - (l - m)) +
sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l &&
y4 <= y <= y5)}, {sige0*
pO2data^(-1/4), (l - m < x <= l && y5 < y <= y1)}}];

(*5) Solve for phi*)
Irxn = -(2*F)*(c*pO2^n);
A = (Irxn - sigi/(4*e)*(D[mu[x, y] /. solmu, x] /. x -> l))/(-sigi);
B = sigi/(4*e)*(D[mu[x, y] /. solmu, x] /.
x -> l)/(sigi + sigedata /. x -> l - m);
bcphi = DirichletCondition[phi[x, y] == phi1, (x == 0 && 0 < y < y1)];
solphi = NDSolve[{Laplacian[phi[x, y], {x, y}] ==
0 + NeumannValue[0,
y == 0 ||
y == y1 || (x == l && 0 <= y < y2) || (x == l &&
y3 < y < y4) || (x == l && y5 < y < y1)] +
NeumannValue[-A[[1]], (x == l && y2 <= y <= y3)] +
NeumannValue[-B[[1]], (x == l && y4 <= y <= y5)], bcphi},
phi, {x, y} \[Element] mesh];

(*6) Print values to check for convergence*)
P[x_, y_] := phi[x, y] /. solphi;
P[l, (y3 - y2)/2 + y2]
P[l, (y5 - y4)/2 + y4]
(*{-0.50881}
{-0.433773}*}


Visualization

{DensityPlot[
Evaluate[Exp[(mu[x, y] /. solmu)/kb/T]], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotPoints -> 200, PlotRange -> All, PlotLabel -> "Exp[mu]/kbT"],
DensityPlot[Evaluate[phi[x, y] /. solphi], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotPoints -> 200, PlotRange -> All, PlotLabel -> "phi"]}


• @ Alex Trounev thanks for your response! Did you change something in the first two sections of my code? When I run the first two sections of mine + the edit you made to mesh and the rest your code I get the error "NDSolve: No places were found on the boundary where x==1/100000 && 0.005<=y<=0.0050001 was True...." Jul 27, 2020 at 15:30
• @kjcole I have removed mesh generator only. I added first part of code to my answer for clearance. What version you are running? Jul 27, 2020 at 15:52