# Streamfunction of a vector field

**URL:** https://community.freefem.org/t/streamfunction-of-a-vector-field/1740
**Category:** General Discussion
**Created:** [May 13, 2022, 6:37pm UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740 "2022-05-13T18:37:00Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![Helene](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/helene/32/1145_2.png) [@Helene](https://community.freefem.org/u/Helene)
#### Post date: [May 13, 2022, 6:37pm UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740/1 "2022-05-13T18:37:00Z")

</div>

Hello,

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

```auto
load "gmsh"
load "iovtk"

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

fespace Ph(Th, P2);

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

cout<<reg<<endl;
cout<<tags<<endl;

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)
;
laplacian;

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);

stream;

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);

```

---

<div class="post-metadata">

### Author: ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [May 16, 2022, 11:44am UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740/2 "2022-05-16T11:44:53Z")

</div>

You can see in example  
[cavityNewton.edp](https://community.freefem.org/uploads/short-url/gcUQYe4igZgesZxIshJOWm8oZX0.edp) (2.8 KB)

TH boundary condition on the `

```auto
// Problem stream-lines (with solve)
solve streamlines (psi, phi)
	= int2d(Th)(
		  dx(psi)*dx(phi)
		+ dy(psi)*dy(phi)
	)
	+ int2d(Th)(
		- phi*(dy(u1) - dx(u2))
	)
	+ on(1, 2, 3, 4, psi=psiBC);

```

where psiBC is the psiBC(x(s),y(s)) = \int\_0^s U.n ds   
where s is a curve abcisse on border and U is the flow field.

---

<div class="post-metadata">

### Author: ![Helene](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/helene/32/1145_2.png) [@Helene](https://community.freefem.org/u/Helene)
#### Post date: [May 17, 2022, 11:39am UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740/3 "2022-05-17T11:39:59Z")

</div>

Hello @frederichecht,

I really want to thank you for the time you took to answer.

- I’ve already showed this example cavityNewton, in that example the velocity field is normal to the borders (slip condition), so that the boundary integral is 0 and psi=0 on the boundaries.

- In my case (a square with 3 holes) that integral is not obvious and i don’t know how to compute it with freefem. Moreover one must choose an arbitrary origin point on the border ?

- I’m not sure to understand why you use a dirichlet bc instead of using the border integral in the weak formulation as follow :  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b886de8be83ee0590b0f5d574450d630bd0343a3.png)  
Could you explain ?

Best regards.

---

<div class="post-metadata">

### Author: ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [May 17, 2022, 11:58am UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740/4 "2022-05-17T11:58:53Z")

</div>

1)Because in this case F.T = 0 so the integral is zero . (T == tangent )  
2) the Dirichlet BC. are more precise , and you can see recirculation zone close to the corner.

So I make a mistake on the previous mail the psiBC(s) = int\_0^1 F.T not int\_0^1 F.N

Best Regards,

Frédéric Hecht.

---

<div class="post-metadata">

### Author: ![Helene](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/helene/32/1145_2.png) [@Helene](https://community.freefem.org/u/Helene)
#### Post date: [May 17, 2022, 12:23pm UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740/5 "2022-05-17T12:23:36Z")

</div>

-Yes F.T=0 !

- What do you mean by recirculation ?  
I cant see how to compute psiBC for each node on the border to prescribe the dirichlet bc on(gamma, psi=psiBC)

Thanks for your support.

---

<div class="post-metadata">

### Author: ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [May 17, 2022, 1:01pm UTC](https://community.freefem.org/t/streamfunction-of-a-vector-field/1740/6 "2022-05-17T13:01:31Z")

</div>

Of coarse, if you do not know the way to compute de Dirichlel BC you can use Neuman BC.