# Solve a PDE over a region defined by a Bézier patch

I am using NDSolve to find the solution to a PDE over an arbitrary domain. The domain is specified by a Bézier patch.

pts = {{{10, 0}, {0, 30}}, {{15, 10}, {8, 30}}, {{25, 37}, {20,
37}}, {{35, 40}, {10, 50}}};
bezierfunc = BezierFunction[pts];
plot = ParametricPlot[bezierfunc[ξ, η], {ξ, 0, 1}, {η, 0, 1}]


The plot is then discretized into a region over which the PDE can be solved.

Ω = BoundaryDiscretizeGraphics[plot]
sol = NDSolveValue[{
D[u[x, y], x, x] + D[u[x, y], y, y] == -3,
DirichletCondition[u[x, y] == 0, True]
},
u,
{x, y} ∈ Ω];
Plot3D[sol[x, y], {x, y} ∈ Ω]


However, the problem comes when trying to specify boundary conditions along some part of the boundary of the region. For example, we can choose the following curve as the border and discretize it.

border = ParametricPlot[bezierfunc[ξ, 0], {ξ, 0, 1}]
borderRegion = DiscretizeGraphics[border]


Since the border and the region are discretized differently, the border can not be used reliably to specify a boundary region.

sol2 = NDSolveValue[{
D[u[x, y], x, x] + D[u[x, y], y, y] == -3,
DirichletCondition[u[x, y] == 0, {x, y} ∈ borderRegion]
},
u,
{x, y} ∈ Ω];
Plot3D[sol2[x, y], {x, y} ∈ Ω]


Clearly, the boundary condition is only applied near the edges of the boundary, since the discretization of the region and the boundary do not coincide everywhere.

So my question is: what would be a better way to prescribe a boundary condition along partial edges of the Bézier patch?

Here is another way, which is more straightforward than my other answer. At first, I got stumped by couple of things, including, it turns out, a Bug in ArcLength?, and I didn't have time to explore the issues.

Instead of using a "BoundaryMarkerFunction" we can list the markers directly in LineElement[elements, markers]. We can make a fairly general function that create a boundary mesh from a list of edges, marking each edge in turn by its index in the list. The edges are specified by a list of points representing a polygonal path that is part of the boundary of the region.

edgesToBoundaryMesh[edges_, opts : OptionsPattern[ToBoundaryMesh]] :=
With[{lengths = Length /@ edges, bcoords = Join @@ Most /@ edges},
ToBoundaryMesh[
"Coordinates" -> bcoords,
"BoundaryElements" -> {LineElement[
Partition[Range@Length[bcoords], 2, 1, 1],
Flatten[MapIndexed[ConstantArray[#2, #1] &, lengths - 1]]]},
opts
]]


Given a list of parametrizations defining the edges, such as

{bezierfunc[#, 0] &, bezierfunc[1, #] &, bezierfunc[1 - #, 1] &, bezierfunc[0, 1 - #] &}


it would be nice to have a function that creates a discretization of the edges. This can be done with ParametricPlot as the OP shows in the question.

edges = Map[First@Cases[Normal@ParametricPlot[#[t], {t, 0, 1}], Line[p_] :> p, Infinity] &,
{bezierfunc[#, 0] &, bezierfunc[1, #] &, bezierfunc[1 - #, 1] &, bezierfunc[0, 1 - #] &}];
mesh = ToElementMesh@edgesToBoundaryMesh[edges];
ElementMeshWireframe@mesh


The result is okay, if a little nonuniform due to the recursive subdivision of ParametricPlot where the boundary has greater curvature. (One could also simply use Table instead of ParametricPlot to generate the points directly.) The size of the cells in the middle (and overall) can be controlled by the ToElementMesh option "MaxCellMeasure". The initial boundary can be controlled by ParametricPlot options PlotPoints and MaxRecursion.

The ElementMarker for the boundary for the DirichletCondition is 1. The solution to the PDE is obtained with

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == -3,
DirichletCondition[u[x, y] == 0, ElementMarker == 1]},
u, {x, y} ∈ mesh]


For some reason, I was initially curious about what seemed to be somewhat excessive refinement by ParametricPlot and also with controlling the mesh size on each segment of the boundary. When the edges have greatly different lengths, this will be desirable. If we assume that the parametrizations do not change their velocity much, then it seemed one might obtain a more uniform mesh with Table or Range without too much trouble. The subdivision of the interval {t, 0, 1} of the parameterizations should depend on the arc length of the edge and maximum cell measure.

parametricEdgesToBoundaryMesh[edgefns_,
opts : OptionsPattern[ToBoundaryMesh]] :=
With[{lengths = Map[ArcLength@DiscretizeRegion@Line@Table[#[t], {t, 0., 1., 1/32}] &, edgefns]},
With[{ninteverals = Map[Ceiling[(1/If[NumericQ[#], #, Max[lengths]/50.] &@
OptionValue["MaxBoundaryCellMeasure"]) #] &, lengths]},
edgesToBoundaryMesh[
MapThread[Function[t, #[t]] /@ Rescale[Range[0, #2], {0., #2}] &, {edgefns, ninteverals}],
opts]
]];

ClearAll[parametricPatch];
SetAttributes[parametricPatch, HoldAll];
parametricPatch[f_, domain : {{u1_, u2_}, {v1_, v2_}} : {{0, 1}, {0, 1}}] :=
Activate[
{Transpose[{{#1, v1}, {u2, #1}, {u2 - #1, v2}, {u1, v2 - #1}}],
{{u1, u2}, {v1, v2}}}]
];


The utility parametricPatch takes a parametrization f over a rectangular domain and returns a list of parametrizations of the edges over the unit interval.

bmesh = parametricEdgesToBoundaryMesh[
parametricPatch[bezierfunc], "MaxBoundaryCellMeasure" -> 2.];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]


This has 254 triangles compared to over 3000 with the default ParametricPlot. This is not a great advantage in this particular case. I think the important consideration here is to divide the edge segments of the boundary into similarly sized line elements.

• Deeply impressed. Thanks for improving your initial answer even further! – cfdguy Feb 8 '15 at 19:10
• I made a little syntax correction to your code, but it seems it was rejected by somebody. There is an error in the function 'edgesToBoundaryMesh': at the end of line 7 of this function, there is an orphan comma. Please delete this comma. Otherwise the code will give errors. – cfdguy Feb 9 '15 at 8:18
• @cfdguy Thanks! I was notified of the rejection, but only saw the dom/mesh fix. – Michael E2 Feb 9 '15 at 10:59

You might create a NearestFunction to help pick the particular boundary you want. You can use it to mark the boundary elements of an ElementMesh (FEM).

plot = ParametricPlot[
bezierfunc[ξ, η], {ξ, 0, 1}, {η, 0, 1}];

edges = Map[
First@Cases[
Normal@ParametricPlot[#, {t, 10^-5, 1 - 10^-5},
PlotPoints -> 100], Line[p_] :> p, Infinity] &,
{bezierfunc[t, 0],
bezierfunc[1, t],
bezierfunc[1 - t, 1],
bezierfunc[0, 1 - t]
}
];
markernf = Nearest[markerRules];

plotregion = DiscretizeGraphics[plot];

mesh = ToElementMesh[
plotregion,
"BoundaryMarkerFunction" :> (First@*markernf@*Mean /@ # &)];
Show[
mesh["Wireframe"],
mesh["Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementMarkerStyle" -> Red]],
Frame -> True, FrameTicks -> All
]


It's a little hard to see, but the η == 0 boundary elements have ElementMarker equal to 1.

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == -3,
DirichletCondition[u[x, y] == 0, ElementMarker == 1]},
u, {x, y} ∈ mesh];
Plot3D[sol[x, y], {x, y} ∈ Ω]


• Thank you! Very clever solution. – cfdguy Feb 6 '15 at 9:54