Run this tutorial

Click here to run this tutorial on mybinder.org: try on mybinder.org
Please note that starting the notebook server may take a couple of minutes.

Example: working with gmsh grids¶

1: creating a gmsh file¶

We use pyMORs PolygonalDomain description and discretize_gmsh to obtain a grid file that gmsh can read. Note that dune-grid can only handle gmsh version 2 files. Depending on the installed gmsh version, the file written below may use a newer format; we convert it to a dune-grid-readable version 2.2 mesh with meshio further down (see also these instructions).

# wurlitzer: display dune's output in the notebook
%load_ext wurlitzer


import numpy as np

from pymor.analyticalproblems.domaindescriptions import PolygonalDomain

L_shaped_domain = PolygonalDomain(points=[
    [0, 0],
    [0, 1],
    [-1, 1],
    [-1, -1],
    [1, -1],
    [1, 0],
], boundary_types={'dirichlet': [1, 2, 3, 4, 5, 6]})

from pymor.discretizers.builtin.domaindiscretizers.gmsh import discretize_gmsh

grid, bi = discretize_gmsh(L_shaped_domain, msh_file_path='L_shaped_domain.msh')

2: discretization using pyMOR¶

We may use pyMOR to solve an elliptic PDE on this grid.

from pymor.analyticalproblems.functions import ConstantFunction
from pymor.analyticalproblems.elliptic import StationaryProblem

problem = StationaryProblem(
    domain=L_shaped_domain,
    rhs=ConstantFunction(1, dim_domain=2),
    diffusion=ConstantFunction(1, dim_domain=2),
)

from pymor.discretizers.builtin import discretize_stationary_cg

fom, fom_data = discretize_stationary_cg(problem, grid=grid, boundary_info=bi)

fom.visualize(fom.solve())

3: using the gmsh grid in dune¶

dune-grid only supports gmsh version 2 files, and only a subset of the specification. Depending on the installed gmsh version, discretize_gmsh above may have written a newer (version 4) mesh file, and the file also contains a boundary type ($PhysicalNames) definition which dune-grid cannot correctly parse (we do not require it, we have our own boundary info). We therefore use meshio to re-write a clean version 2.2 ASCII mesh that contains only the points and the triangle cells (Note that you have to provide the same filename here as in the call to discretize_gmsh):

import meshio

_mesh = meshio.read('L_shaped_domain.msh')
meshio.write_points_cells(
    'L_shaped_domain_dune.msh',
    _mesh.points,
    [('triangle', _mesh.get_cells_type('triangle'))],
    file_format='gmsh22',
    binary=False,
)



Warning: Appending zeros to replace the missing physical tag data.

Warning: Appending zeros to replace the missing geometrical tag data.

from dune.xt.grid import make_gmsh_grid, Dim, Simplex, visualize_grid

grid = make_gmsh_grid('L_shaped_domain_dune.msh', Dim(2), Simplex())

This grid can now be used as any other grid, e.g. for visualization …

_ = visualize_grid(grid)

… or discretization:

from dune.gdt import visualize_function

from discretize_elliptic_cg import discretize_elliptic_cg_dirichlet_zero

u_h = discretize_elliptic_cg_dirichlet_zero(grid, diffusion=1, source=1)
_ = visualize_function(u_h)

Download the code: example__gmsh_grid.md example__gmsh_grid.ipynb