Smoothing the gradient of the potential function from a finite element calculation

Question

The results from a finite element calculation using NDSolve are, in this case, a potential function expressed as an interpolation function. To get values from the potential function it is necessary to calculate the gradient of the potential function and this leads to a poor quality function that needs some smoothing. How can one smooth to extract the best results?

Here is a minimum working example in which I generate a mesh and show the whole mesh and the region of interest. The rather convoluted method is due to a workaround. For my real case the thickness tk is much smaller.

Needs["NDSolve`FEM`"];
Lp = 1; (* length of half plate *)
tk = 0.2;(* plate thickness *)
Ls = 9; (* from plate to boundary *)
Lu = 5; (* upstream distance *)
Ld = 5; (* downstream distance *)
r1 = Rectangle[{-Lp, -Ld}, {Ls, Lu}];
r2 = ImplicitRegion[x^2 + y^2 <= (tk/2)^2, {x, y}];
r3 = Rectangle[{-Lp, -tk/2}, {0, tk/2}];
rd = RegionDifference[r1, RegionUnion[r2, r3]];

area1 = (tk/10)^2 Sqrt[3]/4; area2 = (Ls/10.)^2 Sqrt[3]/4;
cf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
   Block[{r, r0 = tk/2},
    r = Norm@(Total[c]/3);
    If[a > area1 + (area2 - area1) ((r - tk/2)/(Ls - tk/2))^2, True, 
     False]
    ]
   ];
nr = ToNumericalRegion[rd];
m1 = ToElementMesh[rd, MeshRefinementFunction -> cf, "MeshOrder" -> 1,
    "MaxBoundaryCellMeasure" -> tk/10];
m2 = ImproveBoundaryPosition[m1, nr];
m3 = MeshOrderAlteration[m2, 2];
mesh = ImproveBoundaryPosition[m3, nr];
mesh["Wireframe"]
Show[
 m4["Wireframe"],
 m4["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementStyle" -> Red]], PlotRange -> {{-tk, tk}, {-tk, tk}}, 
 Frame -> True]

We now solve Laplace's equation and then calculate the gradient of the solution to get velocities which I plot in a StreamPlot.

sol = NDSolveValue[{
    D[u[x, y], x, x] + D[u[x, y], y, y] ==
     NeumannValue[-1, -Lp <= x <= Ls && y == -Ld] + 
      NeumannValue[1, -Lp <= x <= Ls && y == Lu],
    DirichletCondition[u[x, y] == 0, x == -Lp && y == Lu]
    },
   u, {x, y} ∈ mesh];
ClearAll[vel];
vel[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}]]
subReg = RegionDifference[Rectangle[2 {-tk, -tk}, 2 {tk, tk}], 
   RegionUnion[r2, r3]];
StreamPlot[vel[x, y], {x, y} ∈ subReg]

Now I work out the magnitude of the velocity at various locations close to the curved boundary.

ClearAll[mag];
mag[θ_, eps_] := 
 Norm@vel[((1 + eps) tk)/2 Cos[θ], ((1 + eps) tk)/
    2 Sin[θ]]
Plot[Evaluate@
  Table[mag[θ, 
    eps], {eps, {0, 0.01, 0.02, 0.05, 0.1}}], {θ, -π/
  8, π/8}, PlotLegends -> LineLegend[{0, 0.01, 0.02, 0.05, 0.1}], 
 Frame -> True, FrameLabel -> {"Angular Position", "Velocity"}]

Clearly the velocity values are very ragged. What, for example is the best value at θ = 0? I think this must be expected. We have a second order grid and then take the gradient so we are back to first order interpolation. How should one smooth these results? Should one evaluate the velocities at the mesh nodes and use these values? Or perhaps use the values from the StreamPlot? How do you find the relevant values? Ideas needed. Of course one can go to a finer mesh but that will still have some scatter that needs smoothing. Thanks.

Answer 1

One thing you can do is refine the mesh just behind the curvature by adding "IncludePoints" and changing the geometry factor of 0.2 to 2/10 to get a better mesh.

Set up the region:

Needs["NDSolve`FEM`"];
Lp = 1;(*length of half plate*)tk = 
 2/10;(*plate thickness*)Ls = 9;(*from plate to boundary*)Lu = \
5;(*upstream distance*)Ld = 5;(*downstream distance*)r1 = 
 Rectangle[{-Lp, -Ld}, {Ls, Lu}];
r2 = ImplicitRegion[x^2 + y^2 <= (tk/2)^2, {x, y}];
r3 = Rectangle[{-Lp, -tk/2}, {0, tk/2}];
rd = RegionDifference[r1, RegionUnion[r2, r3]];

This computes the points:

magL[\[Theta]_, eps_] := 
 List[((1 + eps) tk)/2 Cos[\[Theta]], ((1 + eps) tk)/2 Sin[\[Theta]]]

Create the mesh:

mesh = ToElementMesh[rd, "MaxBoundaryCellMeasure" -> tk/15, 
   "IncludePoints" -> 
    Table[magL[\[Theta], 
      10^-2], {\[Theta], -\[Pi]/8, \[Pi]/8, \[Pi]/64}]];

Check the area difference:

Area[rd] - Total[Join @@ mesh["MeshElementMeasure"]]
-4.220800688126136`*^-9

Inspect the mesh:

Show[mesh["Wireframe"], 
 mesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementStyle" -> Red]], PlotRange -> {{-tk, tk}, {-tk, tk}}, 
 Frame -> True]

Very careful: these newly included points can be used as boundary conditions (DirichletConditions) but since none are specified in that part of the region, that's probably fine.

Solve:

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 
     NeumannValue[-1, -Lp <= x <= Ls && y == -Ld] + 
      NeumannValue[1, -Lp <= x <= Ls && y == Lu], 
    DirichletCondition[u[x, y] == 0, x == -Lp && y == Lu]}, 
   u, {x, y} \[Element] mesh];

Visualize:

ClearAll[mag];
mag[\[Theta]_, eps_] := 
 Norm@vel[((1 + eps) tk)/2 Cos[\[Theta]], ((1 + eps) tk)/
     2 Sin[\[Theta]]]
Plot[Evaluate@
  Table[mag[\[Theta], 
    eps], {eps, {0, 0.01, 0.02, 0.05, 0.1}}], {\[Theta], -\[Pi]/
   8, \[Pi]/8}, PlotLegends -> LineLegend[{0, 0.01, 0.02, 0.05, 0.1}],
  Frame -> True, FrameLabel -> {"Angular Position", "Velocity"}]