Streamfunction of a vector field


I just discovered Freefem and use it for my master degree in bioscience.
I’d like to visualize the streamlines of a 2d vector field so that their isolines must be always orthogonal to the vector field. The “streamfunction” S satisfied the “vorticity” equation : div(grad(S))=-curl(F)[z] since grad(S)=[-Fy,Fx].

Unfortunalty i can’t get streamlines orthogonal to the isolines, but not so far, do I miss something ?
My example is very simple, it’s a square with three holes. First I solve the laplacian equation using :

  • dirichlet bc on borders 1,2,3,4,5
  • zero flux bc on border 6
    Then I solve the streamfunction :
  • vort is the z component of the curl(F)
  • borderDerivative is the normal derivative using the definition grad(S) = [-Fy,Fx]

Thanks for the time to help.

Hélène Ségalier.

PS : I’m not a mathematician (a biologist who learn math sometimes)

Here my edp file

load "gmsh"
load "iovtk"

mesh Th = gmshload("square.msh");

fespace Ph(Th, P2);

int[int] reg = regions(Th);
int[int] tags = labels(Th);


Ph u,v;

problem laplacian(u, v,solver=UMFPACK,eps=1e-5)
= int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))
+ on(1, u=0)
+ on(2, u=0)
+ on(3, u=1)
+ on(4, u=1)
+ on(5, u=1)

Ph Fx,Fy;
Fx = -dx(u);
Fy = -dy(u);

Ph S;
func rotZ = -dx(Fy) + dy(Fx);
func borderDerivative = (-Fy*N.x + Fx*N.y);

problem stream(S,v,solver=UMFPACK,eps=1e-5)
= int2d(Th)(dx(S) * dx(v) + dy(S) * dy(v))
- int2d(Th)(rotZ*v)
- int1d(Th)(borderDerivative*v);


real SMin = S[].min;
cout<<"S min = "<<SMin;
S = S-SMin;

real SMax = S[].max;
cout<<"S max = "<<SMax;
S = S/SMax;

int[int] Order = [1, 1];
string DataName = "u S";
savevtk("export.vtu", Th, u, S, dataname = DataName, order = Order, bin = false);