This tutorial is part of multiphenicsx.

Tutorial 07: understanding restrictions

import os
import typing

import basix
import basix.ufl
import dolfinx.fem
import dolfinx.io
import dolfinx.mesh
import gmsh
import matplotlib as mpl
import matplotlib.collections
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import mpi4py.MPI
import mpl_toolkits.axes_grid1
import numpy as np
import numpy.typing

import multiphenicsx.fem

Auxiliary functions for plotting

def plot_mesh(mesh: dolfinx.mesh.Mesh, ax: typing.Optional[plt.Axes] = None) -> plt.Axes:
    """Plot a mesh object on the provied (or, if None, the current) axes object."""
    if ax is None:
        ax = plt.gca()
    ax.set_aspect("equal")
    points = mesh.geometry.x
    cells = mesh.geometry.dofmap
    tria = tri.Triangulation(points[:, 0], points[:, 1], cells)
    ax.triplot(tria, color="k")
    return ax

def plot_mesh_tags(
    mesh: dolfinx.mesh.Mesh, mesh_tags: dolfinx.mesh.MeshTags, ax: typing.Optional[plt.Axes] = None
) -> plt.Axes:
    """Plot a mesh tags object on the provied (or, if None, the current) axes object."""
    if ax is None:
        ax = plt.gca()
    ax.set_aspect("equal")
    points = mesh.geometry.x
    colors = ["b", "r"]
    cmap = mpl.colors.ListedColormap(colors)
    cmap_bounds = [0, 0.5, 1]
    norm = mpl.colors.BoundaryNorm(cmap_bounds, cmap.N)
    assert mesh_tags.dim in (mesh.topology.dim, mesh.topology.dim - 1)
    if mesh_tags.dim == mesh.topology.dim:
        cells = mesh.geometry.dofmap
        tria = tri.Triangulation(points[:, 0], points[:, 1], cells)
        cell_colors = np.zeros((cells.shape[0], ))
        cell_colors[mesh_tags.indices[mesh_tags.values == 1]] = 1
        mappable: mpl.collections.Collection = ax.tripcolor(
            tria, cell_colors, edgecolor="k", cmap=cmap, norm=norm)
    elif mesh_tags.dim == mesh.topology.dim - 1:
        tdim = mesh.topology.dim
        cells_map = mesh.topology.index_map(mesh.topology.dim)
        num_cells = cells_map.size_local + cells_map.num_ghosts
        connectivity_cells_to_facets = mesh.topology.connectivity(tdim, tdim - 1)
        connectivity_cells_to_vertices = mesh.topology.connectivity(tdim, 0)
        connectivity_facets_to_vertices = mesh.topology.connectivity(tdim - 1, 0)
        vertex_map = {
            topology_index: geometry_index for c in range(num_cells) for (topology_index, geometry_index) in zip(
                connectivity_cells_to_vertices.links(c), mesh.geometry.dofmap[c])
        }
        linestyles = [(0, (5, 10)), "solid"]
        lines = list()
        lines_colors_as_int = list()
        lines_colors_as_str = list()
        lines_linestyles = list()
        mesh_tags_1 = mesh_tags.indices[mesh_tags.values == 1]
        for c in range(num_cells):
            facets = connectivity_cells_to_facets.links(c)
            for f in facets:
                if f in mesh_tags_1:
                    value_f = 1
                else:
                    value_f = 0
                vertices = [vertex_map[v] for v in connectivity_facets_to_vertices.links(f)]
                lines.append(points[vertices][:, :2])
                lines_colors_as_int.append(value_f)
                lines_colors_as_str.append(colors[value_f])
                lines_linestyles.append(linestyles[value_f])
        mappable: mpl.collections.Collection = mpl.collections.LineCollection(  # type: ignore[no-redef]
            lines, cmap=cmap, norm=norm, colors=lines_colors_as_str, linestyles=lines_linestyles)
        mappable.set_array(np.array(lines_colors_as_int))
        ax.add_collection(mappable)
        ax.autoscale()
    divider = mpl_toolkits.axes_grid1.make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(mappable, cax=cax, boundaries=cmap_bounds, ticks=cmap_bounds)
    return ax

def _plot_dofmap(coordinates: np.typing.NDArray[np.float64], ax: typing.Optional[plt.Axes] = None) -> plt.Axes:
    if ax is None:
        ax = plt.gca()
    text_offset = [1e-2, 1e-2]
    ax.scatter(coordinates[:, 0], coordinates[:, 1], c="k", s=50)
    for c in np.unique(coordinates, axis=0):
        dofs_c = (coordinates == c).all(axis=1).nonzero()[0]
        text_c = np.array2string(dofs_c, separator=", ", max_line_width=10)
        ax.text(c[0] + text_offset[0], c[1] + text_offset[1], text_c, fontsize=20)
    return ax

def plot_dofmap(V: dolfinx.fem.FunctionSpace, ax: typing.Optional[plt.Axes] = None) -> plt.Axes:
    """Plot the DOFs in a function space object on the provied (or, if None, the current) axes object."""
    coordinates = V.tabulate_dof_coordinates().round(decimals=3)
    return _plot_dofmap(coordinates, ax)

def plot_dofmap_restriction(
    V: dolfinx.fem.FunctionSpace, restriction: multiphenicsx.fem.DofMapRestriction,
    ax: typing.Optional[plt.Axes] = None
) -> plt.Axes:
    """Plot the DOFs in a DofMapRestriction object on the provied (or, if None, the current) axes object."""
    coordinates = V.tabulate_dof_coordinates().round(decimals=3)
    return _plot_dofmap(coordinates[list(restriction.unrestricted_to_restricted.keys())], ax)

Auxiliary functions for asserts

def count_dofs(restriction: multiphenicsx.fem.DofMapRestriction, comm: mpi4py.MPI.Intracomm) -> int:
    """Count the DOFs in a DofMapRestriction object."""
    u2r = restriction.unrestricted_to_restricted
    restricted_local_indices = np.array([r for (_, r) in u2r.items()], dtype=np.int32)
    dofs_V_restriction_global = restriction.index_map.local_to_global(restricted_local_indices)
    return len(set(gdof for gdofs in comm.allgather(dofs_V_restriction_global) for gdof in gdofs))

def locate_dofs_by_polar_coordinates(
    r: typing.Union[int, float, np.float64], theta: typing.Union[int, float, np.float64],
    V: dolfinx.fem.FunctionSpace, restriction: multiphenicsx.fem.DofMapRestriction
) -> set[np.int32]:
    """Determine which DOFs in a DofMapRestriction object are located at a point (written in polar coordinates)."""
    p = [r * np.cos(theta), r * np.sin(theta), 0.]
    dofs_V = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x.T, p).all(axis=1)).reshape(-1)
    u2r = restriction.unrestricted_to_restricted
    restricted_local_indices = np.array([u2r[dof] for dof in dofs_V if dof in u2r], dtype=np.int32)
    dofs_V_restriction_global = restriction.index_map.local_to_global(restricted_local_indices)
    return set(gdof for gdofs in V.mesh.comm.allgather(dofs_V_restriction_global) for gdof in gdofs)

Geometrical parameters

r = 1
mesh_size = 2

Mesh

Generate a simple mesh consisting in an hexagon discretized with six equilateral triangle cells.

gmsh.initialize()
gmsh.model.add("mesh")

points = [
    gmsh.model.geo.addPoint(np.cos(t / 3 * np.pi), np.sin(t / 3 * np.pi), 0.0, mesh_size) for t in range(6)]
lines = [gmsh.model.geo.addLine(points[t], points[(t + 1) % 6]) for t in range(6)]
line_loop = gmsh.model.geo.addCurveLoop(lines)
domain = gmsh.model.geo.addPlaneSurface([line_loop])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, lines, 0)
gmsh.model.addPhysicalGroup(2, [domain], 0)
gmsh.model.mesh.generate(2)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 20%] Meshing curve 2 (Line)
Info    : [ 40%] Meshing curve 3 (Line)
Info    : [ 60%] Meshing curve 4 (Line)
Info    : [ 70%] Meshing curve 5 (Line)
Info    : [ 90%] Meshing curve 6 (Line)
Info    : Done meshing 1D (Wall 0.000424279s, CPU 0.000821s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.000281093s, CPU 0.000282s)
Info    : 7 nodes 18 elements

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

if "PYTEST_CURRENT_TEST" not in os.environ:
    plot_mesh(mesh)

Create connectivities required by the rest of the code.

mesh.topology.create_connectivity(mesh.topology.dim - 1, 0)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, 0)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)

Mesh restrictions on cells

Define mesh tags on cells, which are equal to one on all cells.

cell_entities_all = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim, lambda x: np.full((x.shape[1], ), True))
cell_values_all = np.full(cell_entities_all.shape, 1, dtype=np.int32)
cell_restriction_all = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim, cell_entities_all, cell_values_all)
cell_restriction_all.name = "cell_restriction_all"

if "PYTEST_CURRENT_TEST" not in os.environ:
    plot_mesh_tags(mesh, cell_restriction_all)

Define mesh tags on cells, which are equal to one on one half of the cells

eps = np.finfo(float).eps
cell_entities_subset = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim,
    lambda x: np.logical_or(x[0] < eps, np.logical_and(x[1] < eps, x[0] < 0.5 + eps)))
cell_values_subset = np.full(cell_entities_subset.shape, 1, dtype=np.int32)
cell_restriction_subset = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim, cell_entities_subset, cell_values_subset)
cell_restriction_subset.name = "cell_restriction_subset"

if "PYTEST_CURRENT_TEST" not in os.environ:
    plot_mesh_tags(mesh, cell_restriction_subset)

Mesh restrictions on facets

Define mesh tags on facets, which are equal to one on all facets

facet_entities_all = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim - 1, lambda x: np.full((x.shape[1], ), True))
facet_values_all = np.full(facet_entities_all.shape, 1, dtype=np.int32)
facet_restriction_all = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim - 1, facet_entities_all, facet_values_all)
facet_restriction_all.name = "facet_restriction_all"

if "PYTEST_CURRENT_TEST" not in os.environ:
    plot_mesh_tags(mesh, facet_restriction_all)

Define mesh tags on facets, which are equal to one on two facets

facet_entities_subset = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim - 1, lambda x: np.fabs(x[1] + np.sqrt(3) * x[0]) < 0.01)
facet_values_subset = np.full(facet_entities_subset.shape, 1, dtype=np.int32)
facet_restriction_subset = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim - 1, facet_entities_subset, facet_values_subset)
facet_restriction_subset.name = "facet_restriction_subset"

if "PYTEST_CURRENT_TEST" not in os.environ:
    plot_mesh_tags(mesh, facet_restriction_subset)

Mesh restrictions summary

cell_restrictions = (cell_restriction_all, cell_restriction_subset)
facet_restrictions = (facet_restriction_all, facet_restriction_subset)
all_restrictions = cell_restrictions + facet_restrictions

Lagrange spaces

Define Lagrange FE spaces of order $k=1, 2, 3$, and plot the associated DofMap.

CG_elem = [
    basix.ufl.element(
        "Lagrange", mesh.basix_cell(), k, lagrange_variant=basix.LagrangeVariant.equispaced
    ) for k in (1, 2, 3)
]
CG = [dolfinx.fem.functionspace(mesh, CG_elem_k) for CG_elem_k in CG_elem]

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 1, figsize=(10, 30))
    assert isinstance(ax, np.ndarray)
    for (k, CGk) in enumerate(CG):
        plot_mesh(mesh, ax[k])
        plot_dofmap(CGk, ax[k])
        ax[k].set_title("CG " + str(k + 1) + " DofMap", fontsize=30)

Define DofMapRestriction objects associated to the Lagrange FE spaces, for all four restrictions

dofmap_restriction_CG: dict[
    dolfinx.mesh.MeshTags, list[multiphenicsx.fem.DofMapRestriction]] = dict()
for restriction in all_restrictions:
    dofmap_restriction_CG[restriction] = list()
    for CGk in CG:
        restrict_CGk = dolfinx.fem.locate_dofs_topological(
            CGk, restriction.dim, restriction.indices[restriction.values == 1])
        dofmap_restriction_CG[restriction].append(
            multiphenicsx.fem.DofMapRestriction(CGk.dofmap, restrict_CGk))

Compare DOFs for the case of cell restriction equal to one on the entire domain. There is indeed no difference.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, CGk) in enumerate(CG):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(CGk, ax[k, 0])
        ax[k, 0].set_title("CG " + str(k + 1) + " DofMap", fontsize=30)
    for (k, (CGk, dofmap_restriction_CGk)) in enumerate(zip(CG, dofmap_restriction_CG[cell_restriction_all])):
        plot_mesh_tags(mesh, cell_restriction_all, ax[k, 1])
        plot_dofmap_restriction(CGk, dofmap_restriction_CGk, ax[k, 1])
        ax[k, 1].set_title("CG " + str(k + 1) + " DofMapRestriction", fontsize=30)

Assert that DOFs are at the expected locations

CG_1 = CG[0]
dofmap_restriction_CG_1 = dofmap_restriction_CG[cell_restriction_all][0]
assert count_dofs(dofmap_restriction_CG_1, mesh.comm) == 7
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1

CG_2 = CG[1]
dofmap_restriction_CG_2 = dofmap_restriction_CG[cell_restriction_all][1]
assert count_dofs(dofmap_restriction_CG_2, mesh.comm) == 19
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 3 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 11 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1

CG_3 = CG[2]
dofmap_restriction_CG_3 = dofmap_restriction_CG[cell_restriction_all][2]
assert count_dofs(dofmap_restriction_CG_3, mesh.comm) == 37
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 3 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 5 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 7 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 9 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 11 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 2 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 2 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 4 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 4 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 5 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 5 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1

Compare DOFs for che case of cell restriction equal to one on a subset of the domain. Note how the DofMapRestriction has only a subset of the DOFs of the DofMap, and properly renumbers them.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, CGk) in enumerate(CG):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(CGk, ax[k, 0])
        ax[k, 0].set_title("CG " + str(k + 1) + " DofMap", fontsize=30)
    for (k, (CGk, dofmap_restriction_CGk)) in enumerate(zip(CG, dofmap_restriction_CG[cell_restriction_subset])):
        plot_mesh_tags(mesh, cell_restriction_subset, ax[k, 1])
        plot_dofmap_restriction(CGk, dofmap_restriction_CGk, ax[k, 1])
        ax[k, 1].set_title("CG " + str(k + 1) + " DofMapRestriction", fontsize=30)

Assert that DOFs are at the expected locations

CG_1 = CG[0]
dofmap_restriction_CG_1 = dofmap_restriction_CG[cell_restriction_subset][0]
assert count_dofs(dofmap_restriction_CG_1, mesh.comm) == 5
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1

CG_2 = CG[1]
dofmap_restriction_CG_2 = dofmap_restriction_CG[cell_restriction_subset][1]
assert count_dofs(dofmap_restriction_CG_2, mesh.comm) == 12
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1

CG_3 = CG[2]
dofmap_restriction_CG_3 = dofmap_restriction_CG[cell_restriction_subset][2]
assert count_dofs(dofmap_restriction_CG_3, mesh.comm) == 22
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 5 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 7 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 9 * np.pi / 6, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 2 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 4 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 4 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 5 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1

Compare DOFs for che case of facet restriction equal to one on the entire domain. Note how there is no difference for $k=1, 2$, but the cases $k=3$ differ (the DofMapRestriction does not have the DOF at the cell center).

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, CGk) in enumerate(CG):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(CG[k], ax[k, 0])
        ax[k, 0].set_title("CG " + str(k + 1) + " DofMap", fontsize=30)
    for (k, (CGk, dofmap_restriction_CGk)) in enumerate(zip(CG, dofmap_restriction_CG[facet_restriction_all])):
        plot_mesh_tags(mesh, facet_restriction_all, ax[k, 1])
        plot_dofmap_restriction(CGk, dofmap_restriction_CGk, ax[k, 1])
        ax[k, 1].set_title("CG " + str(k + 1) + " DofMapRestriction", fontsize=30)

Assert that DOFs are at the expected locations

CG_1 = CG[0]
dofmap_restriction_CG_1 = dofmap_restriction_CG[facet_restriction_all][0]
assert count_dofs(dofmap_restriction_CG_1, mesh.comm) == 7
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1

CG_2 = CG[1]
dofmap_restriction_CG_2 = dofmap_restriction_CG[facet_restriction_all][1]
assert count_dofs(dofmap_restriction_CG_2, mesh.comm) == 19
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 3 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 11 * np.pi / 6, CG_2, dofmap_restriction_CG_2)) == 1

CG_3 = CG[2]
dofmap_restriction_CG_3 = dofmap_restriction_CG[facet_restriction_all][2]
assert count_dofs(dofmap_restriction_CG_3, mesh.comm) == 31
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 2 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 2 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 4 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 4 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 5 * np.pi / 3 - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, 5 * np.pi / 3 + np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(7) / 3, np.pi - np.arctan(np.sqrt(3) / 5), CG_3, dofmap_restriction_CG_3)) == 1

Compare DOFs for che case of facet restriction equal to one on a subset of the domain. Note how the DofMapRestriction has only a subset of the DOFs of the DofMap, and properly renumbers them.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, CGk) in enumerate(CG):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(CG[k], ax[k, 0])
        ax[k, 0].set_title("CG " + str(k + 1) + " DofMap", fontsize=30)
    for (k, (CGk, dofmap_restriction_CGk)) in enumerate(zip(CG, dofmap_restriction_CG[facet_restriction_subset])):
        plot_mesh_tags(mesh, facet_restriction_subset, ax[k, 1])
        plot_dofmap_restriction(CGk, dofmap_restriction_CGk, ax[k, 1])
        ax[k, 1].set_title("CG " + str(k + 1) + " DofMapRestriction", fontsize=30)

Assert that DOFs are at the expected locations

CG_1 = CG[0]
dofmap_restriction_CG_1 = dofmap_restriction_CG[facet_restriction_subset][0]
assert count_dofs(dofmap_restriction_CG_1, mesh.comm) == 3
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1

CG_2 = CG[1]
dofmap_restriction_CG_2 = dofmap_restriction_CG[facet_restriction_subset][1]
assert count_dofs(dofmap_restriction_CG_2, mesh.comm) == 5
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1

CG_3 = CG[2]
dofmap_restriction_CG_3 = dofmap_restriction_CG[facet_restriction_subset][2]
assert count_dofs(dofmap_restriction_CG_3, mesh.comm) == 7
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    2 / 3, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_3, dofmap_restriction_CG_3)) == 1

Discontinuous Galerkin spaces

Define Discontinuous Galerkin spaces of order $k = 0, 1, 2$, and plot the associated DofMap. This spaces will be used in combination with a cell restriction, as DG DOFs are associated to cells.

DG = [dolfinx.fem.functionspace(mesh, ("Discontinuous Lagrange", k)) for k in (0, 1, 2)]

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 1, figsize=(10, 30))
    assert isinstance(ax, np.ndarray)
    for (k, DGk) in enumerate(DG):
        plot_mesh(mesh, ax[k])
        plot_dofmap(DGk, ax[k])
        ax[k].set_title("DG " + str(k) + " DofMap", fontsize=30)

Define Discontinuous Galerkin Trace spaces of order $k = 0, 1, 2, 3$, and plot the associated DofMap. This spaces will be used in combination with a facet restriction, as DGT DOFs are associated to facets.

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT = [dolfinx.fem.functionspace(mesh, ("Discontinuous Lagrange Trace", k)) for k in (0, 1, 2)]
"""

'Expect this cell to fail.\n\nDGT = [dolfinx.fem.functionspace(mesh, ("Discontinuous Lagrange Trace", k)) for k in (0, 1, 2)]\n'

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 1, figsize=(10, 30))
    assert isinstance(ax, np.ndarray)
    for (k, DGTk) in enumerate(DGT):
        plot_mesh(mesh, ax[k])
        plot_dofmap(DGTk, ax[k])
        ax[k].set_title("DGT " + str(k) + " DofMap\n", fontsize=30)
"""

'Expect this cell to fail.\n\nif "PYTEST_CURRENT_TEST" not in os.environ:\n    _, ax = plt.subplots(3, 1, figsize=(10, 30))\n    assert isinstance(ax, np.ndarray)\n    for (k, DGTk) in enumerate(DGT):\n        plot_mesh(mesh, ax[k])\n        plot_dofmap(DGTk, ax[k])\n        ax[k].set_title("DGT " + str(k) + " DofMap\n", fontsize=30)\n'

Define DofMapRestriction objects associated to the Discontinuos Galerkin FE spaces, for all cell restrictions

dofmap_restriction_DG: dict[
    dolfinx.mesh.MeshTags, list[multiphenicsx.fem.DofMapRestriction]] = dict()
for restriction in cell_restrictions:
    dofmap_restriction_DG[restriction] = list()
    for DGk in DG:
        restrict_DGk = dolfinx.fem.locate_dofs_topological(
            DGk, restriction.dim, restriction.indices[restriction.values == 1])
        dofmap_restriction_DG[restriction].append(
            multiphenicsx.fem.DofMapRestriction(DGk.dofmap, restrict_DGk))

Define DofMapRestriction objects associated to the Discontinuos Galerkin Trace FE spaces, for all facet restrictions

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

dofmap_restriction_DGT: dict[
    dolfinx.mesh.MeshTags, list[multiphenicsx.fem.DofMapRestriction]] = dict()
for restriction in facet_restrictions:
    dofmap_restriction_DGT[restriction] = list()
    for DGTk in DGT:
        restrict_DGTk = dolfinx.fem.locate_dofs_topological(
            DGTk, restriction.dim, restriction.indices[restriction.values == 1])
        dofmap_restriction_DGT[restriction].append(
            multiphenicsx.fem.DofMapRestriction(DGTk.dofmap, restrict_DGTk))
"""

'Expect this cell to fail.\n\ndofmap_restriction_DGT: dict[\n    dolfinx.mesh.MeshTags, list[multiphenicsx.fem.DofMapRestriction]] = dict()\nfor restriction in facet_restrictions:\n    dofmap_restriction_DGT[restriction] = list()\n    for DGTk in DGT:\n        restrict_DGTk = dolfinx.fem.locate_dofs_topological(\n            DGTk, restriction.dim, restriction.indices[restriction.values == 1])\n        dofmap_restriction_DGT[restriction].append(\n            multiphenicsx.fem.DofMapRestriction(DGTk.dofmap, restrict_DGTk))\n'

Compare DOFs for the case of cell restriction equal to one on the entire domain. There is indeed no difference.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, DGk) in enumerate(DG):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(DGk, ax[k, 0])
        ax[k, 0].set_title("DG " + str(k) + " DofMap", fontsize=30)
    for (k, (DGk, dofmap_restriction_DGk)) in enumerate(zip(DG, dofmap_restriction_DG[cell_restriction_all])):
        plot_mesh_tags(mesh, cell_restriction_all, ax[k, 1])
        plot_dofmap_restriction(DGk, dofmap_restriction_DGk, ax[k, 1])
        ax[k, 1].set_title("DG " + str(k) + " DofMapRestriction", fontsize=30)

Assert that DOFs are at the expected locations

DG_0 = DG[0]
dofmap_restriction_DG_0 = dofmap_restriction_DG[cell_restriction_all][0]
assert count_dofs(dofmap_restriction_DG_0, mesh.comm) == 6
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 3 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 5 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 7 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 9 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 11 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1

DG_1 = DG[1]
dofmap_restriction_DG_1 = dofmap_restriction_DG[cell_restriction_all][1]
assert count_dofs(dofmap_restriction_DG_1, mesh.comm) == 18
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DG_1, dofmap_restriction_DG_1)) == 6
assert len(locate_dofs_by_polar_coordinates(
    1, 0, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 2

DG_2 = DG[2]
dofmap_restriction_DG_2 = dofmap_restriction_DG[cell_restriction_all][2]
assert count_dofs(dofmap_restriction_DG_2, mesh.comm) == 36
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DG_2, dofmap_restriction_DG_2)) == 6
assert len(locate_dofs_by_polar_coordinates(
    0.5, 0, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 0, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 3 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 11 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1

Compare DOFs for che case of cell restriction equal to one on a subset of the domain. Note how the DofMapRestriction has only a subset of the DOFs of the DofMap, and properly renumbers them. Note also that the number of DOFs at the same physical location might be different between DofMap and DofMapRestriction (see e.g. the center of the hexagon).

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, DGk) in enumerate(DG):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(DGk, ax[k, 0])
        ax[k, 0].set_title("DG " + str(k) + " DofMap", fontsize=30)
    for (k, (DGk, dofmap_restriction_DGk)) in enumerate(zip(DG, dofmap_restriction_DG[cell_restriction_subset])):
        plot_mesh_tags(mesh, cell_restriction_subset, ax[k, 1])
        plot_dofmap_restriction(DGk, dofmap_restriction_DGk, ax[k, 1])
        ax[k, 1].set_title("DG " + str(k) + " DofMapRestriction", fontsize=30)

Assert that DOFs are at the expected locations

DG_0 = DG[0]
dofmap_restriction_DG_0 = dofmap_restriction_DG[cell_restriction_subset][0]
assert count_dofs(dofmap_restriction_DG_0, mesh.comm) == 3
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 5 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 7 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 3, 9 * np.pi / 6, DG_0, dofmap_restriction_DG_0)) == 1

DG_1 = DG[1]
dofmap_restriction_DG_1 = dofmap_restriction_DG[cell_restriction_subset][1]
assert count_dofs(dofmap_restriction_DG_1, mesh.comm) == 9
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DG_1, dofmap_restriction_DG_1)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DG_1, dofmap_restriction_DG_1)) == 1

DG_2 = DG[2]
dofmap_restriction_DG_2 = dofmap_restriction_DG[cell_restriction_subset][2]
assert count_dofs(dofmap_restriction_DG_2, mesh.comm) == 18
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DG_2, dofmap_restriction_DG_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, DG_2, dofmap_restriction_DG_2)) == 1

Compare DOFs for che case of facet restriction equal to one on the entire domain. There is indeed no difference.

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, DGTk) in enumerate(DGT):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(DGTk, ax[k, 0])
        ax[k, 0].set_title("DGT " + str(k) + " DofMap\n", fontsize=30)
    for (k, (DGTk, dofmap_restriction_DGTk)) in enumerate(zip(DGT, dofmap_restriction_DGT[facet_restriction_all])):
        plot_mesh_tags(mesh, facet_restriction_all, ax[k, 1])
        plot_dofmap_restriction(DGTk, dofmap_restriction_DGTk, ax[k, 1])
        ax[k, 1].set_title("DGT " + str(k) + " DofMapRestriction\n", fontsize=30)
"""

'Expect this cell to fail.\n\nif "PYTEST_CURRENT_TEST" not in os.environ:\n    _, ax = plt.subplots(3, 2, figsize=(20, 30))\n    assert isinstance(ax, np.ndarray)\n    for (k, DGTk) in enumerate(DGT):\n        plot_mesh(mesh, ax[k, 0])\n        plot_dofmap(DGTk, ax[k, 0])\n        ax[k, 0].set_title("DGT " + str(k) + " DofMap\n", fontsize=30)\n    for (k, (DGTk, dofmap_restriction_DGTk)) in enumerate(zip(DGT, dofmap_restriction_DGT[facet_restriction_all])):\n        plot_mesh_tags(mesh, facet_restriction_all, ax[k, 1])\n        plot_dofmap_restriction(DGTk, dofmap_restriction_DGTk, ax[k, 1])\n        ax[k, 1].set_title("DGT " + str(k) + " DofMapRestriction\n", fontsize=30)\n'

Assert that DOFs are at the expected locations

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT_0 = DGT[0]
dofmap_restriction_DGT_0 = dofmap_restriction_DGT[facet_restriction_all][0]
assert count_dofs(dofmap_restriction_DGT_0, mesh.comm) == 12
assert len(locate_dofs_by_polar_coordinates(
    0.5, 0, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 3 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 11 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1
"""

'Expect this cell to fail.\n\nDGT_0 = DGT[0]\ndofmap_restriction_DGT_0 = dofmap_restriction_DGT[facet_restriction_all][0]\nassert count_dofs(dofmap_restriction_DGT_0, mesh.comm) == 12\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 0, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 2 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, np.pi, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 4 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 5 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 3 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 5 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 7 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 9 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 11 * np.pi / 6, DGT_0, dofmap_restriction_DGT_0)) == 1\n'

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT_1 = DGT[1]
dofmap_restriction_DGT_1 = dofmap_restriction_DGT[facet_restriction_all][1]
assert count_dofs(dofmap_restriction_DGT_1, mesh.comm) == 24
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DGT_1, dofmap_restriction_DGT_1)) == 6
assert len(locate_dofs_by_polar_coordinates(
    1, 0, DGT_1, dofmap_restriction_DGT_1)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, DGT_1, dofmap_restriction_DGT_1)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3
"""

'Expect this cell to fail.\n\nDGT_1 = DGT[1]\ndofmap_restriction_DGT_1 = dofmap_restriction_DGT[facet_restriction_all][1]\nassert count_dofs(dofmap_restriction_DGT_1, mesh.comm) == 24\nassert len(locate_dofs_by_polar_coordinates(\n    0, 0, DGT_1, dofmap_restriction_DGT_1)) == 6\nassert len(locate_dofs_by_polar_coordinates(\n    1, 0, DGT_1, dofmap_restriction_DGT_1)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, 2 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, np.pi, DGT_1, dofmap_restriction_DGT_1)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, 4 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, 5 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 3\n'

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT_2 = DGT[2]
dofmap_restriction_DGT_2 = dofmap_restriction_DGT[facet_restriction_all][2]
assert count_dofs(dofmap_restriction_DGT_2, mesh.comm) == 36
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DGT_2, dofmap_restriction_DGT_2)) == 6
assert len(locate_dofs_by_polar_coordinates(
    0.5, 0, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, np.pi, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 4 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 0, DGT_2, dofmap_restriction_DGT_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, np.pi, DGT_2, dofmap_restriction_DGT_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 4 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 3 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 5 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 7 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 9 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    np.sqrt(3) / 2, 11 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1
"""

'Expect this cell to fail.\n\nDGT_2 = DGT[2]\ndofmap_restriction_DGT_2 = dofmap_restriction_DGT[facet_restriction_all][2]\nassert count_dofs(dofmap_restriction_DGT_2, mesh.comm) == 36\nassert len(locate_dofs_by_polar_coordinates(\n    0, 0, DGT_2, dofmap_restriction_DGT_2)) == 6\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 0, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, np.pi, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 4 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    1, 0, DGT_2, dofmap_restriction_DGT_2)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, np.pi, DGT_2, dofmap_restriction_DGT_2)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, 4 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    1, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 3\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 3 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 5 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 7 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 9 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    np.sqrt(3) / 2, 11 * np.pi / 6, DGT_2, dofmap_restriction_DGT_2)) == 1\n'

Compare DOFs for che case of facet restriction equal to one on a subset of the domain. Note how the DofMapRestriction has only a subset of the DOFs of the DofMap, and properly renumbers them. Note also that the number of DOFs at the same physical location might be different between DofMap and DofMapRestriction (see e.g. the center of the hexagon).

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

if "PYTEST_CURRENT_TEST" not in os.environ:
    _, ax = plt.subplots(3, 2, figsize=(20, 30))
    assert isinstance(ax, np.ndarray)
    for (k, DGTk) in enumerate(DGT):
        plot_mesh(mesh, ax[k, 0])
        plot_dofmap(DGTk, ax[k, 0])
        ax[k, 0].set_title("DGT " + str(k) + " DofMap\n", fontsize=30)
    for (k, (DGTk, dofmap_restriction_DGTk)) in enumerate(zip(DGT, dofmap_restriction_DGT[facet_restriction_subset])):
        plot_mesh_tags(mesh, facet_restriction_subset, ax[k, 1])
        plot_dofmap_restriction(DGTk, dofmap_restriction_DGTk, ax[k, 1])
        ax[k, 1].set_title("DGT " + str(k) + " DofMapRestriction\n", fontsize=30)
"""

'Expect this cell to fail.\n\nif "PYTEST_CURRENT_TEST" not in os.environ:\n    _, ax = plt.subplots(3, 2, figsize=(20, 30))\n    assert isinstance(ax, np.ndarray)\n    for (k, DGTk) in enumerate(DGT):\n        plot_mesh(mesh, ax[k, 0])\n        plot_dofmap(DGTk, ax[k, 0])\n        ax[k, 0].set_title("DGT " + str(k) + " DofMap\n", fontsize=30)\n    for (k, (DGTk, dofmap_restriction_DGTk)) in enumerate(zip(DGT, dofmap_restriction_DGT[facet_restriction_subset])):\n        plot_mesh_tags(mesh, facet_restriction_subset, ax[k, 1])\n        plot_dofmap_restriction(DGTk, dofmap_restriction_DGTk, ax[k, 1])\n        ax[k, 1].set_title("DGT " + str(k) + " DofMapRestriction\n", fontsize=30)\n'

Assert that DOFs are at the expected locations

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT_0 = DGT[0]
dofmap_restriction_DGT_0 = dofmap_restriction_DGT[facet_restriction_subset][0]
assert count_dofs(dofmap_restriction_DGT_0, mesh.comm) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1
"""

'Expect this cell to fail.\n\nDGT_0 = DGT[0]\ndofmap_restriction_DGT_0 = dofmap_restriction_DGT[facet_restriction_subset][0]\nassert count_dofs(dofmap_restriction_DGT_0, mesh.comm) == 2\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 2 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 5 * np.pi / 3, DGT_0, dofmap_restriction_DGT_0)) == 1\n'

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT_1 = DGT[1]
dofmap_restriction_DGT_1 = dofmap_restriction_DGT[facet_restriction_subset][1]
assert count_dofs(dofmap_restriction_DGT_1, mesh.comm) == 4
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DGT_1, dofmap_restriction_DGT_1)) == 2
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 1
"""

'Expect this cell to fail.\n\nDGT_1 = DGT[1]\ndofmap_restriction_DGT_1 = dofmap_restriction_DGT[facet_restriction_subset][1]\nassert count_dofs(dofmap_restriction_DGT_1, mesh.comm) == 4\nassert len(locate_dofs_by_polar_coordinates(\n    0, 0, DGT_1, dofmap_restriction_DGT_1)) == 2\nassert len(locate_dofs_by_polar_coordinates(\n    1, 2 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    1, 5 * np.pi / 3, DGT_1, dofmap_restriction_DGT_1)) == 1\n'

# PYTEST_XFAIL: Temporarily broken due to lack of suport for DGT in basix
"""Expect this cell to fail.

DGT_2 = DGT[2]
dofmap_restriction_DGT_2 = dofmap_restriction_DGT[facet_restriction_subset][2]
assert count_dofs(dofmap_restriction_DGT_2, mesh.comm) == 6
assert len(locate_dofs_by_polar_coordinates(
    0, 0, DGT_2, dofmap_restriction_DGT_2)) == 2
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1
"""

'Expect this cell to fail.\n\nDGT_2 = DGT[2]\ndofmap_restriction_DGT_2 = dofmap_restriction_DGT[facet_restriction_subset][2]\nassert count_dofs(dofmap_restriction_DGT_2, mesh.comm) == 6\nassert len(locate_dofs_by_polar_coordinates(\n    0, 0, DGT_2, dofmap_restriction_DGT_2)) == 2\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    0.5, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    1, 2 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\nassert len(locate_dofs_by_polar_coordinates(\n    1, 5 * np.pi / 3, DGT_2, dofmap_restriction_DGT_2)) == 1\n'