# Streamlines using the Elmer method

**URL:** <https://community.freefem.org/t/streamlines-using-the-elmer-method/3407>\
**Category:** General Discussion\
**Created:** [July 20, 2024, 3:56pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407 "2024-07-20T15:56:32Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![PierreAntonetti](https://avatars.discourse-cdn.com/v4/letter/p/f1d935/32.png) [@PierreAntonetti](https://community.freefem.org/u/PierreAntonetti)\
**Post date:** [July 20, 2024, 3:56pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/1 "2024-07-20T15:56:32Z")

</div>

Hello,

I want to get the streamline function \psi of a 2d incompressible field U=(u,v) cause the isoline of \psi are the streamlines (orthogonal to U).

I’ve seen two related topics on this forum mainly using the vorticity equation, and it seems to no working.

I’m a new Freefem user because I come from ElmerFEM. Based on the Elmer manual, the “Streamline” ElmerFEM module use the equation (46.2). I don’t understand how to implement this using Freefem cause in the weak formulation the space test function of \vec{v} is twice bigger than the space of \psi

Thank you very much for helping

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/2/20cc5040fcb8e8776f9a90f3114169c30e563d7a.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 20, 2024, 4:52pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/2 "2024-07-20T16:52:51Z")

</div>

A method is proposed in

> [@Uzawa scheme for stokes bingham flow](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/16):
>
> In that paper the flow is taken incompressible, whereas in your code you use compressible. This is why the solution is different (but looks correct). You can modify easily your code to do the incompressible case, by using the extra variables p,q as solve bing([ux,uy,p],[vx,vy,q]) and adding -int2d(Th)(p\*Divergence(vx, vy)+q\*Divergence(ux, uy)) But then it happens that the scheme converges extremely slowly (you need thousands of iterations). In order to improve it you should use the Augmente…

from “Another point is that your computation of the stream function…”

The important thing is the \varepsilon.  
If your u is P1, you should take the stream function in P2.  
If your u is P2, you should take the stream function in P3.

---

<div class="post-metadata">

**Author:** ![PierreAntonetti](https://avatars.discourse-cdn.com/v4/letter/p/f1d935/32.png) [@PierreAntonetti](https://community.freefem.org/u/PierreAntonetti)\
**Post date:** [July 20, 2024, 7:39pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/3 "2024-07-20T19:39:04Z")

</div>

Hi @fb77

Thank you for your answer !

I tried to implement your method on a classic thermal conduction problem (laplace and dirichlet bc), and it seems quite good, but not enough to get orthogonal streamlines with respect to the temperature.  
I played with the penalty coefficient from 1e-2 to 1e-12 without any success. The P1 flux field = -grad(T) is derived from a P2 scalar field, so i used a P2 streamline scalar field on Th.  
As you can see on the picture, the streamline (continuous black lines are not aligned with the vectors or not orthogonal to the temperature surface field).

Any idea ?

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/ee339ebd7efde2c114bc5bfec077d66d17bd3ca9.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 21, 2024, 11:52am UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/4 "2024-07-21T11:52:53Z")

</div>

Can you upload your full freefem test code?

---

<div class="post-metadata">

**Author:** ![PierreAntonetti](https://avatars.discourse-cdn.com/v4/letter/p/f1d935/32.png) [@PierreAntonetti](https://community.freefem.org/u/PierreAntonetti)\
**Post date:** [July 21, 2024, 3:17pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/5 "2024-07-21T15:17:58Z")

</div>

Hi,  
Here the zip file containing the mesh and edp file

Very strange, some streamlines are very good (on the top and on left) and others (bottom)^are not aligned with the vector flux…  
Legend : fluxes (on the isostreamline meshes) are back and iso temperature are in rianbow.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/7/791c41d66e76beb7be7ae2e020fafdb5606179ca.jpeg)

Thans you !

[main.zip](https://community.freefem.org/uploads/short-url/5LFJhhkEmNFADR1HC9mResPz2Ii.zip) (143.8 KB)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 21, 2024, 5:12pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/6 "2024-07-21T17:12:14Z")

</div>

I cannot see the result with Paraview, but I think that the following change can improve: integrate by parts to replace two integrals by one (result is shown below).  
Indeed in the `int1d` there could be some internal boundaries that separate the regions. These should not be taken in the formulation.

```auto
// Solve the streamline and normalise
real coef=1e-8;
P2h psi,psitest;
solve streamlines (psi, psitest, solver=UMFPACK)
    = int2d(Th)(
          dx(psi)*dx(psitest)
        + dy(psi)*dy(psitest)
    )
// + int2d(Th)(- psitest*(dy(fluxX) - dx(fluxY)))
// + int1d(Th)(- psitest*(fluxY*N.x-fluxX*N.y))
    + int2d(Th)(dy(psitest)*fluxX - dx(psitest)*fluxY )
    + int1d(Th)(coef*psi*psitest)
    ;

```

About the `int1d` that includes `coef` there could be a similar problem.  
There you can eventually replace `int1d` by `int2d`. This determines the additive constant so that \int\_\Omega\psi=0 instead of \int\_{\partial\Omega}\psi=0.

---

<div class="post-metadata">

**Author:** ![PierreAntonetti](https://avatars.discourse-cdn.com/v4/letter/p/f1d935/32.png) [@PierreAntonetti](https://community.freefem.org/u/PierreAntonetti)\
**Post date:** [July 22, 2024, 9:14pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/7 "2024-07-22T21:14:48Z")

</div>

Hi,  
I just tried your solution and it works a little bit better !

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c004502df9540712a47dd8b0d85dc9570d4b02d0.png)

Few remarks :

- replacing \int\_{\Omega}\psi =0 by \int\_{\partial\Omega}\psi=0 doesn’t change anything
- sometimes for “large” model the streamlines are not aligned with the flux vector (see picture below, example with a ground)

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/f/faa8385c267071a51ae4f0484f486642dc89ffaa.png)

Best regards.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 23, 2024, 10:52am UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/8 "2024-07-23T10:52:54Z")

</div>

It looks rather nice! There is always a discrepancy because the free divergence condition on the vector field (here `[fluxX,fluxY]`), that is satisfied only in the weak sense. This is inherent to the discretization by finite elements of the Laplace equation (I think it is the case for any numerical method).

In your picture below it is quite surprising to have lots of discrepancies. The fact is that there are some holes in the domain, which may lead to a problem of the definition of the stream function psi. In terms of differential forms, the vector field is “closed”, but maybe not “exact”. The condition to be exact writes here  
\int\_C \vec U\cdot\vec n dl=0,  
where \vec U is the free divergence vector field, \vec n is the normal to the closed loop C, and dl is the length element along C. This must hold for any closed loop in the domain, in particular for any loop around each hole. If this is not satisfied, if follows that there exist no function \psi satisfying \mathop{\rm curl}\psi=\vec U in the domain. The variational formulation always gives a solution \psi, but in this case it does not satisfy \mathop{\rm curl}\psi=\vec U.

An explicit example of such pathology is in the annulus 1\<r\<2 (with r=\sqrt{x^2+y^2}) with temperature T=\log r. It satisfies \Delta T=0. We have \vec U=\nabla T. There is no stream function \psi because it should be the argument of (x,y), but we cannot define it as a continuous function in the annulus. We have \int\_C \vec U\cdot\vec n dl=2\pi .The solution to the variational problem “streamlines” is \psi=0.

---

<div class="post-metadata">

**Author:** ![PierreAntonetti](https://avatars.discourse-cdn.com/v4/letter/p/f1d935/32.png) [@PierreAntonetti](https://community.freefem.org/u/PierreAntonetti)\
**Post date:** [July 24, 2024, 11:22am UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/9 "2024-07-24T11:22:53Z")

</div>

Hello,

Yes, I do confirm the discrepancy only occurs when there are holes in the geometry otherwise it works very well (I did the same example with a ground but without holes) and the streamlines are aligned to the vector field

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/8/89d03c75bbc4c32fbb9c8795c1c8abcafc4078dc.png)

So, as far as I understand, it’s not possible to correct the weak formulation with holes ?

Thanl you

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [July 24, 2024, 1:51pm UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/10 "2024-07-24T13:51:22Z")

</div>

The problem with holes arises at the level of PDEs, it is not a problem of discrete formulation. However, a way to bypass the difficulty is to solve the streamlines problem over a modified domain, as you did (but the temperature problem can still be solved on the original mesh). We have just to cut small parts so that there remains no hole.  
In the example of the annulus, one can cut the domain by considering the subdomain where the angle \theta in polar coordinates is such that |\theta|\<\pi-\epsilon. Solving the `streamlines` formulation over this subdomain will give a correct stream function (with a jump across the cut part). The key point is to allow \psi to jump somewhere when making a loop around a hole.

---

<div class="post-metadata">

**Author:** ![PierreAntonetti](https://avatars.discourse-cdn.com/v4/letter/p/f1d935/32.png) [@PierreAntonetti](https://community.freefem.org/u/PierreAntonetti)\
**Post date:** [July 25, 2024, 9:23am UTC](https://community.freefem.org/t/streamlines-using-the-elmer-method/3407/11 "2024-07-25T09:23:29Z")

</div>

Hello,

I do understand the differential form is not exact due to \Omega is not a star domain, and “cutting” into small parts is not an easy task for a “general” geometry.

I just wonder myself how the cfd community manage to display the velocity streamlines when dealing with domains with obstacles (cf. picture below) ?

![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/d5e8cf5550540102720c27573d809879e47776c9.jpeg)

Best regards.