# Solve 3D Electric Field using FEA

The objective is to solve the 3D electric fields for the geometry below and generate the continuous function eField3D[x_,y_,z_]. The method I was planning on using was FEA to solve Laplace PDEs.

I have the following mesh generation code:

Needs["NDSolveFEM"]

cylinOffset = 0.032;
cylinLength = 0.100;

shellLength = 0.100;

cylinA = Cylinder[{{-cylinOffset, 0, -cylinLength /2},
cylinB = Cylinder[{{cylinOffset, 0, -cylinLength /2},

shell = Cylinder[{{0, 0, -shellLength /2}, {0, 0, shellLength /2}}, shellRadius];

bmeshShell = ToBoundaryMesh[shell]["Wireframe"];
bmeshCylinA = ToBoundaryMesh[cylinA]["Wireframe"];
bmeshCylinB = ToBoundaryMesh[cylinB]["Wireframe"];

Show[bmeshShell, bmeshCylinA, bmeshCylinB] I'm quite confused on how to progress further. The following is an attempt to setup the boundary conditions of voltage potential and region where the equations should be evaluated, but I'm far away from a solution. Guidance would be very much appreciated.

reg = {bmeshShell, bmeshCylinA, bmeshCylinB};
sol = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 0,
DirichletCondition[bmeshShell == 0,
bmeshCylinA == bmeshCylinB == 30000]},
u, {x, y, z} \[Element] reg];

• How do your boundary conditions work? Are you saying you want the function u[x, y, z] to be zero at the bmeshShell boundary and 30000 at the other two boundaries? I think the syntax of your DirichletCondition is not correct – Jason B. Jun 25 '16 at 18:08
• Yes. Zero potential at the shell and 30 000 at each cylinder. I'm sure DirichletCondition is wrong and the region isn't setup correctly either. I'm working on figuring out the region definition now. – Young Jun 25 '16 at 18:14
• I think it should be more like this, sol = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 0, DirichletCondition[u[x, y, z] == 0, bmeshShell], DirichletCondition[u[x, y, z] == 30000, bmeshCylinA], DirichletCondition[u[x, y, z] == 30000, bmeshCylinB]}, u, {x, y, z} \[Element] shell] but I still don't get it to give me an answer. Maybe we can get @user21 to take a look at this, he's the expert – Jason B. Jun 25 '16 at 18:18
• Try a 2D example first and the extend that to 3D. Also, reg = {bmeshShell, bmeshCylinA, bmeshCylinB} is not how it works. I think I'd start expressing the region as an ImplicitRegion. – user21 Jun 27 '16 at 3:26

Using user21's advice, I built the 2D model of the system then expanded it to 3D.

Here is the general 3D solution:

ClearAll["Global*"]
Needs["NDSolveFEM"];

cylinLength = 73.0;
cylinBuffer =  5.0;
cylinOffset = 32.0;

cylinVoltage = 30000;

shell = Cylinder[{{0, 0, -cylinLength/2}, {0, 0, cylinLength/2}}, shellRadius];
cylin = RegionUnion[
Cylinder[{{-cylinOffset, 0, -cylinLength/2}, {-cylinOffset, 0,
Cylinder[{{cylinOffset, 0, -cylinLength/2}, {cylinOffset, 0,

region = RegionDifference[shell, cylin];

solMV = 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] == 0, z^2 == (cylinLength/2)^2],
DirichletCondition[u[x, y, z] == 0,
x^2 + y^2 == shellRadius^2 && z^2 < (cylinLength/2)^2 + (cylinBuffer)^2],
DirichletCondition[
u[x, y, z] ==
cylinVoltage, (x + cylinOffset)^2 + y^2 == cylinRadius^2 &&
z^2 < (cylinLength/2)^2 + (cylinBuffer)^2],
DirichletCondition[
u[x, y, z] ==
cylinVoltage, (x - cylinOffset)^2 + y^2 == cylinRadius^2 &&
z^2 < (cylinLength/2)^2 + (cylinBuffer)^2]},
u, {x, y, z} ∈ region, InterpolationOrder -> All,
Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 1}}];

solME[x_, y_, z_] = -Grad[solMV[x, y, z], {x, y, z}];

Show[
SliceContourPlot3D[solMV[x, y, z],
"CenterPlanes", {x, y, z} ∈ sol["ElementMesh"], Contours -> 20,
PlotRange -> All, PlotLegends -> Automatic],
VectorPlot3D[
` 