ToElementMesh Error/General meshing strategies
Mathematica Asked on September 2, 2021
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];
One Answer
We can solve this problem using
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^-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"]}
Answered by Alex Trounev on September 2, 2021
Add your own answers!
Ask a Question
Get help from others!
Recent Questions
- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?
Recent Answers
- Lex on Does Google Analytics track 404 page responses as valid page views?
- Peter Machado on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?