9
$\begingroup$

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]

Mathematica graphics Mathematica graphics

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]

Mathematica graphics

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"}]

Mathematica graphics

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.

$\endgroup$

1 Answer 1

3
$\begingroup$

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]

enter image description here

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"}]

enter image description here

$\endgroup$
9
  • $\begingroup$ Thanks. I will study this in detail as soon as possible. Is this new in Version 12? $\endgroup$
    – Hugh
    Apr 22, 2019 at 7:57
  • $\begingroup$ @Hugh I think it should work in 11.3. But I have not tried it. $\endgroup$
    – user21
    Apr 22, 2019 at 8:02
  • $\begingroup$ I have examined the solution. What is remarkable is how just a few extra points make a difference. My MeshRefinementFunction adds many points and does not give a good solution. Your few extra points makes a huge difference without adding many points. Also, your replaced 0.2 with 2/10 and this removes a big error in the mesh generation. $\endgroup$
    – Hugh
    Apr 23, 2019 at 11:34
  • $\begingroup$ @Hugh, sorry I should have mentioned the 0.2 vs. 2/10 issue. I'll update the post. $\endgroup$
    – user21
    Apr 23, 2019 at 11:38
  • 1
    $\begingroup$ @user21 Some of my current issues are very similar to the OP's problem here. Can anyone comment on why changing 0.2 to 2/10 makes any difference? After all, they are equal to each other in most people's understanding (including mine). How is Mathematica interpreting them differently, and is there a general principle here. I should say I have noticed people using integer fractions instead of floats in SE answers a lot, but I have never seen this referred to explicitly in MMA documentation. I read once somewhere that MMA does not distinguish between different number types. $\endgroup$ Apr 21, 2020 at 10:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.