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 \ diff --git a/docs/source/conf.py b/docs/source/conf.py index 1377b2d2a8..9af63da1af 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/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/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ec7ab631af..8a4b8d2ebb 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,131 @@ 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() + coords_shape = ((cEnd-cStart) * ndofs[0], dim) + + 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 + # 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_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() + 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_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: + assert 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_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() + + 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.getCellCoordinatesLocal().array_r.reshape(coords_shape) + dm_sec = dm.getCellCoordinateSection() + + 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 + + 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..1f7c77e2a6 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -103,6 +103,10 @@ 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 *[]) + PetscErrorCode DMGetSparseLocalize(PETSc.PetscDM,PetscBool *) + PetscErrorCode DMSetSparseLocalize(PETSc.PetscDM,PetscBool) + cdef extern from "petscdmswarm.h" nogil: PetscErrorCode DMSwarmGetLocalSize(PETSc.PetscDM,PetscInt*) PetscErrorCode DMSwarmGetCellDM(PETSc.PetscDM, PETSc.PetscDM*) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 59531066c1..36b38aa0b4 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -37,7 +37,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 @@ -201,14 +201,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): @@ -283,9 +275,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.""" diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index f75cd514ae..399bd020f9 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 + * 3: plane y == 0 + * 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. + """ + 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 + * 3: plane y == 0 + * 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. """ 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) + * 3: plane y == 0 + * 4: plane y == Ly - # 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..e523945651 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], # 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')) - 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], # 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: # 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..470dae4807 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"], @@ -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() 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..8d5eb5fac6 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,9 @@ 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"): + if direction == "x": + bcs = DirichletBC(V, Constant(0), (3, 4)) + elif direction == "y": bcs = DirichletBC(V, Constant(0), (1, 2)) elif direction == "both": bcs = [] @@ -72,12 +56,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..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) @@ -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], # noqa + [0, 0, 0, 0, 1/3, 0, 1/6, 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) + 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], # noqa + [1/3, 0, 1/6, 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) + 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], # 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) @@ -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]]) # noqa + M10ww = np.array([[1/30, 0, 0, 1/20, 0, 0], + [1/20, 0, 0, 1/5, 0, 0]]) # noqa # 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..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 @@ -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: @@ -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() 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