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["NDSolve`FEM`"]
(* 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];
ToBoundaryMeshwon't resolve your problem.NDSolvewill 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. $\endgroup$