# How to specify PDE Boundary condition on a B-spline?

## Context

I would like to solve a PDE on a boundary which is parametrized as a BSpline. I am trying to solve the force-free Grad-Shafranov equation on a boundary whose shape I do not know in advance.

Specifically I need to solve for the toroidal flux of the magnetic field above an accretion disc.

R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] - D[P[R, z], R] == - R/2;


and I am seeking solution satisfying P==0 on a spline, see below.

This question is related to the physical context of that question, where we try in to explain astrophysical jets like this: Eventually I would like to optimize the problem while changing the shape of the spline.

## First attempt

I define my region via a BSpline:

ff0 = BSplineFunction[pts = {{1, 0}, {1.2, 2}, {0, 2}}]


So the upper envelope of the jet looks like this:

pl0 = ParametricPlot[ ff0[t] // Release, {t, 0, 1},
Frame -> False, Axes -> False, PlotPoints -> 15, ImageSize -> Small] and the region like that:

pl = ParametricPlot[r ff0[t] // Release, {t, 0, 1}, {r, 0.01, 1},
Frame -> False, Axes -> False, PlotPoints -> 15, ImageSize -> Small] I can then discretize both the boundary and the region:

Ω = DiscretizeGraphics[pl] δΩ = DiscretizeGraphics[pl0, MaxCellMeasure -> 0.1] and then solve for the PDE

eqn0 = R D[P[R, z], {R, 2}] +  R D[P[R, z], {z, 2}] - D[P[R, z], R] == - R/2;
P0 = NDSolveValue[{eqn0,
DirichletCondition[P[R, z] == 0, R == 0],
DirichletCondition[P[R, z] == 0, {R, z} ∈ δΩ],
DirichletCondition[P[R, z] == E  R^2 Log[1/R^2], z == 0]},
P, {R, z} ∈ Ω, Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 1/10000},
"IntegrationOrder" -> 3}}]


If I then try and plot the resulting PDE solution, P0,

ContourPlot[P0[R, z], {R, z} ∈ Ω,
PlotLegends -> Automatic, PlotPoints -> 30,
ColorFunction -> "LightTemperatureMap", ImageSize -> Small,
PlotRange -> All,
FrameLabel -> {R, z},
AspectRatio -> 1] Even though it seems happy, it satisfies very poorly the boundary on the spine:

Plot[ P0 @@ ff0[t], {t, 0, 1}, ImageSize -> Small] This should be zero…

## Second attempt

Following J. M., I have attempted using explicit splines and ParametricRegion as follows:

pts = {{1, 0}, {1.8, 3}, {0, 2}};
{xu, yu} = Transpose[pts];
n = 2;m = Length[pts];
knots = {ConstantArray[0, n + 1], Range[m - n - 1]/(m - n),
ConstantArray[1, n + 1]} // Flatten;
fx[t_] = xu.Table[ BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];
fy[t_] = yu.Table[ BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];


Indeed

ParametricPlot[{fx[t], fy[t]}, {t, 0, 1}, Axes -> None, Frame -> True,
Epilog -> {Directive[AbsolutePointSize, Red], Point[pts]}] seems to return the same spine; now I can define my region and triangulate it:

pr = ParametricRegion[{{r fx[t], r fy[t]}, 1 <= t <= 1 && 0 <= r <= 1}, {t, r}];
Ω = DiscretizeRegion[pr, MaxCellMeasure -> 0.001]
RegionPlot[Ω] and similarly its boundary:

  dpr = ParametricRegion[{{ fx[t], fy[t]}, 0 <= t <= 1}, t];
δΩ = DiscretizeRegion[dpr, MaxCellMeasure -> 0.001];


But applying the same PDE on these regions/boundary with these newly regions yields the same inaccuracies as before (boundary condition not satisfied properly on δΩ).

The problem might be with the second discretize region: indeed

   Show[δΩ, Axes -> True] presents some defect in the triangulation. Note in particular the two points at the origin and at coordinate (0.9,-0.2).

## Questions

Any suggestion on why it fails to satisfy the boundary?

Any suggestion on how to avoid going through DiscretizeGraphics ?

Any suggestion on how to specify DirichletCondition on BSplineFunction?

I feel I am not using the most straightforward method here but…

Thanks!

• Have you tried decomposing your BSplineFunction[] and forming the corresponding ParametricRegion[]? – J. M. will be back soon Jul 23 '15 at 7:57
• no because i did not know about parametric region. Thanks – chris Jul 23 '15 at 8:06
• @user21 would you please be able to throw some suggestions? It seems it is a problem of general interest from the point of view of solving large classes of PDEs? – chris Jul 26 '15 at 21:19

The best way (as pointed out by J. M.) is to convert splines into implicit functions. The real issue you are having is that you'd need a second order mesh to get a decent solution. Note that DiscretizeGraphics and DiscretizeRegion create first order meshes. So you'd need to use ToElementMesh. We also would like to have a finer boundary resolution, thus use "MaxBoundaryCellMeasure". Another thing to think about is the way the boundary condition is specified on the spline. A better way to specify is to say "all boundary elements where R and z are not 0 instead of the code to rest for region member ship on the boundary with Element.

This then gives:

Needs["NDSolveFEM"]
pts = {{1, 0}, {1.8, 3}, {0, 2}};
{xu, yu} = Transpose[pts];
n = 2; m = Length[pts];
knots = {ConstantArray[0, n + 1], Range[m - n - 1]/(m - n),
ConstantArray[1, n + 1]} // Flatten;
fx[t_] = xu.Table[
BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];
fy[t_] = yu.Table[
BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];
pr = ParametricRegion[{{r fx[t], r fy[t]}, -1 <= t <= 1 &&
0 <= r <= 1}, {t, r}];
mesh = ToElementMesh[pr, "MaxBoundaryCellMeasure" -> 0.01];
mesh["Wireframe"] Note that the mesh order is 2.

mesh["MeshOrder"]
2


From there we go to:

eqn0 = R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] -
D[P[R, z], R] == -R/2;
P0 = NDSolveValue[{eqn0,
DirichletCondition[P[R, z] == 0, R == 0],
DirichletCondition[P[R, z] == 0, R != 0 && z != 0],
DirichletCondition[
P[R, z] == E R^2 Log[1/(R + $MachineEpsilon)^2], z == 0] }, P, {R, z} \[Element] mesh];  Note the $MachineEpsilon to avoid division by zero.

ContourPlot[P0[R, z], {R, z} \[Element] mesh,
PlotLegends -> Automatic, PlotPoints -> 30,
ColorFunction -> "LightTemperatureMap", PlotRange -> All,
FrameLabel -> {R, z}, AspectRatio -> 1] And then this is about 2 order of magnitude better:

ff0 = BSplineFunction[pts];
Plot[P0 @@ ff0[t], {t, 0, 1}] Which I hope is reasonable.

Note, that the boundary conditions are set to zero in the interpolating function:

bmesh = ToBoundaryMesh[mesh];
bc = bmesh["Coordinates"];
nodes = DeleteCases[bc, {x_ /; x < 10^-3, y_} | {x_, y_ /; y < 10^-3}];
MinMax[P0 @@@ nodes]
{-1.3877787807814457*^-17, 2.7755575615628914*^-17}


So what you see above is a an interpolation "limitiation" (it's "only" second order accurate). What I am not sure about is why it does not deteriorate further if the boundary is refined. Nevertheless, I think, it's OK to take the derivative of the interpolating function since doing that (currently V10.2) does not evaluates points that are not on the mesh.

• …so that was the secret sauce! :) – J. M. will be back soon Jul 27 '15 at 9:39
• @chris, added some more info. I am not sure how to get a better interpolation, but the solution on the mesh is accurate. – user21 Jul 27 '15 at 10:38
• @chris, if there is no objections, I'd like to use this as a basis for a FEM Best Practice to show the use of splines in the boundary. – user21 Jul 27 '15 at 10:40
• @user21 I found a workaround: this works: Jet0[{{1, 0}, {18, 18}/10, {0, 2.2}}] ! – chris Jul 27 '15 at 19:33
• @chris, yes try to Rationalize all numbers, that may help. Sometimes you can find a solution with more symbolic processing time e.g.: SetSystemOptions[ "FiniteElementOptions" -> {"SymbolicProcessing" -> 10.}] – user21 Jul 28 '15 at 6:42