From e76eb74facefd16ddeea3af199e9b5b5564f5078 Mon Sep 17 00:00:00 2001 From: ksagiyam Date: Fri, 24 Oct 2025 10:32:15 +0100 Subject: [PATCH] Enable solving multi-domain problems involving codim-1 submeshes Co-authored-by: Connor Ward --- firedrake/assemble.py | 33 ++ firedrake/embedding.py | 2 +- firedrake/functionspaceimpl.py | 4 + firedrake/interpolation.py | 2 +- firedrake/mesh.py | 242 +++++++++++---- firedrake/mg/kernels.py | 9 +- firedrake/pointeval_utils.py | 1 + firedrake/pointquery_utils.py | 1 + firedrake/slate/slac/compiler.py | 1 + .../meshes/mixed_cell_unit_square.msh | 231 ++++++++++++++ .../submesh/test_submesh_assemble.py | 206 ++++++++++++ tests/firedrake/submesh/test_submesh_base.py | 2 +- tests/firedrake/submesh/test_submesh_solve.py | 293 ++++++++++++++++++ tsfc/driver.py | 3 + tsfc/fem.py | 151 +++++---- tsfc/kernel_args.py | 4 + tsfc/kernel_interface/__init__.py | 4 + tsfc/kernel_interface/common.py | 13 +- tsfc/kernel_interface/firedrake_loopy.py | 29 +- 19 files changed, 1102 insertions(+), 129 deletions(-) create mode 100644 tests/firedrake/meshes/mixed_cell_unit_square.msh diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 641125f735..5db7c007ce 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1661,6 +1661,7 @@ def __init__(self, form, local_knl, subdomain_id, all_integer_subdomain_ids, dia self._constants = _FormHandler.iter_constants(form, local_knl.kinfo) self._active_exterior_facets = _FormHandler.iter_active_exterior_facets(form, local_knl.kinfo) self._active_interior_facets = _FormHandler.iter_active_interior_facets(form, local_knl.kinfo) + self._active_orientations_cell = _FormHandler.iter_active_orientations_cell(form, local_knl.kinfo) self._active_orientations_exterior_facet = _FormHandler.iter_active_orientations_exterior_facet(form, local_knl.kinfo) self._active_orientations_interior_facet = _FormHandler.iter_active_orientations_interior_facet(form, local_knl.kinfo) @@ -1683,6 +1684,7 @@ def build(self): assert_empty(self._constants) assert_empty(self._active_exterior_facets) assert_empty(self._active_interior_facets) + assert_empty(self._active_orientations_cell) assert_empty(self._active_orientations_exterior_facet) assert_empty(self._active_orientations_interior_facet) @@ -1880,6 +1882,17 @@ def _as_global_kernel_arg_interior_facet(_, self): return op2.DatKernelArg((2,), m._global_kernel_arg) +@_as_global_kernel_arg.register(kernel_args.OrientationsCellKernelArg) +def _(_, self): + mesh = next(self._active_orientations_cell) + if mesh is self._mesh: + return op2.DatKernelArg((1,)) + else: + m, integral_type = mesh.topology.trans_mesh_entity_map(self._mesh.topology, self._integral_type, self._subdomain_id, self._all_integer_subdomain_ids) + assert integral_type == "cell" + return op2.DatKernelArg((1,), m._global_kernel_arg) + + @_as_global_kernel_arg.register(kernel_args.OrientationsExteriorFacetKernelArg) def _(_, self): mesh = next(self._active_orientations_exterior_facet) @@ -1951,6 +1964,7 @@ def __init__(self, form, bcs, local_knl, subdomain_id, self._constants = _FormHandler.iter_constants(form, local_knl.kinfo) self._active_exterior_facets = _FormHandler.iter_active_exterior_facets(form, local_knl.kinfo) self._active_interior_facets = _FormHandler.iter_active_interior_facets(form, local_knl.kinfo) + self._active_orientations_cell = _FormHandler.iter_active_orientations_cell(form, local_knl.kinfo) self._active_orientations_exterior_facet = _FormHandler.iter_active_orientations_exterior_facet(form, local_knl.kinfo) self._active_orientations_interior_facet = _FormHandler.iter_active_orientations_interior_facet(form, local_knl.kinfo) @@ -2216,6 +2230,17 @@ def _as_parloop_arg_interior_facet(_, self): return op2.DatParloopArg(mesh.interior_facets.local_facet_dat, m) +@_as_parloop_arg.register(kernel_args.OrientationsCellKernelArg) +def _(_, self): + mesh = next(self._active_orientations_cell) + if mesh is self._mesh: + m = None + else: + m, integral_type = mesh.topology.trans_mesh_entity_map(self._mesh.topology, self._integral_type, self._subdomain_id, self._all_integer_subdomain_ids) + assert integral_type == "cell" + return op2.DatParloopArg(mesh.local_cell_orientation_dat, m) + + @_as_parloop_arg.register(kernel_args.OrientationsExteriorFacetKernelArg) def _(_, self): mesh = next(self._active_orientations_exterior_facet) @@ -2312,6 +2337,14 @@ def iter_active_interior_facets(form, kinfo): mesh = all_meshes[i] yield mesh + @staticmethod + def iter_active_orientations_cell(form, kinfo): + """Yield the form cell orientations referenced in ``kinfo``.""" + all_meshes = extract_domains(form) + for i in kinfo.active_domain_numbers.orientations_cell: + mesh = all_meshes[i] + yield mesh + @staticmethod def iter_active_orientations_exterior_facet(form, kinfo): """Yield the form exterior facet orientations referenced in ``kinfo``.""" diff --git a/firedrake/embedding.py b/firedrake/embedding.py index 8bbbb5af42..15f8a09075 100644 --- a/firedrake/embedding.py +++ b/firedrake/embedding.py @@ -5,7 +5,7 @@ def get_embedding_dg_element(element, value_shape, broken_cg=False): - cell = element.cell + cell, = set(element.cell.cells) family = lambda c: "DG" if c.is_simplex else "DQ" if isinstance(cell, ufl.TensorProductCell): degree = element.degree() diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index 1edc0f4fb8..6dd1aebbdc 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -13,6 +13,7 @@ import ufl import finat.ufl +from ufl.cell import CellSequence from ufl.duals import is_dual, is_primal from pyop2 import op2, mpi from pyop2.utils import as_tuple @@ -52,6 +53,9 @@ def check_element(element, top=True): ValueError If the element is illegal. """ + if isinstance(element.cell, CellSequence) and \ + type(element) is not finat.ufl.MixedElement: + raise ValueError("MixedElement modifier must be outermost") if element.cell.cellname == "hexahedron" and \ element.family() not in ["Q", "DQ", "Real"]: raise NotImplementedError("Currently can only use 'Q', 'DQ', and/or 'Real' elements on hexahedral meshes, not", element.family()) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 85f17400cd..b656a47b21 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -1198,7 +1198,7 @@ def get_interp_node_map(source_mesh, target_mesh, fs): else: raise ValueError("Have coefficient with unexpected mesh") else: - m_ = fs.entity_node_map(target_mesh.topology, "cell", None, None) + m_ = fs.entity_node_map(target_mesh.topology, "cell", "everywhere", None) return m_ diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 296ab0e32b..215ea8f93c 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -10,6 +10,7 @@ from collections import OrderedDict, defaultdict from collections.abc import Sequence from ufl.classes import ReferenceGrad +from ufl.cell import CellSequence from ufl.domain import extract_unique_domain import enum import numbers @@ -157,13 +158,13 @@ class _Facets(object): The unique_markers argument **must** be the same on all processes.""" @PETSc.Log.EventDecorator() - def __init__(self, mesh, facets, classes, kind, facet_cell, local_facet_number, + def __init__(self, mesh, facets, classes, set_, kind, facet_cell, local_facet_number, unique_markers=None): self.mesh = mesh self.facets = facets - classes = as_tuple(classes, int, 3) self.classes = classes + self.set = set_ self.kind = kind assert kind in ["interior", "exterior"] @@ -188,17 +189,6 @@ def __init__(self, mesh, facets, classes, kind, facet_cell, local_facet_number, self.unique_markers = [] if unique_markers is None else unique_markers self._subsets = {} - @utils.cached_property - def set(self): - size = self.classes - if isinstance(self.mesh, ExtrudedMeshTopology): - label = "%s_facets" % self.kind - layers = self.mesh.entity_layers(1, label) - base = getattr(self.mesh._base_mesh, label).set - return op2.ExtrudedSet(base, layers=layers) - return op2.Set(size, "%sFacets" % self.kind.capitalize()[:3], - comm=self.mesh.comm) - @utils.cached_property def _null_subset(self): '''Empty subset for the case in which there are no facets with @@ -318,22 +308,36 @@ def local_facet_orientation_dat(self): local_facet_end = offsets[-2] map_from_cell_to_facet_orientations = self.mesh.entity_orientations[:, local_facet_start:local_facet_end] # Make output data; - # this is a map from an exterior/interior facet to the corresponding local facet orientation/orientations. - # Halo data are required by design, but not actually used. - # -- Reshape as (-1, self._rank) to uniformly handle exterior and interior facets. - data = np.empty_like(self.local_facet_dat.data_ro_with_halos).reshape((-1, self._rank)) - data.fill(np.iinfo(dtype).max) - # Set local facet orientations on the block corresponding to the owned facets; i.e., data[:shape[0], :] below. - local_facets = self.local_facet_dat.data_ro # do not need halos. - # -- Reshape as (-1, self._rank) to uniformly handle exterior and interior facets. - local_facets = local_facets.reshape((-1, self._rank)) - shape = local_facets.shape - map_from_owned_facet_to_cells = self.facet_cell[:shape[0], :] - data[:shape[0], :] = np.take_along_axis( - map_from_cell_to_facet_orientations[map_from_owned_facet_to_cells], - local_facets.reshape(shape + (1, )), # reshape as required by take_along_axis. + # this is a map from an exterior/interior facet to the corresponding + # local facet orientation/orientations. + # The local facet orientation/orientations of a halo facet is/are also + # used in some submesh problems. + # + # Example: + # + # +-------+-------+ + # | | | + # meshA | g g o | + # | | | + # +-------+-------+ + # +-------+ + # | | + # meshB o o | o: owned + # | | g: ghost + # +-------+ + # + # form = FacetNormal(meshA)[0] * ds(meshB, interface) + # + # Reshape local_facets as (-1, self._rank) to uniformly handle exterior and interior facets. + local_facets = self.local_facet_dat.data_ro_with_halos.reshape((-1, self._rank)) + # Make slice for masking out rows for which orientations are not needed. + slice_ = (self.facet_cell != -1).all(axis=1) + data = np.full_like(local_facets, np.iinfo(dtype).max) + data[slice_, :] = np.take_along_axis( + map_from_cell_to_facet_orientations[self.facet_cell[slice_, :]], + local_facets.reshape(local_facets.shape + (1, ))[slice_, :, :], # reshape as required by take_along_axis. axis=2, - ).reshape(shape) + ).reshape((-1, self._rank)) return op2.Dat( self.local_facet_dat.dataset, data, @@ -705,6 +709,12 @@ def entity_orientations(self): """ pass + @property + @abc.abstractmethod + def local_cell_orientation_dat(self): + """Local cell orientation dat.""" + pass + @abc.abstractmethod def _facets(self, kind): pass @@ -1223,7 +1233,8 @@ def cell_closure(self): cell = self.ufl_cell() assert tdim == cell.topological_dimension if self.submesh_parent is not None and \ - not (self.submesh_parent.ufl_cell().cellname == "hexahedron" and cell.cellname == "quadrilateral"): + not (self.submesh_parent.ufl_cell().cellname == "hexahedron" and cell.cellname == "quadrilateral") and \ + len(self.submesh_parent.dm_cell_types) == 1: # Codim-1 submesh of a hex mesh (i.e. a quad submesh) can not # inherit cell_closure from the hex mesh as the cell_closure # must follow the special orientation restriction. This means @@ -1282,14 +1293,23 @@ def cell_closure(self): def entity_orientations(self): return dmcommon.entity_orientations(self, self.cell_closure) + @utils.cached_property + def local_cell_orientation_dat(self): + """Local cell orientation dat.""" + return op2.Dat( + op2.DataSet(self.cell_set, 1), + self.entity_orientations[:, [-1]], + gem.uint_type, + f"{self.name}_local_cell_orientation" + ) + @PETSc.Log.EventDecorator() def _facets(self, kind): if kind not in ["interior", "exterior"]: raise ValueError("Unknown facet type '%s'" % kind) dm = self.topology_dm - facets, classes = dmcommon.get_facets_by_class(dm, (kind + "_facets"), - self._facet_ordering) + facets, classes, set_ = getattr(self, "_" + kind + "_facet_numbers_classes_set") label = dmcommon.FACE_SETS_LABEL if dm.hasLabel(label): from mpi4py import MPI @@ -1313,7 +1333,7 @@ def merge_ids(x, y, datatype): point2facetnumber = np.full(facets.max(initial=0)+1, -1, dtype=IntType) point2facetnumber[facets] = np.arange(len(facets), dtype=IntType) - obj = _Facets(self, facets, classes, kind, + obj = _Facets(self, facets, classes, set_, kind, facet_cell, local_facet_number, unique_markers=unique_markers) obj.point2facetnumber = point2facetnumber @@ -1327,6 +1347,25 @@ def exterior_facets(self): def interior_facets(self): return self._facets("interior") + def _facet_numbers_classes_set(self, kind): + if kind not in ["interior", "exterior"]: + raise ValueError("Unknown facet type '%s'" % kind) + # Can not call target.{interior, exterior}_facets.facets + # if target is a mixed cell mesh (cell_closure etc. can not be defined), + # so directly call dmcommon.get_facets_by_class. + _numbers, _classes = dmcommon.get_facets_by_class(self.topology_dm, (kind + "_facets"), self._facet_ordering) + _classes = as_tuple(_classes, int, 3) + _set = op2.Set(_classes, f"{kind.capitalize()[:3]}Facets", comm=self.comm) + return _numbers, _classes, _set + + @utils.cached_property + def _exterior_facet_numbers_classes_set(self): + return self._facet_numbers_classes_set("exterior") + + @utils.cached_property + def _interior_facet_numbers_classes_set(self): + return self._facet_numbers_classes_set("interior") + @utils.cached_property def cell_to_facets(self): """Returns a :class:`pyop2.types.dat.Dat` that maps from a cell index to the local @@ -1493,11 +1532,15 @@ def submesh_child_cell_parent_cell_map(self): @utils.cached_property def submesh_child_exterior_facet_parent_exterior_facet_map(self): - return self._submesh_make_entity_entity_map(self.exterior_facets.set, self.submesh_parent.exterior_facets.set, self.exterior_facets.facets, self.submesh_parent.exterior_facets.facets, True) + _self_numbers, _, _self_set = self._exterior_facet_numbers_classes_set + _parent_numbers, _, _parent_set = self.submesh_parent._exterior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_self_set, _parent_set, _self_numbers, _parent_numbers, True) @utils.cached_property def submesh_child_exterior_facet_parent_interior_facet_map(self): - return self._submesh_make_entity_entity_map(self.exterior_facets.set, self.submesh_parent.interior_facets.set, self.exterior_facets.facets, self.submesh_parent.interior_facets.facets, True) + _self_numbers, _, _self_set = self._exterior_facet_numbers_classes_set + _parent_numbers, _, _parent_set = self.submesh_parent._interior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_self_set, _parent_set, _self_numbers, _parent_numbers, True) @utils.cached_property def submesh_child_interior_facet_parent_exterior_facet_map(self): @@ -1505,7 +1548,19 @@ def submesh_child_interior_facet_parent_exterior_facet_map(self): @utils.cached_property def submesh_child_interior_facet_parent_interior_facet_map(self): - return self._submesh_make_entity_entity_map(self.interior_facets.set, self.submesh_parent.interior_facets.set, self.interior_facets.facets, self.submesh_parent.interior_facets.facets, True) + _self_numbers, _, _self_set = self._interior_facet_numbers_classes_set + _parent_numbers, _, _parent_set = self.submesh_parent._interior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_self_set, _parent_set, _self_numbers, _parent_numbers, True) + + @utils.cached_property + def submesh_child_cell_parent_interior_facet_map(self): + _parent_numbers, _, _parent_set = self.submesh_parent._interior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(self.cell_set, _parent_set, self.cell_closure[:, -1], _parent_numbers, True) + + @utils.cached_property + def submesh_child_cell_parent_exterior_facet_map(self): + _parent_numbers, _, _parent_set = self.submesh_parent._exterior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(self.cell_set, _parent_set, self.cell_closure[:, -1], _parent_numbers, True) @utils.cached_property def submesh_parent_cell_child_cell_map(self): @@ -1513,7 +1568,9 @@ def submesh_parent_cell_child_cell_map(self): @utils.cached_property def submesh_parent_exterior_facet_child_exterior_facet_map(self): - return self._submesh_make_entity_entity_map(self.submesh_parent.exterior_facets.set, self.exterior_facets.set, self.submesh_parent.exterior_facets.facets, self.exterior_facets.facets, False) + _self_numbers, _, _self_set = self._exterior_facet_numbers_classes_set + _parent_numbers, _, _parent_set = self.submesh_parent._exterior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_parent_set, _self_set, _parent_numbers, _self_numbers, False) @utils.cached_property def submesh_parent_exterior_facet_child_interior_facet_map(self): @@ -1521,11 +1578,25 @@ def submesh_parent_exterior_facet_child_interior_facet_map(self): @utils.cached_property def submesh_parent_interior_facet_child_exterior_facet_map(self): - return self._submesh_make_entity_entity_map(self.submesh_parent.interior_facets.set, self.exterior_facets.set, self.submesh_parent.interior_facets.facets, self.exterior_facets.facets, False) + _self_numbers, _, _self_set = self._exterior_facet_numbers_classes_set + _parent_numbers, _, _parent_set = self.submesh_parent._interior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_parent_set, _self_set, _parent_numbers, _self_numbers, False) @utils.cached_property def submesh_parent_interior_facet_child_interior_facet_map(self): - return self._submesh_make_entity_entity_map(self.submesh_parent.interior_facets.set, self.interior_facets.set, self.submesh_parent.interior_facets.facets, self.interior_facets.facets, False) + _self_numbers, _, _self_set = self._interior_facet_numbers_classes_set + _parent_numbers, _, _parent_set = self.submesh_parent._interior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_parent_set, _self_set, _parent_numbers, _self_numbers, False) + + @utils.cached_property + def submesh_parent_exterior_facet_child_cell_map(self): + _parent_numbers, _, _parent_set = self.submesh_parent._exterior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_parent_set, self.cell_set, _parent_numbers, self.cell_closure[:, -1], False) + + @utils.cached_property + def submesh_parent_interior_facet_child_cell_map(self): + _parent_numbers, _, _parent_set = self.submesh_parent._interior_facet_numbers_classes_set + return self._submesh_make_entity_entity_map(_parent_set, self.cell_set, _parent_numbers, self.cell_closure[:, -1], False) def submesh_map_child_parent(self, source_integral_type, source_subset_points, reverse=False): """Return the map from submesh child entities to submesh parent entities or its reverse. @@ -1555,20 +1626,52 @@ def submesh_map_child_parent(self, source_integral_type, source_subset_points, r target = self.submesh_parent target_dim = target.topology_dm.getDimension() source_dim = source.topology_dm.getDimension() - if source_dim != target_dim: - raise NotImplementedError(f"Not implemented for (source_dim, target_dim) == ({source_dim}, {target_dim})") - if source_integral_type == "cell": + if target_dim == source_dim: + if source_integral_type == "cell": + target_integral_type_temp = "cell" + elif source_integral_type in ["interior_facet", "exterior_facet"]: + target_integral_type_temp = "facet" + else: + raise NotImplementedError("Unsupported combination") + elif target_dim - 1 == source_dim: + if source_integral_type == "cell": + target_integral_type_temp = "facet" + else: + raise NotImplementedError("Unsupported combination") + elif target_dim == source_dim - 1: + if source_integral_type in ["interior_facet", "exterior_facet"]: + target_integral_type_temp = "cell" + else: + raise NotImplementedError("Unsupported combination") + else: + raise NotImplementedError("Unsupported combination") + if target_integral_type_temp == "cell": + _cell_numbers = target.cell_closure[:, -1] + with self.topology_dm.getSubpointIS() as subpoints: + if reverse: + _, target_indices_cell, source_indices_cell = np.intersect1d(subpoints[_cell_numbers], source_subset_points, return_indices=True) + else: + target_subset_points = subpoints[source_subset_points] + _, target_indices_cell, source_indices_cell = np.intersect1d(_cell_numbers, target_subset_points, return_indices=True) + n_cell = len(source_indices_cell) + n_cell_max = self._comm.allreduce(n_cell, op=MPI.MAX) + if n_cell_max > 0: + if n_cell > len(source_subset_points): + raise RuntimeError("Found inconsistent data") target_integral_type = "cell" - target_subset_points = None - elif source_integral_type in ["interior_facet", "exterior_facet"]: + if reverse: + target_subset_points = _cell_numbers[target_indices_cell] + elif target_integral_type_temp == "facet": + _exterior_facet_numbers, _, _ = target._exterior_facet_numbers_classes_set + _interior_facet_numbers, _, _ = target._interior_facet_numbers_classes_set with self.topology_dm.getSubpointIS() as subpoints: if reverse: - _, target_indices_int, source_indices_int = np.intersect1d(subpoints[target.interior_facets.facets], source_subset_points, return_indices=True) - _, target_indices_ext, source_indices_ext = np.intersect1d(subpoints[target.exterior_facets.facets], source_subset_points, return_indices=True) + _, target_indices_int, source_indices_int = np.intersect1d(subpoints[_interior_facet_numbers], source_subset_points, return_indices=True) + _, target_indices_ext, source_indices_ext = np.intersect1d(subpoints[_exterior_facet_numbers], source_subset_points, return_indices=True) else: target_subset_points = subpoints[source_subset_points] - _, target_indices_int, source_indices_int = np.intersect1d(target.interior_facets.facets, target_subset_points, return_indices=True) - _, target_indices_ext, source_indices_ext = np.intersect1d(target.exterior_facets.facets, target_subset_points, return_indices=True) + _, target_indices_int, source_indices_int = np.intersect1d(_interior_facet_numbers, target_subset_points, return_indices=True) + _, target_indices_ext, source_indices_ext = np.intersect1d(_exterior_facet_numbers, target_subset_points, return_indices=True) n_int = len(source_indices_int) n_ext = len(source_indices_ext) n_int_max = self._comm.allreduce(n_int, op=MPI.MAX) @@ -1589,11 +1692,11 @@ def submesh_map_child_parent(self, source_integral_type, source_subset_points, r raise RuntimeError("Can not find a map from source to target.") if reverse: if target_integral_type == "interior_facet": - target_subset_points = target.interior_facets.facets[target_indices_int] + target_subset_points = _interior_facet_numbers[target_indices_int] elif target_integral_type == "exterior_facet": - target_subset_points = target.exterior_facets.facets[target_indices_ext] + target_subset_points = _exterior_facet_numbers[target_indices_ext] else: - raise NotImplementedError(f"Not implemented for (source_dim, target_dim, source_integral_type) == ({source_dim}, {target_dim}, {source_integral_type})") + raise NotImplementedError if reverse: map_ = getattr(self, f"submesh_parent_{source_integral_type}_child_{target_integral_type}_map") else: @@ -1629,13 +1732,16 @@ def trans_mesh_entity_map(self, base_mesh, base_integral_type, base_subdomain_id raise NotImplementedError("Currenlty can not return identity map") else: if base_integral_type == "cell": - base_subset_points = None + base_subset = base_mesh.measure_set(base_integral_type, base_subdomain_id, all_integer_subdomain_ids=base_all_integer_subdomain_ids) + base_subset_points = base_mesh.cell_closure[:, -1][base_subset.indices] elif base_integral_type in ["interior_facet", "exterior_facet"]: base_subset = base_mesh.measure_set(base_integral_type, base_subdomain_id, all_integer_subdomain_ids=base_all_integer_subdomain_ids) if base_integral_type == "interior_facet": - base_subset_points = base_mesh.interior_facets.facets[base_subset.indices] + _interior_facet_numbers, _, _ = base_mesh._interior_facet_numbers_classes_set + base_subset_points = _interior_facet_numbers[base_subset.indices] elif base_integral_type == "exterior_facet": - base_subset_points = base_mesh.exterior_facets.facets[base_subset.indices] + _exterior_facet_numbers, _, _ = base_mesh._exterior_facet_numbers_classes_set + base_subset_points = _exterior_facet_numbers[base_subset.indices] else: raise NotImplementedError(f"Unknown integration type : {base_integral_type}") composed_map, integral_type, _ = self.submesh_map_composed(base_mesh, base_integral_type, base_subset_points) @@ -1732,11 +1838,19 @@ def cell_closure(self): def entity_orientations(self): return self._base_mesh.entity_orientations + @utils.cached_property + def local_cell_orientation_dat(self): + """Local cell orientation dat.""" + return self._base_mesh.local_cell_orientation_dat + def _facets(self, kind): if kind not in ["interior", "exterior"]: raise ValueError("Unknown facet type '%s'" % kind) - base = getattr(self._base_mesh, "%s_facets" % kind) - return _Facets(self, base.facets, base.classes, + label = f"{kind}_facets" + base = getattr(self._base_mesh, label) + layers = self.entity_layers(1, label) + set_ = op2.ExtrudedSet(base.set, layers=layers) + return _Facets(self, base.facets, base.classes, set_, kind, base.facet_cell, base.local_facet_dat.data_ro_with_halos, @@ -1991,6 +2105,11 @@ def cell_closure(self): entity_orientations = None + @property + def local_cell_orientation_dat(self): + """Local cell orientation dat.""" + raise NotImplementedError("Not implemented for VertexOnlyMeshTopology") + def _facets(self, kind): """Raises an AttributeError since cells in a `VertexOnlyMeshTopology` have no facets. @@ -2273,6 +2392,8 @@ def _init_topology(self, topology): import firedrake.function as function self._topology = topology + if len(topology.dm_cell_types) > 1: + return coordinates_fs = functionspace.FunctionSpace(self.topology, self.ufl_coordinate_element()) coordinates_data = dmcommon.reordered_coords(topology.topology_dm, coordinates_fs.dm.getDefaultSection(), (self.num_vertices(), self.geometric_dimension)) @@ -4846,13 +4967,14 @@ def topological(self): return self def ufl_cell(self): - cell, = set(m.ufl_cell() for m in self._meshes) - return cell + return CellSequence([m.ufl_cell() for m in self._meshes]) def ufl_mesh(self): - cell = self.ufl_cell() - return ufl.MeshSequence([ufl.Mesh(finat.ufl.VectorElement("Lagrange", cell, 1, dim=cell.topological_dimension)) - for _ in self._meshes]) + dim = self.ufl_cell().topological_dimension + return ufl.MeshSequence( + [ufl.Mesh(finat.ufl.VectorElement("Lagrange", cell, 1, dim=dim)) + for cell in self.ufl_cell().cells] + ) def __eq__(self, other): if type(other) != type(self): diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py index 864422541d..76d68445cc 100644 --- a/firedrake/mg/kernels.py +++ b/firedrake/mg/kernels.py @@ -144,6 +144,7 @@ def compile_element(expression, dual_space=None, parameters=None, # Replace coordinates (if any) builder = firedrake_interface.KernelBuilderBase(scalar_type=ScalarType) builder._domain_integral_type_map = {domain: "cell"} + builder._entity_ids = {domain: (0,)} # Translate to GEM cell = domain.ufl_cell() dim = cell.topological_dimension @@ -527,6 +528,7 @@ def dg_injection_kernel(Vf, Vc, ncell): raise NotImplementedError("In complex mode we are waiting for Slate") macro_builder = MacroKernelBuilder(ScalarType, ncell) macro_builder._domain_integral_type_map = {Vf.mesh(): "cell"} + macro_builder._entity_ids = {Vf.mesh(): (0,)} f = ufl.Coefficient(Vf) macro_builder.set_coefficients([f]) macro_builder.set_coordinates(Vf.mesh()) @@ -541,11 +543,10 @@ def dg_injection_kernel(Vf, Vc, ncell): macro_quadrature_rule = make_quadrature(ref_complex, estimate_total_polynomial_degree(ufl.inner(f, f))) index_cache = {} parameters = default_parameters() - integration_dim, entity_ids = lower_integral_type(Vfe.cell, "cell") + integration_dim, _ = lower_integral_type(Vfe.cell, "cell") macro_cfg = dict(interface=macro_builder, ufl_cell=Vf.ufl_cell(), integration_dim=integration_dim, - entity_ids=entity_ids, index_cache=index_cache, quadrature_rule=macro_quadrature_rule, scalar_type=parameters["scalar_type"]) @@ -572,18 +573,18 @@ def dg_injection_kernel(Vf, Vc, ncell): coarse_builder = firedrake_interface.KernelBuilder(info, parameters["scalar_type"]) coarse_builder.set_coordinates([Vc.mesh()]) + coarse_builder.set_entity_numbers([Vc.mesh()]) argument_multiindices = coarse_builder.argument_multiindices argument_multiindex, = argument_multiindices return_variable, = coarse_builder.return_variables - integration_dim, entity_ids = lower_integral_type(Vce.cell, "cell") + integration_dim, _ = lower_integral_type(Vce.cell, "cell") # Midpoint quadrature for jacobian on coarse cell. quadrature_rule = make_quadrature(Vce.cell, 0) coarse_cfg = dict(interface=coarse_builder, ufl_cell=Vc.ufl_cell(), integration_dim=integration_dim, - entity_ids=entity_ids, index_cache=index_cache, quadrature_rule=quadrature_rule, scalar_type=parameters["scalar_type"]) diff --git a/firedrake/pointeval_utils.py b/firedrake/pointeval_utils.py index 6f01a5c631..d5cf3ab60e 100644 --- a/firedrake/pointeval_utils.py +++ b/firedrake/pointeval_utils.py @@ -55,6 +55,7 @@ def compile_element(expression, coordinates, parameters=None): # Initialise kernel builder builder = firedrake_interface.KernelBuilderBase(utils.ScalarType) builder._domain_integral_type_map = {domain: "cell"} + builder._entity_ids = {domain: (0,)} builder.domain_coordinate[domain] = coordinates builder._coefficient(coordinates, "x") x_arg = builder.generate_arg_from_expression(builder.coefficient_map[coordinates]) diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 5b0a38b1b0..d8c0fe1764 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -144,6 +144,7 @@ def to_reference_coords_newton_step(ufl_coordinate_element, parameters, x0_dtype builder = firedrake_interface.KernelBuilderBase(ScalarType) builder._domain_integral_type_map = {domain: "cell"} + builder._entity_ids = {domain: (0,)} builder.domain_coordinate[domain] = C Cexpr = builder._coefficient(C, "C") x0_expr = builder._coefficient(x0, "x0") diff --git a/firedrake/slate/slac/compiler.py b/firedrake/slate/slac/compiler.py index 8751348b16..9f68648833 100644 --- a/firedrake/slate/slac/compiler.py +++ b/firedrake/slate/slac/compiler.py @@ -200,6 +200,7 @@ def generate_loopy_kernel(slate_expr, compiler_parameters=None): cell_sizes=(0, ) if builder.bag.needs_cell_sizes else (), exterior_facets=(), interior_facets=(), + orientations_cell=(), orientations_exterior_facet=(), orientations_interior_facet=(),), coefficient_numbers=coefficient_numbers, diff --git a/tests/firedrake/meshes/mixed_cell_unit_square.msh b/tests/firedrake/meshes/mixed_cell_unit_square.msh new file mode 100644 index 0000000000..582012ce11 --- /dev/null +++ b/tests/firedrake/meshes/mixed_cell_unit_square.msh @@ -0,0 +1,231 @@ +$MeshFormat +2.200000 0 8 +$EndMeshFormat +$Nodes +121 +1 0.000000 0.000000 0.000000 +2 1.000000 0.000000 0.000000 +3 0.000000 1.000000 0.000000 +4 1.000000 1.000000 0.000000 +5 0.100000 0.000000 0.000000 +6 0.200000 0.000000 0.000000 +7 0.300000 0.000000 0.000000 +8 0.400000 0.000000 0.000000 +9 0.500000 0.000000 0.000000 +10 0.600000 0.000000 0.000000 +11 0.700000 0.000000 0.000000 +12 0.800000 0.000000 0.000000 +13 0.900000 0.000000 0.000000 +14 0.000000 0.100000 0.000000 +15 0.000000 0.200000 0.000000 +16 0.000000 0.300000 0.000000 +17 0.000000 0.400000 0.000000 +18 0.000000 0.500000 0.000000 +19 0.000000 0.600000 0.000000 +20 0.000000 0.700000 0.000000 +21 0.000000 0.800000 0.000000 +22 0.000000 0.900000 0.000000 +23 1.000000 0.100000 0.000000 +24 1.000000 0.200000 0.000000 +25 1.000000 0.300000 0.000000 +26 1.000000 0.400000 0.000000 +27 1.000000 0.500000 0.000000 +28 1.000000 0.600000 0.000000 +29 1.000000 0.700000 0.000000 +30 1.000000 0.800000 0.000000 +31 1.000000 0.900000 0.000000 +32 0.900000 1.000000 0.000000 +33 0.800000 1.000000 0.000000 +34 0.700000 1.000000 0.000000 +35 0.600000 1.000000 0.000000 +36 0.500000 1.000000 0.000000 +37 0.400000 1.000000 0.000000 +38 0.300000 1.000000 0.000000 +39 0.200000 1.000000 0.000000 +40 0.100000 1.000000 0.000000 +41 0.099910 0.099602 0.000000 +42 0.199647 0.099017 0.000000 +43 0.298763 0.097861 0.000000 +44 0.397030 0.096831 0.000000 +45 0.494086 0.096232 0.000000 +46 0.592555 0.098553 0.000000 +47 0.692780 0.101000 0.000000 +48 0.795403 0.102370 0.000000 +49 0.897918 0.101440 0.000000 +50 0.099976 0.199139 0.000000 +51 0.100276 0.298871 0.000000 +52 0.100754 0.398576 0.000000 +53 0.100758 0.498986 0.000000 +54 0.100564 0.599366 0.000000 +55 0.100269 0.699723 0.000000 +56 0.100125 0.799876 0.000000 +57 0.100042 0.899958 0.000000 +58 0.895799 0.203814 0.000000 +59 0.895030 0.307167 0.000000 +60 0.899822 0.405711 0.000000 +61 0.901364 0.503314 0.000000 +62 0.901355 0.601763 0.000000 +63 0.900674 0.700631 0.000000 +64 0.900353 0.800260 0.000000 +65 0.900126 0.900083 0.000000 +66 0.800284 0.900132 0.000000 +67 0.700338 0.900083 0.000000 +68 0.600524 0.900015 0.000000 +69 0.500444 0.899908 0.000000 +70 0.400365 0.899842 0.000000 +71 0.300156 0.899889 0.000000 +72 0.200102 0.899912 0.000000 +73 0.394894 0.192601 0.000000 +74 0.298413 0.195663 0.000000 +75 0.487192 0.189602 0.000000 +76 0.582637 0.196271 0.000000 +77 0.682960 0.202753 0.000000 +78 0.789894 0.207261 0.000000 +79 0.783414 0.322202 0.000000 +80 0.199715 0.198087 0.000000 +81 0.200548 0.296574 0.000000 +82 0.201983 0.396430 0.000000 +83 0.202035 0.497478 0.000000 +84 0.201235 0.598747 0.000000 +85 0.200629 0.699407 0.000000 +86 0.200213 0.799802 0.000000 +87 0.804703 0.506541 0.000000 +88 0.803135 0.412674 0.000000 +89 0.803344 0.603044 0.000000 +90 0.801425 0.700957 0.000000 +91 0.800669 0.800362 0.000000 +92 0.701315 0.800417 0.000000 +93 0.601314 0.800070 0.000000 +94 0.501275 0.799737 0.000000 +95 0.400862 0.799573 0.000000 +96 0.300662 0.799553 0.000000 +97 0.477928 0.273171 0.000000 +98 0.396177 0.285892 0.000000 +99 0.665542 0.306018 0.000000 +100 0.712544 0.509728 0.000000 +101 0.721484 0.426103 0.000000 +102 0.300142 0.292896 0.000000 +103 0.305573 0.390645 0.000000 +104 0.305015 0.494723 0.000000 +105 0.303243 0.597330 0.000000 +106 0.301423 0.698947 0.000000 +107 0.702926 0.701156 0.000000 +108 0.706906 0.603770 0.000000 +109 0.604087 0.700387 0.000000 +110 0.402800 0.698655 0.000000 +111 0.503558 0.699294 0.000000 +112 0.565943 0.292189 0.000000 +113 0.617065 0.502789 0.000000 +114 0.633291 0.406568 0.000000 +115 0.608525 0.601090 0.000000 +116 0.534001 0.392869 0.000000 +117 0.416004 0.377213 0.000000 +118 0.410393 0.491619 0.000000 +119 0.405409 0.596818 0.000000 +120 0.508620 0.598318 0.000000 +121 0.517096 0.496276 0.000000 +$EndNodes +$Elements +101 +1 3 2 1 1 1 5 41 14 +2 3 2 1 1 5 6 42 41 +3 3 2 1 1 6 7 43 42 +4 3 2 1 1 7 8 44 43 +5 3 2 1 1 8 9 45 44 +6 3 2 1 1 9 10 46 45 +7 3 2 1 1 10 11 47 46 +8 3 2 1 1 11 12 48 47 +9 3 2 1 1 12 13 49 48 +10 3 2 1 1 13 2 23 49 +11 3 2 1 1 15 14 41 50 +12 3 2 1 1 16 15 50 51 +13 3 2 1 1 17 16 51 52 +14 3 2 1 1 18 17 52 53 +15 3 2 1 1 19 18 53 54 +16 3 2 1 1 20 19 54 55 +17 3 2 1 1 21 20 55 56 +18 3 2 1 1 22 21 56 57 +19 3 2 1 1 3 22 57 40 +20 3 2 1 1 23 24 58 49 +21 3 2 1 1 24 25 59 58 +22 3 2 1 1 25 26 60 59 +23 3 2 1 1 26 27 61 60 +24 3 2 1 1 27 28 62 61 +25 3 2 1 1 28 29 63 62 +26 3 2 1 1 29 30 64 63 +27 3 2 1 1 30 31 65 64 +28 3 2 1 1 31 4 32 65 +29 3 2 1 1 32 33 66 65 +30 3 2 1 1 33 34 67 66 +31 3 2 1 1 34 35 68 67 +32 3 2 1 1 35 36 69 68 +33 3 2 1 1 36 37 70 69 +34 3 2 1 1 37 38 71 70 +35 3 2 1 1 38 39 72 71 +36 3 2 1 1 39 40 57 72 +37 3 2 1 1 43 44 73 74 +38 3 2 1 1 44 45 75 73 +39 3 2 1 1 45 46 76 75 +40 3 2 1 1 46 47 77 76 +41 3 2 1 1 47 48 78 77 +42 3 2 1 1 48 49 58 78 +43 3 2 1 1 58 59 79 78 +44 3 2 1 1 41 42 80 50 +45 3 2 1 1 51 50 80 81 +46 3 2 1 1 52 51 81 82 +47 3 2 1 1 53 52 82 83 +48 3 2 1 1 54 53 83 84 +49 3 2 1 1 55 54 84 85 +50 3 2 1 1 56 55 85 86 +51 3 2 1 1 57 56 86 72 +52 3 2 1 1 60 61 87 88 +53 3 2 1 1 61 62 89 87 +54 3 2 1 1 62 63 90 89 +55 3 2 1 1 63 64 91 90 +56 3 2 1 1 64 65 66 91 +57 3 2 1 1 66 67 92 91 +58 3 2 1 1 67 68 93 92 +59 3 2 1 1 68 69 94 93 +60 3 2 1 1 69 70 95 94 +61 3 2 1 1 70 71 96 95 +62 3 2 1 1 71 72 86 96 +63 3 2 1 1 59 60 88 79 +64 3 2 1 1 42 43 74 80 +65 3 2 1 1 73 75 97 98 +66 3 2 1 1 77 78 79 99 +67 3 2 1 1 88 87 100 101 +68 3 2 1 1 74 73 98 102 +69 3 2 1 1 81 80 74 102 +70 3 2 1 1 82 81 102 103 +71 3 2 1 1 83 82 103 104 +72 3 2 1 1 84 83 104 105 +73 3 2 1 1 85 84 105 106 +74 3 2 1 1 86 85 106 96 +75 3 2 1 1 90 91 92 107 +76 3 2 1 1 87 89 108 100 +77 3 2 1 1 92 93 109 107 +78 3 2 1 1 94 95 110 111 +79 3 2 1 1 93 94 111 109 +80 3 2 1 1 95 96 106 110 +81 2 2 1 1 79 88 101 +82 3 2 1 1 89 90 107 108 +83 3 2 1 1 75 76 112 97 +84 3 2 1 1 76 77 99 112 +85 3 2 1 1 101 100 113 114 +86 3 2 1 1 107 109 115 108 +87 3 2 1 1 112 99 114 116 +88 3 2 1 1 104 103 117 118 +89 3 2 1 1 105 104 118 119 +90 3 2 1 1 106 105 119 110 +91 3 2 1 1 97 112 116 117 +92 3 2 1 1 102 98 117 103 +93 3 2 1 1 111 110 119 120 +94 3 2 1 1 100 108 115 113 +95 3 2 1 1 109 111 120 115 +96 3 2 1 1 114 113 121 116 +97 3 2 1 1 119 118 121 120 +98 3 2 1 1 117 116 121 118 +99 3 2 1 1 115 120 121 113 +100 2 2 1 1 97 117 98 +101 3 2 1 1 101 114 99 79 +$EndElements diff --git a/tests/firedrake/submesh/test_submesh_assemble.py b/tests/firedrake/submesh/test_submesh_assemble.py index 7a3e0136e3..28f244581f 100644 --- a/tests/firedrake/submesh/test_submesh_assemble.py +++ b/tests/firedrake/submesh/test_submesh_assemble.py @@ -1,5 +1,12 @@ +import os +import pytest import numpy as np from firedrake import * +from firedrake.cython import dmcommon +from petsc4py import PETSc + + +cwd = os.path.abspath(os.path.dirname(__file__)) def test_submesh_assemble_cell_cell_integral_cell(): @@ -326,3 +333,202 @@ def test_submesh_assemble_cell_cell_equation_bc(): [- 1. / 3., - 1. / 6., 0., 0.]]) assert np.allclose(A.M[0][0].values, M00) assert np.allclose(A.M[0][1].values, M01) + + +def test_submesh_assemble_cell_facet_integral_various(): + # CG1 DoF numbers (nprocs = 1): + # + # 5-------1-------2 + # | | | + # | | | mesh + # | | | + # 4-------0-------3 + # + # 0 + # | + # | subm + # | + # 1 + # + distribution_parameters = { + "overlap_type": (DistributedMeshOverlapType.RIDGE, 1), + } + subdomain_id = 777 + mesh = RectangleMesh(2, 1, 2., 1., quadrilateral=True, distribution_parameters=distribution_parameters) + x, y = SpatialCoordinate(mesh) + V1 = FunctionSpace(mesh, "HDiv Trace", 0) + f1 = Function(V1).interpolate(conditional(And(x > 0.9, x < 1.1), 1., 0.)) + mesh = RelabeledMesh(mesh, [f1], [subdomain_id]) + x, y = SpatialCoordinate(mesh) + subm = Submesh(mesh, mesh.topological_dimension - 1, subdomain_id) + subx, suby = SpatialCoordinate(subm) + V0 = FunctionSpace(mesh, "CG", 1) + V1 = FunctionSpace(subm, "CG", 1) + V = V0 * V1 + u = TrialFunction(V) + v = TestFunction(V) + u0, u1 = split(u) + v0, v1 = split(v) + coordV0 = VectorFunctionSpace(mesh, "CG", 1) + coordV1 = VectorFunctionSpace(subm, "CG", 1) + coordV = coordV0 * coordV1 + coords = Function(coordV) + 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.], + ] + ) + # Use subm as primal integration domain. + measure = Measure( + "dx", subm, + intersect_measures=( + 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, + intersect_measures=( + Measure("dx", subm), + ), + ) + a = inner(u0('+'), v1) * measure(subdomain_id) + A = assemble(a, mat_type="nest") + assert np.allclose(A.M[1][0].values, M10) + a = inner(u1, v0('-')) * measure(subdomain_id) + A = assemble(a, mat_type="nest") + assert np.allclose(A.M[0][1].values, np.transpose(M10)) + + +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_assemble_quad_triangle_base(): + dim = 2 + label_ext = 1 + label_interf = 2 + mesh = Mesh(os.path.join(cwd, "..", "meshes", "mixed_cell_unit_square.msh")) + mesh.topology_dm.markBoundaryFaces(dmcommon.FACE_SETS_LABEL, label_ext) + mesh_t = Submesh(mesh, dim, PETSc.DM.PolytopeType.TRIANGLE, label_name="celltype", name="mesh_tri") + x_t, y_t = SpatialCoordinate(mesh_t) + n_t = FacetNormal(mesh_t) + mesh_q = Submesh(mesh, dim, PETSc.DM.PolytopeType.QUADRILATERAL, label_name="celltype", name="mesh_quad") + x_q, y_q = SpatialCoordinate(mesh_q) + n_q = FacetNormal(mesh_q) + # pgfplot(f, "mesh_tri.dat", degree=2) + dx_t = Measure("dx", mesh_t) + dx_q = Measure("dx", mesh_q) + ds_t = Measure("ds", mesh_t, intersect_measures=(Measure("ds", mesh_q),)) + ds_q = Measure("ds", mesh_q, intersect_measures=(Measure("ds", mesh_t),)) + A_t = assemble(Constant(1) * dx_t) + A_q = assemble(Constant(1) * dx_q) + assert abs(A_t + A_q - 1.0) < 1.e-13 + HDiv_t = FunctionSpace(mesh_t, "BDM", 3) + HDiv_q = FunctionSpace(mesh_q, "RTCF", 3) + hdiv_t = Function(HDiv_t).interpolate(as_vector([x_t**2, y_t**2])) + hdiv_q = Function(HDiv_q).project(as_vector([x_q**2, y_q**2]), solver_parameters={"ksp_rtol": 1.e-13}) + v_t = assemble(dot(hdiv_q, as_vector([x_q, y_q])) * ds_t(label_interf)) + v_q = assemble(dot(hdiv_t, as_vector([x_t, y_t])) * ds_q(label_interf)) + assert abs(v_q - v_t) < 1.e-13 + v_t = assemble(dot(hdiv_q, as_vector([x_t, y_t])) * ds_t(label_interf)) + v_q = assemble(dot(hdiv_t, as_vector([x_q, y_q])) * ds_q(label_interf)) + assert abs(v_q - v_t) < 1.e-13 + v_t = assemble(dot(hdiv_q, as_vector([x_q, y_t])) * ds_t(label_interf)) + v_q = assemble(dot(hdiv_t, as_vector([x_t, y_q])) * ds_q(label_interf)) + assert abs(v_q - v_t) < 1.e-13 + v = assemble(inner(n_t, as_vector([888., 999.])) * ds_t(label_interf)) + assert abs(v) < 1.e-13 + v = assemble(inner(n_q, as_vector([888., 999.])) * ds_q(label_interf)) + assert abs(v) < 1.e-13 + v = assemble(inner(n_q, as_vector([888., 999.])) * ds_t(label_interf)) + assert abs(v) < 1.e-13 + v = assemble(inner(n_t, as_vector([888., 999.])) * ds_q(label_interf)) + assert abs(v) < 1.e-13 + v = assemble(dot(n_q + n_t, n_q + n_t) * ds_t(label_interf)) + assert abs(v) < 1.e-30 + v = assemble(dot(n_q + n_t, n_q + n_t) * ds_q(label_interf)) + assert abs(v) < 1.e-30 + + +def test_submesh_assemble_quad_triangle(): + dim = 2 + label_ext = 1 + label_interf = 2 + mesh = Mesh(os.path.join(cwd, "..", "meshes", "mixed_cell_unit_square.msh")) + mesh.topology_dm.markBoundaryFaces(dmcommon.FACE_SETS_LABEL, label_ext) + mesh_t = Submesh(mesh, dim, PETSc.DM.PolytopeType.TRIANGLE, label_name="celltype", name="mesh_tri") + x_t, y_t = SpatialCoordinate(mesh_t) + n_t = FacetNormal(mesh_t) + mesh_q = Submesh(mesh, dim, PETSc.DM.PolytopeType.QUADRILATERAL, label_name="celltype", name="mesh_quad") + x_q, y_q = SpatialCoordinate(mesh_q) + n_q = FacetNormal(mesh_q) + V_t = FunctionSpace(mesh_t, "P", 4) + V_q = FunctionSpace(mesh_q, "Q", 3) + V = V_t * V_q + u = TrialFunction(V) + v = TestFunction(V) + u_t, u_q = split(u) + v_t, v_q = split(v) + ds_t = Measure("ds", mesh_t, intersect_measures=(Measure("ds", mesh_q),)) + ds_q = Measure("ds", mesh_q, intersect_measures=(Measure("ds", mesh_t),)) + # Test against the base cases. + c = x_t**2 * y_t**2 + a = c * inner(u_t, v_q) * ds_t(label_interf) + A = assemble(a) + c_ref = x_q**2 * y_q**2 + a_ref = c_ref * inner(TrialFunction(V_t), TestFunction(V_q)) * ds_t(label_interf) + A_ref = assemble(a_ref) + assert np.allclose(A.M[1][0].values, A_ref.M.values) + c = x_t**2 * y_q**2 + a = c * inner(u_q, v_t) * ds_t(label_interf) + A = assemble(a) + c_ref = x_q**2 * y_t**2 + a_ref = c_ref * inner(TrialFunction(V_q), TestFunction(V_t)) * ds_t(label_interf) + A_ref = assemble(a_ref) + assert np.allclose(A.M[0][1].values, A_ref.M.values) + c = dot(n_t, n_t) + a = c * inner(u_t, v_q) * ds_q(label_interf) + A = assemble(a) + c_ref = dot(n_q, n_q) + a_ref = c_ref * inner(TrialFunction(V_t), TestFunction(V_q)) * ds_q(label_interf) + A_ref = assemble(a_ref) + assert np.allclose(A.M[1][0].values, A_ref.M.values) + c = dot(n_t, n_q) + a = c * inner(u_q, v_t) * ds_q(label_interf) + A = assemble(a) + c_ref = dot(n_q, n_t) + a_ref = c_ref * inner(TrialFunction(V_q), TestFunction(V_t)) * ds_q(label_interf) + A_ref = assemble(a_ref) + assert np.allclose(A.M[0][1].values, A_ref.M.values) diff --git a/tests/firedrake/submesh/test_submesh_base.py b/tests/firedrake/submesh/test_submesh_base.py index 5c4a0bc61b..db877ff23c 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.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() - composed_map, integral_type = mesh.topology.trans_mesh_entity_map(submesh.topology, "cell", None, None) + composed_map, integral_type = mesh.topology.trans_mesh_entity_map(submesh.topology, "cell", "everywhere", None) assert integral_type == "cell" if rank == 0: assert (composed_map.maps_[0].values_with_halo == np.array([0])).all() diff --git a/tests/firedrake/submesh/test_submesh_solve.py b/tests/firedrake/submesh/test_submesh_solve.py index 0691646748..aa5e6725de 100644 --- a/tests/firedrake/submesh/test_submesh_solve.py +++ b/tests/firedrake/submesh/test_submesh_solve.py @@ -1,7 +1,10 @@ +import os import pytest from os.path import abspath, dirname, join import numpy as np from firedrake import * +from firedrake.cython import dmcommon +from petsc4py import PETSc cwd = abspath(dirname(__file__)) @@ -458,3 +461,293 @@ def test_submesh_solve_cell_cell_equation_bc(nref, degree, simplex): solve(a == L, sol, bcs=[dbc, ebc]) assert sqrt(assemble(inner(sol[0] - x * y, sol[0] - x * y) * dx_outer)) < 1.e-12 assert sqrt(assemble(inner(sol[1] - x * y, sol[1] - x * y) * dx_inner)) < 1.e-12 + + +def _test_submesh_solve_quad_triangle_poisson(nref, degree): + dim = 2 + label_ext = 1 + label_interf = 2 + distribution_parameters_noop = { + "partition": True, + "overlap_type": (DistributedMeshOverlapType.NONE, 0), + } + mesh = Mesh(os.path.join(cwd, "..", "meshes", "mixed_cell_unit_square.msh"), distribution_parameters=distribution_parameters_noop) + plex = mesh.topology_dm + for _ in range(nref): + plex = plex.refine() + plex.removeLabel("pyop2_core") + plex.removeLabel("pyop2_owned") + plex.removeLabel("pyop2_ghost") + mesh = Mesh(plex) + h = 0.1 / 2**nref # roughly + mesh.topology_dm.markBoundaryFaces(dmcommon.FACE_SETS_LABEL, label_ext) + mesh_t = Submesh(mesh, dim, PETSc.DM.PolytopeType.TRIANGLE, label_name="celltype", name="mesh_tri") + x_t, y_t = SpatialCoordinate(mesh_t) + n_t = FacetNormal(mesh_t) + mesh_q = Submesh(mesh, dim, PETSc.DM.PolytopeType.QUADRILATERAL, label_name="celltype", name="mesh_quad") + x_q, y_q = SpatialCoordinate(mesh_q) + n_q = FacetNormal(mesh_q) + V_t = FunctionSpace(mesh_t, "P", degree) + V_q = FunctionSpace(mesh_q, "Q", degree) + V = V_t * V_q + u = TrialFunction(V) + v = TestFunction(V) + u_t, u_q = split(u) + v_t, v_q = split(v) + dx_t = Measure("dx", mesh_t) + dx_q = Measure("dx", mesh_q) + ds_t = Measure("ds", mesh_t, intersect_measures=(Measure("ds", mesh_q),)) + ds_q = Measure("ds", mesh_q, intersect_measures=(Measure("ds", mesh_t),)) + g_t = cos(2 * pi * x_t) * cos(2 * pi * y_t) + g_q = cos(2 * pi * x_q) * cos(2 * pi * y_q) + f_t = 8 * pi**2 * g_t + f_q = 8 * pi**2 * g_q + a = ( + inner(grad(u_t), grad(v_t)) * dx_t + inner(grad(u_q), grad(v_q)) * dx_q + - inner( + (grad(u_q) + grad(u_t)) / 2, + (v_q * n_q + v_t * n_t) + ) * ds_q(label_interf) + - inner( + (u_q * n_q + u_t * n_t), + (grad(v_q) + grad(v_t)) / 2 + ) * ds_t(label_interf) + + 100 / h * inner(u_q - u_t, v_q - v_t) * ds_q(label_interf) + ) + L = ( + inner(f_t, v_t) * dx_t + inner(f_q, v_q) * dx_q + ) + sol = Function(V) + bc_q = DirichletBC(V.sub(1), g_q, label_ext) + solve(a == L, sol, bcs=[bc_q]) + sol_t, sol_q = split(sol) + L2Error_t = assemble(inner(sol_t - g_t, sol_t - g_t) * dx_t) + L2Error_q = assemble(inner(sol_q - g_q, sol_q - g_q) * dx_q) + H1Error_t = L2Error_t + assemble(inner(grad(sol_t - g_t), grad(sol_t - g_t)) * dx_t) + H1Error_q = L2Error_q + assemble(inner(grad(sol_q - g_q), grad(sol_q - g_q)) * dx_q) + return sqrt(L2Error_t + L2Error_q), sqrt(H1Error_t + H1Error_q) + + +@pytest.mark.parallel(nprocs=8) +def test_submesh_solve_quad_triangle_poisson_convergence(): + for degree in range(1, 5): + L2Errors = [] + H1Errors = [] + for nref in range(4): + L2Error, H1Error = _test_submesh_solve_quad_triangle_poisson(nref, degree) + L2Errors.append(L2Error) + H1Errors.append(H1Error) + L2Errors = [np.log2(c) - np.log2(f) for c, f in zip(L2Errors[:-1], L2Errors[1:])] + H1Errors = [np.log2(c) - np.log2(f) for c, f in zip(H1Errors[:-1], H1Errors[1:])] + assert (np.array(L2Errors) > (degree + 1) * 0.995).all() + assert (np.array(H1Errors) > (degree) * 0.995).all() + + +def _test_submesh_solve_3d_2d_poisson(simplex, direction, nref, degree): + distribution_parameters_noop = { + "partition": True, + "overlap_type": (DistributedMeshOverlapType.NONE, 0), + } + distribution_parameters = { + "overlap_type": (DistributedMeshOverlapType.RIDGE, 1), + } + dim = 3 + interf_at = 0.499 + if simplex: + nref_simplex = 3 + mesh = BoxMesh(2 ** nref_simplex, 2 ** nref_simplex, 2 ** nref_simplex, 1., 1., 1., hexahedral=False, distribution_parameters=distribution_parameters_noop) + xyz = SpatialCoordinate(mesh) + DG0 = FunctionSpace(mesh, "DG", 0) + c1 = Function(DG0).interpolate(conditional(xyz[direction] < interf_at, 1, 0)) + c2 = Function(DG0).interpolate(conditional(xyz[direction] > interf_at, 1, 0)) + mesh = RelabeledMesh(mesh, [c1, c2], [1, 2]) + family = "P" + else: + mesh = Mesh(join(cwd, "..", "meshes", "cube_hex.msh"), distribution_parameters=distribution_parameters_noop) + xyz = SpatialCoordinate(mesh) + DG0 = FunctionSpace(mesh, "DQ", 0) + c1 = Function(DG0).interpolate(conditional(xyz[direction] < interf_at, 1, 0)) + c2 = Function(DG0).interpolate(conditional(xyz[direction] > interf_at, 1, 0)) + HDivTrace0 = FunctionSpace(mesh, "Q", 2) + f1 = Function(HDivTrace0).interpolate(conditional(xyz[0] < .001, 1, 0)) + f2 = Function(HDivTrace0).interpolate(conditional(xyz[0] > .999, 1, 0)) + f3 = Function(HDivTrace0).interpolate(conditional(xyz[1] < .001, 1, 0)) + f4 = Function(HDivTrace0).interpolate(conditional(xyz[1] > .999, 1, 0)) + f5 = Function(HDivTrace0).interpolate(conditional(xyz[2] < .001, 1, 0)) + f6 = Function(HDivTrace0).interpolate(conditional(xyz[2] > .999, 1, 0)) + mesh = RelabeledMesh(mesh, [c1, c2, f1, f2, f3, f4, f5, f6], [1, 2, 1, 2, 3, 4, 5, 6]) + family = "Q" + plex = mesh.topology_dm + for _ in range(nref): + plex = plex.refine() + plex.removeLabel("pyop2_core") + plex.removeLabel("pyop2_owned") + plex.removeLabel("pyop2_ghost") + mesh = Mesh(plex, distribution_parameters=distribution_parameters) + mesh1 = Submesh(mesh, dim, 1) + x1, y1, z1 = SpatialCoordinate(mesh1) + mesh2 = Submesh(mesh, dim, 2) + x2, y2, z2 = SpatialCoordinate(mesh2) + label_interf = 7 # max + 1 + mesh12 = Submesh(mesh2, dim - 1, label_interf) + dx1 = Measure("dx", mesh1) + dx2 = Measure("dx", mesh2) + ds1_ds2 = Measure("ds", mesh1, intersect_measures=(Measure("ds", mesh2),)) + dx12_ds1_ds2 = Measure( + "dx", mesh12, + intersect_measures=( + Measure("ds", mesh1), + Measure("ds", mesh2), + ) + ) + # Check sanity. + vol1 = assemble(Constant(1) * dx1) + vol2 = assemble(Constant(1) * dx2) + assert abs(vol1 + vol2 - 1.) < 1.e-13 + # Solve Poisson problem. + V1 = FunctionSpace(mesh1, family, degree) + V12 = FunctionSpace(mesh12, family, degree) + V2 = FunctionSpace(mesh2, family, degree) + V = V1 * V12 * V2 + u = TrialFunction(V) + v = TestFunction(V) + u1, u12, u2 = split(u) + v1, v12, v2 = split(v) + g1 = cos(2 * pi * x1) * cos(2 * pi * y1) * cos(2 * pi * z1) + g2 = cos(2 * pi * x2) * cos(2 * pi * y2) * cos(2 * pi * z2) + f1 = 12 * pi**2 * g1 + f2 = 12 * pi**2 * g2 + n1 = FacetNormal(mesh1) + n2 = FacetNormal(mesh2) + h = 0.1 / 2**nref # roughly + a = ( + inner(grad(u1), grad(v1)) * dx1 + inner(grad(u2), grad(v2)) * dx2 + - inner( + u12, + (v1 - v2) + ) * dx12_ds1_ds2 + - inner( + (u1 * n1 + u2 * n2), + (grad(v1) + grad(v2)) / 2 + ) * dx12_ds1_ds2 + + 100 / h * inner(u1 - u2, v1 - v2) * ds1_ds2(label_interf) # Can also use dx12_ds1_ds2. + + inner( + (dot(grad(u1), n1) - dot(grad(u2), n2)) / 2 - u12, + v12 + ) * dx12_ds1_ds2 + ) + L = ( + inner(f1, v1) * dx1 + inner(f2, v2) * dx2 + ) + sol = Function(V) + bc1 = DirichletBC(V.sub(0), g1, [i for i in range(1, 7) if i != 2 * direction + 2]) + bc2 = DirichletBC(V.sub(2), g2, [i for i in range(1, 7) if i != 2 * direction + 1]) + solver_parameters = { + "mat_type": "matfree", + "ksp_type": "preonly", + "pc_type": "fieldsplit", + "pc_fieldsplit_type": "schur", + "pc_fieldsplit_schur_fact_type": "full", + "pc_fieldsplit_0_fields": "1", + "pc_fieldsplit_1_fields": "0, 2", + "fieldsplit_0_ksp_type": "cg", + "fieldsplit_0_ksp_rtol": 1e-14, + "fieldsplit_0_pc_type": "jacobi", + "fieldsplit_1_ksp_type": "cg", + "fieldsplit_1_ksp_rtol": 1e-14, + "fieldsplit_1_pc_type": "jacobi", + } + solve(a == L, sol, bcs=[bc1, bc2], solver_parameters=solver_parameters) + sol1, sol12, sol2 = split(sol) + L2Error1 = assemble(inner(sol1 - g1, sol1 - g1) * dx1) + L2Error2 = assemble(inner(sol2 - g2, sol2 - g2) * dx2) + H1Error1 = L2Error1 + assemble(inner(grad(sol1 - g1), grad(sol1 - g1)) * dx1) + H1Error2 = L2Error2 + assemble(inner(grad(sol2 - g2), grad(sol2 - g2)) * dx2) + return sqrt(L2Error1 + L2Error2), sqrt(H1Error1 + H1Error2) + + +@pytest.mark.parallel(nprocs=6) +@pytest.mark.parametrize('simplex', [True, False]) +@pytest.mark.parametrize('direction', [0, 1, 2]) +def test_submesh_solve_3d_2d_poisson_sanity(simplex, direction): + nref = 0 + degree = 4 + L2Error, H1Error = _test_submesh_solve_3d_2d_poisson(simplex, direction, nref, degree) + assert L2Error < 6.e-5 + assert H1Error < 5.e-3 + + +@pytest.mark.parallel(nprocs=8) +@pytest.mark.parametrize('simplex', [False]) +@pytest.mark.parametrize('direction', [0]) +@pytest.mark.parametrize('degree', [3]) +def test_submesh_solve_3d_2d_poisson_convergence(simplex, direction, degree): + L2Errors = [] + H1Errors = [] + for nref in range(2): + L2Error, H1Error = _test_submesh_solve_3d_2d_poisson(simplex, direction, nref, degree) + L2Errors.append(L2Error) + H1Errors.append(H1Error) + L2Errors = [np.log2(c) - np.log2(f) for c, f in zip(L2Errors[:-1], L2Errors[1:])] + H1Errors = [np.log2(c) - np.log2(f) for c, f in zip(H1Errors[:-1], H1Errors[1:])] + assert (np.array(L2Errors) > (degree + 1) * 0.96).all() + assert (np.array(H1Errors) > (degree) * 0.96).all() + + +@pytest.mark.parallel(nprocs=7) +def test_submesh_solve_2d_1d_poisson_hermite(): + distribution_parameters_noop = { + "partition": True, + "overlap_type": (DistributedMeshOverlapType.NONE, 0), + } + distribution_parameters = { + "overlap_type": (DistributedMeshOverlapType.RIDGE, 1), + } + mesh3d = Mesh(join(cwd, "..", "meshes", "cube_hex.msh"), distribution_parameters=distribution_parameters_noop) + plex = mesh3d.topology_dm + for _ in range(2): + plex = plex.refine() + plex.removeLabel("pyop2_core") + plex.removeLabel("pyop2_owned") + plex.removeLabel("pyop2_ghost") + mesh3d = Mesh(plex, distribution_parameters=distribution_parameters) + xyz = SpatialCoordinate(mesh3d) + HDivTrace0 = FunctionSpace(mesh3d, "Q", 2) + f1 = Function(HDivTrace0).interpolate(conditional(xyz[0] < .001, 1, 0)) + f2 = Function(HDivTrace0).interpolate(conditional(xyz[0] > .999, 1, 0)) + f3 = Function(HDivTrace0).interpolate(conditional(xyz[1] < .001, 1, 0)) + f4 = Function(HDivTrace0).interpolate(conditional(xyz[1] > .999, 1, 0)) + f5 = Function(HDivTrace0).interpolate(conditional(xyz[2] < .001, 1, 0)) + f6 = Function(HDivTrace0).interpolate(conditional(xyz[2] > .999, 1, 0)) + mesh3d = RelabeledMesh(mesh3d, [f1, f2, f3, f4, f5, f6], [1, 2, 3, 4, 5, 6]) + mesh2d = Submesh(mesh3d, mesh3d.topological_dimension - 1, 6) + mesh1d = Submesh(mesh2d, mesh2d.topological_dimension - 1, 4) + x2d = SpatialCoordinate(mesh2d) + x1d = SpatialCoordinate(mesh1d) + g2d = sin(2 * pi * x2d[0]) * sin(2 * pi * x2d[1]) + f2d = 2 * (2 * pi)**2 * g2d + g1d = sin(2 * pi * x1d[0]) + f1d = (2 * pi)**4 * g1d + V2d = FunctionSpace(mesh2d, "Q", 3) + V1d = FunctionSpace(mesh1d, "Hermite", 3) + V = V2d * V1d + u = TrialFunction(V) + v = TestFunction(V) + u2d, u1d = split(u) + v2d, v1d = split(v) + dx2d = Measure("dx", mesh2d) + dx1d = Measure("dx", mesh1d) + ds2d_dx1d = Measure("ds", mesh2d, intersect_measures=(Measure("dx", mesh1d),)) + a = inner(grad(u2d), grad(v2d)) * dx2d - \ + inner(2 * pi * u1d, v2d) * ds2d_dx1d(4) + \ + inner(grad(grad(u1d)), grad(grad(v1d))) * dx1d + L = inner(f2d, v2d) * dx2d + \ + inner(f1d, v1d) * dx1d + sol = Function(V) + bc2d = DirichletBC(V.sub(0), g2d, (1, 2, 3)) + bc1d = DirichletBC(V.sub(1), g1d, (1, 2)) + solve(a == L, sol, bcs=[bc2d, bc1d]) + error2d = assemble(inner(sol[0] - g2d, sol[0] - g2d) * dx2d) + error1d = assemble(inner(sol[1] - g1d, sol[1] - g1d) * dx1d) + assert sqrt(error2d) < 1.e-5 + assert sqrt(error1d) < 5.e-5 diff --git a/tsfc/driver.py b/tsfc/driver.py index eda810bdf7..c82a4eb5ed 100644 --- a/tsfc/driver.py +++ b/tsfc/driver.py @@ -16,6 +16,7 @@ import gem.impero_utils as impero_utils import finat +from finat.element_factory import as_fiat_cell from tsfc import fem, ufl_utils from tsfc.logging import logger @@ -247,6 +248,7 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, domain = extract_unique_domain(expression) assert domain is not None builder._domain_integral_type_map = {domain: "cell"} + builder._entity_ids = {domain: (0,)} # Collect required coefficients and determine numbering coefficients = extract_coefficients(expression) @@ -280,6 +282,7 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, # Set up kernel config for translation of UFL expression to gem kernel_cfg = dict(interface=builder, ufl_cell=domain.ufl_cell(), + integration_dim=as_fiat_cell(domain.ufl_cell()).get_dimension(), # FIXME: change if we ever implement # interpolation on facets. argument_multiindices=argument_multiindices, diff --git a/tsfc/fem.py b/tsfc/fem.py index fc468d8fc7..1c13a37a52 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -9,7 +9,7 @@ import numpy import ufl from FIAT.orientation_utils import Orientation as FIATOrientation -from FIAT.reference_element import UFCHexahedron, UFCSimplex, make_affine_mapping +from FIAT.reference_element import UFCHexahedron, UFCQuadrilateral, UFCSimplex, make_affine_mapping from FIAT.reference_element import TensorProductCell from finat.physically_mapped import (NeedsCoordinateMappingElement, PhysicalGeometry) @@ -35,6 +35,7 @@ from tsfc import ufl2gem from tsfc.kernel_interface import ProxyKernelInterface +from tsfc.kernel_interface.common import lower_integral_type from tsfc.modified_terminals import (analyse_modified_terminal, construct_modified_terminal) from tsfc.parameters import is_complex @@ -52,7 +53,6 @@ class ContextBase(ProxyKernelInterface): 'ufl_cell', 'fiat_cell', 'integration_dim', - 'entity_ids', 'argument_multiindices', 'facetarea', 'index_cache', @@ -74,9 +74,14 @@ def fiat_cell(self): @cached_property def integration_dim(self): - return self.fiat_cell.get_dimension() - - entity_ids = [0] + integration_dims = set() + for domain, integral_type in self.domain_integral_type_map.items(): + cell = domain.ufl_cell() + fiat_cell = as_fiat_cell(cell) + integration_dim, _ = lower_integral_type(fiat_cell, integral_type) + integration_dims.add(integration_dim) + integration_dim, = integration_dims + return integration_dim @cached_property def epsilon(self): @@ -97,11 +102,11 @@ def entity_selector(self, callback, domain, restriction): :arg restriction: Restriction of the modified terminal, used for entity selection. """ - if len(self.entity_ids) == 1: - return callback(self.entity_ids[0]) + if len(self.entity_ids(domain)) == 1: + return callback(self.entity_ids(domain)[0]) else: f = self.entity_number(domain, restriction) - return gem.select_expression(list(map(callback, self.entity_ids)), f) + return gem.select_expression(list(map(callback, self.entity_ids(domain))), f) argument_multiindices = () @@ -119,19 +124,52 @@ def use_canonical_quadrature_point_ordering(self): # Directly set use_canonical_quadrature_point_ordering = False in context # for translation of special nodes, e.g., CellVolume, FacetArea, CellOrigin, and CellVertices, # as quadrature point ordering is not relevant for those node types. - cell_integral_type_map = { - as_fiat_cell(domain.ufl_cell()): integral_type - for domain, integral_type in self.domain_integral_type_map.items() - if integral_type is not None - } - if all(integral_type == 'cell' for integral_type in cell_integral_type_map.values()): + def _any(cell_type, integral_types): + for d, it in self.domain_integral_type_map.items(): + if it is None: + continue + c = as_fiat_cell(d.ufl_cell()) + if isinstance(c, cell_type) and it in integral_types: + return True return False - elif all(integral_type in ['exterior_facet', 'interior_facet'] for integral_type in cell_integral_type_map.values()): - if all(isinstance(cell, UFCHexahedron) for cell in cell_integral_type_map): - return True - elif len(set(cell_integral_type_map)) > 1: # mixed cell types - return True - return False + if _any(UFCHexahedron, ['exterior_facet', 'interior_facet']): + return True + elif _any(UFCQuadrilateral, ['exterior_facet', 'interior_facet']) and _any(UFCSimplex, ['cell', 'exterior_facet', 'interior_facet']): + return True + else: + return False + + +class CellKernelInterface(ProxyKernelInterface): + # For a single-domain cell integration kernel. + + def __init__(self, wrapee, domain): + super().__init__(wrapee) + self._domain = domain + + def entity_ids(self, domain): + if domain is not self._domain: + raise ValueError(f"{domain} != {self._domain}") + return (0,) + + @cached_property + def domain_integral_type_map(self): + return {self._domain: "cell"} + + +class CellVolumeKernelInterface(CellKernelInterface): + # Since CellVolume is evaluated as a cell integral, we must ensure + # that the right restriction is applied when it is used in an + # interior facet integral. This proxy diverts coefficient + # translation to use a specified restriction. + + def __init__(self, wrapee, domain, restriction): + super().__init__(wrapee, domain) + self.restriction = restriction + + def coefficient(self, ufl_coefficient, r): + assert r is None + return self._wrapee.coefficient(ufl_coefficient, self.restriction) class CoordinateMapping(PhysicalGeometry): @@ -168,7 +206,7 @@ def config(self): config["interface"] = self.interface return config - def translate_point_expression(self, expr, point=None): + def translate_point_expression(self, expr, point=None, interface=None): if self.mt.restriction == '+': expr = PositiveRestricted(expr) elif self.mt.restriction == '-': @@ -181,6 +219,8 @@ def translate_point_expression(self, expr, point=None): config = {"point_set": PointSingleton(point)} config.update(self.config) config.update(use_canonical_quadrature_point_ordering=False) # quad point ordering not relevant. + if interface: + config.update(interface=interface) context = PointSetContext(**config) expr = self.preprocess(expr, context) return map_expr_dag(context.translator, expr) @@ -189,12 +229,16 @@ def cell_size(self): return self.interface.cell_size(extract_unique_domain(self.mt.terminal), self.mt.restriction) def jacobian_at(self, point): - expr = Jacobian(extract_unique_domain(self.mt.terminal)) - return self.translate_point_expression(expr, point=point) + domain = extract_unique_domain(self.mt.terminal) + expr = Jacobian(domain) + interface = CellKernelInterface(self.interface, domain) + return self.translate_point_expression(expr, point=point, interface=interface) def detJ_at(self, point): - expr = JacobianDeterminant(extract_unique_domain(self.mt.terminal)) - return self.translate_point_expression(expr, point=point) + domain = extract_unique_domain(self.mt.terminal) + expr = JacobianDeterminant(domain) + interface = CellKernelInterface(self.interface, domain) + return self.translate_point_expression(expr, point=point, interface=interface) def reference_normals(self): cell = self.interface.fiat_cell @@ -257,7 +301,7 @@ def physical_points(self, point_set, entity=None): config.update(self.config) if entity is not None: config.update({name: getattr(self.interface, name) - for name in ["integration_dim", "entity_ids"]}) + for name in ["integration_dim"]}) config.update(use_canonical_quadrature_point_ordering=False) # quad point ordering not relevant. context = PointSetContext(**config) expr = self.preprocess(expr, context) @@ -406,7 +450,7 @@ def facet_avg(self, o): config = {name: getattr(self.context, name) for name in ["ufl_cell", "index_cache", "scalar_type", - "integration_dim", "entity_ids"]} + "integration_dim"]} config.update(quadrature_degree=degree, interface=self.context, argument_multiindices=argument_multiindices) expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) @@ -489,10 +533,13 @@ def make_cell_facet_jacobian(cell, facet_dim, facet_i): @translate.register(ReferenceNormal) def translate_reference_normal(terminal, mt, ctx): + domain = extract_unique_domain(terminal) + fiat_cell = as_fiat_cell(domain.ufl_cell()) + def callback(facet_i): - n = ctx.fiat_cell.compute_reference_normal(ctx.integration_dim, facet_i) + n = fiat_cell.compute_reference_normal(ctx.integration_dim, facet_i) return gem.Literal(n) - return ctx.entity_selector(callback, extract_unique_domain(terminal), mt.restriction) + return ctx.entity_selector(callback, domain, mt.restriction) @translate.register(ReferenceCellEdgeVectors) @@ -547,29 +594,19 @@ def translate_spatialcoordinate(terminal, mt, ctx): return ctx.translator(expr) -class CellVolumeKernelInterface(ProxyKernelInterface): - # Since CellVolume is evaluated as a cell integral, we must ensure - # that the right restriction is applied when it is used in an - # interior facet integral. This proxy diverts coefficient - # translation to use a specified restriction. - - def __init__(self, wrapee, restriction): - ProxyKernelInterface.__init__(self, wrapee) - self.restriction = restriction - - def coefficient(self, ufl_coefficient, r): - assert r is None - return self._wrapee.coefficient(ufl_coefficient, self.restriction) - - @translate.register(CellVolume) def translate_cellvolume(terminal, mt, ctx): - integrand, degree = one_times(ufl.dx(domain=extract_unique_domain(terminal))) - interface = CellVolumeKernelInterface(ctx, mt.restriction) + domain = extract_unique_domain(terminal) + integrand, degree = one_times(ufl.dx(domain=domain)) + interface = CellVolumeKernelInterface(ctx, domain, mt.restriction) config = {name: getattr(ctx, name) for name in ["ufl_cell", "index_cache", "scalar_type"]} - config.update(interface=interface, quadrature_degree=degree, use_canonical_quadrature_point_ordering=False) + config.update( + interface=interface, + quadrature_degree=degree, + use_canonical_quadrature_point_ordering=False, + ) expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) return expr @@ -583,7 +620,7 @@ def translate_facetarea(terminal, mt, ctx): config = {name: getattr(ctx, name) for name in ["ufl_cell", "integration_dim", "scalar_type", - "entity_ids", "index_cache"]} + "index_cache"]} config.update(interface=ctx, quadrature_degree=degree, use_canonical_quadrature_point_ordering=False) expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) return expr @@ -605,13 +642,18 @@ def translate_cellorigin(terminal, mt, ctx): @translate.register(CellVertices) def translate_cell_vertices(terminal, mt, ctx): - coords = SpatialCoordinate(extract_unique_domain(terminal)) + domain = extract_unique_domain(terminal) + coords = SpatialCoordinate(domain) ufl_expr = construct_modified_terminal(mt, coords) ps = PointSet(numpy.array(ctx.fiat_cell.get_vertices())) - + interface = CellKernelInterface(ctx, domain) config = {name: getattr(ctx, name) for name in ["ufl_cell", "index_cache", "scalar_type"]} - config.update(interface=ctx, point_set=ps, use_canonical_quadrature_point_ordering=False) + config.update( + interface=interface, + point_set=ps, + use_canonical_quadrature_point_ordering=False, + ) context = PointSetContext(**config) expr = context.translator(ufl_expr) @@ -694,6 +736,7 @@ def translate_constant_value(terminal, mt, ctx): @translate.register(Coefficient) def translate_coefficient(terminal, mt, ctx): + domain = extract_unique_domain(terminal) vec = ctx.coefficient(terminal, mt.restriction) if terminal.ufl_element().family() == 'Real': @@ -704,7 +747,7 @@ def translate_coefficient(terminal, mt, ctx): # Collect FInAT tabulation for all entities per_derivative = collections.defaultdict(list) - for entity_id in ctx.entity_ids: + for entity_id in ctx.entity_ids(domain): finat_dict = ctx.basis_evaluation(element, mt, entity_id) for alpha, table in finat_dict.items(): # Filter out irrelevant derivatives @@ -716,14 +759,14 @@ def translate_coefficient(terminal, mt, ctx): per_derivative[alpha].append(table) # Merge entity tabulations for each derivative - if len(ctx.entity_ids) == 1: + if len(ctx.entity_ids(domain)) == 1: def take_singleton(xs): x, = xs # asserts singleton return x per_derivative = {alpha: take_singleton(tables) for alpha, tables in per_derivative.items()} else: - f = ctx.entity_number(extract_unique_domain(terminal), mt.restriction) + f = ctx.entity_number(domain, mt.restriction) per_derivative = {alpha: gem.select_expression(tables, f) for alpha, tables in per_derivative.items()} diff --git a/tsfc/kernel_args.py b/tsfc/kernel_args.py index 80b1bed77f..aa5e5472b5 100644 --- a/tsfc/kernel_args.py +++ b/tsfc/kernel_args.py @@ -54,6 +54,10 @@ class InteriorFacetKernelArg(KernelArg): ... +class OrientationsCellKernelArg(KernelArg): + ... + + class OrientationsExteriorFacetKernelArg(KernelArg): ... diff --git a/tsfc/kernel_interface/__init__.py b/tsfc/kernel_interface/__init__.py index 9b7419e35d..3c20720c33 100644 --- a/tsfc/kernel_interface/__init__.py +++ b/tsfc/kernel_interface/__init__.py @@ -29,6 +29,10 @@ def cell_orientation(self, domain, restriction): def cell_size(self, domain, restriction): """Mesh cell size as a GEM expression. Shape (nvertex, ) in FIAT vertex ordering.""" + @abstractmethod + def entity_ids(self, domain): + """Target indices of entity_number.""" + @abstractmethod def entity_number(self, domain, restriction): """Facet or vertex number as a GEM index.""" diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index d1bf7653db..c371490c69 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -87,6 +87,12 @@ def cell_size(self, domain, restriction): else: return self._cell_sizes[domain] + def entity_ids(self, domain): + """Target indices of entity_number.""" + if not hasattr(self, "_entity_ids"): + raise RuntimeError("Haven't called set_entity_numbers") + return self._entity_ids[domain] + def entity_number(self, domain, restriction): """Facet or vertex number as a GEM index.""" if not hasattr(self, "_entity_numbers"): @@ -244,11 +250,10 @@ def fem_config(self): integral_type = info.integral_type cell = info.domain.ufl_cell() fiat_cell = as_fiat_cell(cell) - integration_dim, entity_ids = lower_integral_type(fiat_cell, integral_type) + integration_dim, _ = lower_integral_type(fiat_cell, integral_type) return dict(interface=self, ufl_cell=cell, integration_dim=integration_dim, - entity_ids=entity_ids, scalar_type=self.fem_scalar_type) def create_context(self): @@ -332,8 +337,8 @@ def set_quad_rule(params, cell, integral_type, functions): fiat_cell = as_fiat_cell(cell) finat_elements = set(create_element(e) for e in elements if e.family() != "Real") fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] - fiat_cell = max_complex(fiat_cells) - + if any(c.is_macrocell() for c in fiat_cells): + fiat_cell = max_complex(fiat_cells) integration_dim, _ = lower_integral_type(fiat_cell, integral_type) quad_rule = fem.get_quadrature_rule(fiat_cell, integration_dim, quadrature_degree, scheme) params["quadrature_rule"] = quad_rule diff --git a/tsfc/kernel_interface/firedrake_loopy.py b/tsfc/kernel_interface/firedrake_loopy.py index 84ef5dd939..bb282e1d25 100644 --- a/tsfc/kernel_interface/firedrake_loopy.py +++ b/tsfc/kernel_interface/firedrake_loopy.py @@ -12,8 +12,8 @@ import loopy as lp from tsfc import kernel_args -from finat.element_factory import create_element -from tsfc.kernel_interface.common import KernelBuilderBase as _KernelBuilderBase, KernelBuilderMixin, get_index_names, check_requirements, prepare_coefficient, prepare_arguments, prepare_constant +from finat.element_factory import as_fiat_cell, create_element +from tsfc.kernel_interface.common import KernelBuilderBase as _KernelBuilderBase, KernelBuilderMixin, get_index_names, check_requirements, prepare_coefficient, prepare_arguments, prepare_constant, lower_integral_type from tsfc.loopy import generate as generate_loopy @@ -30,6 +30,7 @@ 'cell_sizes', 'exterior_facets', 'interior_facets', + 'orientations_cell', 'orientations_exterior_facet', 'orientations_interior_facet']) ActiveDomainNumbers.__doc__ = """ @@ -323,9 +324,16 @@ def set_entity_numbers(self, domains): """ self._entity_numbers = {} + self._entity_ids = {} for i, domain in enumerate(domains): - # Facet number + fiat_cell = as_fiat_cell(domain.ufl_cell()) integral_type = self.integral_data_info.domain_integral_type_map[domain] + if integral_type is None: + # Set placeholder for unused domain. + entity_ids = None + else: + _, entity_ids = lower_integral_type(fiat_cell, integral_type) + self._entity_ids[domain] = entity_ids if integral_type in ['exterior_facet', 'exterior_facet_vert']: facet = gem.Variable(f'facet_{i}', (1,), dtype=gem.uint_type) self._entity_numbers[domain] = {None: gem.VariableIndex(gem.Indexed(facet, (0,))), } @@ -369,7 +377,8 @@ def set_entity_orientations(self, domains): '-': gem.OrientationVariableIndex(gem.Indexed(o, (0,))) } else: - self._entity_orientations[domain] = {None: None} + o = gem.Variable(variable_name, (1,), dtype=gem.uint_type) + self._entity_orientations[domain] = {None: gem.OrientationVariableIndex(gem.Indexed(o, (0,))), } def set_coefficients(self): """Prepare the coefficients of the form.""" @@ -475,6 +484,17 @@ def construct_kernel(self, name, ctx, log=False): dtype=numpy.uint32, ) args.extend(args_) + cell_dict = {} + for domain, expr in self._entity_orientations.items(): + integral_type = info.domain_integral_type_map[domain] + cell_dict[domain] = expr[None].expression if integral_type == "cell" else None + active_domain_numbers_orientations_cell, args_ = self.make_active_domain_numbers( + cell_dict, + active_variables, + kernel_args.OrientationsCellKernelArg, + dtype=gem.uint_type, + ) + args.extend(args_) ext_dict = {} for domain, expr in self._entity_orientations.items(): integral_type = info.domain_integral_type_map[domain] @@ -515,6 +535,7 @@ def construct_kernel(self, name, ctx, log=False): cell_sizes=tuple(active_domain_numbers_cell_sizes), exterior_facets=tuple(active_domain_numbers_exterior_facets), interior_facets=tuple(active_domain_numbers_interior_facets), + orientations_cell=tuple(active_domain_numbers_orientations_cell), orientations_exterior_facet=tuple(active_domain_numbers_orientations_exterior_facet), orientations_interior_facet=tuple(active_domain_numbers_orientations_interior_facet), ),