Run this tutorial

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

Example: mixed-element and prism grids

Most grid managers hold a single geometry type: YaspGrid is all cubes, ALUGrid is all cubes or all simplices. UGGrid is the exception – it can hold a mesh that mixes cubes and simplices (and, in 3d, prisms). The dune-gdt C++ test suite exercises those element types through the make_mixed_grid / make_prism_grid fixtures (dune/gdt/test/spaces/base.hh); this example uses the matching Python factories dune.xt.grid.make_mixed_grid / make_prism_grid, which build the same meshes.

UGGrid is an optional dependency (the uggrid vcpkg feature). The cell below detects whether the installed wheel was built with it; if not, the demonstration cells further down are skipped with a note instead of failing, so this page still builds. (The property tests in python/gdt/test/test_hypothesis_mixed_grids.py exercise the same factories whenever the build provides them.)

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

import math

import numpy as np

from dune.xt.grid import Dim, visualize_grid

try:
    from dune.xt.grid import make_mixed_grid, make_prism_grid

    HAVE_UGGRID = True
except ImportError:
    HAVE_UGGRID = False
    print(
        "This wheel was built without the 'uggrid' vcpkg feature, so make_mixed_grid /\n"
        "make_prism_grid are not available and the cells below are skipped. Rebuild with the\n"
        "'uggrid' feature enabled (see CMakePresets.json) to run the full example."
    )

This wheel was built without the 'uggrid' vcpkg feature, so make_mixed_grid /
make_prism_grid are not available and the cells below are skipped. Rebuild with the
'uggrid' feature enabled (see CMakePresets.json) to run the full example.

A small helper classifies each leaf element by its number of corners (in 2d: 3 -> triangle/simplex, 4 -> quadrilateral/cube; in 3d a prism has 6) and counts how many of each geometry type the grid holds – this is what makes a mixed grid visibly different from the structured ones:

_GEOMETRY_BY_CORNERS = {
    (2, 3): "simplex",
    (2, 4): "cube",
    (3, 4): "simplex",
    (3, 6): "prism",
    (3, 8): "cube",
}


def element_geometry_counts(grid):
    dim = grid.dimension
    counts = {}

    def visit(element):
        geometry_type = _GEOMETRY_BY_CORNERS.get((dim, element.corners.shape[0]), "unknown")
        counts[geometry_type] = counts.get(geometry_type, 0) + 1

    grid.apply_on_each_element(visit)
    return counts

1: a mixed cube/simplex grid in 2d

make_mixed_grid(Dim(2)) inserts two quadrilaterals and two triangles that together tile a small polygonal domain:

if HAVE_UGGRID:
    grid = make_mixed_grid(Dim(2))
    print(f"{grid.size(0)} elements, {grid.size(grid.dimension)} vertices")
    print(element_geometry_counts(grid))

We can look at the mesh directly; the two triangles and two quadrilaterals share a common edge skeleton:

if HAVE_UGGRID:
    _ = visualize_grid(grid)

2: assembling and solving on the mixed grid

A discontinuous Galerkin (interior penalty) discretization is the natural choice here: it assembles element by element and couples neighbours across intersections, so it does not care whether two adjacent elements are the same geometry type. We refine the grid a couple of times for a finer solution – refinement multiplies the element counts but keeps both geometry types present:

if HAVE_UGGRID:
    grid.global_refine(2)
    print(element_geometry_counts(grid))

if HAVE_UGGRID:
    from discretize_elliptic_ipdg import discretize_elliptic_ipdg_dirichlet_zero

    u_h = discretize_elliptic_ipdg_dirichlet_zero(grid, diffusion=1, source=1)
    _ = visualize_grid(grid)

if HAVE_UGGRID:
    from dune.gdt import visualize_function

    _ = visualize_function(u_h)

The size of the DG space is bookkept per element and per geometry type: an order-k cube contributes (k + 1)^d (tensor-product Q_k) local degrees of freedom, a simplex contributes binomial(k + d, d) (P_k). Summing over the elements actually present must reproduce the space’s global DoF count exactly – the mixed-grid generalization of the per-element DG DoF property tested in python/gdt/test/test_hypothesis_mixed_grids.py:

def dg_dofs_per_element(geometry_type, order, dim):
    if geometry_type == "cube":
        return (order + 1) ** dim  # Q_k
    if geometry_type == "simplex":
        return math.comb(order + dim, dim)  # P_k
    raise NotImplementedError(geometry_type)


if HAVE_UGGRID:
    from dune.gdt import DiscontinuousLagrangeSpace

    order = 1
    space = DiscontinuousLagrangeSpace(grid, order=order)
    expected = sum(
        count * dg_dofs_per_element(geometry_type, order, grid.dimension)
        for geometry_type, count in element_geometry_counts(grid).items()
    )
    print(f"space.num_DoFs = {space.num_DoFs}, expected = {expected}")
    assert space.num_DoFs == expected

3: a prism grid in 3d

In three dimensions UGGrid additionally supports prisms (a triangle extruded along a line). make_prism_grid(Dim(3)) builds a single prism, which we refine once to obtain several:

if HAVE_UGGRID:
    prism_grid = make_prism_grid(Dim(3), num_refinements=1)
    print(element_geometry_counts(prism_grid))

All elements are prisms (6 corners each), a geometry type that is unreachable from Python with any of the other bound grid managers.

Download the code: example__mixed_and_prism_grids.md example__mixed_and_prism_grids.ipynb