I'm trying to do the same thing I did in this thread (ignore the actual question I was asking about, I just linked it to show the geometry and equation solving I'm doing), but this time for 3D instead of 2D. So, instead of having two rectangles at different potentials, I now have two rectangular prisms.
Here is the code, I apologize for its lengthiness but it's mostly boundary conditions:
interfaceheight = 5;
interfacewidth = 20;
separation = 5;
Clear[x, y, z];
simxbds = {0, 100};
simybds = {0, 100};
simzbds = {0, 100};
xandbds = Flatten@({x, simxbds});
yandbds = Flatten@({y, simybds});
zandbds = Flatten@({z, simzbds});
xmidpt = (Last@simxbds)/2;
ymidpt = (Last@simybds)/2;
zmidpt = (Last@simzbds)/2;
lrectxbd = xmidpt - separation/2;
rrectxbd = xmidpt + separation/2;
ybdbot = ymidpt - interfaceheight/2;
ybdtop = ymidpt + interfaceheight/2;
zbdlow = zmidpt - interfacewidth/2;
zbdhigh = zmidpt + interfacewidth/2;
xbdslrect = First@simxbds <= x <= lrectxbd;
xbdsrrect = rrectxbd <= x <= Last@simxbds;
ybds = ybdbot <= y <= ybdtop;
zbds = zbdlow <= z <= zbdhigh;
(*Left rectangle BC's.*)
(*XY faces.*)
lrectzlowXYBC = (xbdslrect && ybds && z == zbdlow);
lrectzhighXYBC = (xbdslrect && ybds && z == zbdhigh);
(*XZ faces.*)
lrectybotXZBC = (xbdslrect && zbds && y == ybdbot);
lrectytopXZBC = (xbdslrect && zbds && y == ybdtop);
(*YZ faces.*)
lrectxlowYZBC = (ybds && zbds && x == First@simxbds);
lrectxhighYZBC = (ybds && zbds && x == lrectxbd);
lrectinterior = (xbdslrect && ybds && zbds);
lrectBCs =
lrectzlowXYBC || lrectzhighXYBC || lrectybotXZBC || lrectytopXZBC ||
lrectxlowYZBC || lrectxhighYZBC || lrectinterior;
(*Right rectangle BC's.*)
(*XY faces.*)
rrectzlowXYBC = (xbdsrrect && ybds && z == zbdlow);
rrectzhighXYBC = (xbdsrrect && ybds && z == zbdhigh);
(*XZ faces.*)
rrectybotXZBC = (xbdsrrect && zbds && y == ybdbot);
rrectytopXZBC = (xbdsrrect && zbds && y == ybdtop);
(*YZ faces.*)
rrectxlowYZBC = (ybds && zbds && x == rrectxbd);
rrectxhighYZBC = (ybds && zbds && x == Last@simxbds);
rrectinterior = (xbdsrrect && ybds && zbds);
rrectBCs =
rrectzlowXYBC || rrectzhighXYBC || rrectybotXZBC || rrectytopXZBC ||
rrectxlowYZBC || rrectxhighYZBC || rrectinterior;
Print@lrectBCs;
Print@rrectBCs;
lrectregion =
Cuboid[{First@simxbds, ybdbot, zbdlow}, {lrectxbd, ybdtop, zbdhigh}];
rrectregion =
Cuboid[{rrectxbd, ybdbot, zbdlow}, {Last@simxbds, ybdtop, zbdhigh}];
wholeregion =
Cuboid[{First@simxbds, First@simybds, First@simzbds}, {Last@simxbds,
Last@simybds, Last@simzbds}];
\[CapitalOmega] =
RegionDifference[wholeregion, RegionUnion[lrectregion, rrectregion]];
sol = NDSolveValue[{D[u[x, y, z], x, x] + D[u[x, y, z], y, y] +
D[u[x, y, z], z, z] == 0,
DirichletCondition[u[x, y, z] == Vapp, lrectBCs],
DirichletCondition[u[x, y, z] == 0.0, rrectBCs]},
u, {x, y, z} \[Element] \[CapitalOmega]];
Print@DensityPlot3D[
sol[x, y, z], {x, y, z} \[Element] \[CapitalOmega],
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic, ImageSize -> Large,
AxesLabel -> {"x", "y", "z"}];
gradField =
ComplexExpand[{D[sol[x, y, z], x], D[sol[x, y, z], y],
D[sol[x, y, z], z]}];
Print@Show[
Graphics3D[{Opacity[0], wholeregion, Red, Opacity[.5], lrectregion,
Blue, Opacity[.5], rrectregion}, ImageSize -> Large],
SliceVectorPlot3D[
gradField, {{y == ymidpt}, {z == zmidpt}}, {x, First@simxbds,
Last@simxbds}, {y, First@simybds, Last@simybds}, {z,
First@simzbds, Last@simzbds}, VectorPoints -> 30,
VectorScale -> Small]];
To give an idea of the geometry produced by this (and to show that I think it's working), here is what the Graphics3D part produces:
The first issue is the following. In the 2D example before, I specified the boundary conditions (BCs) using the DirichletConditions with just the boundaries of the rectangles, i.e., the edges, and it worked. Adding the interior of them as a boundary condition changed nothing as far as I could tell.
However, in 3D, if I only specify the walls of the rectangles (without the interior), it gives me this warning/error:
NDSolveValue::bcnop: "No places were found on the boundary where (105/2<=x<=100&&95/2<=y<=105/2&&z==40)||(105/2<=x<=100&&95/2<=y<=105/2&&z==60)||(105/2<=x<=100&&40<=z<=60&&y==95/2)||(105/2<=x<=100&&40<=z<=60&&y==105/2)||(95/2<=y<=105/2&&40<=z<=60&&x==105/2)||(95/2<=y<=105/2&&40<=z<=60&&x==100) was True, so DirichletCondition[u==0.,(105/2<=x<=100&&95/2<=y<=105/2&&z==40)||(105/2<=x<=100&&95/2<=y<=105/2&&z==60)||(105/2<=x<=100&&40<=z<=60&&y==95/2)||(105/2<=x<=100&&40<=z<=60&&y==105/2)||(95/2<=y<=105/2&&40<=z<=60&&x==105/2)||(95/2<=y<=105/2&&40<=z<=60&&x==100)] will effectively be ignored. "
If I add the interiors (the BCs called lrectinterior and rrectinterior), it gets rid of the error. So my question is, why does the interior seem to be required in 3D but not 2D?
I found this thread that seems related, but I can't tell what the answer actually changed.
My 2nd related question is regarding DensityPlot3D. When I use it to plot the solution of NDSolveValue (with the interiors), it gives this:
(the rectangles are turned on their side from in the first picture.)
There are several issues. It's seemingly displaying just two parts of the density spectrum. Also, if you look very close around the rectangles, you don't see the colors corresponding to the values I assigned to them with DirichletConditions. You can actually see a little (of the red in this case) if you turn it to the side and view the back of one of the rectangles, but it's still nothing like for the 2D version:
What could be causing it behave so differently in 3D than in 2D? Am I missing something obvious? thank you!
