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: structured 1d grids (YASP_1D vs OneDGrid)

Until now the Python bindings exposed exactly one one-dimensional grid, Dune::OneDGrid (make_cube_grid(Dim(1), ...)). WP2 of issue #320 adds the structured YaspGrid<1, EquidistantOffsetCoordinates> next to it – the 1d sibling of the YASP_2D / YASP_3D cube grids that already back the 2d/3d examples, and the grid the C++ SimplicialGrids typed-test list uses in 1d. This notebook

  • builds the same 1d domain as both an OneDGrid and a YASP_1D grid,

  • solves the stationary heat equation on each and checks the two discretizations agree, and

  • demonstrates the equidistant-offset feature (a structured grid on an arbitrary interval [a, b], a != 0) that motivates having a structured 1d grid at all.

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


import numpy as np

Two ways to mesh the unit interval

OneDGrid is selected through the dimension-only make_cube_grid overload, while the structured YASP_1D grid is selected through the same (Dimension, Cube) overload as its 2d/3d counterparts – so both one-dimensional grids can coexist without an ambiguous factory signature.

from dune.xt.grid import Cube, Dim, make_cube_grid

n = 8  # number of intervals
lower_left, upper_right = [0.0], [1.0]

oned_grid = make_cube_grid(Dim(1), lower_left=lower_left, upper_right=upper_right,
                           num_elements=[n])
yasp_grid = make_cube_grid(Dim(1), Cube(), lower_left=lower_left, upper_right=upper_right,
                           num_elements=[n])

for name, grid in (("OneDGrid", oned_grid), ("YASP_1D", yasp_grid)):
    print(f'{name:9s}: {grid.size(0)} elements, {grid.size(1)} vertices')

OneDGrid : 8 elements, 9 vertices
YASP_1D  : 8 elements, 9 vertices

Both provide the identical structured skeleton (num_elements equal intervals of width h = 1/num_elements); the element centers coincide up to ordering:

oned_centers = np.sort(np.array(oned_grid.centers(0), copy=False).ravel())
yasp_centers = np.sort(np.array(yasp_grid.centers(0), copy=False).ravel())

assert np.allclose(oned_centers, yasp_centers)
print('element centers:', yasp_centers)

element centers: [0.0625 0.1875 0.3125 0.4375 0.5625 0.6875 0.8125 0.9375]

The stationary heat equation on both grids

We reuse the continuous-Lagrange discretization from discretize_elliptic_cg.py (the same helper the 2d CG tutorial uses) to solve

\[-u''(x) = 1 \quad\text{on } (0, 1), \qquad u(0) = u(1) = 0,\]

whose exact solution is \(u(x) = \tfrac{1}{2}\,x\,(1 - x)\).

from discretize_elliptic_cg import discretize_elliptic_cg_dirichlet_zero

def solve(grid):
    # constant unit diffusion and unit source, homogeneous Dirichlet everywhere
    u_h = discretize_elliptic_cg_dirichlet_zero(grid, diffusion=1.0, source=1.0)
    return np.array(u_h.dofs.vector, copy=False)

oned_dofs = solve(oned_grid)
yasp_dofs = solve(yasp_grid)

The two grids give rise to the identical P1 linear system (up to a DoF permutation), so the sorted DoF vectors must match to solver accuracy:

assert oned_dofs.shape == yasp_dofs.shape == (n + 1,)
grid_to_grid = np.max(np.abs(np.sort(oned_dofs) - np.sort(yasp_dofs)))
print(f'max |OneDGrid - YASP_1D| over the DoFs: {grid_to_grid:.2e}')
assert grid_to_grid < 1e-6

max |OneDGrid - YASP_1D| over the DoFs: 0.00e+00

In one dimension the P1 Galerkin solution of the Poisson problem is nodally exact, so the sorted DoF values reproduce the exact solution sampled on the nodes \(x_i = i/N\):

nodes = np.linspace(0.0, 1.0, n + 1)
exact_nodal = np.sort(0.5 * nodes * (1.0 - nodes))
nodal_error = np.max(np.abs(np.sort(yasp_dofs) - exact_nodal))
print(f'max nodal error vs. exact solution: {nodal_error:.2e}')
assert nodal_error < 1e-6

max nodal error vs. exact solution: 4.73e-15

The equidistant-offset feature: a structured grid on [a, b]

The whole point of EquidistantOffsetCoordinates is that the structured grid need not start at the origin. We put a YASP_1D grid on the offset interval [2, 5] and solve the same problem there,

\[-u''(x) = 1 \quad\text{on } (2, 5), \qquad u(2) = u(5) = 0,\]

with exact solution \(u(x) = \tfrac{1}{2}\,(x - 2)\,(5 - x)\).

a, b = 2.0, 5.0
offset_grid = make_cube_grid(Dim(1), Cube(), lower_left=[a], upper_right=[b],
                             num_elements=[n])

centers = np.sort(np.array(offset_grid.centers(0), copy=False).ravel())
h = (b - a) / n
print(f'offset grid spans [{centers.min() - 0.5 * h:.1f}, {centers.max() + 0.5 * h:.1f}]')

offset_dofs = solve(offset_grid)
offset_nodes = np.linspace(a, b, n + 1)
offset_exact = np.sort(0.5 * (offset_nodes - a) * (b - offset_nodes))
offset_error = np.max(np.abs(np.sort(offset_dofs) - offset_exact))
print(f'max nodal error on [{a}, {b}]: {offset_error:.2e}')
assert offset_error < 1e-6

offset grid spans [2.0, 5.0]

max nodal error on [2.0, 5.0]: 4.26e-14

The structured 1d grid solves the shifted problem exactly, without any special handling of the non-zero offset – the same behaviour the 2d/3d YASP cube grids already provide, now available in one dimension for property testing and structured-grid experiments alongside OneDGrid.

Download the code: example__structured_1d_grids.md example__structured_1d_grids.ipynb