5
$\begingroup$

I am trying to plot some streamlines, where the domain is given by a mesh. The results are good, but sometimes there are streamlines, that do not end at the boundary of the domain and are plotted a little bit outside. Here is MVE:

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

\[Alpha] = Pi/4;
Q = {{Cos[\[Alpha]], -Sin[\[Alpha]]}, {Sin[\[Alpha]], Cos[\[Alpha]]}}; 
d = Q.{1, 0};

\[CapitalOmega] = RegionDifference[Rectangle[{0, 0}, {1, 1}], Ellipsoid[{1/2, 1/2}, {1/6, 1/50}]];
mesh = ToElementMesh[\[CapitalOmega]];
StreamPlot[d, {x, y} \[Element] mesh, StreamPoints -> Coarse, StreamMarkers -> "Segment"]

This produces the following plot: MVE As you can see, there is a small line inside the ellipse. Is there a way, how to remove it?

EDIT Another example:

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

\[Alpha] = Pi/4*x;
Q = {{Cos[\[Alpha]], -Sin[\[Alpha]]}, {Sin[\[Alpha]], Cos[\[Alpha]]}}; 
vars = {{u[x, y], v[x, y]}, {x, y}};
pars = <|"MaterialModel" -> "NeoHookeanIsotropic", "LameParameter" -> 10^9, "ShearModulus" -> 5*10^8, "Thickness" -> 1|>;

\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];
mesh = ToElementMesh[\[CapitalOmega]];
pde := {SolidMechanicsPDEComponent[vars, pars] == SolidBoundaryLoadValue[x == 1, vars, pars, <|"Pressure" -> {p, 0}|>], SolidFixedCondition[x == 0, vars, pars]};
displacement = NDSolveValue[pde /. p -> 3*10^8, {u[x, y], v[x, y]}, {x, y} \[Element] \[CapitalOmega]];

deformedMesh = ElementMeshDeformation[mesh, displacement, "ScalingFactor" -> 1];
dispGrad = Map[Grad[#, {x, y}] &, displacement];
direction = Partition[ElementMeshInterpolation[deformedMesh, #["ValuesOnGrid"]][x, y] & /@ Flatten[dispGrad][[All, 0]], 2] . Q . {1, 0} + Q.{1, 0};
StreamPlot[direction, {x, y} \[Element] deformedMesh, StreamPoints -> Coarse, StreamMarkers -> "Segment"]

This produces the following plot. One streamline in the lower right corner starts outside of the domain. another example

$\endgroup$
5
  • $\begingroup$ The second example doesn't work for me -- I get FindRoot::jsing1, FindRoot::sszero and Set::shape errors and it shows a square region. $\endgroup$
    – Chris K
    Apr 15, 2023 at 23:22
  • $\begingroup$ ^^ is the second example working for everyone else? $Version = "13.2.1 for Mac OS X ARM (64-bit) (January 27, 2023)" $\endgroup$
    – Chris K
    Apr 16, 2023 at 3:50
  • $\begingroup$ I think you should report these issues to WRI. $\endgroup$
    – user21
    Apr 17, 2023 at 4:43
  • $\begingroup$ @ChrisK, MKL ( hence Paradiso) does not work on ARM. Use Rosetta . $\endgroup$
    – user21
    Apr 17, 2023 at 4:44
  • $\begingroup$ @ChrisK, even though it's a known limitation, you could report it to create some traction for this issue. $\endgroup$
    – user21
    Apr 17, 2023 at 4:48

3 Answers 3

3
$\begingroup$
  • At first, we use reg=MeshRegion@deformedMesh,not only deformedMesh.
  • We use SignedRegionDistance to define a RegionFunction to mask the outside of the region MeshRegion@deformedMesh.
reg = MeshRegion@deformedMesh;
streamplot = 
  StreamPlot[direction, {x, y} ∈ reg, StreamPoints -> Coarse,
    StreamMarkers -> "Segment"];
dist = SignedRegionDistance[reg];
masking = 
  RegionPlot[dist@{x, y} >= 10^-3, {x, 0, 1.25}, {y, 0, 1.2}, 
   PlotStyle -> White, BoundaryStyle -> None];
boundary = RegionPlot[RegionBoundary[reg], PlotStyle -> Blue];
Show[streamplot, masking, boundary]

enter image description here

* enter image description here

$\endgroup$
1
  • $\begingroup$ Looks like this works. Thank you. $\endgroup$
    – lemurman
    Apr 16, 2023 at 11:11
5
$\begingroup$

Just change your mesh definition to

mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> .1];

enter image description here

$\endgroup$
6
  • $\begingroup$ Thank you for your answer, but the different mesh causes me some problems in other parts of the code (mainly with ElementMeshDeformation). Do you think it is possible to achieve the same result without changing the mesh? $\endgroup$
    – lemurman
    Apr 15, 2023 at 14:42
  • $\begingroup$ Is option Coarse necessary? $\endgroup$ Apr 15, 2023 at 16:04
  • $\begingroup$ It is not, but it is appreciated. I use the StreamPlot to plot direction of fibers in material and it looks much better, if the number of streamlines is lower. $\endgroup$
    – lemurman
    Apr 15, 2023 at 16:14
  • $\begingroup$ Aren't there easier ways to draw a group of straight lines? $\endgroup$ Apr 15, 2023 at 16:34
  • $\begingroup$ Straight lines in a square? Probably yes, but this is just a MVE. In the more complex cases, the lines are not necessarily straight and the domain is more complicated. I use Mathematica's solid mechanics framework to model deformations o fiber-reinforced materials (only in 2D). After I solve the equations, I would like to plot the fibers, for which I use the StreamPlot. The fibers are not alway straight and the domain is given by the deformed material (I use ElementMeshDeformation to deform the mesh). $\endgroup$
    – lemurman
    Apr 15, 2023 at 18:36
3
$\begingroup$

Instead of making StreamPlot deal with the Region, how about just covering up the excluded region afterwards?

Show[
  StreamPlot[d, {x, 0, 1}, {y, 0, 1}, StreamPoints -> Coarse, StreamMarkers -> "Segment"],
  Graphics[{White, EdgeForm[Blue], Ellipsoid[{1/2, 1/2}, {1/6, 1/50}]}]
]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you for your answer. This works only in this specific situation. It doesn't solve more complicated problems. I have added another example into the post. $\endgroup$
    – lemurman
    Apr 15, 2023 at 21:57

Your Answer

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

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