# Setting Up Boundary Conditions for Magnetostatic PDE

Bug introduced in 11.2.0 and fixed in 11.3.0

The system is a hollow cylinder (thin solenoid) with a current density $\text{J}$ and I'm looking to solve the magnetic potential ($\text{A}$) inside the radius of the cylinder using FEA. I'm struggling with the boundary conditions for the open ends (and in general). Suggestions are appreciated.

u0 = 4 Pi 1*^-7;

j = 100; (*current density parallel to z-axis*)

solution =
NDSolveValue[{
-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"],
{r, th, z}, "Cylindrical"] == NeumannValue[j, r == 1],
DirichletCondition[u[r, z] == 0, z^2 == 1]},
u, {r, 0, 1}, {z, -1, 1},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxBoundaryCellMeasure" -> 0.01}}}]

ContourPlot[solution[r, z], {r, z} ∈ solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium]


To compute the magnetic density field ($\text{B}$):

bField = Curl[solution[r, z], {r, th, z}, "Cylindrical"]


Update:

I've done additional research and the following is an approach by FEMM to solve such problems. The cross section is still in cylindrical coordinates but the open boundary is spherical.

bound = Table[{boundRadius Sin[theta],
boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}];
boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1];

coil = {{1, -1}, {1, 1}};
coilIndex = {{Length[bound] + 1, Length[bound] + 2}};

reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}];

bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil],
"BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}];
mesh = ToElementMesh[bmesh];

Show[
RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2,
ImageSize -> Medium],
bmesh["Wireframe"]
]


I need to apply the following mixed boundary condition to the spherical surface:

$\frac{1}{\mu_0}\frac{\partial A}{\partial r}+\frac{A}{\mu_0 R} = 0$

I have the following failed attempt (NDSolveValue spits out bcnop warning and the internal b.c. is ignored):

u0 = 4 Pi 1*^-7;
j = 100;

solution =
NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"],
{r, th, z}, "Cylindrical"] - j ==
NeumannValue[0, r == 0] -
DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]},
u, {r, z} ∈ mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement"}}]

ContourPlot[solution[r, z], {r, z} ∈ solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium,
AspectRatio -> 2]

• Just a precison for beginners in magnetostatic : this problem can be solved anatically. See here Sep 22, 2017 at 18:27
• That's what I was thinking. My message was for "visitors" Sep 22, 2017 at 18:30
• The problem of open ends (=boundaries at infinity) arrises very often in magnetostatics and electrostatics. Example : In electrostatics in 2 dimensions, on may use a conformal mapping to map infinity to finite boundaries, and then use FEM. Otherwise, FEM solvers have special complex methods for dealing with some "infinite mesh elements" (?). I'm afraid that's not a simple affair. Would be happy to be wrong. Sep 22, 2017 at 20:03
• OK, @xzczd there is a regression in V11.2 that triggered the message and that will need to get fixed. Never the less I think the wire should have some extension. Oct 1, 2017 at 12:07
• @xzczd, I just merged a fix for this problem, the next release will work as expected. Thanks for digging this out. Oct 17, 2017 at 0:57

This is fixed in version 11.3:

Needs["NDSolveFEM"]
boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}];
boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1];

coil = {{1, -1}, {1, 1}};
coilIndex = {{Length[bound] + 1, Length[bound] + 2}};

reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}];

bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil],
"BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}];
mesh = ToElementMesh[bmesh];

Show[RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2,
ImageSize -> Medium], bmesh["Wireframe"]]


Show[Graphics[
GraphicsComplex[mesh["Coordinates"],
Point[Flatten[ElementIncidents[mesh["PointElements"]]]]]]]


Note, how the interior points are now (again) part of the boundary mesh. And the solution:

u0 = 4 Pi 1*^-7;
j = 100;

solution =
NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z},
"Cylindrical"], {r, th, z}, "Cylindrical"] - j ==
NeumannValue[0, r == 0] -
DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]},
u, {r, z} \[Element] mesh];

ContourPlot[solution[r, z], {r, z} \[Element] solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium,
AspectRatio -> 2, PlotPoints -> 50]


Well, I'm not quite familiar with electromagnetism and it's not immediately clear to me how to compare the numeric solution with the analytic solution, but I've circumvented the bcnop warning anyway, so let me post an answer.

The idea is creating a very narrow slit in the domain to approximate the wire i.e. the LineElement in your failed attempt. Notice I've also enlarge boundRadius because after some trial I noticed boundRadius = 2 is too small:

boundRadius = 4;
u0 = 4 Pi 1*^-7;
j = 100;

eps = 10^-4;

reg = ImplicitRegion[
r^2 + z^2 < boundRadius^2 && r >= 0 && ! (1 - eps < r < 1 + eps && z^2 < 1), {r, z}];

solution = NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th,
z}, "Cylindrical"] -
DirichletCondition[u[r, z] == u0 j, 1 - eps < r < 1 + eps && z^2 < 1]},
u, {r, z} ∈ reg,
Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.001}]

interestedreg =
ImplicitRegion[
r^2 + z^2 < interestedRadius^2 && ! (1 - eps < r < 1 + eps && z^2 < 1), {r, z}];

ContourPlot[
Evaluate@PiecewiseExpand@If[r > 0, solution[r, z], solution[-r, z]], {r, z} ∈
interestedreg, PlotRange -> All, PlotPoints -> 100, Contours -> 19]


• I think this is the way to go - also from a physics point of view. You may want to create the mesh manually like in the question for very 'thin' wires. Oct 1, 2017 at 8:41
• I eventually obtained a result close to your graph but when you take the Curl of solution to obtain the magnetic field it doesn't match the analytical solution. Oct 1, 2017 at 17:32
• I think there is something wrong with how I formed the differential equation. Oct 1, 2017 at 17:39
• @Young Er…sorry, but the $A$ i.e. u in your equation is a scalar, how can you calculate its Curl? Oct 2, 2017 at 2:50
• The use of PiecewiseExpand seems unnecessary; replace Evaluate@PiecewiseExpand@If[r > 0, solution[r, z], solution[-r, z]] by Evaluate@solution[Abs@r, z] Oct 7, 2021 at 11:09

Extended comment

The result of the OP's code are not the same on Mathematica 11.0 and Mathematica 11.2

Here is the code I have used :

<<NDSolveFEM

boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}];
boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1];

coil = {{1, -1}, {1, 1}};
coilIndex = {{Length[bound] + 1, Length[bound] + 2}};

reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}];

bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil],
"BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}];
mesh = ToElementMesh[bmesh];

(*Show[
RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2,
ImageSize -> Medium],
bmesh["Wireframe"]
]*)

u0 = 4 Pi 1*^-7;
j = 100;

solution =
NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"],
{r, th, z}, "Cylindrical"] - j ==
NeumannValue[0, r == 0] -
DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]},
u, {r, z} \[Element] mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement"}}]

ContourPlot[solution[r, z], {r, z} \[Element] solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium,
AspectRatio -> 2]


Mathematica 11.2

Mathematica 11.0 - first evaluation (on a fresh kernel)

It doesn't work :

Mathematica 11.0 - second evaluation (= reevaluation of the same code)

Edit

same as before, but with PlotPoints->50added to ContourPlot:

• Will the last result be improved if you add e.g. PlotPoints->50 to ContourPlot? Oct 1, 2017 at 10:55
• @xzczd see my edit. Nice ! Oct 1, 2017 at 11:07
• It helps if you define the boundRadius = 2; before the bound. Oct 1, 2017 at 11:54
• OK, thanks Andre (@xzczd) this is a regression in 11.2 and I'll see to get it fixed for the next release. Oct 1, 2017 at 12:05
• A precision for visitors : the comment of @user21 above (boundRadius =2 before the bound) is destinated to fix the problem of the failing first evaluation on Mathematica version 11.0 Oct 1, 2017 at 12:08

Extended comment:

Something like this

epsi = 10^-1;
coil = {{1 - epsi, -1}, {1 - epsi, 1}, {1 + epsi,
1}, {1 + epsi, -1}};
coilIndex = {{1, 2}, {2, 3}, {3, 4}, {4, 1}} + Length[bound];