diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 447345d563..d3f4cfb8ba 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1749,7 +1749,7 @@ def _as_global_kernel_arg_coefficient(_, self): ufl_element = V.ufl_element() if ufl_element.family() == "Real": - return op2.GlobalKernelArg((ufl_element.value_size,)) + return op2.GlobalKernelArg((V.value_size,)) else: return self._make_dat_global_kernel_arg(V, index=index) diff --git a/firedrake/bcs.py b/firedrake/bcs.py index ee38b524e4..272b74ddbf 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -342,20 +342,21 @@ def function_arg(self, g): del self._function_arg_update except AttributeError: pass + V = self.function_space() if isinstance(g, firedrake.Function) and g.ufl_element().family() != "Real": - if g.function_space() != self.function_space(): + if g.function_space() != V: raise RuntimeError("%r is defined on incompatible FunctionSpace!" % g) self._function_arg = g elif isinstance(g, ufl.classes.Zero): - if g.ufl_shape and g.ufl_shape != self.function_space().ufl_element().value_shape: + if g.ufl_shape and g.ufl_shape != V.value_shape: raise ValueError(f"Provided boundary value {g} does not match shape of space") # Special case. Scalar zero for direct Function.assign. self._function_arg = g elif isinstance(g, ufl.classes.Expr): - if g.ufl_shape != self.function_space().ufl_element().value_shape: + if g.ufl_shape != V.value_shape: raise RuntimeError(f"Provided boundary value {g} does not match shape of space") try: - self._function_arg = firedrake.Function(self.function_space()) + self._function_arg = firedrake.Function(V) # Use `Interpolator` instead of assembling an `Interpolate` form # as the expression compilation needs to happen at this stage to # determine if we should use interpolation or projection @@ -363,7 +364,7 @@ def function_arg(self, g): self._function_arg_update = firedrake.Interpolator(g, self._function_arg)._interpolate except (NotImplementedError, AttributeError): # Element doesn't implement interpolation - self._function_arg = firedrake.Function(self.function_space()).project(g) + self._function_arg = firedrake.Function(V).project(g) self._function_arg_update = firedrake.Projector(g, self._function_arg).project else: try: diff --git a/firedrake/checkpointing.py b/firedrake/checkpointing.py index 28aaf04a94..5070a96315 100644 --- a/firedrake/checkpointing.py +++ b/firedrake/checkpointing.py @@ -935,7 +935,7 @@ def save_function(self, f, idx=None, name=None, timestepping_info={}): self._update_function_name_function_space_name_map(tmesh.name, mesh.name, {f.name(): V_name}) # Embed if necessary element = V.ufl_element() - _element = get_embedding_element_for_checkpointing(element) + _element = get_embedding_element_for_checkpointing(element, V.value_shape) if _element != element: path = self._path_to_function_embedded(tmesh.name, mesh.name, V_name, f.name()) self.require_group(path) @@ -1337,7 +1337,7 @@ def load_function(self, mesh, name, idx=None): _name = self.get_attr(path, PREFIX_EMBEDDED + "_function") _f = self.load_function(mesh, _name, idx=idx) element = V.ufl_element() - _element = get_embedding_element_for_checkpointing(element) + _element = get_embedding_element_for_checkpointing(element, V.value_shape) method = get_embedding_method_for_checkpointing(element) assert _element == _f.function_space().ufl_element() f = Function(V, name=name) @@ -1436,8 +1436,7 @@ def _get_shared_data_key_for_checkpointing(self, mesh, ufl_element): shape = ufl_element.reference_value_shape block_size = np.prod(shape) elif isinstance(ufl_element, finat.ufl.VectorElement): - shape = ufl_element.value_shape[:1] - block_size = np.prod(shape) + block_size = ufl_element.reference_value_shape[0] else: block_size = 1 return (nodes_per_entity, real_tensorproduct, block_size) diff --git a/firedrake/constant.py b/firedrake/constant.py index 9363e32f91..145edd5b5f 100644 --- a/firedrake/constant.py +++ b/firedrake/constant.py @@ -79,7 +79,7 @@ def __new__(cls, value, domain=None, name=None, count=None): if not isinstance(domain, ufl.AbstractDomain): cell = ufl.as_cell(domain) - coordinate_element = finat.ufl.VectorElement("Lagrange", cell, 1, gdim=cell.geometric_dimension) + coordinate_element = finat.ufl.VectorElement("Lagrange", cell, 1, dim=cell.topological_dimension()) domain = ufl.Mesh(coordinate_element) cell = domain.ufl_cell() diff --git a/firedrake/embedding.py b/firedrake/embedding.py index 3a618a44f1..e16d799902 100644 --- a/firedrake/embedding.py +++ b/firedrake/embedding.py @@ -4,7 +4,7 @@ import ufl -def get_embedding_dg_element(element, broken_cg=False): +def get_embedding_dg_element(element, value_shape, broken_cg=False): cell = element.cell family = lambda c: "DG" if c.is_simplex() else "DQ" if isinstance(cell, ufl.TensorProductCell): @@ -19,7 +19,7 @@ def get_embedding_dg_element(element, broken_cg=False): scalar_element = finat.ufl.FiniteElement(family(cell), cell=cell, degree=degree) if broken_cg: scalar_element = finat.ufl.BrokenElement(scalar_element.reconstruct(family="Lagrange")) - shape = element.value_shape + shape = value_shape if len(shape) == 0: DG = scalar_element elif len(shape) == 1: @@ -37,12 +37,12 @@ def get_embedding_dg_element(element, broken_cg=False): native_elements_for_checkpointing = {"Lagrange", "Discontinuous Lagrange", "Q", "DQ", "Real"} -def get_embedding_element_for_checkpointing(element): +def get_embedding_element_for_checkpointing(element, value_shape): """Convert the given UFL element to an element that :class:`~.CheckpointFile` can handle.""" if element.family() in native_elements_for_checkpointing: return element else: - return get_embedding_dg_element(element) + return get_embedding_dg_element(element, value_shape) def get_embedding_method_for_checkpointing(element): diff --git a/firedrake/external_operators/point_expr_operator.py b/firedrake/external_operators/point_expr_operator.py index 8d6b8fd992..3aa40e1d5b 100644 --- a/firedrake/external_operators/point_expr_operator.py +++ b/firedrake/external_operators/point_expr_operator.py @@ -44,7 +44,7 @@ def __init__(self, *operands, function_space, derivatives=None, argument_slots=( if not isinstance(operator_data["func"], types.FunctionType): raise TypeError("Expecting a FunctionType pointwise expression") expr_shape = operator_data["func"](*operands).ufl_shape - if expr_shape != function_space.ufl_element().value_shape: + if expr_shape != function_space.value_shape: raise ValueError("The dimension does not match with the dimension of the function space %s" % function_space) @property diff --git a/firedrake/formmanipulation.py b/firedrake/formmanipulation.py index eda7b730a1..35a6789107 100644 --- a/firedrake/formmanipulation.py +++ b/firedrake/formmanipulation.py @@ -124,8 +124,7 @@ def argument(self, o): args += [a_[j] for j in numpy.ndindex(a_.ufl_shape)] else: args += [Zero() - for j in numpy.ndindex( - V_is[i].ufl_element().value_shape)] + for j in numpy.ndindex(V_is[i].value_shape)] return self._arg_cache.setdefault(o, as_vector(args)) diff --git a/firedrake/function.py b/firedrake/function.py index ddedebad10..83136407ea 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -633,7 +633,7 @@ def at(self, arg, *args, **kwargs): raise NotImplementedError("Point evaluation not implemented for variable layers") # Validate geometric dimension - gdim = mesh.ufl_cell().geometric_dimension() + gdim = mesh.geometric_dimension() if arg.shape[-1] == gdim: pass elif len(arg.shape) == 1 and gdim == 1: diff --git a/firedrake/functionspace.py b/firedrake/functionspace.py index 2d6726b2a1..74dab04703 100644 --- a/firedrake/functionspace.py +++ b/firedrake/functionspace.py @@ -184,7 +184,7 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, """ sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) if dim is None: - dim = mesh.ufl_cell().geometric_dimension() + dim = mesh.geometric_dimension() if not isinstance(dim, numbers.Integral) and dim > 0: raise ValueError(f"Can't make VectorFunctionSpace with dim={dim}") element = finat.ufl.VectorElement(sub_element, dim=dim) @@ -237,7 +237,7 @@ def TensorFunctionSpace(mesh, family, degree=None, shape=None, """ sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) - shape = shape or (mesh.ufl_cell().geometric_dimension(),) * 2 + shape = shape or (mesh.geometric_dimension(),) * 2 element = finat.ufl.TensorElement(sub_element, shape=shape, symmetry=symmetry) return FunctionSpace(mesh, element, name=name) diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index a0cb6accc3..53f69e92ce 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -176,16 +176,13 @@ def split(self): def _components(self): if len(self) == 1: return tuple(type(self).create(self.topological.sub(i), self.mesh()) - for i in range(self.value_size)) + for i in range(self.block_size)) else: return self.subfunctions @PETSc.Log.EventDecorator() def sub(self, i): - if len(self) == 1: - bound = self.value_size - else: - bound = len(self) + bound = len(self._components) if i < 0 or i >= bound: raise IndexError("Invalid component %d, not in [0, %d)" % (i, bound)) return self._components[i] @@ -489,7 +486,7 @@ def __init__(self, mesh, element, name=None): shape_element = element if isinstance(element, finat.ufl.WithMapping): shape_element = element.wrapee - sub = shape_element.sub_elements[0].value_shape + sub = shape_element.sub_elements[0].reference_value_shape self.shape = rvs[:len(rvs) - len(sub)] else: self.shape = () @@ -497,6 +494,9 @@ def __init__(self, mesh, element, name=None): self._ufl_function_space = ufl.FunctionSpace(mesh.ufl_mesh(), element, label=self._label) self._mesh = mesh + self.value_size = self._ufl_function_space.value_size + r"""The number of scalar components of this :class:`FunctionSpace`.""" + self.rank = len(self.shape) r"""The rank of this :class:`FunctionSpace`. Spaces where the element is scalar-valued (or intrinsically vector-valued) have @@ -505,7 +505,7 @@ def __init__(self, mesh, element, name=None): the number of components of their :attr:`finat.ufl.finiteelementbase.FiniteElementBase.value_shape`.""" - self.value_size = int(numpy.prod(self.shape, dtype=int)) + self.block_size = int(numpy.prod(self.shape, dtype=int)) r"""The total number of degrees of freedom at each function space node.""" self.name = name @@ -654,7 +654,7 @@ def __getitem__(self, i): @utils.cached_property def _components(self): - return tuple(ComponentFunctionSpace(self, i) for i in range(self.value_size)) + return tuple(ComponentFunctionSpace(self, i) for i in range(self.block_size)) def sub(self, i): r"""Return a view into the ith component.""" @@ -684,7 +684,7 @@ def node_count(self): def dof_count(self): r"""The number of degrees of freedom (includes halo dofs) of this function space on this process. Cf. :attr:`FunctionSpace.node_count` .""" - return self.node_count*self.value_size + return self.node_count*self.block_size def dim(self): r"""The global number of degrees of freedom for this function space. @@ -821,7 +821,7 @@ def local_to_global_map(self, bcs, lgmap=None): else: indices = lgmap.block_indices.copy() bsize = lgmap.getBlockSize() - assert bsize == self.value_size + assert bsize == self.block_size else: # MatBlock case, LGMap is already unrolled. indices = lgmap.block_indices.copy() @@ -830,11 +830,11 @@ def local_to_global_map(self, bcs, lgmap=None): nodes = [] for bc in bcs: if bc.function_space().component is not None: - nodes.append(bc.nodes * self.value_size + nodes.append(bc.nodes * self.block_size + bc.function_space().component) elif unblocked: - tmp = bc.nodes * self.value_size - for i in range(self.value_size): + tmp = bc.nodes * self.block_size + for i in range(self.block_size): nodes.append(tmp + i) else: nodes.append(bc.nodes) @@ -1300,9 +1300,9 @@ def ComponentFunctionSpace(parent, component): """ element = parent.ufl_element() assert type(element) in frozenset([finat.ufl.VectorElement, finat.ufl.TensorElement]) - if not (0 <= component < parent.value_size): + if not (0 <= component < parent.block_size): raise IndexError("Invalid component %d. not in [0, %d)" % - (component, parent.value_size)) + (component, parent.block_size)) new = ProxyFunctionSpace(parent.mesh(), element.sub_elements[0], name=parent.name) new.identifier = "component" new.component = component @@ -1346,7 +1346,7 @@ def make_dof_dset(self): def make_dat(self, val=None, valuetype=None, name=None): r"""Return a newly allocated :class:`pyop2.types.glob.Global` representing the data for a :class:`.Function` on this space.""" - return op2.Global(self.value_size, val, valuetype, name, self._comm) + return op2.Global(self.block_size, val, valuetype, name, self._comm) def entity_node_map(self, source_mesh, source_integral_type, source_subdomain_id, source_all_integer_subdomain_ids): return None diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 70f6c31f86..1a26285904 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -547,7 +547,7 @@ def __init__( # VectorFunctionSpace equivalent is built from the scalar # sub-element. ufl_scalar_element = ufl_scalar_element.sub_elements[0] - if ufl_scalar_element.value_shape != (): + if ufl_scalar_element.reference_value_shape != (): raise NotImplementedError( "Can't yet cross-mesh interpolate onto function spaces made from VectorElements or TensorElements made from sub elements with value shape other than ()." ) @@ -614,7 +614,7 @@ def __init__( # I first point evaluate my expression at these locations, giving a # P0DG function on the VOM. As described in the manual, this is an # interpolation operation. - shape = V_dest.ufl_element().value_shape + shape = V_dest.ufl_function_space().value_shape if len(shape) == 0: fs_type = firedrake.FunctionSpace elif len(shape) == 1: @@ -988,18 +988,16 @@ def callable(): else: # Make sure we have an expression of the right length i.e. a value for # each component in the value shape of each function space - dims = [numpy.prod(fs.ufl_element().value_shape, dtype=int) - for fs in V] loops = [] - if numpy.prod(expr.ufl_shape, dtype=int) != sum(dims): + if numpy.prod(expr.ufl_shape, dtype=int) != V.value_size: raise RuntimeError('Expression of length %d required, got length %d' - % (sum(dims), numpy.prod(expr.ufl_shape, dtype=int))) + % (V.value_size, numpy.prod(expr.ufl_shape, dtype=int))) if len(V) > 1: raise NotImplementedError( "UFL expressions for mixed functions are not yet supported.") loops.extend(_interpolator(V, tensor, expr, subset, arguments, access, bcs=bcs)) if bcs and len(arguments) == 0: - loops.extend([partial(bc.apply, f) for bc in bcs]) + loops.extend(partial(bc.apply, f) for bc in bcs) def callable(loops, f): for l in loops: @@ -1024,13 +1022,13 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None): if access is op2.READ: raise ValueError("Can't have READ access for output function") - if len(expr.ufl_shape) != len(V.ufl_element().value_shape): + if len(expr.ufl_shape) != len(V.value_shape): raise RuntimeError('Rank mismatch: Expression rank %d, FunctionSpace rank %d' - % (len(expr.ufl_shape), len(V.ufl_element().value_shape))) + % (len(expr.ufl_shape), len(V.value_shape))) - if expr.ufl_shape != V.ufl_element().value_shape: + if expr.ufl_shape != V.value_shape: raise RuntimeError('Shape mismatch: Expression shape %r, FunctionSpace shape %r' - % (expr.ufl_shape, V.ufl_element().value_shape)) + % (expr.ufl_shape, V.value_shape)) # NOTE: The par_loop is always over the target mesh cells. target_mesh = as_domain(V) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 05d7a4a77f..fd37e7b575 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -66,10 +66,10 @@ } -_supported_embedded_cell_types = [ufl.Cell('interval', 2), - ufl.Cell('triangle', 3), - ufl.Cell("quadrilateral", 3), - ufl.TensorProductCell(ufl.Cell('interval'), ufl.Cell('interval'), geometric_dimension=3)] +_supported_embedded_cell_types_and_gdims = [(ufl.Cell('interval'), 2), + (ufl.Cell('triangle'), 3), + (ufl.Cell("quadrilateral"), 3), + (ufl.TensorProductCell(ufl.Cell('interval'), ufl.Cell('interval')), 3)] unmarked = -1 @@ -1477,8 +1477,8 @@ def mark_entities(self, tf, label_value, label_name=None): elem = tV.ufl_element() if tV.mesh() is not self: raise RuntimeError(f"tf must be defined on {self}: {tf.mesh()} is not {self}") - if elem.value_shape != (): - raise RuntimeError(f"tf must be scalar: {elem.value_shape} != ()") + if elem.reference_value_shape != (): + raise RuntimeError(f"tf must be scalar: {elem.reference_value_shape} != ()") if elem.family() in {"Discontinuous Lagrange", "DQ"} and elem.degree() == 0: # cells height = 0 @@ -2303,7 +2303,7 @@ def callback(self): self.topology.init() 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.ufl_coordinate_element().cell.geometric_dimension())) + (self.num_vertices(), self.geometric_dimension())) coordinates = function.CoordinatelessFunction(coordinates_fs, val=coordinates_data, name=_generate_default_mesh_coordinates_name(self.name)) @@ -2470,7 +2470,7 @@ def spatial_index(self): from firedrake import function, functionspace from firedrake.parloops import par_loop, READ, MIN, MAX - gdim = self.ufl_cell().geometric_dimension() + gdim = self.geometric_dimension() if gdim <= 1: info_red("libspatialindex does not support 1-dimension, falling back on brute force.") return None @@ -2765,7 +2765,7 @@ def init_cell_orientations(self, expr): import firedrake.function as function import firedrake.functionspace as functionspace - if self.ufl_cell() not in _supported_embedded_cell_types: + if (self.ufl_cell(), self.geometric_dimension()) not in _supported_embedded_cell_types_and_gdims: raise NotImplementedError('Only implemented for intervals embedded in 2d and triangles and quadrilaterals embedded in 3d') if hasattr(self, '_cell_orientations'): @@ -2774,8 +2774,8 @@ def init_cell_orientations(self, expr): if not isinstance(expr, ufl.classes.Expr): raise TypeError("UFL expression expected!") - if expr.ufl_shape != (self.ufl_cell().geometric_dimension(), ): - raise ValueError(f"Mismatching shapes: expr.ufl_shape ({expr.ufl_shape}) != (self.ufl_cell().geometric_dimension(), ) (({self.ufl_cell().geometric_dimension}, ))") + if expr.ufl_shape != (self.geometric_dimension(), ): + raise ValueError(f"Mismatching shapes: expr.ufl_shape ({expr.ufl_shape}) != (self.geometric_dimension(), ) (({self.geometric_dimension}, ))") fs = functionspace.FunctionSpace(self, 'DG', 0) x = ufl.SpatialCoordinate(self) @@ -2848,12 +2848,9 @@ def make_mesh_from_coordinates(coordinates, name, tolerance=0.5): V = coordinates.function_space() element = coordinates.ufl_element() - if V.rank != 1 or len(element.value_shape) != 1: + if V.rank != 1 or len(element.reference_value_shape) != 1: raise ValueError("Coordinates must be from a rank-1 FunctionSpace with rank-1 value_shape.") assert V.mesh().ufl_cell().topological_dimension() <= V.value_size - # Build coordinate element - cell = element.cell.reconstruct(geometric_dimension=V.value_size) - element = element.reconstruct(cell=cell) mesh = MeshGeometry.__new__(MeshGeometry, element, coordinates.comm) mesh.__init__(coordinates) @@ -2886,11 +2883,10 @@ def make_mesh_from_mesh_topology(topology, name, tolerance=0.5): # TODO: meshfile might indicates higher-order coordinate element cell = topology.ufl_cell() geometric_dim = topology.topology_dm.getCoordinateDim() - cell = cell.reconstruct(geometric_dimension=geometric_dim) if not topology.topology_dm.getCoordinatesLocalized(): - element = finat.ufl.VectorElement("Lagrange", cell, 1) + element = finat.ufl.VectorElement("Lagrange", cell, 1, dim=geometric_dim) else: - element = finat.ufl.VectorElement("DQ" if cell in [ufl.quadrilateral, ufl.hexahedron] else "DG", cell, 1, variant="equispaced") + element = finat.ufl.VectorElement("DQ" if cell in [ufl.quadrilateral, ufl.hexahedron] else "DG", cell, 1, dim=geometric_dim, variant="equispaced") # Create mesh object mesh = MeshGeometry.__new__(MeshGeometry, element, topology.comm) mesh._init_topology(topology) @@ -2922,9 +2918,8 @@ def make_vom_from_vom_topology(topology, name, tolerance=0.5): import firedrake.function as function gdim = topology.topology_dm.getCoordinateDim() - tcell = topology.ufl_cell() - cell = tcell.reconstruct(geometric_dimension=gdim) - element = finat.ufl.VectorElement("DG", cell, 0) + cell = topology.ufl_cell() + element = finat.ufl.VectorElement("DG", cell, 0, dim=gdim) vmesh = MeshGeometry.__new__(MeshGeometry, element, topology.comm) vmesh._init_topology(topology) # Save vertex reference coordinate (within reference cell) in function @@ -3214,7 +3209,7 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', peri pass elif extrusion_type in ("radial", "radial_hedgehog"): # do not allow radial extrusion if tdim = gdim - if mesh.ufl_cell().geometric_dimension() == mesh.ufl_cell().topological_dimension(): + if mesh.geometric_dimension() == mesh.topological_dimension(): raise RuntimeError("Cannot radially-extrude a mesh with equal geometric and topological dimension") else: # check for kernel @@ -3234,7 +3229,7 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', peri element = finat.ufl.TensorProductElement(helement, velement) if gdim is None: - gdim = mesh.ufl_cell().geometric_dimension() + (extrusion_type == "uniform") + gdim = mesh.geometric_dimension() + (extrusion_type == "uniform") coordinates_fs = functionspace.VectorFunctionSpace(topology, element, dim=gdim) coordinates = function.CoordinatelessFunction(coordinates_fs, name=_generate_default_mesh_coordinates_name(name)) @@ -4537,8 +4532,8 @@ def RelabeledMesh(mesh, indicator_functions, subdomain_ids, **kwargs): plex1.createLabel(label_name) for f, subid in zip(indicator_functions, subdomain_ids): elem = f.topological.function_space().ufl_element() - if elem.value_shape != (): - raise RuntimeError(f"indicator functions must be scalar: got {elem.value_shape} != ()") + if elem.reference_value_shape != (): + raise RuntimeError(f"indicator functions must be scalar: got {elem.reference_value_shape} != ()") if elem.family() in {"Discontinuous Lagrange", "DQ"} and elem.degree() == 0: # cells height = 0 diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py index 752ec6ff13..551c09a7e2 100644 --- a/firedrake/mg/embedded.py +++ b/firedrake/mg/embedded.py @@ -17,9 +17,9 @@ non_native_variants = frozenset(["integral", "fdm", "alfeld"]) -def get_embedding_element(element): +def get_embedding_element(element, value_shape): broken_cg = element.sobolev_space in {ufl.H1, ufl.H2} - dg_element = get_embedding_dg_element(element, broken_cg=broken_cg) + dg_element = get_embedding_dg_element(element, value_shape, broken_cg=broken_cg) variant = element.variant() or "default" family = element.family() # Elements on Alfeld splits are embedded onto DG Powell-Sabin. @@ -40,8 +40,8 @@ class Cache(object): """A caching object for work vectors and matrices. :arg element: The element to use for the caching.""" - def __init__(self, element): - self.embedding_element = get_embedding_element(element) + def __init__(self, ufl_element, value_shape): + self.embedding_element = get_embedding_dg_element(ufl_element, value_shape) self._dat_versions = {} self._V_DG_mass = {} self._DG_inv_mass = {} @@ -83,11 +83,12 @@ def _native_transfer(self, element, op): return self.native_transfers.setdefault(element, ops)[op] return None - def cache(self, element): + def cache(self, V): + key = (V.ufl_element(), V.value_shape) try: - return self.caches[element] + return self.caches[key] except KeyError: - return self.caches.setdefault(element, TransferManager.Cache(element)) + return self.caches.setdefault(key, TransferManager.Cache(*key)) def V_dof_weights(self, V): """Dof weights for averaging projection. @@ -95,7 +96,7 @@ def V_dof_weights(self, V): :arg V: function space to compute weights for. :returns: A PETSc Vec. """ - cache = self.cache(V.ufl_element()) + cache = self.cache(V) key = V.dim() try: return cache._V_dof_weights[key] @@ -106,7 +107,7 @@ def V_dof_weights(self, V): # global Vector counting the number of cells that see each # dof. f = firedrake.Function(V) - firedrake.par_loop(("{[i, j]: 0 <= i < A.dofs and 0 <= j < %d}" % V.value_size, + firedrake.par_loop(("{[i, j]: 0 <= i < A.dofs and 0 <= j < %d}" % V.block_size, "A[i, j] = A[i, j] + 1"), firedrake.dx, {"A": (f, firedrake.INC)}) @@ -120,7 +121,7 @@ def V_DG_mass(self, V, DG): :arg DG: the DG space :returns: A PETSc Mat mapping from V -> DG """ - cache = self.cache(V.ufl_element()) + cache = self.cache(V) key = V.dim() try: return cache._V_DG_mass[key] @@ -135,7 +136,7 @@ def DG_inv_mass(self, DG): :arg DG: the DG space :returns: A PETSc Mat. """ - cache = self.caches[DG.ufl_element()] + cache = self.cache(DG) key = DG.dim() try: return cache._DG_inv_mass[key] @@ -151,7 +152,7 @@ def V_approx_inv_mass(self, V, DG): :arg DG: the DG space :returns: A PETSc Mat mapping from V -> DG. """ - cache = self.cache(V.ufl_element()) + cache = self.cache(V) key = V.dim() try: return cache._V_approx_inv_mass[key] @@ -169,7 +170,7 @@ def V_inv_mass_ksp(self, V): :arg V: a function space. :returns: A PETSc KSP for inverting (V, V). """ - cache = self.cache(V.ufl_element()) + cache = self.cache(V) key = V.dim() try: return cache._V_inv_mass_ksp[key] @@ -191,7 +192,7 @@ def DG_work(self, V): :returns: A Function in the embedding DG space. """ needs_dual = ufl.duals.is_dual(V) - cache = self.cache(V.ufl_element()) + cache = self.cache(V) key = (V.dim(), needs_dual) try: return cache._DG_work[key] @@ -208,28 +209,28 @@ def work_vec(self, V): :arg V: a function space. :returns: A PETSc Vec for V. """ - cache = self.cache(V.ufl_element()) + cache = self.cache(V) key = V.dim() try: return cache._work_vec[key] except KeyError: return cache._work_vec.setdefault(key, V.dof_dset.layout_vec.duplicate()) - def requires_transfer(self, element, transfer_op, source, target): + def requires_transfer(self, V, transfer_op, source, target): """Determine whether either the source or target have been modified since the last time a grid transfer was executed with them.""" key = (transfer_op, weakref.ref(source.dat), weakref.ref(target.dat)) dat_versions = (source.dat.dat_version, target.dat.dat_version) try: - return self.cache(element)._dat_versions[key] != dat_versions + return self.cache(V)._dat_versions[key] != dat_versions except KeyError: return True - def cache_dat_versions(self, element, transfer_op, source, target): + def cache_dat_versions(self, V, transfer_op, source, target): """Record the returned dat_versions of the source and target.""" key = (transfer_op, weakref.ref(source.dat), weakref.ref(target.dat)) dat_versions = (source.dat.dat_version, target.dat.dat_version) - self.cache(element)._dat_versions[key] = dat_versions + self.cache(V)._dat_versions[key] = dat_versions @PETSc.Log.EventDecorator() def op(self, source, target, transfer_op): @@ -243,7 +244,7 @@ def op(self, source, target, transfer_op): Vt = target.function_space() source_element = Vs.ufl_element() target_element = Vt.ufl_element() - if not self.requires_transfer(source_element, transfer_op, source, target): + if not self.requires_transfer(Vs, transfer_op, source, target): return if all(self.is_native(e, transfer_op) for e in (source_element, target_element)): @@ -280,7 +281,7 @@ def op(self, source, target, transfer_op): work = self.work_vec(Vt) self.V_DG_mass(Vt, VDGt).multTranspose(dgv, work) self.V_inv_mass_ksp(Vt).solve(work, t) - self.cache_dat_versions(source_element, transfer_op, source, target) + self.cache_dat_versions(Vs, transfer_op, source, target) def prolong(self, uc, uf): """Prolong a function. @@ -308,7 +309,7 @@ def restrict(self, source, target): Vt_star = target.function_space() source_element = Vs_star.ufl_element() target_element = Vt_star.ufl_element() - if not self.requires_transfer(source_element, Op.RESTRICT, source, target): + if not self.requires_transfer(Vs_star, Op.RESTRICT, source, target): return if all(self.is_native(e, Op.RESTRICT) for e in (source_element, target_element)): @@ -344,4 +345,4 @@ def restrict(self, source, target): with dgtarget.dat.vec_ro as dgv, target.dat.vec_wo as t: self.DG_inv_mass(VDGt).mult(dgv, dgwork) self.V_DG_mass(Vt, VDGt).multTranspose(dgwork, t) - self.cache_dat_versions(source_element, Op.RESTRICT, source, target) + self.cache_dat_versions(Vs_star, Op.RESTRICT, source, target) diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py index 2c82092f51..5405b6d726 100644 --- a/firedrake/mg/kernels.py +++ b/firedrake/mg/kernels.py @@ -46,11 +46,11 @@ def to_reference_coordinates(ufl_coordinate_element, parameters=None): # Create FInAT element element = tsfc.finatinterface.create_element(ufl_coordinate_element) - + gdim, = ufl_coordinate_element.reference_value_shape cell = ufl_coordinate_element.cell code = { - "geometric_dimension": cell.geometric_dimension(), + "geometric_dimension": gdim, "topological_dimension": cell.topological_dimension(), "to_reference_coords_newton_step": to_reference_coords_newton_step_body(ufl_coordinate_element, parameters, x0_dtype=ScalarType, dX_dtype="double"), "init_X": init_X(element.cell, parameters), @@ -177,11 +177,10 @@ def compile_element(expression, dual_space=None, parameters=None, return_variable = gem.Indexed(gem.Variable('R', finat_elem.index_shape), argument_multiindex) result = gem.Indexed(result, tensor_indices) if dual_space: - elem = create_element(dual_space.ufl_element()) - if elem.value_shape: - var = gem.Indexed(gem.Variable("b", elem.value_shape), - tensor_indices) - b_arg = [lp.GlobalArg("b", dtype=ScalarType, shape=elem.value_shape)] + value_shape = dual_space.value_shape + if value_shape: + var = gem.Indexed(gem.Variable("b", value_shape), tensor_indices) + b_arg = [lp.GlobalArg("b", dtype=ScalarType, shape=value_shape)] else: var = gem.Indexed(gem.Variable("b", (1, )), (0, )) b_arg = [lp.GlobalArg("b", dtype=ScalarType, shape=(1,))] @@ -220,7 +219,7 @@ def prolong_kernel(expression): assert hierarchy._meshes[int(idx)].cell_set._extruded V = expression.function_space() key = (("prolong",) - + expression.ufl_element().value_shape + + V.value_shape + entity_dofs_key(V.finat_element.complex.get_topology()) + entity_dofs_key(V.finat_element.entity_dofs()) + entity_dofs_key(coordinates.function_space().finat_element.entity_dofs())) @@ -284,7 +283,7 @@ def prolong_kernel(expression): "evaluate": eval_code, "spacedim": element.cell.get_spatial_dimension(), "ncandidate": hierarchy.fine_to_coarse_cells[levelf].shape[1], - "Rdim": numpy.prod(element.value_shape), + "Rdim": V.value_size, "inside_cell": inside_check(element.cell, eps=1e-8, X="Xref"), "celldist_l1_c_expr": celldist_l1_c_expr(element.cell, X="Xref"), "Xc_cell_inc": coords_element.space_dimension(), @@ -302,7 +301,7 @@ def restrict_kernel(Vf, Vc): if Vf.extruded: assert Vc.extruded key = (("restrict",) - + Vf.ufl_element().value_shape + + Vf.value_shape + entity_dofs_key(Vf.finat_element.complex.get_topology()) + entity_dofs_key(Vc.finat_element.complex.get_topology()) + entity_dofs_key(Vf.finat_element.entity_dofs()) @@ -390,7 +389,7 @@ def inject_kernel(Vf, Vc): else: level_ratio = 1 key = (("inject", level_ratio) - + Vf.ufl_element().value_shape + + Vf.value_shape + entity_dofs_key(Vc.finat_element.complex.get_topology()) + entity_dofs_key(Vf.finat_element.complex.get_topology()) + entity_dofs_key(Vc.finat_element.entity_dofs()) @@ -465,7 +464,7 @@ def inject_kernel(Vf, Vc): "celldist_l1_c_expr": celldist_l1_c_expr(Vc.finat_element.cell, X="Xref"), "tdim": Vc.mesh().topological_dimension(), "ncandidate": ncandidate, - "Rdim": numpy.prod(Vf_element.value_shape), + "Rdim": numpy.prod(Vf.value_shape), "Xf_cell_inc": coords_element.space_dimension(), "f_cell_inc": Vf_element.space_dimension() } diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index a40fad62f7..42ea985666 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -173,7 +173,7 @@ def MeshHierarchy(mesh, refinement_levels, scale = mesh._radius / np.linalg.norm(coords, axis=1).reshape(-1, 1) coords *= scale - meshes = [mesh] + [mesh_builder(dm, dim=mesh.ufl_cell().geometric_dimension(), + meshes = [mesh] + [mesh_builder(dm, dim=mesh.geometric_dimension(), distribution_parameters=distribution_parameters, reorder=reorder, comm=mesh.comm) for dm in dms] diff --git a/firedrake/mg/utils.py b/firedrake/mg/utils.py index f3eb475b50..80e5e820f6 100644 --- a/firedrake/mg/utils.py +++ b/firedrake/mg/utils.py @@ -135,7 +135,7 @@ def coarse_cell_to_fine_node_map(Vc, Vf): def physical_node_locations(V): element = V.ufl_element() - if element.value_shape: + if V.value_shape: assert isinstance(element, (finat.ufl.VectorElement, finat.ufl.TensorElement)) element = element.sub_elements[0] mesh = V.mesh() @@ -146,7 +146,7 @@ def physical_node_locations(V): try: return cache[key] except KeyError: - Vc = firedrake.FunctionSpace(mesh, finat.ufl.VectorElement(element)) + Vc = firedrake.VectorFunctionSpace(mesh, element) # FIXME: This is unsafe for DG coordinates and CG target spaces. locations = firedrake.assemble(firedrake.Interpolate(firedrake.SpatialCoordinate(mesh), Vc)) return cache.setdefault(key, locations) diff --git a/firedrake/pointeval_utils.py b/firedrake/pointeval_utils.py index f8a985273d..75ed2b4564 100644 --- a/firedrake/pointeval_utils.py +++ b/firedrake/pointeval_utils.py @@ -116,7 +116,7 @@ def predicate(index): extruded = isinstance(cell, TensorProductCell) code = { - "geometric_dimension": cell.geometric_dimension(), + "geometric_dimension": domain.geometric_dimension(), "layers_arg": ", int const *__restrict__ layers" if extruded else "", "layers": ", layers" if extruded else "", "extruded_define": "1" if extruded else "0", diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 82059def13..4793af0bfe 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -7,6 +7,7 @@ from pyop2 import op2 from pyop2.parloop import generate_single_cell_wrapper +from firedrake.mesh import MeshGeometry from firedrake.petsc import PETSc from firedrake.utils import IntType, as_cstr, ScalarType, ScalarType_c, complex_mode, RealType_c @@ -34,7 +35,7 @@ def make_wrapper(function, **kwargs): def src_locate_cell(mesh, tolerance=None): src = ['#include '] - src.append(compile_coordinate_element(mesh.ufl_coordinate_element(), tolerance)) + src.append(compile_coordinate_element(mesh, tolerance)) src.append(make_wrapper(mesh.coordinates, forward_args=["void*", "double*", RealType_c+"*"], kernel_name="to_reference_coords_kernel", @@ -129,9 +130,10 @@ def to_reference_coords_newton_step(ufl_coordinate_element, parameters, x0_dtype # Set up UFL form cell = ufl_coordinate_element.cell domain = ufl.Mesh(ufl_coordinate_element) + gdim = domain.geometric_dimension() K = ufl.JacobianInverse(domain) x = ufl.SpatialCoordinate(domain) - x0_element = finat.ufl.VectorElement("Real", cell, 0) + x0_element = finat.ufl.VectorElement("Real", cell, 0, dim=gdim) x0 = ufl.Coefficient(ufl.FunctionSpace(domain, x0_element)) expr = ufl.dot(K, x - x0) @@ -192,11 +194,23 @@ def predicate(index): @PETSc.Log.EventDecorator() -def compile_coordinate_element(ufl_coordinate_element, contains_eps, parameters=None): +def compile_coordinate_element(mesh: MeshGeometry, contains_eps: float, parameters: dict | None = None): """Generates C code for changing to reference coordinates. - :arg ufl_coordinate_element: UFL element of the coordinates - :returns: C code as string + Parameters + ---------- + mesh : + The mesh. + + contains_eps : + The tolerance used to verify that a point is contained by a cell. + parameters : + Form compiler parameters, defaults to whatever TSFC defaults to. + + Returns + ------- + str + A string of C code. """ if parameters is None: parameters = tsfc.default_parameters() @@ -204,25 +218,24 @@ def compile_coordinate_element(ufl_coordinate_element, contains_eps, parameters= _ = tsfc.default_parameters() _.update(parameters) parameters = _ + + ufl_coordinate_element = mesh.ufl_coordinate_element() # Create FInAT element element = tsfc.finatinterface.create_element(ufl_coordinate_element) - cell = ufl_coordinate_element.cell - extruded = isinstance(cell, ufl.TensorProductCell) - code = { - "geometric_dimension": cell.geometric_dimension(), - "topological_dimension": cell.topological_dimension(), + "geometric_dimension": mesh.geometric_dimension(), + "topological_dimension": mesh.topological_dimension(), "celldist_l1_c_expr": celldist_l1_c_expr(element.cell, "X"), "to_reference_coords_newton_step": to_reference_coords_newton_step(ufl_coordinate_element, parameters), "init_X": init_X(element.cell, parameters), "max_iteration_count": 1 if is_affine(ufl_coordinate_element) else 16, "convergence_epsilon": 1e-12, - "dX_norm_square": dX_norm_square(cell.topological_dimension()), - "X_isub_dX": X_isub_dX(cell.topological_dimension()), - "extruded_arg": ", int const *__restrict__ layers" if extruded else "", - "extr_comment_out": "//" if extruded else "", - "non_extr_comment_out": "//" if not extruded else "", + "dX_norm_square": dX_norm_square(mesh.topological_dimension()), + "X_isub_dX": X_isub_dX(mesh.topological_dimension()), + "extruded_arg": ", int const *__restrict__ layers" if mesh.extruded else "", + "extr_comment_out": "//" if mesh.extruded else "", + "non_extr_comment_out": "//" if not mesh.extruded else "", "IntType": as_cstr(IntType), "ScalarType": ScalarType_c, "RealType": RealType_c, diff --git a/firedrake/preconditioners/asm.py b/firedrake/preconditioners/asm.py index 86100138e6..19738fe438 100644 --- a/firedrake/preconditioners/asm.py +++ b/firedrake/preconditioners/asm.py @@ -81,8 +81,8 @@ def initialize(self, pc): tinyasm.SetASMLocalSubdomains( asmpc, ises, [W.dm.getDefaultSF() for W in V], - [W.value_size for W in V], - sum(W.value_size * W.dof_dset.total_size for W in V)) + [W.block_size for W in V], + sum(W.block_size * W.dof_dset.total_size for W in V)) asmpc.setUp() else: raise ValueError(f"Unknown backend type {backend}") @@ -168,7 +168,7 @@ def get_patches(self, V): continue off = section.getOffset(p) # Local indices within W - W_indices = slice(off*W.value_size, W.value_size * (off + dof)) + W_indices = slice(off*W.block_size, W.block_size * (off + dof)) indices.extend(V_local_ises_indices[i][W_indices]) iset = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) ises.append(iset) @@ -247,7 +247,7 @@ def get_patches(self, V): continue off = section.getOffset(p) # Local indices within W - W_indices = slice(off*W.value_size, W.value_size * (off + dof)) + W_indices = slice(off*W.block_size, W.block_size * (off + dof)) indices.extend(V_local_ises_indices[i][W_indices]) iset = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) ises.append(iset) @@ -296,7 +296,7 @@ def get_patches(self, V): if dof <= 0: continue off = section.getOffset(p) - indices = numpy.arange(off*V.value_size, V.value_size * (off + dof), dtype=IntType) + indices = numpy.arange(off*V.block_size, V.block_size * (off + dof), dtype=IntType) iset = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) ises.append(iset) @@ -471,7 +471,7 @@ def get_patches(self, V): else: begin = off + min(k, k+plane) * blayer_offset + dof end = off + max(k, k+plane) * blayer_offset - zlice = slice(W.value_size * begin, W.value_size * end) + zlice = slice(W.block_size * begin, W.block_size * end) indices.extend(iset[zlice]) iset = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF) diff --git a/firedrake/preconditioners/facet_split.py b/firedrake/preconditioners/facet_split.py index a1c43dca7c..4c8e654ce4 100644 --- a/firedrake/preconditioners/facet_split.py +++ b/firedrake/preconditioners/facet_split.py @@ -258,7 +258,7 @@ def get_restriction_indices(V, W): eperm = numpy.concatenate((eperm, numpy.setdiff1d(numpy.arange(vsize, dtype=PETSc.IntType), eperm))) pmap = PermutedMap(V.cell_node_map(), eperm) - wsize = sum(Vsub.finat_element.space_dimension() * Vsub.value_size for Vsub in W) + wsize = sum(Vsub.finat_element.space_dimension() * Vsub.block_size for Vsub in W) kernel_code = f""" void copy(PetscInt *restrict w, const PetscInt *restrict v) {{ for (PetscInt i=0; i<{wsize}; i++) w[i] = v[i]; diff --git a/firedrake/preconditioners/fdm.py b/firedrake/preconditioners/fdm.py index 53d03c2a7f..82e56057e7 100644 --- a/firedrake/preconditioners/fdm.py +++ b/firedrake/preconditioners/fdm.py @@ -225,10 +225,10 @@ def allocate_matrix(self, Amat, V, J, bcs, fcp, pmat_type, use_static_condensati raise ValueError("Expecting at most one FunctionSpace restricted onto facets.") self.embedding_element = ebig - if Vbig.value_size == 1: + if Vbig.block_size == 1: self.fises = PETSc.IS().createGeneral(fdofs, comm=COMM_SELF) else: - self.fises = PETSc.IS().createBlock(Vbig.value_size, fdofs, comm=COMM_SELF) + self.fises = PETSc.IS().createBlock(Vbig.block_size, fdofs, comm=COMM_SELF) # Create data structures needed for assembly self.lgmaps = {Vsub: Vsub.local_to_global_map([bc for bc in bcs if bc.function_space() == Vsub]) for Vsub in V} @@ -245,7 +245,7 @@ def allocate_matrix(self, Amat, V, J, bcs, fcp, pmat_type, use_static_condensati self.schur_kernel[V] = SchurComplementPattern elif Vfacet and use_static_condensation: # If we are in a facet space, we build the Schur complement on its diagonal block - if Vfacet.finat_element.formdegree == 0 and Vfacet.value_size == 1: + if Vfacet.finat_element.formdegree == 0 and Vfacet.block_size == 1: self.schur_kernel[Vfacet] = SchurComplementDiagonal interior_pc_type = PETSc.PC.Type.JACOBI elif symmetric: @@ -560,7 +560,7 @@ def assemble_reference_tensor(self, V, transpose=False, sort_interior=False): :returns: a :class:`PETSc.Mat` interpolating V^k * d(V^k) onto broken(V^k) * broken(V^{k+1}) on the reference element. """ - value_size = V.value_size + bsize = V.block_size fe = V.finat_element tdim = fe.cell.get_spatial_dimension() formdegree = fe.formdegree @@ -570,7 +570,7 @@ def assemble_reference_tensor(self, V, transpose=False, sort_interior=False): if formdegree == tdim: degree = degree + 1 is_interior, is_facet = is_restricted(fe) - key = (value_size, tdim, degree, formdegree, is_interior, is_facet, transpose, sort_interior) + key = (bsize, tdim, degree, formdegree, is_interior, is_facet, transpose, sort_interior) cache = self._cache.setdefault("reference_tensor", {}) try: return cache[key] @@ -637,8 +637,8 @@ def assemble_reference_tensor(self, V, transpose=False, sort_interior=False): A00.destroy() A11.destroy() A10.destroy() - if value_size != 1: - eye = petsc_sparse(numpy.eye(value_size), comm=result.getComm()) + if bsize != 1: + eye = petsc_sparse(numpy.eye(bsize), comm=result.getComm()) temp = result result = temp.kron(eye) temp.destroy() @@ -1352,10 +1352,10 @@ def __init__(self, kernel, name=None): V = FunctionSpace(Q.mesh(), unrestrict_element(Q.ufl_element())) V0 = FunctionSpace(Q.mesh(), restrict_element(V.ufl_element(), "interior")) V1 = FunctionSpace(Q.mesh(), restrict_element(V.ufl_element(), "facet")) - idofs = PETSc.IS().createBlock(V.value_size, restricted_dofs(V0.finat_element, V.finat_element), comm=comm) - fdofs = PETSc.IS().createBlock(V.value_size, restricted_dofs(V1.finat_element, V.finat_element), comm=comm) + idofs = PETSc.IS().createBlock(V.block_size, restricted_dofs(V0.finat_element, V.finat_element), comm=comm) + fdofs = PETSc.IS().createBlock(V.block_size, restricted_dofs(V1.finat_element, V.finat_element), comm=comm) size = idofs.size + fdofs.size - assert size == V.finat_element.space_dimension() * V.value_size + assert size == V.finat_element.space_dimension() * V.block_size # Bilinear form on the space with interior and interface a = form if Q == V else form(*(t.reconstruct(function_space=V) for t in args)) # Generate code to apply the action of A within the Schur complement action @@ -1621,15 +1621,14 @@ def tabulate_exterior_derivative(Vc, Vf, cbcs=[], fbcs=[], comm=None): fises.destroy() cises.destroy() - if Vf.value_size > 1: + if Vf.block_size > 1: temp = Dhat - eye = petsc_sparse(numpy.eye(Vf.value_size, dtype=PETSc.RealType), comm=temp.getComm()) + eye = petsc_sparse(numpy.eye(Vf.block_size, dtype=PETSc.RealType), comm=temp.getComm()) Dhat = temp.kron(eye) temp.destroy() eye.destroy() sizes = tuple(V.dof_dset.layout_vec.getSizes() for V in (Vf, Vc)) - block_size = Vf.dof_dset.layout_vec.getBlockSize() preallocator = PETSc.Mat().create(comm=comm) preallocator.setType(PETSc.Mat.Type.PREALLOCATOR) preallocator.setSizes(sizes) @@ -1647,7 +1646,7 @@ def tabulate_exterior_derivative(Vc, Vf, cbcs=[], fbcs=[], comm=None): nnz = get_preallocation(preallocator, sizes[0][0]) preallocator.destroy() - Dmat = PETSc.Mat().createAIJ(sizes, block_size, nnz=nnz, comm=comm) + Dmat = PETSc.Mat().createAIJ(sizes, Vf.block_size, nnz=nnz, comm=comm) Dmat.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, True) assembler.arguments[0].data = Dmat.handle assembler() @@ -1865,8 +1864,7 @@ def cell_to_global(lgmap, cell_to_local, cell_index, result=None): result = cell_to_local(cell_index, result=result) return lgmap.apply(result, result=result) - bsize = Vrow.dof_dset.layout_vec.getBlockSize() - cell_to_local, nel = extrude_node_map(Vrow.cell_node_map(), bsize=bsize) + cell_to_local, nel = extrude_node_map(Vrow.cell_node_map(), bsize=Vrow.block_size) get_rindices = partial(cell_to_global, self.lgmaps[Vrow], cell_to_local) Afdm, Dfdm, bdof, axes_shifts = self.assemble_reference_tensor(Vrow) @@ -1877,7 +1875,7 @@ def cell_to_global(lgmap, cell_to_local, cell_index, result=None): PT_facet = self.coefficients.get("PT_facet") V = Vrow - bsize = V.value_size + bsize = V.block_size ncomp = V.ufl_element().reference_value_size sdim = (V.finat_element.space_dimension() * bsize) // ncomp # dimension of a single component tdim = V.mesh().topological_dimension() diff --git a/firedrake/preconditioners/hiptmair.py b/firedrake/preconditioners/hiptmair.py index 51adace67e..858a07d0d5 100644 --- a/firedrake/preconditioners/hiptmair.py +++ b/firedrake/preconditioners/hiptmair.py @@ -215,7 +215,7 @@ def curl_to_grad(ele): if isinstance(ele, finat.ufl.VectorElement): return type(ele)(curl_to_grad(ele._sub_element), dim=ele.num_sub_elements) elif isinstance(ele, finat.ufl.TensorElement): - return type(ele)(curl_to_grad(ele._sub_element), shape=ele.value_shape, symmetry=ele.symmetry()) + return type(ele)(curl_to_grad(ele._sub_element), shape=ele._shape, symmetry=ele.symmetry()) elif isinstance(ele, finat.ufl.MixedElement): return type(ele)(*(curl_to_grad(e) for e in ele.sub_elements)) elif isinstance(ele, finat.ufl.RestrictedElement): @@ -242,7 +242,7 @@ def div_to_curl(ele): if isinstance(ele, finat.ufl.VectorElement): return type(ele)(div_to_curl(ele._sub_element), dim=ele.num_sub_elements) elif isinstance(ele, finat.ufl.TensorElement): - return type(ele)(div_to_curl(ele._sub_element), shape=ele.value_shape, symmetry=ele.symmetry()) + return type(ele)(div_to_curl(ele._sub_element), shape=ele._shape, symmetry=ele.symmetry()) elif isinstance(ele, finat.ufl.MixedElement): return type(ele)(*(div_to_curl(e) for e in ele.sub_elements)) elif isinstance(ele, finat.ufl.RestrictedElement): diff --git a/firedrake/preconditioners/patch.py b/firedrake/preconditioners/patch.py index b74897f574..0a7bad5575 100644 --- a/firedrake/preconditioners/patch.py +++ b/firedrake/preconditioners/patch.py @@ -580,23 +580,23 @@ def bcdofs(bc, ghost=True): if isinstance(Z.ufl_element(), VectorElement): offset += idx assert i == len(indices)-1 # assert we're at the end of the chain - assert Z.sub(idx).value_size == 1 + assert Z.sub(idx).block_size == 1 elif isinstance(Z.ufl_element(), MixedElement): if ghost: offset += sum(Z.sub(j).dof_count for j in range(idx)) else: - offset += sum(Z.sub(j).dof_dset.size * Z.sub(j).value_size for j in range(idx)) + offset += sum(Z.sub(j).dof_dset.size * Z.sub(j).block_size for j in range(idx)) else: raise NotImplementedError("How are you taking a .sub?") Z = Z.sub(idx) if Z.parent is not None and isinstance(Z.parent.ufl_element(), VectorElement): - bs = Z.parent.value_size + bs = Z.parent.block_size start = 0 stop = 1 else: - bs = Z.value_size + bs = Z.block_size start = 0 stop = bs nodes = bc.nodes @@ -868,7 +868,7 @@ def initialize(self, obj): offsets = numpy.append([0], numpy.cumsum([W.dof_count for W in V])).astype(PETSc.IntType) patch.setPatchDiscretisationInfo([W.dm for W in V], - numpy.array([W.value_size for + numpy.array([W.block_size for W in V], dtype=PETSc.IntType), [W.cell_node_list for W in V], offsets, diff --git a/firedrake/preconditioners/pmg.py b/firedrake/preconditioners/pmg.py index 7775219e60..251f585cbe 100644 --- a/firedrake/preconditioners/pmg.py +++ b/firedrake/preconditioners/pmg.py @@ -131,9 +131,7 @@ def initialize(self, obj): elements = [ele] while True: try: - ele_ = self.coarsen_element(ele) - assert ele_.value_shape == ele.value_shape - ele = ele_ + ele = self.coarsen_element(ele) except ValueError: break elements.append(ele) @@ -906,19 +904,19 @@ def make_kron_code(Vc, Vf, t_in, t_out, mat_name, scratch): shifts = fshifts in_place = False if len(felems) == len(celems): - in_place = all((len(fs)*Vf.value_size == len(cs)*Vc.value_size) for fs, cs in zip(fshifts, cshifts)) - psize = Vf.value_size + in_place = all((len(fs)*Vf.block_size == len(cs)*Vc.block_size) for fs, cs in zip(fshifts, cshifts)) + psize = Vf.block_size if not in_place: if len(celems) == 1: - psize = Vc.value_size + psize = Vc.block_size pelem = celems[0] perm_name = "perm_%s" % t_in celems = celems*len(felems) elif len(felems) == 1: shifts = cshifts - psize = Vf.value_size + psize = Vf.block_size pelem = felems[0] perm_name = "perm_%s" % t_out felems = felems*len(celems) @@ -926,7 +924,7 @@ def make_kron_code(Vc, Vf, t_in, t_out, mat_name, scratch): raise ValueError("Cannot assign fine to coarse DOFs") if set(cshifts) == set(fshifts): - csize = Vc.value_size * Vc.finat_element.space_dimension() + csize = Vc.block_size * Vc.finat_element.space_dimension() prolong_code.append(f""" for({IntType_c} j=1; j<{len(fshifts)}; j++) for({IntType_c} i=0; i<{csize}; i++) @@ -940,8 +938,8 @@ def make_kron_code(Vc, Vf, t_in, t_out, mat_name, scratch): elif pelem == celems[0]: for k in range(len(shifts)): - if Vc.value_size*len(shifts[k]) < Vf.value_size: - shifts[k] = shifts[k]*(Vf.value_size//Vc.value_size) + if Vc.block_size*len(shifts[k]) < Vf.block_size: + shifts[k] = shifts[k]*(Vf.block_size//Vc.block_size) pshape = [e.space_dimension() for e in pelem] pargs = ", ".join(map(str, pshape+[1]*(3-len(pshape)))) @@ -1098,7 +1096,7 @@ def make_mapping_code(Q, cmapping, fmapping, t_in, t_out): if B: tensor = ufl.dot(B, tensor) if tensor else B if tensor is None: - tensor = ufl.Identity(Q.ufl_element().value_shape[0]) + tensor = ufl.Identity(Q.value_shape[0]) u = ufl.Coefficient(Q) expr = ufl.dot(tensor, u) @@ -1117,7 +1115,7 @@ def make_mapping_code(Q, cmapping, fmapping, t_in, t_out): coef_args = "".join([", c%d" % i for i in range(len(coefficients))]) coef_decl = "".join([", PetscScalar const *restrict c%d" % i for i in range(len(coefficients))]) - qlen = Q.value_size * Q.finat_element.space_dimension() + qlen = Q.block_size * Q.finat_element.space_dimension() prolong_code = f""" for({IntType_c} i=0; i<{qlen}; i++) {t_out}[i] = 0.0E0; @@ -1223,7 +1221,7 @@ def work_function(self, V): @cached_property def _weight(self): weight = firedrake.Function(self.Vf) - size = self.Vf.finat_element.space_dimension() * self.Vf.value_size + size = self.Vf.finat_element.space_dimension() * self.Vf.block_size kernel_code = f""" void weight(PetscScalar *restrict w){{ for(PetscInt i=0; i<{size}; i++) w[i] += 1.0; @@ -1347,12 +1345,12 @@ def make_blas_kernels(self, Vf, Vc): in_place_mapping = True except Exception: qelem = finat.ufl.FiniteElement("DQ", cell=felem.cell, degree=PMGBase.max_degree(felem)) - if felem.value_shape: - qelem = finat.ufl.TensorElement(qelem, shape=felem.value_shape, symmetry=felem.symmetry()) + if Vf.value_shape: + qelem = finat.ufl.TensorElement(qelem, shape=Vf.value_shape, symmetry=felem.symmetry()) Qf = firedrake.FunctionSpace(Vf.mesh(), qelem) mapping_output = make_mapping_code(Qf, cmapping, fmapping, "t0", "t1") - qshape = (Qf.value_size, Qf.finat_element.space_dimension()) + qshape = (Qf.block_size, Qf.finat_element.space_dimension()) # interpolate to embedding fine space decl[0], prolong[0], restrict[0], shapes = make_kron_code(Vc, Qf, "t0", "t1", "J0", "t2") @@ -1377,8 +1375,8 @@ def make_blas_kernels(self, Vf, Vc): # We could benefit from loop tiling for the transpose, but that makes the code # more complicated. - fshape = (Vf.value_size, Vf.finat_element.space_dimension()) - cshape = (Vc.value_size, Vc.finat_element.space_dimension()) + fshape = (Vf.block_size, Vf.finat_element.space_dimension()) + cshape = (Vc.block_size, Vc.finat_element.space_dimension()) lwork = numpy.prod([max(*dims) for dims in zip(*shapes)]) lwork = max(lwork, max(numpy.prod(fshape), numpy.prod(cshape))) @@ -1470,8 +1468,8 @@ def make_kernels(self, Vf, Vc): element_kernel = element_kernel.replace("void expression_kernel", "static void expression_kernel") coef_args = "".join([", c%d" % i for i in range(len(coefficients))]) coef_decl = "".join([", const %s *restrict c%d" % (ScalarType_c, i) for i in range(len(coefficients))]) - dimc = Vc.finat_element.space_dimension() * Vc.value_size - dimf = Vf.finat_element.space_dimension() * Vf.value_size + dimc = Vc.finat_element.space_dimension() * Vc.block_size + dimf = Vf.finat_element.space_dimension() * Vf.block_size restrict_code = f""" {element_kernel} diff --git a/firedrake/pyplot/pgf.py b/firedrake/pyplot/pgf.py index 4f4862b722..6640454679 100644 --- a/firedrake/pyplot/pgf.py +++ b/firedrake/pyplot/pgf.py @@ -217,12 +217,12 @@ def pgfplot(f, filename, degree=1, complex_component='real', print_latex_example raise NotImplementedError(f"Not yet implemented for functions in spatial dimension {dim}") if mesh.extruded: raise NotImplementedError("Not yet implemented for functions on extruded meshes") - if elem.value_shape(): + if V.value_shape: raise NotImplementedError("Currently only implemeted for scalar functions") - coordelem = get_embedding_dg_element(mesh.coordinates.function_space().ufl_element()).reconstruct(degree=degree, variant="equispaced") + coordelem = get_embedding_dg_element(mesh.coordinates.function_space().ufl_element(), (dim, )).reconstruct(degree=degree, variant="equispaced") coordV = FunctionSpace(mesh, coordelem) coords = Function(coordV).interpolate(SpatialCoordinate(mesh)) - elemdg = get_embedding_dg_element(elem).reconstruct(degree=degree, variant="equispaced") + elemdg = get_embedding_dg_element(elem, V.value_shape).reconstruct(degree=degree, variant="equispaced") Vdg = FunctionSpace(mesh, elemdg) fdg = Function(Vdg) method = get_embedding_method_for_checkpointing(elem) diff --git a/firedrake/slate/slac/compiler.py b/firedrake/slate/slac/compiler.py index 1333a04039..7e1b14281c 100644 --- a/firedrake/slate/slac/compiler.py +++ b/firedrake/slate/slac/compiler.py @@ -82,7 +82,8 @@ def _compile_expression_hashkey(slate_expr, compiler_parameters=None): def _compile_expression_comm(*args, **kwargs): # args[0] is a slate_expr - return args[0].ufl_domains()[0].comm + domain, = args[0].ufl_domains() + return domain.comm @memory_and_disk_cache( diff --git a/firedrake/slate/slate.py b/firedrake/slate/slate.py index 79b65ba79a..85b7af3635 100644 --- a/firedrake/slate/slate.py +++ b/firedrake/slate/slate.py @@ -31,7 +31,7 @@ from ufl.algorithms.map_integrands import map_integrand_dags from ufl.corealg.multifunction import MultiFunction from ufl.classes import Zero -from ufl.domain import join_domains +from ufl.domain import join_domains, sort_domains from ufl.form import Form import hashlib @@ -198,7 +198,7 @@ def shapes(self): """ shapes = OrderedDict() for i, fs in enumerate(self.arg_function_spaces): - shapes[i] = tuple(int(V.finat_element.space_dimension() * V.value_size) + shapes[i] = tuple(int(V.finat_element.space_dimension() * V.block_size) for V in fs) return shapes @@ -247,11 +247,11 @@ def ufl_domain(self): The function will fail if multiple domains are found. """ - domains = self.ufl_domains() - assert all(domain == domains[0] for domain in domains), ( - "All integrals must share the same domain of integration." - ) - return domains[0] + try: + domain, = self.ufl_domains() + except ValueError: + raise ValueError("All integrals must share the same domain of integration.") + return domain @abstractmethod def ufl_domains(self): @@ -983,7 +983,7 @@ def ufl_domains(self): the tensor. """ collected_domains = [op.ufl_domains() for op in self.operands] - return join_domains(chain(*collected_domains)) + return sort_domains(join_domains(chain(*collected_domains))) def subdomain_data(self): """Returns a mapping on the tensor: diff --git a/firedrake/slate/static_condensation/hybridization.py b/firedrake/slate/static_condensation/hybridization.py index 5dd592c44e..997f969684 100644 --- a/firedrake/slate/static_condensation/hybridization.py +++ b/firedrake/slate/static_condensation/hybridization.py @@ -61,7 +61,7 @@ def initialize(self, pc): if len(V) != 2: raise ValueError("Expecting two function spaces.") - if all(Vi.ufl_element().value_shape for Vi in V): + if all(Vi.value_shape for Vi in V): raise ValueError("Expecting an H(div) x L2 pair of spaces.") # Automagically determine which spaces are vector and scalar diff --git a/firedrake/ufl_expr.py b/firedrake/ufl_expr.py index e129c26037..1f9a66df75 100644 --- a/firedrake/ufl_expr.py +++ b/firedrake/ufl_expr.py @@ -76,7 +76,7 @@ def reconstruct(self, function_space=None, return self if not isinstance(number, int): raise TypeError(f"Expecting an int, not {number}") - if function_space.ufl_element().value_shape != self.ufl_element().value_shape: + if function_space.value_shape != self.function_space().value_shape: raise ValueError("Cannot reconstruct an Argument with a different value shape.") return Argument(function_space, number, part=part) @@ -140,7 +140,7 @@ def reconstruct(self, function_space=None, return self if not isinstance(number, int): raise TypeError(f"Expecting an int, not {number}") - if function_space.ufl_element().value_shape != self.ufl_element().value_shape: + if function_space.value_shape != self.function_space().value_shape: raise ValueError("Cannot reconstruct an Coargument with a different value shape.") return Coargument(function_space, number, part=part) diff --git a/tests/macro/test_macro_multigrid.py b/tests/macro/test_macro_multigrid.py index c45b085555..8f880bd48d 100644 --- a/tests/macro/test_macro_multigrid.py +++ b/tests/macro/test_macro_multigrid.py @@ -136,6 +136,7 @@ def test_macro_grid_transfer(hierarchy, space, degrees, variant, transfer_type): "mg_levels_pc_type": "jacobi", "mg_coarse_pc_type": "python", "mg_coarse_pc_python_type": "firedrake.AssembledPC", + "mg_coarse_assembled_pc_type": "cholesky", } diff --git a/tests/multigrid/test_non_nested.py b/tests/multigrid/test_non_nested.py index 8a68b18b7a..5dfb63ff8b 100644 --- a/tests/multigrid/test_non_nested.py +++ b/tests/multigrid/test_non_nested.py @@ -59,7 +59,7 @@ def test_sphere_mg(): u = TrialFunction(V) v = TestFunction(V) - a = (inner(u, v) + inner(u, v))*dx + a = (inner(grad(u), grad(v)) + inner(u, v))*dx f1 = exp((x+y+z)/R)*x*y*z/R**3 F = inner(f1, v)*dx @@ -91,4 +91,4 @@ def test_sphere_mg(): prob = LinearVariationalProblem(a, F, w) solver = LinearVariationalSolver(prob, solver_parameters=mg_params) solver.solve() - assert solver.snes.ksp.getIterationNumber() < 5 + assert solver.snes.ksp.getIterationNumber() < 7 diff --git a/tests/output/test_io_backward_compat.py b/tests/output/test_io_backward_compat.py index 67cf0d6dff..0e43326a2a 100644 --- a/tests/output/test_io_backward_compat.py +++ b/tests/output/test_io_backward_compat.py @@ -150,7 +150,7 @@ def _get_mesh_and_V(params): def _get_expr(V): mesh = V.mesh() dim = mesh.geometric_dimension() - shape = V.ufl_element().value_shape + shape = V.value_shape if dim == 2: x, y = SpatialCoordinate(mesh) z = x + y @@ -240,9 +240,8 @@ def test_io_backward_compat_base_load(version, params): def _get_expr_timestepping(V, i): mesh = V.mesh() - element = V.ufl_element() x, y = SpatialCoordinate(mesh) - shape = element.value_shape + shape = V.value_shape if shape == (4, ): return as_vector([x + i, y + i, x * y + i, i * i]) else: diff --git a/tests/output/test_io_function.py b/tests/output/test_io_function.py index fc48a220b4..64be091afe 100644 --- a/tests/output/test_io_function.py +++ b/tests/output/test_io_function.py @@ -67,7 +67,7 @@ def _get_mesh(cell_type, comm): def _get_expr(V): mesh = V.mesh() dim = mesh.geometric_dimension() - shape = V.ufl_element().value_shape + shape = V.value_shape if dim == 2: x, y = SpatialCoordinate(mesh) z = x * y diff --git a/tests/output/test_io_timestepping.py b/tests/output/test_io_timestepping.py index c12c30d74a..b2a66648bb 100644 --- a/tests/output/test_io_timestepping.py +++ b/tests/output/test_io_timestepping.py @@ -16,7 +16,7 @@ def _get_expr(V, i): mesh = V.mesh() element = V.ufl_element() x, y = SpatialCoordinate(mesh) - shape = element.value_shape + shape = V.value_shape if element.family() == "Real": return 7. + i * i elif shape == (): diff --git a/tests/regression/test_2dcohomology.py b/tests/regression/test_2dcohomology.py index 73bf409538..85ad069924 100644 --- a/tests/regression/test_2dcohomology.py +++ b/tests/regression/test_2dcohomology.py @@ -51,7 +51,7 @@ def test_betti0(space, mesh): L = assemble(inner(nabla_grad(u), nabla_grad(v))*dx) - bc0 = DirichletBC(V0, Constant(0.0), 9) + bc0 = DirichletBC(V0, 0, 9) L0 = assemble(inner(nabla_grad(u), nabla_grad(v))*dx, bcs=[bc0]) u, s, v = linalg.svd(L.M.values) @@ -94,8 +94,8 @@ def test_betti1(space, mesh): L = assemble((inner(sigma, tau) - inner(u, rot(tau)) + inner(rot(sigma), v) + inner(div(u), div(v))) * dx) - bc0 = DirichletBC(W.sub(0), Constant(0.0), 9) - bc1 = DirichletBC(W.sub(1), Constant((0.0, 0.0)), 9) + bc0 = DirichletBC(W.sub(0), 0, 9) + bc1 = DirichletBC(W.sub(1), 0, 9) L0 = assemble((inner(sigma, tau) - inner(u, rot(tau)) + inner(rot(sigma), v) + inner(div(u), div(v))) * dx, bcs=[bc0, bc1]) @@ -155,7 +155,7 @@ def test_betti2(space, mesh): L = assemble((inner(sigma, tau) - inner(u, div(tau)) + inner(div(sigma), v))*dx) - bc1 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), 9) + bc1 = DirichletBC(W.sub(0), 0, 9) L0 = assemble((inner(sigma, tau) - inner(u, div(tau)) + inner(div(sigma), v))*dx, bcs=[bc1]) dV1 = V1.dof_count diff --git a/tests/regression/test_change_coordinates.py b/tests/regression/test_change_coordinates.py index 05e46bd8d7..a712ffd6d5 100644 --- a/tests/regression/test_change_coordinates.py +++ b/tests/regression/test_change_coordinates.py @@ -13,7 +13,7 @@ def test_immerse_1d(dim): m = Mesh(new_coords) - assert m.ufl_cell().geometric_dimension() == dim + assert m.geometric_dimension() == dim def test_immerse_2d(): @@ -23,7 +23,7 @@ def test_immerse_2d(): m = Mesh(new_coords) - assert m.ufl_cell().geometric_dimension() == 3 + assert m.geometric_dimension() == 3 def test_project_2d(): @@ -33,7 +33,7 @@ def test_project_2d(): m = Mesh(new_coords) - assert m.ufl_cell().geometric_dimension() == 1 + assert m.geometric_dimension() == 1 def test_immerse_extruded(): @@ -44,4 +44,4 @@ def test_immerse_extruded(): m = Mesh(new_coords) - assert m.ufl_cell().geometric_dimension() == 3 + assert m.geometric_dimension() == 3 diff --git a/tests/regression/test_ensembleparallelism.py b/tests/regression/test_ensembleparallelism.py index faa3db99dc..49db8e7d3f 100644 --- a/tests/regression/test_ensembleparallelism.py +++ b/tests/regression/test_ensembleparallelism.py @@ -205,13 +205,13 @@ def test_ensemble_reduce(ensemble, mesh, W, urank, urank_sum, root, blocking): parallel_assert( lambda: error < 1e-12, subset=root_ranks, - msg=f"{error = :.5f}" + msg=f"{error=:.5f}" ) error = errornorm(Function(W).assign(10), u_reduce) parallel_assert( lambda: error < 1e-12, subset={range(COMM_WORLD.size)} - root_ranks, - msg=f"{error = :.5f}" + msg=f"{error=:.5f}" ) # check that u_reduce dat vector is still synchronised @@ -347,7 +347,7 @@ def test_send_and_recv(ensemble, mesh, W, blocking): parallel_assert( lambda: error < 1e-12, subset=root_ranks, - msg=f"{error = :.5f}" + msg=f"{error=:.5f}" ) diff --git a/tests/regression/test_expressions.py b/tests/regression/test_expressions.py index 9387e52684..b1fd8feb64 100644 --- a/tests/regression/test_expressions.py +++ b/tests/regression/test_expressions.py @@ -295,7 +295,7 @@ def test_assign_with_different_meshes_fails(): def test_assign_vector_const_to_vfs(vcg1): f = Function(vcg1) - c = Constant(range(1, f.ufl_element().value_shape[0]+1)) + c = Constant(range(1, f.function_space().value_shape[0]+1)) f.assign(c) assert np.allclose(f.dat.data_ro, c.dat.data_ro) diff --git a/tests/regression/test_fas_snespatch.py b/tests/regression/test_fas_snespatch.py index b140f285ae..63deb3ab5d 100644 --- a/tests/regression/test_fas_snespatch.py +++ b/tests/regression/test_fas_snespatch.py @@ -170,8 +170,7 @@ def test_snespatch(mesh, CG1, solver_params): f = Constant(1, domain=mesh) F = inner(grad(u), grad(v))*dx - inner(f, v)*dx + inner(u**3 - u, v)*dx - z = zero(CG1.ufl_element().value_shape) - bcs = DirichletBC(CG1, z, "on_boundary") + bcs = DirichletBC(CG1, 0, "on_boundary") nvproblem = NonlinearVariationalProblem(F, u, bcs=bcs) solver = NonlinearVariationalSolver(nvproblem, solver_parameters=solver_params) diff --git a/tests/regression/test_fdm.py b/tests/regression/test_fdm.py index fef477190e..947ae43fbc 100644 --- a/tests/regression/test_fdm.py +++ b/tests/regression/test_fdm.py @@ -90,7 +90,7 @@ def build_riesz_map(V, d): x = SpatialCoordinate(V.mesh()) x -= Constant([0.5]*len(x)) - if V.ufl_element().value_shape == (): + if V.value_shape == (): u_exact = exp(-10*dot(x, x)) u_bc = u_exact else: @@ -195,7 +195,7 @@ def test_variable_coefficient(mesh): subs = ("on_boundary",) if mesh.cell_set._extruded: subs += ("top", "bottom") - bcs = [DirichletBC(V, zero(V.ufl_element().value_shape), sub) for sub in subs] + bcs = [DirichletBC(V, 0, sub) for sub in subs] uh = Function(V) problem = LinearVariationalProblem(a, L, uh, bcs=bcs) @@ -228,7 +228,7 @@ def test_ipdg_direct_solver(fs): mesh = fs.mesh() x = SpatialCoordinate(mesh) gdim = mesh.geometric_dimension() - ncomp = fs.ufl_element().value_size + ncomp = fs.value_size homogenize = gdim > 2 if homogenize: diff --git a/tests/regression/test_function.py b/tests/regression/test_function.py index 506d9bcb8e..b42b71762a 100644 --- a/tests/regression/test_function.py +++ b/tests/regression/test_function.py @@ -97,7 +97,7 @@ def test_mismatching_shape_interpolation(V): VV = VectorFunctionSpace(V.mesh(), 'CG', 1) f = Function(VV) with pytest.raises(RuntimeError): - f.interpolate(Constant([1] * (VV.ufl_element().value_shape[0] + 1))) + f.interpolate(Constant([1] * (VV.value_shape[0] + 1))) def test_function_val(V): diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 220c02bfde..dbb81ec700 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -276,7 +276,7 @@ def test_reconstruct_component(space, dg0, rt1, mesh, mesh2, dual): Z = {"dg0": dg0, "rt1": rt1}[space] if dual: Z = Z.dual() - for component in range(Z.value_size): + for component in range(len(Z)): V1 = Z.sub(component) V2 = V1.reconstruct(mesh=mesh2) assert is_dual(V1) == is_dual(V2) == dual @@ -293,7 +293,7 @@ def test_reconstruct_sub_component(dg0, rt1, mesh, mesh2, dual): if dual: Z = Z.dual() for index, Vsub in enumerate(Z): - for component in range(Vsub.value_size): + for component in range(len(Vsub._components)): V1 = Z.sub(index).sub(component) V2 = V1.reconstruct(mesh=mesh2) assert is_dual(V1) == is_dual(V2) == dual diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index b375a6306d..4a963aa36e 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -471,7 +471,7 @@ def test_interpolation_tensor_convergence(): V = TensorFunctionSpace(mesh, "RT", 1) x, y = SpatialCoordinate(mesh) - vs = V.ufl_element().value_shape + vs = V.value_shape expr = as_tensor(np.asarray([ sin(2*pi*x*(i+1))*cos(4*pi*y*i) for i in range(np.prod(vs, dtype=int)) diff --git a/tests/regression/test_interpolate_cross_mesh.py b/tests/regression/test_interpolate_cross_mesh.py index 3f53287885..94252c3df6 100644 --- a/tests/regression/test_interpolate_cross_mesh.py +++ b/tests/regression/test_interpolate_cross_mesh.py @@ -361,7 +361,7 @@ def test_interpolate_unitsquare_mixed(): # Can't go from non-mixed to mixed V_src_2 = VectorFunctionSpace(m_src, "CG", 1) - assert V_src_2.ufl_element().value_shape == V_src.ufl_element().value_shape + assert V_src_2.value_shape == V_src.value_shape f_src_2 = Function(V_src_2) with pytest.raises(NotImplementedError): assemble(interpolate(f_src_2, V_dest)) diff --git a/tests/regression/test_interpolate_vs_project.py b/tests/regression/test_interpolate_vs_project.py index 1a6b38969e..84fabd6a45 100644 --- a/tests/regression/test_interpolate_vs_project.py +++ b/tests/regression/test_interpolate_vs_project.py @@ -37,7 +37,7 @@ def test_interpolate_vs_project(V): elif dim == 3: x, y, z = SpatialCoordinate(mesh) - shape = V.ufl_element().value_shape + shape = V.value_shape if dim == 2: if len(shape) == 0: expression = x + y diff --git a/tests/regression/test_p1pc.py b/tests/regression/test_p1pc.py index 1bdbfe3181..7f1d27e2d2 100644 --- a/tests/regression/test_p1pc.py +++ b/tests/regression/test_p1pc.py @@ -36,7 +36,7 @@ def test_p_independence(mesh, expected): L = inner(Constant(1), v)*dx - bcs = DirichletBC(V, zero(V.ufl_element().value_shape), "on_boundary") + bcs = DirichletBC(V, 0, "on_boundary") uh = Function(V) problem = LinearVariationalProblem(a, L, uh, bcs=bcs) diff --git a/tests/regression/test_patch_pc.py b/tests/regression/test_patch_pc.py index 5a70438089..766da8b611 100644 --- a/tests/regression/test_patch_pc.py +++ b/tests/regression/test_patch_pc.py @@ -38,7 +38,7 @@ def test_jacobi_sor_equivalence(mesh, problem_type, multiplicative): R = TensorFunctionSpace(mesh, "CG", 1) V = P*Q*R - shape = V.ufl_element().value_shape + shape = V.value_shape rhs = numpy.full(shape, 1, dtype=float) u = TrialFunction(V) @@ -48,16 +48,16 @@ def test_jacobi_sor_equivalence(mesh, problem_type, multiplicative): # We also test patch pc with kernel argument compression. i = 1 # only active index f = Function(V) - fval = numpy.full(V.sub(i).ufl_element().value_shape, 1.0, dtype=float) + fval = numpy.full(V.sub(i).value_shape, 1.0, dtype=float) f.sub(i).interpolate(Constant(fval)) a = (inner(f[i], f[i]) * inner(grad(u), grad(v)))*dx L = inner(Constant(rhs), v)*dx - bcs = [DirichletBC(Q, zero(Q.ufl_element().value_shape), "on_boundary") + bcs = [DirichletBC(Q, 0, "on_boundary") for Q in V.subfunctions] else: a = inner(grad(u), grad(v))*dx L = inner(Constant(rhs), v)*dx - bcs = DirichletBC(V, zero(V.ufl_element().value_shape), "on_boundary") + bcs = DirichletBC(V, 0, "on_boundary") uh = Function(V) problem = LinearVariationalProblem(a, L, uh, bcs=bcs) diff --git a/tests/regression/test_patch_precompute_element_tensors.py b/tests/regression/test_patch_precompute_element_tensors.py index 4907dcef56..33bdf64c19 100644 --- a/tests/regression/test_patch_precompute_element_tensors.py +++ b/tests/regression/test_patch_precompute_element_tensors.py @@ -27,8 +27,7 @@ def test_patch_precompute_element_tensors(mesh, V): f = Constant((1, 1)) F = inner(grad(u), grad(v))*dx + gamma*inner(div(u), div(v))*dx - inner(f, v)*dx + avg(inner(u, v))*dS - z = zero(V.ufl_element().value_shape) - bcs = DirichletBC(V, z, "on_boundary") + bcs = DirichletBC(V, 0, "on_boundary") sp = { "mat_type": "matfree", diff --git a/tests/regression/test_tensor_algebra.py b/tests/regression/test_tensor_algebra.py index 8ab3192087..61871d4ef2 100644 --- a/tests/regression/test_tensor_algebra.py +++ b/tests/regression/test_tensor_algebra.py @@ -34,7 +34,7 @@ def mesh(request): "-2*mu_s*inner(grad(u), outer(conj(v), n)) * ds")], ids=lambda x: x[0]) def form_expect(request, mesh): - dim = mesh.ufl_cell().geometric_dimension() + dim = mesh.geometric_dimension() if mesh.ufl_cell().cellname() == "quadrilateral": V = FunctionSpace(mesh, "RTCF", 1) else: diff --git a/tests/submesh/test_submesh_interpolate.py b/tests/submesh/test_submesh_interpolate.py index f2d94459e3..27c7414d4a 100644 --- a/tests/submesh/test_submesh_interpolate.py +++ b/tests/submesh/test_submesh_interpolate.py @@ -21,9 +21,9 @@ def _get_expr(V): x, y, z = SpatialCoordinate(m) else: raise NotImplementedError("Not implemented") - if V.ufl_element().value_shape == (): + if V.value_shape == (): return cos(x) + x * exp(y) + sin(z) - elif V.ufl_element().value_shape == (2, ): + elif V.value_shape == (2, ): return as_vector([cos(x), sin(y)]) @@ -50,7 +50,7 @@ def _test_submesh_interpolate_cell_cell(mesh, subdomain_cond, fe_fesub): # if there is no ambiguity on the subdomain boundary. # For testing, the following suffices. g.interpolate(f) - temp = Constant(999.) if V.ufl_element().value_shape == () else as_vector([Constant(999.)] * np.prod(V.ufl_element().value_shape, dtype=int)) + temp = Constant(999.*np.ones(V.value_shape)) g.interpolate(temp, subset=mesh.topology.cell_subset(label_value)) # pollute the data g.interpolate(gsub, subset=mesh.topology.cell_subset(label_value)) assert assemble(inner(g - f, g - f) * dx(label_value)).real < 1e-14 diff --git a/tests/vertexonly/test_swarm.py b/tests/vertexonly/test_swarm.py index e6e0848fb9..21c0b12eca 100644 --- a/tests/vertexonly/test_swarm.py +++ b/tests/vertexonly/test_swarm.py @@ -27,11 +27,11 @@ def cell_midpoints(m): num_cells_local = len(f.dat.data_ro) num_cells = MPI.COMM_WORLD.allreduce(num_cells_local, op=MPI.SUM) # reshape is for 1D case where f.dat.data_ro has shape (num_cells_local,) - local_midpoints = f.dat.data_ro.reshape(num_cells_local, m.ufl_cell().geometric_dimension()) + local_midpoints = f.dat.data_ro.reshape(num_cells_local, m.geometric_dimension()) local_midpoints_size = np.array(local_midpoints.size) local_midpoints_sizes = np.empty(MPI.COMM_WORLD.size, dtype=int) MPI.COMM_WORLD.Allgatherv(local_midpoints_size, local_midpoints_sizes) - midpoints = np.empty((num_cells, m.ufl_cell().geometric_dimension()), dtype=local_midpoints.dtype) + midpoints = np.empty((num_cells, m.geometric_dimension()), dtype=local_midpoints.dtype) MPI.COMM_WORLD.Allgatherv(local_midpoints, (midpoints, local_midpoints_sizes)) assert len(np.unique(midpoints, axis=0)) == len(midpoints) return midpoints, local_midpoints diff --git a/tests/vertexonly/test_vertex_only_mesh_generation.py b/tests/vertexonly/test_vertex_only_mesh_generation.py index 559d46fef6..ac33b461de 100644 --- a/tests/vertexonly/test_vertex_only_mesh_generation.py +++ b/tests/vertexonly/test_vertex_only_mesh_generation.py @@ -26,11 +26,11 @@ def cell_midpoints(m): num_cells_local = len(f.dat.data_ro) num_cells = MPI.COMM_WORLD.allreduce(num_cells_local, op=MPI.SUM) # reshape is for 1D case where f.dat.data_ro has shape (num_cells_local,) - local_midpoints = f.dat.data_ro.reshape(num_cells_local, m.ufl_cell().geometric_dimension()) + local_midpoints = f.dat.data_ro.reshape(num_cells_local, m.geometric_dimension()) local_midpoints_size = np.array(local_midpoints.size) local_midpoints_sizes = np.empty(MPI.COMM_WORLD.size, dtype=int) MPI.COMM_WORLD.Allgatherv(local_midpoints_size, local_midpoints_sizes) - midpoints = np.empty((num_cells, m.ufl_cell().geometric_dimension()), dtype=local_midpoints.dtype) + midpoints = np.empty((num_cells, m.geometric_dimension()), dtype=local_midpoints.dtype) MPI.COMM_WORLD.Allgatherv(local_midpoints, (midpoints, local_midpoints_sizes)) assert len(np.unique(midpoints, axis=0)) == len(midpoints) return midpoints, local_midpoints