From fcf088c68feaefdd7e6a226b41ba8557d2ee695a Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Tue, 3 Mar 2026 13:58:24 +0000 Subject: [PATCH 01/10] Use DMPlex.createBoxMesh to create our periodic meshes This is necessary because the previous approach to generating periodic meshes resulted in a DMPlex with a slightly invalid state that prevents subsequent transformations (in particular extrusion). We also avoid manually labelling the boundaries in favour of just renumbering the existing 'Face Sets' label produced by the DMPlex. --- .github/workflows/core.yml | 4 +- docs/source/conf.py | 2 + firedrake/cython/dmcommon.pyx | 193 +++- firedrake/cython/petschdr.pxi | 2 + firedrake/utility_meshes.py | 997 +++++++++--------- .../test_equation_bcs_assemble.py | 22 +- .../test_variable_layers_numbering.py | 8 +- .../firedrake/multigrid/test_grid_transfer.py | 2 +- tests/firedrake/regression/test_bcs.py | 8 +- .../regression/test_mark_entities.py | 2 +- .../regression/test_mesh_generation.py | 42 +- .../firedrake/regression/test_periodic_2d.py | 50 +- .../submesh/test_submesh_assemble.py | 78 +- tests/firedrake/submesh/test_submesh_base.py | 2 +- 14 files changed, 773 insertions(+), 639 deletions(-) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index c0682518a4..e54f0f83aa 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -158,7 +158,9 @@ jobs: --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ https://gitlab.com/petsc/petsc.git else - git clone --depth 1 https://gitlab.com/petsc/petsc.git + : # UNDO ME + : # git clone --depth 1 https://gitlab.com/petsc/petsc.git + git clone --depth 1 https://gitlab.com/connorjward/petsc.git --branch connorjward/periodic-firedrake fi cd petsc python3 ../firedrake-repo/scripts/firedrake-configure \ diff --git a/docs/source/conf.py b/docs/source/conf.py index 76087d111d..024bec43cd 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -157,6 +157,8 @@ # Cofunction.ufl_domains references FormArgument but it isn't picked # up by Sphinx (see https://github.com/sphinx-doc/sphinx/issues/11225) ('py:class', 'FormArgument'), + # Some complex type hints confuse Sphinx (https://github.com/sphinx-doc/sphinx/issues/14159) + ("py:obj", r"typing\.Literal\[.*"), ] # Dodgy links diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ec7ab631af..ea76b9673a 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -788,6 +788,7 @@ def quadrilateral_closure_ordering(PETSc.DM plex, PetscInt nclosure, p, vi, v, fi, i PetscInt start_v, off PetscInt *closure = NULL + PetscInt closure_tmp[2*9] PetscInt c_vertices[4] PetscInt c_facets[4] PetscInt g_vertices[4] @@ -804,13 +805,13 @@ def quadrilateral_closure_ordering(PETSc.DM plex, ncells = cEnd - cStart entity_per_cell = 4 + 4 + 1 + CHKERR(PetscMalloc1(2*9, &closure)) + cell_closure = np.empty((ncells, entity_per_cell), dtype=IntType) for c in range(cStart, cEnd): CHKERR(PetscSectionGetOffset(cell_numbering.sec, c, &cell)) get_transitive_closure(plex.dm, c, PETSC_TRUE, &nclosure, &closure) - # First extract the facets (edges) and the vertices - # from the transitive closure into c_facets and c_vertices. # Here we assume that DMPlex gives entities in the order: # # 8--3--7 @@ -821,7 +822,65 @@ def quadrilateral_closure_ordering(PETSc.DM plex, # # where the starting vertex and order of traversal is arbitrary. # (We fix that later.) + + # If we have a periodic mesh with only a single cell in the periodic + # direction then the closure will look like + # + # 4--1--5 + # | | + # 3 0 2 (vertical periodicity) + # | | + # 4--1--5 + # + # or # + # 5--3--5 + # | | + # 2 0 2 (horizontal periodicity) + # | | + # 4--1--4 + # + # and only have 6 entries instead of 9. For the following to work we have + # to blow this out to a 9 entry array including the repeats. + if nclosure == 4: + raise NotImplementedError("Single-cell periodic quad meshes are " + "not supported") + elif nclosure == 6: + horiz_periodicity, vert_periodicity = _get_periodicity(plex) + (_, horiz_unit_periodic) = horiz_periodicity + (_, vert_unit_periodic) = vert_periodicity + if vert_unit_periodic: + assert not horiz_unit_periodic + closure_tmp[2*0] = closure[2*0] + closure_tmp[2*1] = closure[2*1] + closure_tmp[2*2] = closure[2*2] + closure_tmp[2*3] = closure[2*1] + closure_tmp[2*4] = closure[2*3] + closure_tmp[2*5] = closure[2*4] + closure_tmp[2*6] = closure[2*5] + closure_tmp[2*7] = closure[2*5] + closure_tmp[2*8] = closure[2*4] + else: + assert horiz_unit_periodic + assert not vert_unit_periodic + closure_tmp[2*0] = closure[2*0] + closure_tmp[2*1] = closure[2*1] + closure_tmp[2*2] = closure[2*2] + closure_tmp[2*3] = closure[2*3] + closure_tmp[2*4] = closure[2*2] + closure_tmp[2*5] = closure[2*4] + closure_tmp[2*6] = closure[2*4] + closure_tmp[2*7] = closure[2*5] + closure_tmp[2*8] = closure[2*5] + + nclosure = 9 + for i in range(9): + closure[2*i] = closure_tmp[2*i] + else: + assert nclosure == 9 + + # Extract the facets (edges) and the vertices + # from the transitive closure into c_facets and c_vertices. # For the vertices, we also retrieve the global numbers into g_vertices. vi = 0 fi = 0 @@ -923,8 +982,7 @@ def quadrilateral_closure_ordering(PETSc.DM plex, cell_closure[cell, 4 + 3] = facets[1] cell_closure[cell, 8] = c - if closure != NULL: - restore_transitive_closure(plex.dm, 0, PETSC_TRUE, &nclosure, &closure) + CHKERR(PetscFree(closure)) return cell_closure @@ -1987,7 +2045,7 @@ def reordered_coords(PETSc.DM dm, PETSc.Section global_numbering, shape, referen get_depth_stratum(dm.dm, 0, &vStart, &vEnd) if isinstance(dm, PETSc.DMPlex): if not dm.getCoordinatesLocalized(): - # Use CG coordiantes. + # Use CG coordinates. dm_sec = dm.getCoordinateSection() dm_coords = dm.getCoordinatesLocal().array.reshape(shape) coords = np.empty_like(dm_coords) @@ -1998,12 +2056,11 @@ def reordered_coords(PETSc.DM dm, PETSc.Section global_numbering, shape, referen for i in range(dim): coords[offset, i] = dm_coords[dm_offset, i] else: - # Use DG coordiantes. + # Use DG coordinates. get_height_stratum(dm.dm, 0, &cStart, &cEnd) dim = dm.getCoordinateDim() ndofs, perm, perm_offsets = _get_firedrake_plex_permutation_dg_transitive_closure(dm) - dm_sec = dm.getCellCoordinateSection() - dm_coords = dm.getCellCoordinatesLocal().array.reshape(((cEnd - cStart) * ndofs[0], dim)) + dm_coords, dm_sec = _get_expanded_dm_dg_coords(dm, ndofs) coords = np.empty_like(dm_coords) for c in range(cStart, cEnd): CHKERR(PetscSectionGetOffset(global_numbering.sec, c, &offset)) # scalar offset @@ -2031,6 +2088,126 @@ def reordered_coords(PETSc.DM dm, PETSc.Section global_numbering, shape, referen raise ValueError("Only DMPlex and DMSwarm are supported.") return coords + +def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): + cdef: + const PetscReal *L + + PETSc.Section dm_sec_expanded + + cStart, cEnd = dm.getHeightStratum(0) + dim = dm.getCoordinateDim() + + dm_coords = dm.getCellCoordinatesLocal().array_r + coords_shape = ((cEnd-cStart) * ndofs[0], dim) + + if dm_coords.size < np.prod(coords_shape, dtype=IntType): + # Fewer cell coordinates available, we must be single-cell periodic + if dm.getCellType(cStart) == PETSc.DM.PolytopeType.QUADRILATERAL: + # If we have a periodic mesh with only a single cell in the periodic + # direction then the cell coordinates will be + # + # 1-----2 + # | | + # | | (vertical periodicity) + # | | + # 1-----2 + # + # or + # + # 2-----2 + # | | + # | | (horizontal periodicity) + # | | + # 1-----1 + # + # when the standard layout is + # + # 4-----3 + # | | + # | | + # | | + # 1-----2 + assert ndofs[0] == 4, "Not expecting high order coords here" + dm_coords = dm_coords.reshape(((cEnd-cStart) * 2, dim)) + dm_coords_expanded = np.empty(coords_shape, dtype=dm_coords.dtype) + + # Create a new cell coordinate section + dm_sec_orig = dm.getCellCoordinateSection() + dm_sec_expanded = PETSc.Section().create(comm=dm_sec_orig.comm) + dm_sec_expanded.setChart(*dm_sec_orig.getChart()) + dm_sec_expanded.setPermutation(dm_sec_orig.getPermutation()) + + horiz_periodicity, vert_periodicity = _get_periodicity(dm) + (_, horiz_unit_periodic) = horiz_periodicity + (_, vert_unit_periodic) = vert_periodicity + + # Find the domain sizes + CHKERR(DMGetPeriodicity(dm.dm, NULL, NULL, &L)) + + if horiz_unit_periodic: + if vert_unit_periodic: + raise NotImplementedError("Single-cell periodic quad meshes are " + "not supported") + else: + cell_width = L[0] + + for c in range(cStart, cEnd): + CHKERR(PetscSectionSetDof(dm_sec_expanded.sec, c, 8)) + + dm_coords_expanded[4*c+0, 0] = dm_coords[2*c+0, 0] + dm_coords_expanded[4*c+1, 0] = dm_coords[2*c+0, 0] + cell_width + dm_coords_expanded[4*c+2, 0] = dm_coords[2*c+1, 0] + cell_width + dm_coords_expanded[4*c+3, 0] = dm_coords[2*c+1, 0] + dm_coords_expanded[4*c+0, 1] = dm_coords[2*c+0, 1] + dm_coords_expanded[4*c+1, 1] = dm_coords[2*c+0, 1] + dm_coords_expanded[4*c+2, 1] = dm_coords[2*c+1, 1] + dm_coords_expanded[4*c+3, 1] = dm_coords[2*c+1, 1] + + else: + if vert_unit_periodic: + cell_height = L[1] + + for c in range(cStart, cEnd): + CHKERR(PetscSectionSetDof(dm_sec_expanded.sec, c, 8)) + + dm_coords_expanded[4*c+0, 0] = dm_coords[2*c+0, 0] + dm_coords_expanded[4*c+1, 0] = dm_coords[2*c+1, 0] + dm_coords_expanded[4*c+2, 0] = dm_coords[2*c+1, 0] + dm_coords_expanded[4*c+3, 0] = dm_coords[2*c+0, 0] + dm_coords_expanded[4*c+0, 1] = dm_coords[2*c+0, 1] + dm_coords_expanded[4*c+1, 1] = dm_coords[2*c+1, 1] + dm_coords_expanded[4*c+2, 1] = dm_coords[2*c+1, 1] + cell_height + dm_coords_expanded[4*c+3, 1] = dm_coords[2*c+0, 1] + cell_height + else: + raise AssertionError("This case should not be hit") + + dm_sec_expanded.setUp() + + dm_coords = dm_coords_expanded + dm_sec = dm_sec_expanded + + else: + raise NotImplementedError("Single cell periodicity for cell type " + f"{dm.getCellType(cStart)} is not supported") + + else: + dm_coords = dm_coords.reshape(coords_shape) + dm_sec = dm.getCellCoordinateSection() + + return dm_coords, dm_sec + +def _get_periodicity(dm: PETSc.DM) -> tuple[tuple[bool, bool], ...]: + cdef: + const PetscReal *maxCell, *L + + dim = dm.getCoordinateDim() + CHKERR(DMGetPeriodicity(dm.dm, &maxCell, NULL, &L)) + return tuple( + (L[d] >= 0, maxCell[d] >= L[d]) + for d in range(dim) + ) + @cython.boundscheck(False) @cython.wraparound(False) def mark_entity_classes(PETSc.DM dm): diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index d7b1478ca0..f50cc06811 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -103,6 +103,8 @@ cdef extern from "petscdm.h" nogil: PetscErrorCode DMSetLabelValue(PETSc.PetscDM,char[],PetscInt,PetscInt) PetscErrorCode DMGetLabelValue(PETSc.PetscDM,char[],PetscInt,PetscInt*) + PetscErrorCode DMGetPeriodicity(PETSc.PetscDM,PetscReal *[], PetscReal *[], PetscReal *[]) + cdef extern from "petscdmswarm.h" nogil: PetscErrorCode DMSwarmGetLocalSize(PETSc.PetscDM,PetscInt*) PetscErrorCode DMSwarmGetCellDM(PETSc.PetscDM, PETSc.PetscDM*) diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index f75cd514ae..9294868c5e 100644 --- a/firedrake/utility_meshes.py +++ b/firedrake/utility_meshes.py @@ -1,40 +1,31 @@ +import numbers import numpy as np +import warnings +from typing import Literal +import petsctools import ufl +from mpi4py import MPI from pyop2.mpi import COMM_WORLD from firedrake.utils import IntType, ScalarType from firedrake import ( VectorFunctionSpace, - FunctionSpace, Function, Constant, assemble, interpolate, FiniteElement, - interval, tetrahedron, - atan2, - pi, - as_vector, - SpatialCoordinate, - conditional, - gt, - as_tensor, - dot, - And, - Or, - sin, - cos, real ) from firedrake.cython import dmcommon from firedrake.mesh import ( - Mesh, DistributedMeshOverlapType, DEFAULT_MESH_NAME, + Mesh, DistributedMeshOverlapType, DEFAULT_MESH_NAME, MeshGeometry, plex_from_cell_list, _generate_default_mesh_topology_name, _generate_default_mesh_coordinates_name, MeshTopology, - make_mesh_from_mesh_topology, RelabeledMesh, + make_mesh_from_mesh_topology, make_mesh_from_coordinates, ExtrudedMesh ) from firedrake.parameters import parameters @@ -150,35 +141,22 @@ def IntervalMesh( left = length_or_left if ncells <= 0 or ncells % 1: - raise ValueError("Number of cells must be a postive integer") - length = right - left - if length < 0: + raise ValueError("Number of cells must be a positive integer") + if right - left < 0: raise ValueError("Requested mesh has negative length") - dx = length / ncells - # This ensures the rightmost point is actually present. - coords = np.arange(left, right + 0.01 * dx, dx, dtype=np.double).reshape(-1, 1) - cells = np.dstack( - ( - np.arange(0, len(coords) - 1, dtype=np.int32), - np.arange(1, len(coords), dtype=np.int32), - ) - ).reshape(-1, 2) - plex = plex_from_cell_list( - 1, cells, coords, comm, _generate_default_mesh_topology_name(name) + + plex = PETSc.DMPlex().createBoxMesh( + (ncells,), + lower=(left,), + upper=(right,), + simplex=False, + periodic=False, + interpolate=True, + comm=comm ) - # Apply boundary IDs - plex.createLabel(dmcommon.FACE_SETS_LABEL) - coordinates = plex.getCoordinates() - coord_sec = plex.getCoordinateSection() - vStart, vEnd = plex.getDepthStratum(0) # vertices - for v in range(vStart, vEnd): - vcoord = plex.vecGetClosure(coord_sec, coordinates, v) - if vcoord[0] == coords[0]: - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, v, 1) - if vcoord[0] == coords[-1]: - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, v, 2) + _mark_mesh_boundaries(plex) - m = Mesh( + return Mesh( plex, reorder=reorder, distribution_parameters=distribution_parameters, @@ -187,7 +165,6 @@ def IntervalMesh( permutation_name=permutation_name, comm=comm, ) - return m @PETSc.Log.EventDecorator() @@ -219,7 +196,6 @@ def UnitIntervalMesh( The left hand (:math:`x=0`) boundary point has boundary marker 1, while the right hand (:math:`x=1`) point has marker 2. """ - return IntervalMesh( ncells, length_or_left=1.0, @@ -259,46 +235,27 @@ def PeriodicIntervalMesh( when checkpointing; if `None`, the name is automatically generated. """ + plex = PETSc.DMPlex().createBoxMesh( + (ncells,), + lower=(0.,), + upper=(length,), + simplex=False, + periodic=True, + interpolate=True, + sparseLocalize=False, + comm=comm + ) + _mark_mesh_boundaries(plex) - if ncells < 3: - raise ValueError( - "1D periodic meshes with fewer than 3 \ -cells are not currently supported" - ) - m = CircleManifoldMesh( - ncells, - distribution_parameters=distribution_parameters_no_overlap, - reorder=reorder_noop, - comm=comm, + return Mesh( + plex, + reorder=reorder, + distribution_parameters=distribution_parameters, name=name, distribution_name=distribution_name, permutation_name=permutation_name, + comm=comm, ) - indicator = Function(FunctionSpace(m, "DG", 0)) - coord_fs = VectorFunctionSpace( - m, FiniteElement("DG", interval, 1, variant="equispaced"), dim=1 - ) - new_coordinates = Function( - coord_fs, name=_generate_default_mesh_coordinates_name(name) - ) - x, y = SpatialCoordinate(m) - eps = 1.e-14 - indicator.interpolate(conditional(gt(real(y), 0), 0., 1.)) - new_coordinates.interpolate( - as_vector((conditional( - gt(real(x), real(1. - eps)), indicator, # Periodic break. - # Unwrap rest of circle. - atan2(real(-y), real(-x))/(2 * pi) + 0.5 - ) * length,)) - ) - - return _postprocess_periodic_mesh(new_coordinates, - comm, - distribution_parameters, - reorder, - name, - distribution_name, - permutation_name) @PETSc.Log.EventDecorator() @@ -597,7 +554,7 @@ def RectangleMesh( originY=0., quadrilateral=False, reorder=None, - diagonal="left", + diagonal=None, distribution_parameters=None, comm=COMM_WORLD, name=DEFAULT_MESH_NAME, @@ -635,24 +592,34 @@ def RectangleMesh( * 3: plane y == originY * 4: plane y == Ly """ + if any(n <= 0 or not isinstance(n, numbers.Integral) for n in {nx, ny}): + raise ValueError("Number of cells must be a positive integer") + if quadrilateral and diagonal is not None: + raise ValueError("Cannot specify slope of diagonal on quad meshes") + if not quadrilateral and diagonal is None: + diagonal = "left" - for n in (nx, ny): - if n <= 0 or n % 1: - raise ValueError("Number of cells must be a postive integer") + plex = PETSc.DMPlex().createBoxMesh( + (nx, ny), + lower=(originX, originY), + upper=(Lx, Ly), + simplex=False, + interpolate=True, + comm=comm + ) + _mark_mesh_boundaries(plex) - xcoords = np.linspace(originX, Lx, nx + 1, dtype=np.double) - ycoords = np.linspace(originY, Ly, ny + 1, dtype=np.double) - return TensorRectangleMesh( - xcoords, - ycoords, - quadrilateral=quadrilateral, + if not quadrilateral: + plex = _refine_quads_to_triangles(plex, diagonal) + + return Mesh( + plex, reorder=reorder, - diagonal=diagonal, distribution_parameters=distribution_parameters, - comm=comm, name=name, distribution_name=distribution_name, permutation_name=permutation_name, + comm=comm, ) @@ -661,7 +628,7 @@ def TensorRectangleMesh( ycoords, quadrilateral=False, reorder=None, - diagonal="left", + diagonal=None, distribution_parameters=None, comm=COMM_WORLD, name=DEFAULT_MESH_NAME, @@ -680,6 +647,7 @@ def TensorRectangleMesh( :kwarg diagonal: For triangular meshes, should the diagonal got from bottom left to top right (``"right"``), or top left to bottom right (``"left"``), or put in both diagonals (``"crossed"``). + Default is ``"left"``. The boundary edges in this mesh are numbered as follows: @@ -688,61 +656,33 @@ def TensorRectangleMesh( * 3: plane y == ycoords[0] * 4: plane y == ycoords[-1] """ + if quadrilateral and diagonal is not None: + raise ValueError("Cannot specify slope of diagonal on quad meshes") + if not quadrilateral and diagonal is None: + diagonal = "left" + xcoords = np.unique(xcoords) ycoords = np.unique(ycoords) nx = np.size(xcoords) - 1 ny = np.size(ycoords) - 1 - for n in (nx, ny): - if n <= 0: - raise ValueError("Number of cells must be a postive integer") + if any(n <= 0 for n in {nx, ny}): + raise ValueError("Number of cells must be a positive integer") coords = np.asarray(np.meshgrid(xcoords, ycoords)).swapaxes(0, 2).reshape(-1, 2) # cell vertices i, j = np.meshgrid(np.arange(nx, dtype=np.int32), np.arange(ny, dtype=np.int32)) - if not quadrilateral and diagonal == "crossed": - xs = 0.5 * (xcoords[1:] + xcoords[:-1]) - ys = 0.5 * (ycoords[1:] + ycoords[:-1]) - extra = np.asarray(np.meshgrid(xs, ys)).swapaxes(0, 2).reshape(-1, 2) - coords = np.vstack([coords, extra]) - # - # 2-----3 - # | \ / | - # | 4 | - # | / \ | - # 0-----1 - cells = [ - i * (ny + 1) + j, - i * (ny + 1) + j + 1, - (i + 1) * (ny + 1) + j, - (i + 1) * (ny + 1) + j + 1, - (nx + 1) * (ny + 1) + i * ny + j, - ] - cells = np.asarray(cells).swapaxes(0, 2).reshape(-1, 5) - idx = [0, 1, 4, 0, 2, 4, 2, 3, 4, 3, 1, 4] - cells = cells[:, idx].reshape(-1, 3) - else: - cells = [ - i * (ny + 1) + j, - i * (ny + 1) + j + 1, - (i + 1) * (ny + 1) + j + 1, - (i + 1) * (ny + 1) + j, - ] - cells = np.asarray(cells).swapaxes(0, 2).reshape(-1, 4) - if not quadrilateral: - if diagonal == "left": - idx = [0, 1, 3, 1, 2, 3] - elif diagonal == "right": - idx = [0, 1, 2, 0, 2, 3] - else: - raise ValueError("Unrecognised value for diagonal '%r'", diagonal) - # two cells per cell above... - cells = cells[:, idx].reshape(-1, 3) + cells = [ + i * (ny + 1) + j, + i * (ny + 1) + j + 1, + (i + 1) * (ny + 1) + j + 1, + (i + 1) * (ny + 1) + j, + ] + cells = np.asarray(cells).swapaxes(0, 2).reshape(-1, 4) plex = plex_from_cell_list( 2, cells, coords, comm, _generate_default_mesh_topology_name(name) ) - # mark boundary facets plex.createLabel(dmcommon.FACE_SETS_LABEL) plex.markBoundaryFaces("boundary_faces") @@ -767,7 +707,11 @@ def TensorRectangleMesh( if abs(face_coords[1] - y1) < ytol and abs(face_coords[3] - y1) < ytol: plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 4) plex.removeLabel("boundary_faces") - m = Mesh( + + if not quadrilateral: + plex = _refine_quads_to_triangles(plex, diagonal) + + return Mesh( plex, reorder=reorder, distribution_parameters=distribution_parameters, @@ -776,40 +720,59 @@ def TensorRectangleMesh( permutation_name=permutation_name, comm=comm, ) - return m @PETSc.Log.EventDecorator() def SquareMesh( - nx, - ny, - L, - reorder=None, - quadrilateral=False, - diagonal="left", - distribution_parameters=None, - comm=COMM_WORLD, - name=DEFAULT_MESH_NAME, - distribution_name=None, - permutation_name=None, + nx: numbers.Integral, + ny: numbers.Integral, + L: numbers.Real, + reorder: bool | None = None, + quadrilateral: bool = False, + diagonal: Literal['crossed', 'left', 'right'] | None = None, + distribution_parameters: dict | None = None, + comm: MPI.Comm = COMM_WORLD, + name: str = DEFAULT_MESH_NAME, + distribution_name: str | None = None, + permutation_name: str | None = None, ): - """Generate a square mesh + """Generate a square mesh. - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :arg L: The extent in the x and y directions - :kwarg quadrilateral: (optional), creates quadrilateral mesh. - :kwarg reorder: (optional), should the mesh be reordered - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx + The number of cells in the x direction. + ny + The number of cells in the y direction. + L + The extent in the x and y directions. + reorder + Flag indicating whether to reorder the mesh. + quadrilateral + Flag indicating whether to create a quadrilateral mesh. + diagonal + The refinement strategy used for non-quadrilateral meshes. Must be + one of ``"crossed"``, ``"left"``, ``"right"``. + distribution_parameters + Options controlling mesh distribution, see :func:`.Mesh` for details. + comm + Optional communicator to build the mesh on. + name + Optional name of the mesh. + distribution_name + The name of parallel distribution used when checkpointing; if `None`, + the name is automatically generated. + permutation_name + The name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. + + Returns + ------- + MeshGeometry + The new mesh. + + Notes + ----- The boundary edges in this mesh are numbered as follows: @@ -817,6 +780,7 @@ def SquareMesh( * 2: plane x == L * 3: plane y == 0 * 4: plane y == L + """ return RectangleMesh( nx, @@ -836,33 +800,52 @@ def SquareMesh( @PETSc.Log.EventDecorator() def UnitSquareMesh( - nx, - ny, - reorder=None, - diagonal="left", - quadrilateral=False, - distribution_parameters=None, - comm=COMM_WORLD, - name=DEFAULT_MESH_NAME, - distribution_name=None, - permutation_name=None, + nx: numbers.Integral, + ny: numbers.Integral, + reorder: bool | None = None, + diagonal: Literal["crossed", "left", "right"] | None = None, + quadrilateral: bool = False, + distribution_parameters: dict | None = None, + comm: MPI.Comm = COMM_WORLD, + name: str = DEFAULT_MESH_NAME, + distribution_name: str | None = None, + permutation_name: str | None = None, ): - """Generate a unit square mesh + """Generate a unit square mesh. - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :kwarg quadrilateral: (optional), creates quadrilateral mesh. - :kwarg reorder: (optional), should the mesh be reordered - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx + The number of cells in the x direction. + ny + The number of cells in the y direction. + reorder + Flag indicating whether to reorder the mesh. + diagonal + The refinement strategy used for non-quadrilateral meshes. Must be + one of ``"crossed"``, ``"left"``, ``"right"``. + quadrilateral + Flag indicating whether to create a quadrilateral mesh. + distribution_parameters + Options controlling mesh distribution, see :func:`.Mesh` for details. + comm + Optional communicator to build the mesh on. + name + Optional name of the mesh. + distribution_name + The name of parallel distribution used when checkpointing; if `None`, + the name is automatically generated. + permutation_name + The name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. + + Returns + ------- + MeshGeometry + The new mesh. + + Notes + ----- The boundary edges in this mesh are numbered as follows: @@ -870,6 +853,7 @@ def UnitSquareMesh( * 2: plane x == 1 * 3: plane y == 0 * 4: plane y == 1 + """ return SquareMesh( nx, @@ -888,191 +872,179 @@ def UnitSquareMesh( @PETSc.Log.EventDecorator() def PeriodicRectangleMesh( - nx, - ny, - Lx, - Ly, - direction="both", - quadrilateral=False, - reorder=None, - distribution_parameters=None, - diagonal=None, - comm=COMM_WORLD, - name=DEFAULT_MESH_NAME, - distribution_name=None, - permutation_name=None, -): - """Generate a periodic rectangular mesh + nx: numbers.Integral, + ny: numbers.Integral, + Lx: numbers.Real, + Ly: numbers.Real, + direction: Literal["both", "x", "y"] = "both", + quadrilateral: bool = False, + reorder: bool | None = None, + distribution_parameters: dict | None = None, + diagonal: Literal["crossed", "left", "right"] | None = None, + comm: MPI.Comm = COMM_WORLD, + name: str = DEFAULT_MESH_NAME, + distribution_name: str | None = None, + permutation_name: str | None = None, +) -> MeshGeometry: + """Generate a periodic rectangular mesh. - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :arg Lx: The extent in the x direction - :arg Ly: The extent in the y direction - :arg direction: The direction of the periodicity, one of - ``"both"``, ``"x"`` or ``"y"``. - :kwarg quadrilateral: (optional), creates quadrilateral mesh. - :kwarg reorder: (optional), should the mesh be reordered - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg diagonal: (optional), one of ``"crossed"``, ``"left"``, ``"right"``. - Not valid for quad meshes. Only used for direction ``"x"`` or direction ``"y"``. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx + The number of cells in the x direction. + ny + The number of cells in the y direction. + Lx + The extent in the x direction. + Ly + The extent in the y direction. + direction + The direction of the periodicity, one of ``"both"``, ``"x"`` or ``"y"``. + quadrilateral + Flag indicating whether to create a quadrilateral mesh. + reorder + Flag indicating whether to reorder the mesh. + distribution_parameters + Options controlling mesh distribution, see :func:`.Mesh` for details. + diagonal + The refinement strategy used for non-quadrilateral meshes. Must be + one of ``"crossed"``, ``"left"``, ``"right"``. + comm + Optional communicator to build the mesh on. + name + Optional name of the mesh. + distribution_name + The name of parallel distribution used when checkpointing; if `None`, + the name is automatically generated. + permutation_name + The name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. - If direction == "x" the boundary edges in this mesh are numbered as follows: + Returns + ------- + MeshGeometry + The new mesh. - * 1: plane y == 0 - * 2: plane y == Ly + Notes + ----- - If direction == "y" the boundary edges are: + The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 - * 2: plane x == Lx + * 2: plane x == 1 + * 3: plane y == 0 + * 4: plane y == 1 + + If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and + if periodic in 'y' then 3 and 4 are empty. + """ + if quadrilateral and diagonal is not None: + raise ValueError("Cannot specify slope of diagonal on quad meshes") + if not quadrilateral and diagonal is None: + diagonal = "left" - if direction == "both" and ny == 1 and quadrilateral: - return OneElementThickMesh( - nx, - Lx, - Ly, - distribution_parameters=distribution_parameters, - name=name, - distribution_name=distribution_name, - permutation_name=permutation_name, - comm=comm, - ) + match direction: + case "both": + periodic = (True, True) + case "x": + periodic = (True, False) + case "y": + periodic = (False, True) + case _: + raise ValueError( + f"Cannot have a periodic mesh with periodicity '{direction}'" + ) - if direction not in ("both", "x", "y"): - raise ValueError( - "Cannot have a periodic mesh with periodicity '%s'" % direction - ) - if direction != "both": - return PartiallyPeriodicRectangleMesh( - nx, - ny, - Lx, - Ly, - direction=direction, - quadrilateral=quadrilateral, - reorder=reorder, - distribution_parameters=distribution_parameters, - diagonal=diagonal, - comm=comm, - name=name, - distribution_name=distribution_name, - permutation_name=permutation_name, - ) - if nx < 3 or ny < 3: - raise ValueError( - "2D periodic meshes with fewer than 3 cells in each direction are not currently supported" - ) + plex = PETSc.DMPlex().createBoxMesh( + (nx, ny), + lower=(0., 0.), + upper=(Lx, Ly), + simplex=False, + periodic=periodic, + interpolate=True, + sparseLocalize=False, + comm=comm + ) + _mark_mesh_boundaries(plex) + if not quadrilateral: + plex = _refine_quads_to_triangles(plex, diagonal) - m = TorusMesh( - nx, - ny, - 1.0, - 0.5, - quadrilateral=quadrilateral, - reorder=reorder_noop, - distribution_parameters=distribution_parameters_no_overlap, - comm=comm, + return Mesh( + plex, + reorder=reorder, + distribution_parameters=distribution_parameters, name=name, distribution_name=distribution_name, permutation_name=permutation_name, + comm=comm, ) - coord_family = "DQ" if quadrilateral else "DG" - cell = "quadrilateral" if quadrilateral else "triangle" - - coord_fs = VectorFunctionSpace( - m, FiniteElement(coord_family, cell, 1, variant="equispaced"), dim=2 - ) - new_coordinates = Function( - coord_fs, name=_generate_default_mesh_coordinates_name(name) - ) - x, y, z = SpatialCoordinate(m) - eps = 1.e-14 - indicator_y = Function(FunctionSpace(m, coord_family, 0)) - indicator_y.interpolate(conditional(gt(real(y), 0), 0., 1.)) - x_coord = Function(FunctionSpace(m, coord_family, 1, variant="equispaced")) - x_coord.interpolate( - # Periodic break. - conditional(And(gt(real(eps), real(abs(y))), gt(real(x), 0.)), indicator_y, - # Unwrap rest of circle. - atan2(real(-y), real(-x))/(2*pi)+0.5) - ) - phi_coord = as_vector([cos(2*pi*x_coord), sin(2*pi*x_coord)]) - dr = dot(as_vector((x, y))-phi_coord, phi_coord) - indicator_z = Function(FunctionSpace(m, coord_family, 0)) - indicator_z.interpolate(conditional(gt(real(z), 0), 0., 1.)) - new_coordinates.interpolate(as_vector(( - x_coord * Lx, - # Periodic break. - conditional(And(gt(real(eps), real(abs(z))), gt(real(dr), 0.)), indicator_z, - # Unwrap rest of circle. - atan2(real(-z), real(-dr))/(2*pi)+0.5) * Ly - ))) - - return _postprocess_periodic_mesh(new_coordinates, - comm, - distribution_parameters, - reorder, - name, - distribution_name, - permutation_name) @PETSc.Log.EventDecorator() def PeriodicSquareMesh( - nx, - ny, - L, - direction="both", - quadrilateral=False, - reorder=None, - distribution_parameters=None, - diagonal=None, - comm=COMM_WORLD, - name=DEFAULT_MESH_NAME, - distribution_name=None, - permutation_name=None, + nx: numbers.Integral, + ny: numbers.Integral, + L: numbers.Real, + direction: Literal["both", "x", "y"] = "both", + quadrilateral: bool = False, + reorder: bool | None = None, + distribution_parameters: dict | None = None, + diagonal: Literal["crossed", "left", "right"] | None = None, + comm: MPI.Comm = COMM_WORLD, + name: str = DEFAULT_MESH_NAME, + distribution_name: str | None = None, + permutation_name: str | None = None, ): - """Generate a periodic square mesh + """Generate a periodic square mesh. - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :arg L: The extent in the x and y directions - :arg direction: The direction of the periodicity, one of - ``"both"``, ``"x"`` or ``"y"``. - :kwarg quadrilateral: (optional), creates quadrilateral mesh. - :kwarg reorder: (optional), should the mesh be reordered - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg diagonal: (optional), one of ``"crossed"``, ``"left"``, ``"right"``. - Not valid for quad meshes. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx + The number of cells in the x direction. + ny + The number of cells in the y direction. + L + The extent in the x and y directions. + direction + The direction of the periodicity, one of ``"both"``, ``"x"`` or ``"y"``. + quadrilateral + Flag indicating whether to create a quadrilateral mesh. + reorder + Flag indicating whether to reorder the mesh. + distribution_parameters + Options controlling mesh distribution, see :func:`.Mesh` for details. + diagonal + The refinement strategy used for non-quadrilateral meshes. Must be + one of ``"crossed"``, ``"left"``, ``"right"``. + comm + Optional communicator to build the mesh on. + name + Optional name of the mesh. + distribution_name + The name of parallel distribution used when checkpointing; if `None`, + the name is automatically generated. + permutation_name + The name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. - If direction == "x" the boundary edges in this mesh are numbered as follows: + Returns + ------- + MeshGeometry + The new mesh. - * 1: plane y == 0 - * 2: plane y == L + Notes + ----- - If direction == "y" the boundary edges are: + The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 - * 2: plane x == L + * 2: plane x == 1 + * 3: plane y == 0 + * 4: plane y == 1 + + If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and + if periodic in 'y' then 3 and 4 are empty. """ return PeriodicRectangleMesh( nx, @@ -1093,48 +1065,65 @@ def PeriodicSquareMesh( @PETSc.Log.EventDecorator() def PeriodicUnitSquareMesh( - nx, - ny, - direction="both", - reorder=None, - quadrilateral=False, - distribution_parameters=None, - diagonal=None, - comm=COMM_WORLD, - name=DEFAULT_MESH_NAME, - distribution_name=None, - permutation_name=None, + nx: numbers.Integral, + ny: numbers.Integral, + direction: Literal["both", "x", "y"] = "both", + quadrilateral: bool = False, + reorder: bool | None = None, + distribution_parameters: dict | None = None, + diagonal: Literal["crossed", "left", "right"] | None = None, + comm: MPI.Comm = COMM_WORLD, + name: str = DEFAULT_MESH_NAME, + distribution_name: str | None = None, + permutation_name: str | None = None, ): - """Generate a periodic unit square mesh + """Generate a periodic unit square mesh. - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :arg direction: The direction of the periodicity, one of - ``"both"``, ``"x"`` or ``"y"``. - :kwarg quadrilateral: (optional), creates quadrilateral mesh. - :kwarg reorder: (optional), should the mesh be reordered - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg diagonal: (optional), one of ``"crossed"``, ``"left"``, ``"right"``. - Not valid for quad meshes. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx + The number of cells in the x direction. + ny + The number of cells in the y direction. + direction + The direction of the periodicity, one of ``"both"``, ``"x"`` or ``"y"``. + quadrilateral + Flag indicating whether to create a quadrilateral mesh. + reorder + Flag indicating whether to reorder the mesh. + distribution_parameters + Options controlling mesh distribution, see :func:`.Mesh` for details. + diagonal + The refinement strategy used for non-quadrilateral meshes. Must be + one of ``"crossed"``, ``"left"``, ``"right"``. + comm + Optional communicator to build the mesh on. + name + Optional name of the mesh. + distribution_name + The name of parallel distribution used when checkpointing; if `None`, + the name is automatically generated. + permutation_name + The name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. - If direction == "x" the boundary edges in this mesh are numbered as follows: + Returns + ------- + MeshGeometry + The new mesh. - * 1: plane y == 0 - * 2: plane y == 1 + Notes + ----- - If direction == "y" the boundary edges are: + The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 * 2: plane x == 1 + * 3: plane y == 0 + * 4: plane y == 1 + + If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and + if periodic in 'y' then 3 and 4 are empty. """ return PeriodicSquareMesh( nx, @@ -1626,38 +1615,17 @@ def BoxMesh( if n <= 0 or n % 1: raise ValueError("Number of cells must be a postive integer") if hexahedral: - plex = PETSc.DMPlex().createBoxMesh((nx, ny, nz), lower=(0., 0., 0.), upper=(Lx, Ly, Lz), simplex=False, periodic=False, interpolate=True, comm=comm) - plex.removeLabel(dmcommon.FACE_SETS_LABEL) - nvert = 4 # num. vertices on faect - - # Apply boundary IDs - plex.createLabel(dmcommon.FACE_SETS_LABEL) - plex.markBoundaryFaces("boundary_faces") - coords = plex.getCoordinates() - coord_sec = plex.getCoordinateSection() - cdim = plex.getCoordinateDim() - assert cdim == 3 - if plex.getStratumSize("boundary_faces", 1) > 0: - boundary_faces = plex.getStratumIS("boundary_faces", 1).getIndices() - xtol = Lx / (2 * nx) - ytol = Ly / (2 * ny) - ztol = Lz / (2 * nz) - for face in boundary_faces: - face_coords = plex.vecGetClosure(coord_sec, coords, face) - if all([abs(face_coords[0 + cdim * i]) < xtol for i in range(nvert)]): - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 1) - if all([abs(face_coords[0 + cdim * i] - Lx) < xtol for i in range(nvert)]): - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 2) - if all([abs(face_coords[1 + cdim * i]) < ytol for i in range(nvert)]): - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 3) - if all([abs(face_coords[1 + cdim * i] - Ly) < ytol for i in range(nvert)]): - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 4) - if all([abs(face_coords[2 + cdim * i]) < ztol for i in range(nvert)]): - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 5) - if all([abs(face_coords[2 + cdim * i] - Lz) < ztol for i in range(nvert)]): - plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 6) - plex.removeLabel("boundary_faces") - m = Mesh( + plex = PETSc.DMPlex().createBoxMesh( + (nx, ny, nz), + lower=(0., 0., 0.), + upper=(Lx, Ly, Lz), + simplex=False, + periodic=False, + interpolate=True, + comm=comm, + ) + _mark_mesh_boundaries(plex) + return Mesh( plex, reorder=reorder, distribution_parameters=distribution_parameters, @@ -1666,7 +1634,6 @@ def BoxMesh( permutation_name=permutation_name, comm=comm, ) - return m else: xcoords = np.linspace(0, Lx, nx + 1, dtype=np.double) ycoords = np.linspace(0, Ly, ny + 1, dtype=np.double) @@ -1868,14 +1835,9 @@ def PeriodicBoxMesh( * 5: plane z == 0 * 6: plane z == 1 - where periodic surfaces are regarded as interior, for which dS integral is to be used. + where periodic surfaces are ignored. """ - for n in (nx, ny, nz): - if n < 3: - raise ValueError( - "3D periodic meshes with fewer than 3 cells are not currently supported" - ) if hexahedral: if len(directions) != 3: raise ValueError(f"directions must have exactly dim (=3) elements : Got {directions}") @@ -1889,7 +1851,8 @@ def PeriodicBoxMesh( sparseLocalize=False, comm=comm, ) - m = Mesh( + _mark_mesh_boundaries(plex) + return Mesh( plex, reorder=reorder, distribution_parameters=distribution_parameters, @@ -1897,29 +1860,8 @@ def PeriodicBoxMesh( distribution_name=distribution_name, permutation_name=permutation_name, comm=comm) - x, y, z = SpatialCoordinate(m) - V = FunctionSpace(m, "Q", 2) - eps = min([Lx / nx, Ly / ny, Lz / nz]) / 1000. - if directions[0]: # x - fx0 = Function(V).interpolate(conditional(Or(x < eps, x > Lx - eps), 1., 0.)) - fx1 = fx0 - else: - fx0 = Function(V).interpolate(conditional(x < eps, 1., 0.)) - fx1 = Function(V).interpolate(conditional(x > Lx - eps, 1., 0.)) - if directions[1]: # y - fy0 = Function(V).interpolate(conditional(Or(y < eps, y > Ly - eps), 1., 0.)) - fy1 = fy0 - else: - fy0 = Function(V).interpolate(conditional(y < eps, 1., 0.)) - fy1 = Function(V).interpolate(conditional(y > Ly - eps, 1., 0.)) - if directions[2]: # z - fz0 = Function(V).interpolate(conditional(Or(z < eps, z > Lz - eps), 1., 0.)) - fz1 = fz0 - else: - fz0 = Function(V).interpolate(conditional(z < eps, 1., 0.)) - fz1 = Function(V).interpolate(conditional(z > Lz - eps, 1., 0.)) - return RelabeledMesh(m, [fx0, fx1, fy0, fy1, fz0, fz1], [1, 2, 3, 4, 5, 6], name=name) else: + # TODO: When hexahedra -> simplex refinement is implemented this can go away. if tuple(directions) != (True, True, True): raise NotImplementedError("Can only specify directions with hexahedral = True") xcoords = np.arange(0.0, Lx, Lx / nx, dtype=np.double) @@ -2059,7 +2001,7 @@ def PeriodicUnitCubeMesh( * 5: plane z == 0 * 6: plane z == 1 - where periodic surfaces are regarded as interior, for which dS integral is to be used. + where periodic surfaces are ignored. """ return PeriodicBoxMesh( @@ -3106,74 +3048,123 @@ def PartiallyPeriodicRectangleMesh( when checkpointing; if `None`, the name is automatically generated. - If direction == "x" the boundary edges in this mesh are numbered as follows: - - * 1: plane y == 0 - * 2: plane y == Ly - - If direction == "y" the boundary edges are: + The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 - * 2: plane x == Lx - """ - if direction not in ("x", "y"): - raise ValueError("Unsupported periodic direction '%s'" % direction) + * 2: plane x == 1 + * 3: plane y == 0 + * 4: plane y == 1 - # handle x/y directions: na, La are for the periodic axis - na, nb = nx, ny - if direction == "y": - na, nb = ny, nx + If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and + if periodic in 'y' then 3 and 4 are empty. - if na < 3: - raise ValueError( - "2D periodic meshes with fewer than 3 cells in each direction are not currently supported" - ) + """ + warnings.warn( + "'PartiallyPeriodicRectangleMesh' is deprecated. Please use " + "'PeriodicRectangleMesh' instead, passing 'direction=\"x\"' or " + "'direction=\"y\"'.", + FutureWarning, + ) - m = CylinderMesh( - na, - nb, - 1.0, - 1.0, - longitudinal_direction="z", + return PeriodicRectangleMesh( + nx, + ny, + Lx, + Ly, + direction=direction, quadrilateral=quadrilateral, - reorder=reorder_noop, - distribution_parameters=distribution_parameters_no_overlap, + reorder=reorder, + distribution_parameters=distribution_parameters, diagonal=diagonal, comm=comm, name=name, distribution_name=distribution_name, permutation_name=permutation_name, ) - coord_family = "DQ" if quadrilateral else "DG" - cell = "quadrilateral" if quadrilateral else "triangle" - indicator = Function(FunctionSpace(m, coord_family, 0)) - coord_fs = VectorFunctionSpace( - m, FiniteElement(coord_family, cell, 1, variant="equispaced"), dim=2 - ) - new_coordinates = Function( - coord_fs, name=_generate_default_mesh_coordinates_name(name) - ) - x, y, z = SpatialCoordinate(m) - eps = 1.e-14 - indicator.interpolate(conditional(gt(real(y), 0), 0., 1.)) - if direction == "x": - transform = as_tensor([[Lx, 0.], [0., Ly]]) - else: - transform = as_tensor([[0., Lx], [Ly, 0]]) - new_coordinates.interpolate(dot( - transform, - as_vector(( - conditional(gt(real(x), real(1. - eps)), indicator, # Periodic break. - # Unwrap rest of circle. - atan2(real(-y), real(-x))/(2 * pi) + 0.5), - z - )) - )) - return _postprocess_periodic_mesh(new_coordinates, - comm, - distribution_parameters, - reorder, - name, - distribution_name, - permutation_name) + +def _mark_mesh_boundaries(plex: PETSc.DMPlex) -> None: + """Reorder the 'Face Sets' label of the DMPlex to match Firedrake ordering.""" + match plex.getDimension(): + case 0: + plex_to_firedrake_boundary_labels = {} + case 1: + # Firedrake and DMPlex conventions agree (left is 1, right is 2) + plex_to_firedrake_boundary_labels = {1: 1, 2: 2} + case 2: + # DMPlex Firedrake + # + # 3 4 + # +-----+ +-----+ + # | | | | + # 4| |2 1| |2 + # | | | | + # +-----+ +-----+ + # 1 3 + plex_to_firedrake_boundary_labels = {1: 3, 2: 2, 3: 4, 4: 1} + case 3: + # DMPlex Firedrake + # + # +-------+ +-------+ + # / /| / /| + # / / | / / | + # / 4/3 / | / 4/3 / | + # / / | / / | + # / / | / / | + # +-------+ 5/6 + +-------+ 2/1 + + # | | / | | / + # | | / | | / + # y | 1/2 | / z y | 5/6 | / z + # | | / | | / + # | |/ | |/ + # +-------+ +-------+ + # O x O x + # + # 'x/y' means that 'x' is the label for the visible face and 'y' + # the label for the opposite hidden face. + plex_to_firedrake_boundary_labels = {1: 5, 2: 6, 3: 3, 4: 4, 5: 2, 6: 1} + case _: + raise AssertionError + + # Get the original label + plex_boundary_label = plex.getLabel(dmcommon.FACE_SETS_LABEL) + plex.removeLabel(dmcommon.FACE_SETS_LABEL) + + # Now create the new one + plex.createLabel(dmcommon.FACE_SETS_LABEL) + firedrake_boundary_label = plex.getLabel(dmcommon.FACE_SETS_LABEL) + for plex_value, firedrake_value in plex_to_firedrake_boundary_labels.items(): + points = plex_boundary_label.getStratumIS(plex_value) + if points: + firedrake_boundary_label.setStratumIS(firedrake_value, points) + else: + # create an empty stratum + firedrake_boundary_label.addStratum(firedrake_value) + + +def _refine_quads_to_triangles( + plex: PETSc.DMPlex, + diagonal: Literal["crossed", "left", "right"], +) -> PETSc.DMPlex: + match diagonal: + case "crossed": + options = {"dm_refine": 1, "dm_plex_transform_type": "refine_alfeld"} + case "left": + options = { + "dm_refine": 1, + "dm_plex_transform_type": "refine_tosimplex", + "dm_plex_transform_tosimplex_reflect": True, + } + case "right": + options = {"dm_refine": 1, "dm_plex_transform_type": "refine_tosimplex"} + case _: + raise AssertionError(f"'diagonal' type '{diagonal}' is not recognised") + + tr = PETSc.DMPlexTransform().create(comm=plex.comm) + petsctools.set_from_options(tr, options) + tr.setDM(plex) + tr.setUp() + with petsctools.inserted_options(tr): + refined_plex = tr.apply(plex) + tr.destroy() + return refined_plex diff --git a/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py b/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py index 402148edfd..cc2d67fe65 100644 --- a/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py +++ b/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py @@ -14,7 +14,7 @@ def test_equation_bcs_direct_assemble_one_form(): bc = EquationBC(F1 == 0, u, 1) g = assemble(F, bcs=bc.extract_form('F')) - assert np.allclose(g.dat.data, [0.5, 0.5, 0, 0]) + assert np.allclose(g.dat.data, [0.5, 0, 0, 0.5]) def test_equation_bcs_direct_assemble_two_form(): @@ -31,15 +31,19 @@ def test_equation_bcs_direct_assemble_two_form(): # Must preprocess bc to extract appropriate # `EquationBCSplit` object. A = assemble(a, bcs=bc.extract_form('J')) - assert np.allclose(A.M.values, [[1 / 3, 1 / 6, 0, 0], - [1 / 6, 1 / 3, 0, 0], - [-1 / 3, -1 / 6, 2 / 3, -1 / 6], - [-1 / 6, -1 / 3, -1 / 6, 2 / 3]]) + expected = [[ 1/3, 0, 0, 1/6], + [-1/6, 2/3, -1/6, -1/3], + [-1/3, -1/6, 2/3, -1/6], + [ 1/6, 0, 0, 1/3]] + assert np.allclose(A.M.values, expected) + A = assemble(a, bcs=bc.extract_form('Jp')) - assert np.allclose(A.M.values, [[2 / 3, 2 / 6, 0, 0], - [2 / 6, 2 / 3, 0, 0], - [-1 / 3, -1 / 6, 2 / 3, -1 / 6], - [-1 / 6, -1 / 3, -1 / 6, 2 / 3]]) + expected = [[ 2/3, 0, 0, 1/3], + [-1/6, 2/3, -1/6, -1/3], + [-1/3, -1/6, 2/3, -1/6], + [ 1/3, 0, 0, 2/3]] + assert np.allclose(A.M.values, expected) + with pytest.raises(TypeError) as excinfo: # Unable to use raw `EquationBC` object, as # assembler can not infer merely from the rank diff --git a/tests/firedrake/extrusion/test_variable_layers_numbering.py b/tests/firedrake/extrusion/test_variable_layers_numbering.py index e13358f57d..9eadb6d5b2 100644 --- a/tests/firedrake/extrusion/test_variable_layers_numbering.py +++ b/tests/firedrake/extrusion/test_variable_layers_numbering.py @@ -381,16 +381,16 @@ def test_numbering_quad(): [2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 18, 20, 21, 24]).all() assert numpy.equal(DirichletBC(V, 0, 1).nodes, - [0, 1, 2, 3, 4, 5, 17, 18]).all() + [0, 1, 2, 9, 10, 11, 15, 16]).all() assert numpy.equal(DirichletBC(V, 0, 2).nodes, - [12, 13, 14, 15, 16, 22, 23, 24]).all() + [17, 18, 19, 20, 21, 22, 23, 24]).all() assert numpy.equal(DirichletBC(V, 0, 3).nodes, - [0, 1, 2, 9, 10, 11, 15, 16]).all() + [0, 1, 2, 3, 4, 5, 17, 18]).all() assert numpy.equal(DirichletBC(V, 0, 4).nodes, - [17, 18, 19, 20, 21, 22, 23, 24]).all() + [12, 13, 14, 15, 16, 22, 23, 24]).all() @pytest.mark.parametrize(["domain", "expected"], diff --git a/tests/firedrake/multigrid/test_grid_transfer.py b/tests/firedrake/multigrid/test_grid_transfer.py index 86af28b455..53a848fb59 100644 --- a/tests/firedrake/multigrid/test_grid_transfer.py +++ b/tests/firedrake/multigrid/test_grid_transfer.py @@ -329,7 +329,7 @@ def exact_primal_periodic(mesh, shape, degree): return expr -@pytest.mark.parallel(nprocs=3) +@pytest.mark.parallel def test_grid_transfer_periodic(periodic_hierarchy, periodic_space): degrees = [4] shape = "scalar" diff --git a/tests/firedrake/regression/test_bcs.py b/tests/firedrake/regression/test_bcs.py index 4544b95cc0..9e43ba805b 100644 --- a/tests/firedrake/regression/test_bcs.py +++ b/tests/firedrake/regression/test_bcs.py @@ -423,8 +423,8 @@ def test_bcs_mixed_real(): bc = DirichletBC(V.sub(0), 0.0, 1) a = inner(u1, v0) * dx + inner(u0, v1) * dx A = assemble(a, bcs=[bc, ]) - assert np.allclose(A.M[0][1].values, [[0.00], [0.00], [0.25], [0.25]]) - assert np.allclose(A.M[1][0].values, [[0.00, 0.00, 0.25, 0.25]]) + assert np.allclose(A.M[0][1].values, [[0.00], [0.25], [0.25], [0.00]]) + assert np.allclose(A.M[1][0].values, [[0.00, 0.25, 0.25, 0.00]]) def test_bcs_mixed_real_vector(): @@ -437,8 +437,8 @@ def test_bcs_mixed_real_vector(): bc = DirichletBC(V.sub(0).sub(1), 0.0, 1) a = inner(as_vector([u1, u1]), v0) * dx + inner(u0, as_vector([v1, v1])) * dx A = assemble(a, bcs=[bc, ]) - assert np.allclose(A.M[0][1].values, [[[0.25], [0.], [0.25], [0.], [0.25], [0.25], [0.25], [0.25]]]) - assert np.allclose(A.M[1][0].values, [[0.25, 0., 0.25, 0., 0.25, 0.25, 0.25, 0.25]]) + assert np.allclose(A.M[0][1].values, [[[0.25], [0.], [0.25], [0.25], [0.25], [0.25], [0.25], [0.]]]) + assert np.allclose(A.M[1][0].values, [[0.25, 0., 0.25, 0.25, 0.25, 0.25, 0.25, 0.]]) def test_homogeneous_bc_residual(): diff --git a/tests/firedrake/regression/test_mark_entities.py b/tests/firedrake/regression/test_mark_entities.py index a9b824fe35..c525f6be59 100644 --- a/tests/firedrake/regression/test_mark_entities.py +++ b/tests/firedrake/regression/test_mark_entities.py @@ -103,7 +103,7 @@ def test_mark_entities_mesh_mark_entities_2d(): plex = mesh.topology.topology_dm label = plex.getLabel(label_name) assert label.getStratumIS(label_value).getSize() == 2 - assert all(label.getStratumIS(label_value).getIndices() == [20, 30]) + assert all(label.getStratumIS(label_value).getIndices() == [20, 23]) def test_mark_entities_mesh_mark_entities_1d(): diff --git a/tests/firedrake/regression/test_mesh_generation.py b/tests/firedrake/regression/test_mesh_generation.py index 4cb7cc109a..824cd09277 100644 --- a/tests/firedrake/regression/test_mesh_generation.py +++ b/tests/firedrake/regression/test_mesh_generation.py @@ -74,7 +74,8 @@ def test_tensor_box(): assert abs(integrate_one(TensorBoxMesh(xcoords, ycoords, zcoords)) - 0.6) < 1e-3 -def run_one_element_advection(): +@pytest.mark.parallel([1, 2]) +def test_one_element_advection(): nx = 20 m = PeriodicRectangleMesh(nx, 1, 1.0, 1.0, quadrilateral=True) nlayers = 20 @@ -122,16 +123,8 @@ def run_one_element_advection(): assert assemble(inner(q0-q_init, q0-q_init)*dx)**0.5 < 0.005 -def test_one_element_advection(): - run_one_element_advection() - - -@pytest.mark.parallel(nprocs=2) -def test_one_element_advection_parallel(): - run_one_element_advection() - - -def run_one_element_mesh(): +@pytest.mark.parallel([1, 3]) +def test_one_element_mesh(): mesh = PeriodicRectangleMesh(20, 1, Lx=1.0, Ly=1.0, quadrilateral=True) x = SpatialCoordinate(mesh) V = FunctionSpace(mesh, "CG", 1) @@ -160,19 +153,12 @@ def run_one_element_mesh(): assert err > 1.0e-3 -def test_one_element_mesh(): - run_one_element_mesh() - - -@pytest.mark.parallel(nprocs=3) -def test_one_element_mesh_parallel(): - run_one_element_mesh() - - -def test_box(): - assert abs(integrate_one(BoxMesh(3, 3, 3, 1, 2, 3)) - 6) < 1e-3 +@pytest.mark.parametrize("hexahedral", [False, True]) +def test_box(hexahedral): + assert abs(integrate_one(BoxMesh(3, 3, 3, 1, 2, 3, hexahedral=hexahedral)) - 6) < 1e-3 +@pytest.mark.parallel([1, 3]) def test_periodic_unit_cube(): assert abs(integrate_one(PeriodicUnitCubeMesh(3, 3, 3)) - 1) < 1e-3 @@ -239,11 +225,6 @@ def test_tensor_box_parallel(): assert abs(integrate_one(TensorBoxMesh(xcoords, ycoords, zcoords)) - 0.6) < 1e-3 -@pytest.mark.parallel -def test_periodic_unit_cube_parallel(): - assert abs(integrate_one(PeriodicUnitCubeMesh(3, 3, 3)) - 1) < 1e-3 - - def assert_num_exterior_facets_equals_zero(m): # Need to initialise the mesh so that exterior facets have been # built. @@ -471,7 +452,7 @@ def test_boxmesh_kind(kind, num_cells): assert m.num_cells() == num_cells -@pytest.mark.parallel(nprocs=2) +@pytest.mark.parallel(2) def test_periodic_unit_cube_hex_cell(): mesh = PeriodicUnitCubeMesh(3, 3, 3, directions=[True, True, False], hexahedral=True) x, y, z = SpatialCoordinate(mesh) @@ -482,12 +463,9 @@ def test_periodic_unit_cube_hex_cell(): assert error < 1.e-30 -@pytest.mark.parallel(nprocs=4) +@pytest.mark.parallel(4) def test_periodic_unit_cube_hex_facet(): mesh = PeriodicUnitCubeMesh(3, 3, 3, directions=[True, False, False], hexahedral=True) - for subdomain_id in [1, 2]: - area = assemble(Constant(1.) * dS(domain=mesh, subdomain_id=subdomain_id)) - assert abs(area - 1.0) < 1.e-15 for subdomain_id in [3, 4, 5, 6]: area = assemble(Constant(1.) * ds(domain=mesh, subdomain_id=subdomain_id)) assert abs(area - 1.0) < 1.e-15 diff --git a/tests/firedrake/regression/test_periodic_2d.py b/tests/firedrake/regression/test_periodic_2d.py index 8e403479da..b3568a6be5 100644 --- a/tests/firedrake/regression/test_periodic_2d.py +++ b/tests/firedrake/regression/test_periodic_2d.py @@ -20,33 +20,15 @@ from firedrake import * -@pytest.fixture(params=["x", "y", "both"]) -def direction(request): - return request.param - - -@pytest.fixture(params=[False, True], - ids=["tri", "quad"]) -def quadrilateral(request): - return request.param - - -@pytest.fixture(params=["left", "right", "crossed"]) -def diagonal(request): - return request.param - - -def run_periodic_helmholtz(direction, quadrilateral, diagonal): - if quadrilateral: - if diagonal == "left": - # run the test - diagonal = None - else: - # don't run the test - return - - mesh = PeriodicRectangleMesh(100, 60, 5, 3, quadrilateral=quadrilateral, - diagonal=diagonal, direction=direction) +@pytest.mark.parallel([1, 3]) +@pytest.mark.parametrize("direction", ["x", "y", "both"]) +@pytest.mark.parametrize("cell_options", + [{"quadrilateral": True}, + {"quadrilateral": False, "diagonal": "left"}, + {"quadrilateral": False, "diagonal": "right"}, + {"quadrilateral": False, "diagonal": "crossed"}]) +def test_periodic_helmholtz(direction, cell_options): + mesh = PeriodicRectangleMesh(100, 60, 5, 3, **cell_options, direction=direction) x = SpatialCoordinate(mesh) V = FunctionSpace(mesh, "CG", 1) @@ -56,7 +38,10 @@ def run_periodic_helmholtz(direction, quadrilateral, diagonal): f = Function(V).assign((244.0*pi*pi/225.0 + 1.0)*u_exact) - if direction in ("x", "y"): + # FIXME: This suggests some really weird behaviour! + if direction == "x": + bcs = DirichletBC(V, Constant(0), (3, 4)) + if direction == "y": bcs = DirichletBC(V, Constant(0), (1, 2)) elif direction == "both": bcs = [] @@ -72,12 +57,3 @@ def run_periodic_helmholtz(direction, quadrilateral, diagonal): l2err = sqrt(assemble(inner((out-u_exact), (out-u_exact))*dx)) l2norm = sqrt(assemble(inner(u_exact, u_exact)*dx)) assert l2err/l2norm < 0.004 - - -def test_periodic_helmholtz(direction, quadrilateral, diagonal): - run_periodic_helmholtz(direction, quadrilateral, diagonal) - - -@pytest.mark.parallel(nprocs=3) -def test_periodic_helmholtz_parallel(direction, quadrilateral, diagonal): - run_periodic_helmholtz(direction, quadrilateral, diagonal) diff --git a/tests/firedrake/submesh/test_submesh_assemble.py b/tests/firedrake/submesh/test_submesh_assemble.py index fabf9ca411..05d243e8b5 100644 --- a/tests/firedrake/submesh/test_submesh_assemble.py +++ b/tests/firedrake/submesh/test_submesh_assemble.py @@ -63,12 +63,14 @@ def test_submesh_assemble_cell_cell_integral_facet(): assert np.allclose(A.M.sparsity[0][1].nnz, [4, 4, 4, 4, 4, 4, 4, 4]) assert np.allclose(A.M.sparsity[1][0].nnz, [8, 8, 8, 8]) assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M10 = np.array([[0., 0., 0., 0., 0., 0., 1. / 3., 1. / 6.], - [0., 0., 0., 0., 0., 0., 1. / 6., 1. / 3.], - [0., 0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0., 0.]]) + + M10 = [[0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1/3, 0, 1/6, 0], + [0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1/6, 0, 1/3, 0]] assert np.allclose(A.M[0][1].values, np.transpose(M10)) assert np.allclose(A.M[1][0].values, M10) + b = inner(u1, v0('+')) * ds1(5) + inner(u0('+'), v1) * dS0 B = assemble(b, mat_type="nest") assert np.allclose(B.M.sparsity[0][0].nnz, [1, 1, 1, 1, 1, 1, 1, 1]) # bc nodes @@ -167,12 +169,14 @@ def test_submesh_assemble_cell_cell_cell_cell_integral_various(): assert np.allclose(A.M.sparsity[0][1].nnz, [4, 4, 4, 4, 0, 0, 0, 0]) assert np.allclose(A.M.sparsity[1][0].nnz, [4, 4, 4, 4]) assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M10 = np.array([[0., 0., 1. / 3., 1. / 6., 0., 0., 0., 0.], - [0., 0., 1. / 6., 1. / 3., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0., 0., 0.]]) + + M10 = [[ 0, 0, 0, 0, 0, 0, 0, 0], + [1/3, 0, 1/6, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0], + [1/6, 0, 1/3, 0, 0, 0, 0, 0]] assert np.allclose(A.M[0][1].values, np.transpose(M10)) assert np.allclose(A.M[1][0].values, M10) + b = inner(u_rl, v_l) * dS(label_int) + inner(u_l, v_rl) * dS(label_int) B = assemble(b, mat_type="nest") assert np.allclose(B.M.sparsity[0][0].nnz, [1, 1, 1, 1, 1, 1, 1, 1]) # bc nodes @@ -315,22 +319,24 @@ def test_submesh_assemble_cell_cell_equation_bc(): sol = Function(V) bc = EquationBC(a_int == L_int, sol, label_int, V=V.sub(0)) A = assemble(a, bcs=bc.extract_form('J'), mat_type="nest") - assert np.allclose(Function(V_l).interpolate(SpatialCoordinate(mesh_l)[0]).dat.data, [0., 0., 1., 1.]) - assert np.allclose(Function(V_l).interpolate(SpatialCoordinate(mesh_l)[1]).dat.data, [0., 1., 1., 0.]) - assert np.allclose(Function(V_r).interpolate(SpatialCoordinate(mesh_r)[0]).dat.data, [1., 1., 2., 2.]) - assert np.allclose(Function(V_r).interpolate(SpatialCoordinate(mesh_r)[1]).dat.data, [0., 1., 1., 0.]) + assert np.allclose(Function(V_l).interpolate(SpatialCoordinate(mesh_l)[0]).dat.data, [0., 1., 1., 0.]) + assert np.allclose(Function(V_l).interpolate(SpatialCoordinate(mesh_l)[1]).dat.data, [0., 0., 1., 1.]) + assert np.allclose(Function(V_r).interpolate(SpatialCoordinate(mesh_r)[0]).dat.data, [1., 2., 2., 1.]) + assert np.allclose(Function(V_r).interpolate(SpatialCoordinate(mesh_r)[1]).dat.data, [0., 0., 1., 1.]) assert np.allclose(A.M.sparsity[0][0].nnz, [4, 4, 4, 4]) assert np.allclose(A.M.sparsity[0][1].nnz, [4, 4, 4, 4]) assert np.allclose(A.M.sparsity[1][0].nnz, [0, 0, 0, 0]) assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M00 = np.array([[1. / 9. , 1. / 18., 1. / 36., 1. / 18.], # noqa: E203 - [1. / 18., 1. / 9. , 1. / 18., 1. / 36.], # noqa: E203 - [0., 0., 1. / 3., 1. / 6.], - [0., 0., 1. / 6., 1. / 3.]]) - M01 = np.array([[0., 0., 0., 0.], - [0., 0., 0., 0.], - [- 1. / 6., - 1. / 3., 0., 0.], - [- 1. / 3., - 1. / 6., 0., 0.]]) + + + M00 = np.array([[ 1/9, 1/18, 1/36, 1/18], + [ 0, 1/3, 1/6, 0], + [ 0, 1/6, 1/3, 0], + [1/18, 1/36, 1/18, 1/9]]) + M01 = np.array([[ 0, 0, 0, 0], + [-1/3, 0, 0, -1/6], + [-1/6, 0, 0, -1/3], + [ 0, 0, 0, 0]]) assert np.allclose(A.M[0][0].values, M00) assert np.allclose(A.M[0][1].values, M01) @@ -376,24 +382,13 @@ def test_submesh_assemble_cell_facet_integral_various(): coords.sub(0).assign(mesh.coordinates) coords.sub(1).assign(subm.coordinates) coords0, coords1 = split(coords) - M10 = np.array( - [ - [1. / 6., 1. / 3., 0., 0., 0., 0.], - [1. / 3., 1. / 6., 0., 0., 0., 0.], - ] - ) - M10w = np.array( - [ - [1. / 12., 1. / 4., 0., 0., 0., 0.], - [1. / 12., 1. / 12., 0., 0., 0., 0.], - ] - ) - M10ww = np.array( - [ - [1. / 20., 1. / 5., 0., 0., 0., 0.], - [1. / 30., 1. / 20., 0., 0., 0., 0.], - ] - ) + + M10 = np.array([[1/3, 0, 0, 1/6, 0, 0], + [1/6, 0, 0, 1/3, 0, 0]]) + M10w = np.array([[1/12, 0, 0, 1/12, 0, 0], + [1/12, 0, 0, 1/4, 0, 0]]) + M10ww = np.array([[1/30, 0, 0, 1/20, 0, 0], + [1/20, 0, 0, 1/5, 0, 0]]) # Use subm as primal integration domain. measure = Measure( "dx", subm, @@ -401,24 +396,31 @@ def test_submesh_assemble_cell_facet_integral_various(): Measure("dS", mesh), ), ) + a = inner(u0('-'), v1) * measure A = assemble(a, mat_type="nest") assert np.allclose(A.M[1][0].values, M10) + a = inner(u1, v0('+')) * measure A = assemble(a, mat_type="nest") assert np.allclose(A.M[0][1].values, np.transpose(M10)) + a = y * inner(u0('-'), v1) * measure A = assemble(a, mat_type="nest") assert np.allclose(A.M[1][0].values, M10w) + a = y * suby * inner(u0('-'), v1) * measure A = assemble(a, mat_type="nest") assert np.allclose(A.M[1][0].values, M10ww) + a = coords0[1] * inner(u0('-'), v1) * measure A = assemble(a, mat_type="nest") assert np.allclose(A.M[1][0].values, M10w) + a = coords0[1] * coords1[1] * inner(u0('-'), v1) * measure A = assemble(a, mat_type="nest") assert np.allclose(A.M[1][0].values, M10ww) + # Use mesh as primal integration domain. measure = Measure( "dS", mesh, diff --git a/tests/firedrake/submesh/test_submesh_base.py b/tests/firedrake/submesh/test_submesh_base.py index e2c93fdad3..05ccd3de91 100644 --- a/tests/firedrake/submesh/test_submesh_base.py +++ b/tests/firedrake/submesh/test_submesh_base.py @@ -249,7 +249,7 @@ def test_submesh_base_entity_maps(): assert (mesh.interior_facets.facets == np.array([8])).all assert (mesh.exterior_facets.facets == np.array([9, 10, 11, 12, 13, 14])).all assert (submesh.interior_facets.facets == np.array([])).all - assert (submesh.exterior_facets.facets == np.array([6, 7, 5, 8])).all() + assert (submesh.exterior_facets.facets == np.array([6, 5, 7, 8])).all() composed_map, integral_type = mesh.topology.trans_mesh_entity_map(submesh.topology, "cell", "everywhere", None) assert integral_type == "cell" if rank == 0: From 6a4fd9a80ca20d560afd6e8b2646c64b7fff11c7 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Fri, 6 Mar 2026 15:46:26 +0000 Subject: [PATCH 02/10] linting --- .../test_equation_bcs_assemble.py | 16 ++++----- .../submesh/test_submesh_assemble.py | 36 +++++++++---------- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py b/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py index cc2d67fe65..e523945651 100644 --- a/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py +++ b/tests/firedrake/equation_bcs/test_equation_bcs_assemble.py @@ -31,17 +31,17 @@ def test_equation_bcs_direct_assemble_two_form(): # Must preprocess bc to extract appropriate # `EquationBCSplit` object. A = assemble(a, bcs=bc.extract_form('J')) - expected = [[ 1/3, 0, 0, 1/6], - [-1/6, 2/3, -1/6, -1/3], - [-1/3, -1/6, 2/3, -1/6], - [ 1/6, 0, 0, 1/3]] + expected = [[ 1/3, 0, 0, 1/6], # noqa + [-1/6, 2/3, -1/6, -1/3], # noqa + [-1/3, -1/6, 2/3, -1/6], # noqa + [ 1/6, 0, 0, 1/3]] # noqa assert np.allclose(A.M.values, expected) A = assemble(a, bcs=bc.extract_form('Jp')) - expected = [[ 2/3, 0, 0, 1/3], - [-1/6, 2/3, -1/6, -1/3], - [-1/3, -1/6, 2/3, -1/6], - [ 1/3, 0, 0, 2/3]] + expected = [[ 2/3, 0, 0, 1/3], # noqa + [-1/6, 2/3, -1/6, -1/3], # noqa + [-1/3, -1/6, 2/3, -1/6], # noqa + [ 1/3, 0, 0, 2/3]] # noqa assert np.allclose(A.M.values, expected) with pytest.raises(TypeError) as excinfo: diff --git a/tests/firedrake/submesh/test_submesh_assemble.py b/tests/firedrake/submesh/test_submesh_assemble.py index 05d243e8b5..af317a8ec9 100644 --- a/tests/firedrake/submesh/test_submesh_assemble.py +++ b/tests/firedrake/submesh/test_submesh_assemble.py @@ -32,10 +32,10 @@ def test_submesh_assemble_cell_cell_integral_cell(): assert np.allclose(A.M.sparsity[0][1].nnz, [4, 4, 4, 4, 0, 0]) assert np.allclose(A.M.sparsity[1][0].nnz, [4, 4, 4, 4]) assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M10 = np.array([[1./9. , 1./18., 1./36., 1./18., 0., 0.], # noqa: E203 - [1./18., 1./9. , 1./18., 1./36., 0., 0.], # noqa: E203 - [1./36., 1./18., 1./9. , 1./18., 0., 0.], # noqa: E203 - [1./18., 1./36., 1./18., 1./9. , 0., 0.]]) # noqa: E203 + M10 = np.array([[1./9. , 1./18., 1./36., 1./18., 0., 0.], # noqa + [1./18., 1./9. , 1./18., 1./36., 0., 0.], # noqa + [1./36., 1./18., 1./9. , 1./18., 0., 0.], # noqa + [1./18., 1./36., 1./18., 1./9. , 0., 0.]]) # noqa assert np.allclose(A.M[0][1].values, np.transpose(M10)) assert np.allclose(A.M[1][0].values, M10) @@ -64,9 +64,9 @@ def test_submesh_assemble_cell_cell_integral_facet(): assert np.allclose(A.M.sparsity[1][0].nnz, [8, 8, 8, 8]) assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M10 = [[0, 0, 0, 0, 0, 0, 0, 0], + M10 = [[0, 0, 0, 0, 0, 0, 0, 0], # noqa [0, 0, 0, 0, 1/3, 0, 1/6, 0], - [0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0], # noqa [0, 0, 0, 0, 1/6, 0, 1/3, 0]] assert np.allclose(A.M[0][1].values, np.transpose(M10)) assert np.allclose(A.M[1][0].values, M10) @@ -170,9 +170,9 @@ def test_submesh_assemble_cell_cell_cell_cell_integral_various(): assert np.allclose(A.M.sparsity[1][0].nnz, [4, 4, 4, 4]) assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M10 = [[ 0, 0, 0, 0, 0, 0, 0, 0], + M10 = [[ 0, 0, 0, 0, 0, 0, 0, 0], # noqa [1/3, 0, 1/6, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0], # noqa [1/6, 0, 1/3, 0, 0, 0, 0, 0]] assert np.allclose(A.M[0][1].values, np.transpose(M10)) assert np.allclose(A.M[1][0].values, M10) @@ -329,14 +329,14 @@ def test_submesh_assemble_cell_cell_equation_bc(): assert np.allclose(A.M.sparsity[1][1].nnz, [1, 1, 1, 1]) # bc nodes - M00 = np.array([[ 1/9, 1/18, 1/36, 1/18], - [ 0, 1/3, 1/6, 0], - [ 0, 1/6, 1/3, 0], - [1/18, 1/36, 1/18, 1/9]]) - M01 = np.array([[ 0, 0, 0, 0], - [-1/3, 0, 0, -1/6], - [-1/6, 0, 0, -1/3], - [ 0, 0, 0, 0]]) + M00 = np.array([[ 1/9, 1/18, 1/36, 1/18], # noqa + [ 0, 1/3, 1/6, 0], # noqa + [ 0, 1/6, 1/3, 0], # noqa + [1/18, 1/36, 1/18, 1/9]]) # noqa + M01 = np.array([[ 0, 0, 0, 0], # noqa + [-1/3, 0, 0, -1/6], # noqa + [-1/6, 0, 0, -1/3], # noqa + [ 0, 0, 0, 0]]) # noqa assert np.allclose(A.M[0][0].values, M00) assert np.allclose(A.M[0][1].values, M01) @@ -386,9 +386,9 @@ def test_submesh_assemble_cell_facet_integral_various(): M10 = np.array([[1/3, 0, 0, 1/6, 0, 0], [1/6, 0, 0, 1/3, 0, 0]]) M10w = np.array([[1/12, 0, 0, 1/12, 0, 0], - [1/12, 0, 0, 1/4, 0, 0]]) + [1/12, 0, 0, 1/4, 0, 0]]) # noqa M10ww = np.array([[1/30, 0, 0, 1/20, 0, 0], - [1/20, 0, 0, 1/5, 0, 0]]) + [1/20, 0, 0, 1/5, 0, 0]]) # noqa # Use subm as primal integration domain. measure = Measure( "dx", subm, From f4af99958180f760983c5537a2f069cde6b0ff09 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Fri, 6 Mar 2026 15:48:41 +0000 Subject: [PATCH 03/10] vom fix --- tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py b/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py index 2ca7acda3f..4361733dcd 100644 --- a/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py +++ b/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py @@ -354,7 +354,7 @@ def test_point_tolerance(): m = UnitSquareMesh(1, 1) assert m.tolerance == 0.5 # Make the mesh non-axis-aligned. - m.coordinates.dat.data[1, :] = [1.1, 1] + m.coordinates.dat.data[2, :] = [1.1, 1] coords = [[1.0501, 0.5]] vm = VertexOnlyMesh(m, coords, tolerance=0.1) assert vm.cell_set.size == 1 From 8bf3ca2582ad9db3ccb196804ef3a0eae5f859e5 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Fri, 6 Mar 2026 15:53:52 +0000 Subject: [PATCH 04/10] submesh fix --- tests/firedrake/submesh/test_submesh_base.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/firedrake/submesh/test_submesh_base.py b/tests/firedrake/submesh/test_submesh_base.py index 05ccd3de91..58d379bf5a 100644 --- a/tests/firedrake/submesh/test_submesh_base.py +++ b/tests/firedrake/submesh/test_submesh_base.py @@ -240,7 +240,7 @@ def test_submesh_base_entity_maps(): assert (mesh.interior_facets.facets == np.array([11])).all assert (mesh.exterior_facets.facets == np.array([8, 9, 10, 12, 13, 14])).all assert (submesh.interior_facets.facets == np.array([])).all - assert (submesh.exterior_facets.facets == np.array([5, 6, 8, 7])).all() + assert (submesh.exterior_facets.facets == np.array([5, 8, 6, 7])).all() else: assert subdm.getLabel("pyop2_core").getStratumSize(1) == 0 assert subdm.getLabel("pyop2_owned").getStratumSize(1) == 0 @@ -259,18 +259,18 @@ def test_submesh_base_entity_maps(): composed_map, integral_type = mesh.topology.trans_mesh_entity_map(submesh.topology, "exterior_facet", 5, None) assert integral_type == "interior_facet" if rank == 0: - assert (composed_map.maps_[0].values_with_halo == np.array([-1, -1, 0, -1]).reshape((-1, 1))).all() # entire exterior-interior map + assert (composed_map.maps_[0].values_with_halo == np.array([-1, 0, -1, -1]).reshape((-1, 1))).all() # entire exterior-interior map else: - assert (composed_map.maps_[0].values_with_halo == np.array([-1, -1, 0, -1]).reshape((-1, 1))).all() # entire exterior-interior map + assert (composed_map.maps_[0].values_with_halo == np.array([-1, 0, -1, -1]).reshape((-1, 1))).all() # entire exterior-interior map composed_map, integral_type = mesh.topology.trans_mesh_entity_map(submesh.topology, "exterior_facet", 4, None) assert integral_type == "exterior_facet" if rank == 0: - assert (composed_map.maps_[0].values_with_halo == np.array([0, 1, -1, 2]).reshape((-1, 1))).all() # entire exterior-exterior map + assert (composed_map.maps_[0].values_with_halo == np.array([0, -1, 1, 2]).reshape((-1, 1))).all() # entire exterior-exterior map else: - assert (composed_map.maps_[0].values_with_halo == np.array([3, 4, -1, 5]).reshape((-1, 1))).all() # entire exterior-exterior map + assert (composed_map.maps_[0].values_with_halo == np.array([3, -1, 4, 5]).reshape((-1, 1))).all() # entire exterior-exterior map composed_map, integral_type = submesh.topology.trans_mesh_entity_map(mesh.topology, "exterior_facet", 1, None) assert integral_type == "exterior_facet" if rank == 0: - assert (composed_map.maps_[0].values_with_halo == np.array([0, 1, 3, -1, -1, -1]).reshape((-1, 1))).all() + assert (composed_map.maps_[0].values_with_halo == np.array([0, 2, 3, -1, -1, -1]).reshape((-1, 1))).all() else: - assert (composed_map.maps_[0].values_with_halo == np.array([-1, -1, -1, 0, 1, 3]).reshape((-1, 1))).all() + assert (composed_map.maps_[0].values_with_halo == np.array([-1, -1, -1, 0, 2, 3]).reshape((-1, 1))).all() From 077d1d4399ca754f83c942d1aa6eaf7918ea7e2a Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Fri, 6 Mar 2026 16:40:08 +0000 Subject: [PATCH 05/10] More fixes, all of them? --- firedrake/cython/dmcommon.pyx | 42 +++++++++---------- firedrake/cython/petschdr.pxi | 2 + .../test_variable_layers_numbering.py | 6 +-- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ea76b9673a..1cfe46188d 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -2097,11 +2097,9 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): cStart, cEnd = dm.getHeightStratum(0) dim = dm.getCoordinateDim() - - dm_coords = dm.getCellCoordinatesLocal().array_r coords_shape = ((cEnd-cStart) * ndofs[0], dim) - if dm_coords.size < np.prod(coords_shape, dtype=IntType): + if dm.getCellCoordinateSection().getDof(cStart) < ndofs[0] * dim: # Fewer cell coordinates available, we must be single-cell periodic if dm.getCellType(cStart) == PETSc.DM.PolytopeType.QUADRILATERAL: # If we have a periodic mesh with only a single cell in the periodic @@ -2129,8 +2127,8 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): # | | # 1-----2 assert ndofs[0] == 4, "Not expecting high order coords here" - dm_coords = dm_coords.reshape(((cEnd-cStart) * 2, dim)) - dm_coords_expanded = np.empty(coords_shape, dtype=dm_coords.dtype) + dm_coords_orig = dm.getCellCoordinatesLocal().array_r.reshape(((cEnd-cStart) * 2, dim)) + dm_coords_expanded = np.empty(coords_shape, dtype=dm_coords_orig.dtype) # Create a new cell coordinate section dm_sec_orig = dm.getCellCoordinateSection() @@ -2155,14 +2153,14 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): for c in range(cStart, cEnd): CHKERR(PetscSectionSetDof(dm_sec_expanded.sec, c, 8)) - dm_coords_expanded[4*c+0, 0] = dm_coords[2*c+0, 0] - dm_coords_expanded[4*c+1, 0] = dm_coords[2*c+0, 0] + cell_width - dm_coords_expanded[4*c+2, 0] = dm_coords[2*c+1, 0] + cell_width - dm_coords_expanded[4*c+3, 0] = dm_coords[2*c+1, 0] - dm_coords_expanded[4*c+0, 1] = dm_coords[2*c+0, 1] - dm_coords_expanded[4*c+1, 1] = dm_coords[2*c+0, 1] - dm_coords_expanded[4*c+2, 1] = dm_coords[2*c+1, 1] - dm_coords_expanded[4*c+3, 1] = dm_coords[2*c+1, 1] + dm_coords_expanded[4*c+0, 0] = dm_coords_orig[2*c+0, 0] + dm_coords_expanded[4*c+1, 0] = dm_coords_orig[2*c+0, 0] + cell_width + dm_coords_expanded[4*c+2, 0] = dm_coords_orig[2*c+1, 0] + cell_width + dm_coords_expanded[4*c+3, 0] = dm_coords_orig[2*c+1, 0] + dm_coords_expanded[4*c+0, 1] = dm_coords_orig[2*c+0, 1] + dm_coords_expanded[4*c+1, 1] = dm_coords_orig[2*c+0, 1] + dm_coords_expanded[4*c+2, 1] = dm_coords_orig[2*c+1, 1] + dm_coords_expanded[4*c+3, 1] = dm_coords_orig[2*c+1, 1] else: if vert_unit_periodic: @@ -2171,14 +2169,14 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): for c in range(cStart, cEnd): CHKERR(PetscSectionSetDof(dm_sec_expanded.sec, c, 8)) - dm_coords_expanded[4*c+0, 0] = dm_coords[2*c+0, 0] - dm_coords_expanded[4*c+1, 0] = dm_coords[2*c+1, 0] - dm_coords_expanded[4*c+2, 0] = dm_coords[2*c+1, 0] - dm_coords_expanded[4*c+3, 0] = dm_coords[2*c+0, 0] - dm_coords_expanded[4*c+0, 1] = dm_coords[2*c+0, 1] - dm_coords_expanded[4*c+1, 1] = dm_coords[2*c+1, 1] - dm_coords_expanded[4*c+2, 1] = dm_coords[2*c+1, 1] + cell_height - dm_coords_expanded[4*c+3, 1] = dm_coords[2*c+0, 1] + cell_height + dm_coords_expanded[4*c+0, 0] = dm_coords_orig[2*c+0, 0] + dm_coords_expanded[4*c+1, 0] = dm_coords_orig[2*c+1, 0] + dm_coords_expanded[4*c+2, 0] = dm_coords_orig[2*c+1, 0] + dm_coords_expanded[4*c+3, 0] = dm_coords_orig[2*c+0, 0] + dm_coords_expanded[4*c+0, 1] = dm_coords_orig[2*c+0, 1] + dm_coords_expanded[4*c+1, 1] = dm_coords_orig[2*c+1, 1] + dm_coords_expanded[4*c+2, 1] = dm_coords_orig[2*c+1, 1] + cell_height + dm_coords_expanded[4*c+3, 1] = dm_coords_orig[2*c+0, 1] + cell_height else: raise AssertionError("This case should not be hit") @@ -2192,7 +2190,7 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): f"{dm.getCellType(cStart)} is not supported") else: - dm_coords = dm_coords.reshape(coords_shape) + dm_coords = dm.getCellCoordinatesLocal().array_r.reshape(coords_shape) dm_sec = dm.getCellCoordinateSection() return dm_coords, dm_sec diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index f50cc06811..1f7c77e2a6 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -104,6 +104,8 @@ cdef extern from "petscdm.h" nogil: PetscErrorCode DMGetLabelValue(PETSc.PetscDM,char[],PetscInt,PetscInt*) PetscErrorCode DMGetPeriodicity(PETSc.PetscDM,PetscReal *[], PetscReal *[], PetscReal *[]) + PetscErrorCode DMGetSparseLocalize(PETSc.PetscDM,PetscBool *) + PetscErrorCode DMSetSparseLocalize(PETSc.PetscDM,PetscBool) cdef extern from "petscdmswarm.h" nogil: PetscErrorCode DMSwarmGetLocalSize(PETSc.PetscDM,PetscInt*) diff --git a/tests/firedrake/extrusion/test_variable_layers_numbering.py b/tests/firedrake/extrusion/test_variable_layers_numbering.py index 9eadb6d5b2..470dae4807 100644 --- a/tests/firedrake/extrusion/test_variable_layers_numbering.py +++ b/tests/firedrake/extrusion/test_variable_layers_numbering.py @@ -497,8 +497,8 @@ def test_layer_extents_parallel(): [0, 2, 0, 2], [0, 2, 0, 2], [0, 2, 0, 2], - [0, 2, 0, 2], [0, 3, 0, 2], + [0, 2, 0, 2], [0, 2, 0, 2]], dtype=IntType) elif mesh.comm.rank == 1: # Top view, plex points @@ -545,8 +545,8 @@ def test_layer_extents_parallel(): # edges [0, 3, 0, 3], [0, 3, 0, 3], - [0, 2, 0, 2], [0, 3, 0, 2], + [0, 2, 0, 2], [0, 2, 0, 2]], dtype=IntType) elif mesh.comm.rank == 3: # Top view, plex points @@ -569,11 +569,11 @@ def test_layer_extents_parallel(): [0, 2, 0, 2], [0, 3, 0, 3], # edges - [0, 2, 0, 2], [0, 3, 0, 2], [0, 2, 0, 2], [0, 2, 0, 2], [0, 2, 0, 2], + [0, 2, 0, 2], [0, 3, 0, 3], [0, 3, 0, 3]], dtype=IntType) assert numpy.equal(extmesh.layer_extents, expected).all() From 76ab04f3efe572af051a0565ce52ec4d245cf17d Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Thu, 12 Mar 2026 13:01:40 +0000 Subject: [PATCH 06/10] Apply suggestions from code review Co-authored-by: Connor Ward --- firedrake/cython/dmcommon.pyx | 35 +++++++++++-------- firedrake/utility_meshes.py | 12 +++---- .../firedrake/regression/test_periodic_2d.py | 3 +- 3 files changed, 28 insertions(+), 22 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 1cfe46188d..4b76bac971 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -2163,22 +2163,20 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): dm_coords_expanded[4*c+3, 1] = dm_coords_orig[2*c+1, 1] else: - if vert_unit_periodic: - cell_height = L[1] + assert vert_unit_periodic: + cell_height = L[1] - for c in range(cStart, cEnd): - CHKERR(PetscSectionSetDof(dm_sec_expanded.sec, c, 8)) + for c in range(cStart, cEnd): + CHKERR(PetscSectionSetDof(dm_sec_expanded.sec, c, 8)) - dm_coords_expanded[4*c+0, 0] = dm_coords_orig[2*c+0, 0] - dm_coords_expanded[4*c+1, 0] = dm_coords_orig[2*c+1, 0] - dm_coords_expanded[4*c+2, 0] = dm_coords_orig[2*c+1, 0] - dm_coords_expanded[4*c+3, 0] = dm_coords_orig[2*c+0, 0] - dm_coords_expanded[4*c+0, 1] = dm_coords_orig[2*c+0, 1] - dm_coords_expanded[4*c+1, 1] = dm_coords_orig[2*c+1, 1] - dm_coords_expanded[4*c+2, 1] = dm_coords_orig[2*c+1, 1] + cell_height - dm_coords_expanded[4*c+3, 1] = dm_coords_orig[2*c+0, 1] + cell_height - else: - raise AssertionError("This case should not be hit") + dm_coords_expanded[4*c+0, 0] = dm_coords_orig[2*c+0, 0] + dm_coords_expanded[4*c+1, 0] = dm_coords_orig[2*c+1, 0] + dm_coords_expanded[4*c+2, 0] = dm_coords_orig[2*c+1, 0] + dm_coords_expanded[4*c+3, 0] = dm_coords_orig[2*c+0, 0] + dm_coords_expanded[4*c+0, 1] = dm_coords_orig[2*c+0, 1] + dm_coords_expanded[4*c+1, 1] = dm_coords_orig[2*c+1, 1] + dm_coords_expanded[4*c+2, 1] = dm_coords_orig[2*c+1, 1] + cell_height + dm_coords_expanded[4*c+3, 1] = dm_coords_orig[2*c+0, 1] + cell_height dm_sec_expanded.setUp() @@ -2195,7 +2193,15 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): return dm_coords, dm_sec + def _get_periodicity(dm: PETSc.DM) -> tuple[tuple[bool, bool], ...]: + """Return mesh periodicity information. + + This function returns a 2-tuple of bools per dimension where the first entry indicates + whether the mesh is periodic in that dimension, and the second indicates whether the + mesh is single-cell periodic in that dimension. + + """ cdef: const PetscReal *maxCell, *L @@ -2206,6 +2212,7 @@ def _get_periodicity(dm: PETSc.DM) -> tuple[tuple[bool, bool], ...]: for d in range(dim) ) + @cython.boundscheck(False) @cython.wraparound(False) def mark_entity_classes(PETSc.DM dm): diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index 9294868c5e..399bd020f9 100644 --- a/firedrake/utility_meshes.py +++ b/firedrake/utility_meshes.py @@ -931,9 +931,9 @@ def PeriodicRectangleMesh( The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 - * 2: plane x == 1 + * 2: plane x == Lx * 3: plane y == 0 - * 4: plane y == 1 + * 4: plane y == Ly If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and if periodic in 'y' then 3 and 4 are empty. @@ -1039,9 +1039,9 @@ def PeriodicSquareMesh( The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 - * 2: plane x == 1 + * 2: plane x == L * 3: plane y == 0 - * 4: plane y == 1 + * 4: plane y == L If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and if periodic in 'y' then 3 and 4 are empty. @@ -3051,9 +3051,9 @@ def PartiallyPeriodicRectangleMesh( The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 - * 2: plane x == 1 + * 2: plane x == Lx * 3: plane y == 0 - * 4: plane y == 1 + * 4: plane y == Ly If periodic in the 'x' direction then boundary edges 1 and 2 are empty, and if periodic in 'y' then 3 and 4 are empty. diff --git a/tests/firedrake/regression/test_periodic_2d.py b/tests/firedrake/regression/test_periodic_2d.py index b3568a6be5..8d5eb5fac6 100644 --- a/tests/firedrake/regression/test_periodic_2d.py +++ b/tests/firedrake/regression/test_periodic_2d.py @@ -38,10 +38,9 @@ def test_periodic_helmholtz(direction, cell_options): f = Function(V).assign((244.0*pi*pi/225.0 + 1.0)*u_exact) - # FIXME: This suggests some really weird behaviour! if direction == "x": bcs = DirichletBC(V, Constant(0), (3, 4)) - if direction == "y": + elif direction == "y": bcs = DirichletBC(V, Constant(0), (1, 2)) elif direction == "both": bcs = [] From 1f80b72ed6b5c7da63378f6eb1f02509f4678103 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Thu, 12 Mar 2026 13:18:44 +0000 Subject: [PATCH 07/10] fixup --- firedrake/cython/dmcommon.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 4b76bac971..8a4b8d2ebb 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -2163,7 +2163,7 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): dm_coords_expanded[4*c+3, 1] = dm_coords_orig[2*c+1, 1] else: - assert vert_unit_periodic: + assert vert_unit_periodic cell_height = L[1] for c in range(cStart, cEnd): From 9ce5532282650226296c64093b91d82fba0d63ac Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Thu, 12 Mar 2026 17:33:58 +0000 Subject: [PATCH 08/10] Add warning for empty subdomains --- firedrake/bcs.py | 18 +++++++++++++++--- firedrake/mesh.py | 21 ++++++++++----------- 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/firedrake/bcs.py b/firedrake/bcs.py index 921f879de4..2203f2b3d5 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -1,9 +1,11 @@ # A module implementing strong (Dirichlet) boundary conditions. -import numpy as np from functools import partial, reduce, cached_property import itertools +import numpy as np +from mpi4py import MPI + import ufl from ufl import as_ufl, as_tensor from finat.ufl import VectorElement @@ -11,6 +13,7 @@ import pyop2 as op2 from pyop2 import exceptions +from pyop2.mpi import temp_internal_comm from pyop2.utils import as_tuple import firedrake @@ -19,6 +22,7 @@ from firedrake import slate from firedrake import solving from firedrake.formmanipulation import ExtractSubBlock +from firedrake.logging import logger from firedrake.adjoint_utils.dirichletbc import DirichletBCMixin from firedrake.petsc import PETSc @@ -147,7 +151,7 @@ def hermite_stride(bcnodes): bcnodes = np.setdiff1d(bcnodes, deriv_ids) return bcnodes - sub_d = (self.sub_domain, ) if isinstance(self.sub_domain, str) else as_tuple(self.sub_domain) + sub_d = (self.sub_domain,) if isinstance(self.sub_domain, str) else as_tuple(self.sub_domain) sub_d = [s if isinstance(s, str) else as_tuple(s) for s in sub_d] bcnodes = [] for s in sub_d: @@ -168,7 +172,15 @@ def hermite_stride(bcnodes): bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss))) bcnodes1 = reduce(np.intersect1d, bcnodes1) bcnodes.append(bcnodes1) - return np.concatenate(bcnodes) + bcnodes = np.concatenate(bcnodes) + + with temp_internal_comm(self._function_space.mesh().comm) as icomm: + num_global_nodes = icomm.reduce(len(bcnodes), MPI.SUM, root=0) + if num_global_nodes == 0 and icomm.rank == 0: + logger.warn(f"Subdomain {self.sub_domain} is empty. This is likely an error. " + "Did you choose the right label?") + + return bcnodes @cached_property def node_set(self): diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 6bb990f1cd..5b382c89c7 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -35,7 +35,7 @@ import firedrake.cython.spatialindex as spatialindex import firedrake.utils as utils from firedrake.utils import as_cstr, IntType, RealType -from firedrake.logging import info_red +from firedrake.logging import info_red, logger from firedrake.parameters import parameters from firedrake.petsc import PETSc, DEFAULT_PARTITIONER from firedrake.adjoint_utils import MeshGeometryMixin @@ -189,14 +189,6 @@ def __init__(self, mesh, facets, classes, set_, kind, facet_cell, local_facet_nu self.unique_markers = [] if unique_markers is None else unique_markers self._subsets = {} - @cached_property - def _null_subset(self): - '''Empty subset for the case in which there are no facets with - a given marker value. This is required because not all - markers need be represented on all processors.''' - - return op2.Subset(self.set, []) - @PETSc.Log.EventDecorator() def measure_set(self, integral_type, subdomain_id, all_integer_subdomain_ids=None): @@ -271,9 +263,16 @@ def subset(self, markers): marked_points_list.append(self.mesh.topology_dm.getStratumIS(dmcommon.FACE_SETS_LABEL, i).indices) if marked_points_list: _, indices, _ = np.intersect1d(self.facets, np.concatenate(marked_points_list), return_indices=True) - return self._subsets.setdefault(markers, op2.Subset(self.set, indices)) else: - return self._subsets.setdefault(markers, self._null_subset) + indices = np.empty(0, dtype=IntType) + + with temp_internal_comm(self.mesh.comm) as icomm: + num_global_indices = icomm.reduce(len(indices), MPI.SUM, root=0) + if num_global_indices == 0 and icomm.rank == 0: + logger.warn(f"Subdomain {markers} is empty. This is likely an error. " + "Did you choose the right label?") + + return self._subsets.setdefault(markers, op2.Subset(self.set, indices)) def _collect_unmarked_points(self, markers): """Collect points that are not marked by markers.""" From 3b27dbc7c2ac6190d7fe6198b3dc4c1e05be976e Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Wed, 18 Mar 2026 14:45:54 +0000 Subject: [PATCH 09/10] Apply suggestions from code review Co-authored-by: Connor Ward --- .github/workflows/core.yml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index e54f0f83aa..c0682518a4 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -158,9 +158,7 @@ jobs: --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ https://gitlab.com/petsc/petsc.git else - : # UNDO ME - : # git clone --depth 1 https://gitlab.com/petsc/petsc.git - git clone --depth 1 https://gitlab.com/connorjward/petsc.git --branch connorjward/periodic-firedrake + git clone --depth 1 https://gitlab.com/petsc/petsc.git fi cd petsc python3 ../firedrake-repo/scripts/firedrake-configure \ From 94835133496956670dee7529aa3ede5f60be9c56 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Tue, 24 Mar 2026 09:12:53 +0000 Subject: [PATCH 10/10] Try final PETSc fix --- .github/workflows/core.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index e738561b9c..08be16cb31 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -158,7 +158,9 @@ jobs: --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ https://gitlab.com/petsc/petsc.git else - git clone --depth 1 https://gitlab.com/petsc/petsc.git + : # UNDO ME + : # git clone --depth 1 https://gitlab.com/petsc/petsc.git + git clone --depth 1 https://gitlab.com/connorjward/petsc.git --branch connorjward/tosimplex-periodic-fix fi cd petsc python3 ../firedrake-repo/scripts/firedrake-configure \