Download :download:`this testcase <2d_poiseuille_tfem_convergence.py>`.

Bidimensional Poiseuille Flow — TFEM Convergence Test
=====================================================
This example evaluates the convergence of a Two-Field Finite Element Method
(TFEM) applied to a 2D Poiseuille flow. The domain is meshed using Gmsh and
solved with MigFlow’s fluid module. Two mesh configurations are compared:
-  a regular structured mesh
-  a perturbed mesh containing isolated triangular "caps".

Keywords
--------
TFEM, FEM, Gmsh, Convergence, Poiseuille flow, Mesh perturbation

Description
-----------
This benchmark demonstrates:
- How to generate a 2D rectangular channel mesh with named boundaries in Gmsh.
- How to introduce local mesh perturbations ("caps") to assess solver robustness.
- How to impose an analytical Poiseuille velocity profile at the inlet.
- How to solve a steady-state viscous flow using MigFlow’s TFEM formulation.
- How to measure the numerical error with respect to the analytical solution.

Physical Setup
--------------
The flow corresponds to classical 2D Poiseuille flow driven by a pressure
difference. The analytical velocity profile is:

  u_x(y) = (1 / (20 μ)) * y * (1 − y)

Boundary conditions:
- Left boundary: imposed parabolic inflow velocity profile.
- Right boundary: zero relative pressure (open outlet).
- Top/Bottom: no-slip solid walls (velocity = 0).

Numerical parameters:
- ρ = 1000 kg/m³  (fluid density)
- ν = 1e−3 m²/s   (kinematic viscosity)
- μ = ρν          (dynamic viscosity)
- dt = 10         (pseudo–time step)
- tEnd = 1000     (final pseudo–time)

Convergence Study
-----------------
The test loops over several mesh resolutions N = {10, 20, 40, 80, 160}.
For each resolution, two simulations are run:
- with caps (perturbed mesh)
- without caps (regular mesh)

At the end of each simulation, the L² error norm between numerical and
analytical velocity is computed:

  error = || u_num − u_exact ||₂

The code then plots the error in log–log scale to assess convergence order
and to compare robustness with/without mesh perturbations.

-----------------------------------------------------------------------------

.. code-block:: python
  
  from migflow import fluid, time_integration
  import numpy as np
  import os
  import time
  import shutil
  import gmsh
  
  # Geometry
  # --------
  # The computational domain is a 2D channel of height H = 1 and length L = 1.
  # Depending on configuration:
  # - A single rectangle is generated for the regular mesh.
  # - Two adjacent rectangles are generated and remeshed to insert triangular caps
  #   in the perturbed case.
  #
  # Geometrical parameters
  H = 1  # domain height
  L = 1  # domain width
  Ox = -5
  Oy = 0
  
  
  # Mesh Generation
  # ---------------
  # Mesh generation relies on Gmsh with:
  # - transfinite curves along top/bottom boundaries to control refinement;
  # - automatic Delaunay meshing (Mesh.Algorithm = 6);
  # - optional internal triangular elements (“caps”) that locally distort the mesh
  #   for convergence robustness testing.
  # Meshes are exported to .msh files and loaded by MigFlow.
  def mesh(
      N=100, L0=0.5, D=0, withcaps=True, needles=False, one_side=False, special=False
  ):
      L = 2 * L0
      gmsh.model.mesh.set_size_callback(lambda *x: L / (N))
  
      def make_square(x0, y0, x1, y1):
          xs = [(x0, y0), (x0, y1), (x1, y1), (x1, y0)]
          p = list(gmsh.model.geo.add_point(x[0], x[1], 0) for x in xs)
          l = list(gmsh.model.geo.add_line(p0, p1) for (p0, p1) in zip(p, np.roll(p, 1)))
          s = gmsh.model.geo.add_plane_surface([gmsh.model.geo.add_curve_loop(l)])
          return s, l
  
      if not withcaps:
          s0, ls0 = make_square(0, 0, 4 * L, L)
          gmsh.model.geo.synchronize()
          gmsh.model.add_physical_group(1, [ls0[1]], -1, "Left")
          gmsh.model.add_physical_group(1, [ls0[2]], -1, "Top")
          gmsh.model.add_physical_group(1, [ls0[3]], -1, "Right")
          gmsh.model.add_physical_group(1, [ls0[0]], -1, "Bottom")
          gmsh.model.add_physical_group(2, [s0], -1, "domain")
          gmsh.option.set_number("Mesh.Algorithm", 6)
          gmsh.model.mesh.generate(2)
          # gmsh.fltk.run()
          return
      s0, ls0 = make_square(0, 0, 4 * L0, L)
      s1, ls1 = make_square(4 * L0 + D, 0, 4 * L, L)
      gmsh.model.geo.synchronize()
      l0 = gmsh.model.get_entities_in_bounding_box(L, 0, 0, L, L, 0, 1)
      gmsh.model.mesh.set_transfinite_curve(ls0[3], N + 1)
      if needles:
          gmsh.model.mesh.set_transfinite_curve(ls1[1], N + 1)
      elif one_side:
          gmsh.model.mesh.set_transfinite_curve(ls1[1], 2 * N + 1)
      elif special:
          gmsh.model.mesh.set_transfinite_curve(ls1[1], 3 * N + 2)
      else:
          gmsh.model.mesh.set_transfinite_curve(ls1[1], N + 2)
      gmsh.option.set_number("Mesh.Algorithm", 6)
      gmsh.model.mesh.generate(1)
      l1tag, l1, _ = gmsh.model.mesh.get_nodes(1, ls1[1], False, False)
      l0tag, _, _ = gmsh.model.mesh.get_nodes(1, ls0[3], False, False)
      l1 = l1.reshape(-1, 3)
      if needles:
          l1[:, 1] = np.linspace(L / N, L * (1 - 1 / N), N - 1)[::-1]
      elif one_side:
          l1[:, 1] = np.linspace(L / (2 * N), L - (L / (2 * N)), 2 * N - 1)[::-1]
      elif special:
          l1[1 : 3 * N : 3, 1] = np.linspace(L / (2 * N), L - (L / (2 * N)), N)[::-1]
          l1[0 : 3 * N - 1 : 3, 1] = l1[1 : 3 * N : 3, 1] + D
          l1[2 : 3 * N + 1 : 3, 1] = l1[1 : 3 * N : 3, 1] - D
      else:
          l1[:, 1] = np.linspace(L / (2 * N), L - (L / (2 * N)), N)[::-1]  #
      for t, x in zip(l1tag, l1):
          gmsh.model.mesh.set_node(t, x, [])
      # gmsh.model.mesh.synchronize()
      gmsh.model.mesh.generate(2)
      con = gmsh.model.add_discrete_entity(2)
      l1, _, l1param = gmsh.model.mesh.get_nodes(1, ls1[1], True, True)
      l0, _, l0param = gmsh.model.mesh.get_nodes(1, ls0[3], True, True)
      l1 = l1[np.argsort(l1param)]
      l0 = l0[np.argsort(-l0param)]
      if needles:
          tri = np.zeros((N, 3), np.int32)
          tri[:, 1] = l1[:-1]
          tri[:, 0] = l1[1:]
          tri[:, 2] = l0[:-1]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
          tri2 = np.zeros((N, 3), np.int32)
          tri2[:, 1] = l0[:-1]
          tri2[:, 0] = l1[1:]
          tri2[:, 2] = l0[1:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
      elif one_side:
          tri = np.zeros((N, 3), np.int32)
          tri[:, 0] = l1[0 : 2 * N - 1 : 2]
          tri[:, 1] = l1[1 : 2 * N : 2]
          tri[:, 2] = l0[:-1]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
          tri1 = np.zeros((N, 3), np.int32)
          tri1[:, 0] = l1[1 : 2 * N : 2]
          tri1[:, 1] = l1[2 : 2 * N + 1 : 2]
          tri1[:, 2] = l0[1:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri1.flat)
          tri2 = np.zeros((N, 3), np.int32)
          tri2[:, 0] = l0[:-1]
          tri2[:, 1] = l1[1 : 2 * N : 2]
          tri2[:, 2] = l0[1:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
      elif special:
          tri = np.zeros((N + 1, 3), np.int32)
          tri[:, 0] = l1[0 : 3 * N + 1 : 3]
          tri[:, 1] = l1[1 : 3 * N + 2 : 3]
          tri[:, 2] = l0[:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
          tri1 = np.zeros((N, 3), np.int32)
          tri1[:, 0] = l1[1 : 3 * N - 1 : 3]
          tri1[:, 1] = l1[2 : 3 * N : 3]
          tri1[:, 2] = l0[:-1]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri1.flat)
          tri2 = np.zeros((N, 3), np.int32)
          tri2[:, 0] = l0[:-1]
          tri2[:, 1] = l1[2 : 3 * N : 3]
          tri2[:, 2] = l0[1:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
          tri3 = np.zeros((N, 3), np.int32)
          tri3[:, 0] = l1[2 : 3 * N : 3]
          tri3[:, 1] = l1[3 : 3 * N + 1 : 3]
          tri3[:, 2] = l0[1:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri3.flat)
      else:
          tri = np.zeros((N + 1, 3), np.int32)
          tri[:, 0] = l1[:-1]
          tri[:, 2] = l1[1:]
          tri[:, 1] = l0
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
          tri2 = np.zeros((N, 3), np.int32)
          tri2[:, 0] = l0[:-1]
          tri2[:, 2] = l1[1:-1]
          tri2[:, 1] = l0[1:]
          gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
      conl_top = gmsh.model.add_discrete_entity(1)
      conl_bottom = gmsh.model.add_discrete_entity(1)
      gmsh.model.mesh.add_elements_by_type(conl_top, 1, [], [l1[0], l0[0]])
      gmsh.model.mesh.add_elements_by_type(conl_bottom, 1, [], [l0[-1], l1[-1]])
      gmsh.model.add_physical_group(1, [ls0[1]], -1, "Left")
      gmsh.model.add_physical_group(1, [ls0[2], ls1[2], conl_top], -1, "Top")
      gmsh.model.add_physical_group(1, [ls1[3]], -1, "Right")
      gmsh.model.add_physical_group(1, [ls0[0], ls1[0], conl_bottom], -1, "Bottom")
      gmsh.model.add_physical_group(2, [s0, s1, con], -1, "domain")
      # gmsh.fltk.run()
  
  
  def mesh_isolated_caps(N=100, L=1, D=0, N_caps=-1):
      gmsh.model.mesh.set_size_callback(lambda *x: L / (N))
  
      def make_square(x0, y0, x1, y1):
          xs = [(x0, y0), (x0, y1), (x1, y1), (x1, y0)]
          p = list(gmsh.model.geo.add_point(x[0], x[1], 0) for x in xs)
          l = list(gmsh.model.geo.add_line(p0, p1) for (p0, p1) in zip(p, np.roll(p, 1)))
          s = gmsh.model.geo.add_plane_surface([gmsh.model.geo.add_curve_loop(l)])
          return s, l
  
      s0, ls0 = make_square(0, 0, 4 * L, L)
      gmsh.model.geo.synchronize()
      gmsh.option.set_number("Mesh.Algorithm", 6)
      gmsh.model.mesh.generate(2)
      if N_caps < 0:
          N_caps = (N // 2) * 5
  
      nodestag, nodes, nodesparam = gmsh.model.mesh.get_nodes()
      nodesorder = np.argsort(nodestag)
      nodes = nodes.reshape(-1, 3)[:, :]
  
      [etype], [etags], [elements] = gmsh.model.mesh.get_elements(2)
      elements = elements.reshape(etags.shape[0], -1)
      n = elements.shape[1]
      elements = np.searchsorted(nodestag, elements, sorter=nodesorder)
  
      _, bnd_edges_top = gmsh.model.mesh.get_elements_by_type(1, ls0[2])
      _, bnd_edges_bottom = gmsh.model.mesh.get_elements_by_type(1, ls0[0])
      _, bnd_edges_left = gmsh.model.mesh.get_elements_by_type(1, ls0[1])
      _, bnd_edges_right = gmsh.model.mesh.get_elements_by_type(1, ls0[3])
  
      random_iel = np.random.choice(elements.shape[0], N_caps, replace=False)
      random_el = elements[random_iel]
      new_nodes = np.mean(nodes[random_el[:, :2]], axis=1)
      dists = np.linalg.norm(nodes[random_el[:, 2]] - new_nodes, axis=1)[:, None]
      Ds = np.minimum(D, dists / 2)
      new_nodes += Ds * (nodes[random_el[:, 2]] - new_nodes) / dists
      final_nodes = np.vstack([nodes, new_nodes])
  
      new_ids = nodes.shape[0] + np.arange(N_caps)
      new_elements = np.array([random_el[:, 0], random_el[:, 1], new_ids])
      new_elements2 = np.array([random_el[:, 1], random_el[:, 2], new_ids])
      new_elements3 = np.array([random_el[:, 2], random_el[:, 0], new_ids])
      new_elements = np.vstack([new_elements.T, new_elements2.T, new_elements3.T])
  
      elements = np.delete(elements, random_iel, axis=0)
      final_elements = np.vstack([elements, new_elements]) + 1
  
      gmsh.model.mesh.clear()
      # Inject transformed mesh
      gmsh.model.mesh.add_nodes(
          2, 1, np.arange(final_nodes.shape[0]) + 1, final_nodes.flatten()
      )
      # modified
      dom_entity = gmsh.model.add_discrete_entity(2)
      gmsh.model.mesh.add_elements_by_type(
          dom_entity, 2, np.arange(final_elements.shape[0]) + 1, final_elements.flatten()
      )
  
      bnd_entity_top = gmsh.model.add_discrete_entity(1)
      gmsh.model.mesh.add_elements_by_type(bnd_entity_top, 1, [], bnd_edges_top.flatten())
      gmsh.model.add_physical_group(1, [bnd_entity_top], 10, "Top")
      bnd_entity_bottom = gmsh.model.add_discrete_entity(1)
      bnd_entity_left = gmsh.model.add_discrete_entity(1)
      bnd_entity_right = gmsh.model.add_discrete_entity(1)
  
      gmsh.model.mesh.add_elements_by_type(
          bnd_entity_bottom, 1, [], bnd_edges_bottom.flatten()
      )
      gmsh.model.add_physical_group(1, [bnd_entity_bottom], 11, "Bottom")
  
      gmsh.model.mesh.add_elements_by_type(
          bnd_entity_left, 1, [], bnd_edges_left.flatten()
      )
      gmsh.model.add_physical_group(1, [bnd_entity_left], 12, "Left")
  
      gmsh.model.mesh.add_elements_by_type(
          bnd_entity_right, 1, [], bnd_edges_right.flatten()
      )
      gmsh.model.add_physical_group(1, [bnd_entity_right], 13, "Right")
  
      gmsh.model.add_physical_group(2, [dom_entity], -1, "domain")
  
  

Run the simulation for different mesh configurations
----------------------------------------------------

.. code-block:: python
  
  def run(N, withcaps):

  # Output Directory
  # ----------------
  # Create a clean output directory for simulation results.

.. code-block:: python
  
      outputdir = "output"
      shutil.rmtree(outputdir, True)
      if not os.path.isdir(outputdir):
          os.makedirs(outputdir)

  # Mesh Generation
  # ---------------
  # Depending on the `withcaps` flag, either a regular mesh or a perturbed
  # mesh with triangular caps is generated and saved to .msh files.

.. code-block:: python
  
      if withcaps:
          mesh_isolated_caps(N=N)
          filename = outputdir + f"/mesh_{N}_caps.msh"
      else:
          mesh(N=N, D=0, withcaps=False)
          filename = outputdir + f"/mesh_{N}_nocaps.msh"
      gmsh.write(filename)

  # Physical parameters
  # -------------------

.. code-block:: python
  
      g = 0  # gravity
      rho = 1000  # fluid density
      nu = 1e-3  # kinematic viscosity
      mu = nu * rho  # dynamic viscosity
  

  # Numerical parameters
  # --------------------

.. code-block:: python
  
  
      #
      # FLUID PROBLEM
      #
      # Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
      f = fluid.FluidProblem(2, g, mu, rho)
      f.load_msh(filename)
      f.set_open_boundary(
          "Left", velocity=[lambda x: 1 / (20 * mu) * x[:, 1] * (1 - x[:, 1]), 0]
      )
      f.set_wall_boundary("Bottom", velocity=[0, 0])
      f.set_wall_boundary("Top", velocity=[0, 0])
      f.set_strong_boundary("Bottom", velocity=[0, 0])
      f.set_strong_boundary("Top", velocity=[0, 0])
      f.set_open_boundary("Right", pressure=0)
  

  # Simulation Parameters
  # ---------------------

.. code-block:: python
  
      t = 0  # initial time
      i = 0  # iteration counter
      outf = 1  # number of iterations between output files
      dt = 10  # time step
      tEnd = 500  # final time
  

  # Simulation Loop
  # ---------------

.. code-block:: python
  
      while t < tEnd:
          print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
          time_integration.iterate(f, None, dt)
          t += dt
          # Output files writting
          if i % outf == 0:
              f.write_mig(outputdir, t)
          i += 1

  # Error Computation
  # -----------------
  # After reaching the final time, the error at the nodes is computed between the numerical
  # and analytical velocity profiles to assess convergence.

.. code-block:: python
  
      s = f.velocity()
      x = f.coordinates_fields()[f.velocity_index()[:, 0]]
      vel = (s[:, 0] - 1 / (20 * nu * rho) * x[:, 1] * (1 - x[:, 1])) ** 2
      vS = np.sum((1 / (20 * nu * rho) * x[:, 1] * (1 - x[:, 1])) ** 2)
      print(f"Error criteria : {(vel.sum())**.5} \t < {(vS**0.5)/50}")
      print(f"Succeed ? {(vel.sum())**.5 < (vS**0.5)/50}")
      gmsh.model.remove()
      return vel.sum() ** 0.5
  
  

Main Convergence Loop
---------------------
The main script loops over different mesh resolutions and both mesh
configurations (with/without caps), collecting error norms for convergence
analysis.

.. code-block:: python
  
  if __name__ == "__main__":
      Ns = np.array([10, 20, 40, 80, 160])
      data_caps = []
      data_nocaps = []
      plot = False
      for N in Ns:
          for withcaps in [False, True]:
              print(f"N = {N}, withcaps = {withcaps}")
              sol = run(N, withcaps)
              if withcaps:
                  data_caps.append(sol)
              else:
                  data_nocaps.append(sol)
      if plot:
          import matplotlib.pyplot as plt
  
          plt.figure()
          # log log plot
          plt.loglog(1.0 / Ns, data_caps, label="with caps", marker="o")
          plt.loglog(1.0 / Ns, data_nocaps, label="without caps", marker="o")
          plt.legend()
          plt.show()
  

Plot
----
.. code-block:: shell

 python3 -m migflow.plot.migplot output --actors fluid --fluid-field velocity --fluid-vmin 0 --fluid-vmax 0.0125

