TransWikia.com

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"]}

Figure 1

Answered by Alex Trounev on September 2, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP