1
$\begingroup$

I am trying to create two density plots of values but one plot has blank spots and the other is entirely blank. I do not get any errors when I run the code so I am not sure why the plots are turning out like this. Below is my code. I've separated it into 3 sections; defining constants that will be used in proceeding calculations and creating the mesh, solving a differential equation using NDSolve, and using the results of NDSolve to solve for the values I want to plot. Throughout the domain, the value of sige is a function of mu and sigi is constant except in two sections of the right boundary. There is a discontinuity in sige and sigi at the boundary, so to account for this I have created a very small transition region next to the boundary where sigi and sige linearly transition from their values within the domain to the values at the boundary. This is the purpose of the Piecewise functions in the third section.

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
(* 1) Define Constants and create mesh*)
e = 1.60217662*10^-19;
sigi0 = 18; 
sigini = 0;
sigeni = 1.5*10^7;
T = 1000; 
R = 8.314; 
kb = 1.381*10^-23; 
mu2 = -5.984*^-19;
l = 10*10^-6;
y1 = 0.01; 
y2 = 0.0025; 
y3 = 0.005; 
y4 = 0.006; 
y5 = 0.007; 
pts = {{0, 0}, {l, 0}, {l, y2}, {l, y3}, {l, y4}, {l, y5}, {l, 
    y1}, {0, y1}, {0, y5}, {0, y4}, {0, y3}, {0, y2}};
incidents = Partition[FindShortestTour[pts][[2]], 2, 1];
markers = {1, 2, 3, 4, 5, 6, 1, 7, 7, 7, 7, 7};
bcEle = {LineElement[incidents, markers]};
bmesh = ToBoundaryMesh["Coordinates" -> pts, 
   "BoundaryElements" -> bcEle];
mrf = With[{rmf = 
     RegionMember[
      Region@RegionUnion[Disk[{l, y2}, 0.0025], Disk[{l, y3}, 0.0025],
         Disk[{l, y4}, 0.0025], Disk[{l, y5}, 0.0025]]]}, 
   Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
     Which[rmf[{x, y}], 
      area > (2.5*10^-10)/258, (x > 0.9*l) && (y2*0.9 <= y <= y3*1.1),
       area > 2.5*10^-10, (x > 0.9*l) && (y4*0.9 <= y <= y5*1.1), 
      area > 2.5*10^-10, True, area > 2.5*10^-7]]]];
mesh = ToElementMesh[bmesh, MeshRefinementFunction -> mrf];

(* 2) Solve for mu profile*)
mu1 = 0;
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]; 

(* 3) Solve for sigi and sige everywhere and plot results*)
pO2data = Exp[(mu[x, y] /. solmu)/kb/T];
sige0 = 2.77*10^-7;
sigedata = 
  Piecewise[{{sige0*pO2data^(-1/4), 
     0 <= x <= 
      0.99999*l}, {(sigeni - 
          sige0*(pO2data /. x -> 0.99999*l)^(-1/4))/(0.99999*l)*x + 
      sige0*(pO2data /. x -> 0.99999*l)^(-1/4), 0.99999*l < x <= l}}];
sigidata = 
  Piecewise[{{sigi0, 
     0 <= x <= 0.99999*l}, {(sigini - sigi0)/(0.99999*l)*x + 
      sigi0, (0.99999*l < x <= l && 
       y2 <=  y <= y3)}, {(sigini - sigi0)/(0.99999*l)*x + 
      sigi0, (0.99999*l < x <= l && 
       y4 <= y <= y5)}, {sigi0, (0.99999*l < x <= l && 
       0 <= y < y2)}, {sigi0, (0.99999*l < x <= l && 
       y3 < y < y4)}, {sigi0, (0.99999*l < x <= l && y5 < y <= y1)}}];
DensityPlot[sigedata, {x, 0, l}, {y, 0, y1}, 
 PlotLegends -> Automatic, PlotLabel -> "Sige Results"]
DensityPlot[sigidata, {x, 0, l}, {y, 0, y1}, PlotLegends -> Automatic,
  PlotLabel -> "Sigi Results"]

This is what the resulting plots look like:

enter image description here

enter image description here

Does anyone know why this is happening or have suggested edits to my code to improve the image resolution? Thanks!

*Update: the suggested edits in the comments fixed the blank plot issue with sigi and improved the resolution of both plots so the last two lines of the code are

DensityPlot[sigedata, {x, 0, l}, {y, 0, y1}, 
 PlotLegends -> Automatic, PlotLabel -> "Sige Results", 
 Exclusions -> None, PlotPoints -> 50]
DensityPlot[sigidata, {x, 0, l}, {y, 0, y1}, PlotLegends -> Automatic,
  PlotLabel -> "Sigi Results", Exclusions -> None, PlotPoints -> 50]

and the plots look like this

enter image description here enter image description here

The plot of sige still has blank spots, and I am not sure why. I know there should be values where the blank spots are because if I plot a slice of the sigedata at x=8*10^-6

Plot[sigedata /. x -> 8*10^-6, {y, 0, y1}, PlotRange -> All, 
 PlotLabel -> "Sige Results at x = 8*10^-6", 
 AxesLabel -> {"y-Position", "sige"}]

this is the result

enter image description here

Why isn't the DensityPlot for sige showing these values? I added Exclusions->None to the function but it is still excluding values.

$\endgroup$
3
  • 2
    $\begingroup$ Table[sigidata, {y, 0, y1, y1/30}, {x, 0, l, l/30}] // MatrixPlot your data is almost all 18 over the ranges for y1 and l. If you add Exclusions -> None to your DensityPlot then you'll see it. $\endgroup$
    – flinty
    Jul 13 '20 at 15:48
  • $\begingroup$ @flinty thanks for your help! Do you know how to fix the white patches in the sige plot? I added Exclusions -> None to that DensityPlot as well but it didn't help $\endgroup$
    – kjcole
    Jul 13 '20 at 18:49
  • $\begingroup$ Maybe DensityPlot the Log of sige instead. The resolution could be improved with PlotPoints -> 50. $\endgroup$
    – flinty
    Jul 13 '20 at 19:31
1
$\begingroup$

After doing further research into the options for DensityPlot I was able to answer my own question. Adding a Log scale to the function with ScalingFunctions-> "Log" produces a full plot of all sige values

DensityPlot[sigedata, {x, 0, l}, {y, 0, y1}, 
 PlotLegends -> Automatic, PlotLabel -> "Sige Results", 
 ScalingFunctions -> "Log", Exclusions -> None, PlotPoints -> 100]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.