Allow mesh to have multiple coordinate elements

.github/workflows/fenicsx-tests.yml

@@ -45,7 +45,7 @@ jobs:
         with:
           path: ./ffcx
           repository: FEniCS/ffcx
-          ref: main
+          ref: mscroggs/ufl_coordinate_elements
 
       - name: Install FFCx
         run: |
@@ -76,14 +76,14 @@ jobs:
       - name: Install Basix and FFCx
         run: |
           python3 -m pip install --break-system-packages  git+https://github.com/FEniCS/basix.git
-          python3 -m pip install --break-system-packages  git+https://github.com/FEniCS/ffcx.git
+          python3 -m pip install --break-system-packages  git+https://github.com/FEniCS/ffcx.git@mscroggs/ufl_coordinate_elements
 
       - name: Clone DOLFINx
         uses: actions/checkout@v4
         with:
           path: ./dolfinx
           repository: FEniCS/dolfinx
-          ref: main
+          ref: mscroggs/ufl_coordinate_elements
       - name: Install DOLFINx
         run: |
           cmake -G Ninja -DCMAKE_BUILD_TYPE=Developer -B build -S dolfinx/cpp/

test/test_derivative.py

@@ -107,11 +107,11 @@ def make_value(c):
     acell = adomain.ufl_cell()
     bcell = bdomain.ufl_cell()
     assert acell == bcell
-    if adomain.geometric_dimension() == 1:
+    if adomain.geometric_dimension == 1:
         x = (0.3,)
-    elif adomain.geometric_dimension() == 2:
+    elif adomain.geometric_dimension == 2:
         x = (0.3, 0.4)
-    elif adomain.geometric_dimension() == 3:
+    elif adomain.geometric_dimension == 3:
         x = (0.3, 0.4, 0.5)
     av = a(x, amapping)
     bv = b(x, bmapping)

test/test_domains.py

@@ -68,18 +68,19 @@ def test_domains_sort_by_name():
         )
         for cell in sorted(all_cells)
     ]
-    sdomains = sorted(domains1, key=lambda D: (D.topological_dimension(), D.ufl_cell(), D.ufl_id()))
+    sdomains = sorted(
+        domains1, key=lambda D: (D.ufl_cell().topological_dimension(), D.ufl_cell(), D.ufl_id())
+    )
     assert sdomains != domains1
     assert sdomains == domains2
 
 
 def test_topdomain_creation():
     D = Mesh(LagrangeElement(interval, 1, (1,)))
-    assert D.geometric_dimension() == 1
+    assert D.geometric_dimension == 1
     D = Mesh(LagrangeElement(triangle, 1, (2,)))
-    assert D.geometric_dimension() == 2
-    D = Mesh(LagrangeElement(tetrahedron, 1, (3,)))
-    assert D.geometric_dimension() == 3
+    assert D.geometric_dimension == 2
+    D = Mesh(LagrangeElement(triangle, 1, (3,)))
 
 
 def test_cell_legacy_case():
@@ -377,7 +378,7 @@ def test_merge_sort_integral_data():
 
 
 def test_extract_domains():
-    "Test that the domains are extracted properly from a mixed-domain expression"
+    """Test that the domains are extracted properly from a mixed-domain expression."""
 
     # Create domains of different topological dimensions
     gdim = 2
@@ -404,5 +405,29 @@ def test_extract_domains():
 
     domains = extract_domains(expr)
 
-    assert domains[0] == dom_1
-    assert domains[1] == dom_0
+    assert domains[0] == dom_0
+    assert domains[1] == dom_1
+
+
+def test_hybrid_mesh():
+    """Test meshes with multiple cell types."""
+    domain = Mesh(
+        [
+            LagrangeElement(triangle, 1, (2,)),
+            LagrangeElement(quadrilateral, 1, (2,)),
+        ]
+    )
+
+    elements = [LagrangeElement(triangle, 1), LagrangeElement(quadrilateral, 1)]
+
+    space = FunctionSpace(domain, elements)
+
+    # Create test and trial functions
+    u = TrialFunction(space)
+    v = TestFunction(space)
+
+    # Create an expression
+    expr = u * v
+
+    # Create an integral
+    expr * dx

test/test_functionspace.py

@@ -0,0 +1,27 @@
+"""Tests of function spaces."""
+
+import pytest
+from utils import LagrangeElement
+
+from ufl import FunctionSpace, Mesh, quadrilateral, triangle
+
+
+def test_cell_mismatch():
+    domain = Mesh(LagrangeElement(triangle, 1, (2,)))
+    elements = LagrangeElement(quadrilateral, 1)
+
+    with pytest.raises(ValueError):
+        FunctionSpace(domain, elements)
+
+
+def test_wrong_order():
+    domain = Mesh(
+        [
+            LagrangeElement(quadrilateral, 1, (2,)),
+            LagrangeElement(triangle, 1, (2,)),
+        ]
+    )
+    elements = [LagrangeElement(triangle, 1), LagrangeElement(quadrilateral, 1)]
+
+    with pytest.raises(ValueError):
+        FunctionSpace(domain, elements)

test/test_measures.py

@@ -75,8 +75,8 @@ def test_foo():
 
     assert cell.topological_dimension() == tdim
     assert cell.cellname() == "triangle"
-    assert mydomain.topological_dimension() == tdim
-    assert mydomain.geometric_dimension() == gdim
+    assert mydomain.ufl_cell().topological_dimension() == tdim
+    assert mydomain.geometric_dimension == gdim
     assert mydomain.ufl_cell() == cell
     assert mydomain.ufl_id() == 9
     assert mydomain.ufl_cargo() == mymesh

test/test_piecewise_checks.py

@@ -280,14 +280,14 @@ def mappings_are_cellwise_constant(domain, test):
     assert is_cellwise_constant(e) == test
     e = JacobianInverse(domain)
     assert is_cellwise_constant(e) == test
-    if domain.topological_dimension() != 1:
+    if domain.ufl_cell().topological_dimension() != 1:
         e = FacetJacobian(domain)
         assert is_cellwise_constant(e) == test
         e = FacetJacobianDeterminant(domain)
         assert is_cellwise_constant(e) == test
         e = FacetJacobianInverse(domain)
         assert is_cellwise_constant(e) == test
-    if domain.topological_dimension() > 2:
+    if domain.ufl_cell().topological_dimension() > 2:
         e = RidgeJacobian(domain)
         assert is_cellwise_constant(e) == test
         e = RidgeJacobianDeterminant(domain)

ufl/__init__.py

@@ -63,7 +63,6 @@
     -AbstractDomain
     -Mesh
     -MeshSequence
-    -MeshView
 
 * Sobolev spaces::
 
@@ -266,7 +265,7 @@
 from ufl.core.external_operator import ExternalOperator
 from ufl.core.interpolate import Interpolate, interpolate
 from ufl.core.multiindex import Index, indices
-from ufl.domain import AbstractDomain, Mesh, MeshSequence, MeshView
+from ufl.domain import AbstractDomain, Mesh, MeshSequence
 from ufl.finiteelement import AbstractFiniteElement
 from ufl.form import BaseForm, Form, FormSum, ZeroBaseForm
 from ufl.formoperators import (
@@ -494,7 +493,6 @@
     "Measure",
     "Mesh",
     "MeshSequence",
-    "MeshView",
     "MinCellEdgeLength",
     "MinFacetEdgeLength",
     "MixedFunctionSpace",

ufl/algorithms/apply_geometry_lowering.py

@@ -118,7 +118,7 @@ def jacobian_determinant(self, o):
         # TODO: Is "signing" the determinant for manifolds the
         #       cleanest approach?  The alternative is to have a
         #       specific type for the unsigned pseudo-determinant.
-        if domain.topological_dimension() < domain.geometric_dimension():
+        if domain.ufl_cell().topological_dimension() < domain.geometric_dimension:
             co = CellOrientation(domain)
             detJ = co * detJ
         return detJ
@@ -265,7 +265,7 @@ def facet_area(self, o):
             return o
 
         domain = extract_unique_domain(o)
-        tdim = domain.topological_dimension()
+        tdim = domain.ufl_cell().topological_dimension()
         if not domain.is_piecewise_linear_simplex_domain():
             # Don't lower for non-affine cells, instead leave it to
             # form compiler
@@ -420,8 +420,8 @@ def cell_normal(self, o):
             return o
 
         domain = extract_unique_domain(o)
-        gdim = domain.geometric_dimension()
-        tdim = domain.topological_dimension()
+        gdim = domain.geometric_dimension
+        tdim = domain.ufl_cell().topological_dimension()
 
         if tdim == gdim - 1:  # n-manifold embedded in n-1 space
             i = Index()
@@ -453,15 +453,15 @@ def facet_normal(self, o):
             return o
 
         domain = extract_unique_domain(o)
-        tdim = domain.topological_dimension()
+        tdim = domain.ufl_cell().topological_dimension()
 
         if tdim == 1:
             # Special-case 1D (possibly immersed), for which we say
             # that n is just in the direction of J.
             J = self.jacobian(Jacobian(domain))  # dx/dX
             ndir = J[:, 0]
 
-            gdim = domain.geometric_dimension()
+            gdim = domain.geometric_dimension
             if gdim == 1:
                 nlen = abs(ndir[0])
             else:

ufl/algorithms/apply_integral_scaling.py

@@ -24,10 +24,9 @@ def compute_integrand_scaling_factor(integral):
     """Change integrand geometry to the right representations."""
     domain = integral.ufl_domain()
     integral_type = integral.integral_type()
-    # co = CellOrientation(domain)
     weight = QuadratureWeight(domain)
-    tdim = domain.topological_dimension()
-    # gdim = domain.geometric_dimension()
+    assert len(domain.ufl_coordinate_elements()) == 1  # TODO: remove this assumption
+    tdim = domain.ufl_coordinate_elements()[0].cell.topological_dimension()
 
     # Polynomial degree of integrand scaling
     degree = 0

ufl/algorithms/apply_restrictions.py

@@ -175,8 +175,8 @@ def facet_normal(self, o):
         """Restrict a facet_normal."""
         D = extract_unique_domain(o)
         e = D.ufl_coordinate_element()
-        gd = D.geometric_dimension()
-        td = D.topological_dimension()
+        gd = D.geometric_dimension
+        td = D.ufl_cell().topological_dimension()
 
         if e.embedded_superdegree <= 1 and e in H1 and gd == td:
             # For meshes with a continuous linear non-manifold

ufl/algorithms/compute_form_data.py

@@ -397,7 +397,7 @@ def compute_form_data(
     self.rank = len(self.original_form.arguments())
 
     # Extract common geometric dimension (topological is not common!)
-    self.geometric_dimension = self.original_form.integrals()[0].ufl_domain().geometric_dimension()
+    self.geometric_dimension = self.original_form.integrals()[0].ufl_domain().geometric_dimension
 
     # --- Build mapping from old incomplete element objects to new
     # well defined elements.  This is to support the Expression

ufl/algorithms/strip_terminal_data.py

@@ -12,7 +12,6 @@
     FunctionSpace,
     Integral,
     Mesh,
-    MeshView,
     MixedFunctionSpace,
     TensorProductFunctionSpace,
 )
@@ -126,9 +125,5 @@ def strip_domain(domain):
     """Return a new domain with all non-UFL information removed."""
     if isinstance(domain, Mesh):
         return Mesh(domain.ufl_coordinate_element(), domain.ufl_id())
-    elif isinstance(domain, MeshView):
-        return MeshView(
-            strip_domain(domain.ufl_mesh()), domain.topological_dimension(), domain.ufl_id()
-        )
     else:
         raise NotImplementedError(f"{type(domain)} cannot be stripped")

ufl/constant.py

@@ -78,12 +78,12 @@ def _ufl_signature_data_(self, renumbering):
 def VectorConstant(domain, count=None):
     """Vector constant."""
     domain = as_domain(domain)
-    return Constant(domain, shape=(domain.geometric_dimension(),), count=count)
+    return Constant(domain, shape=(domain.geometric_dimension,), count=count)
 
 
 def TensorConstant(domain, count=None):
     """Tensor constant."""
     domain = as_domain(domain)
     return Constant(
-        domain, shape=(domain.geometric_dimension(), domain.geometric_dimension()), count=count
+        domain, shape=(domain.geometric_dimension, domain.geometric_dimension), count=count
     )

ufl/differentiation.py

@@ -318,15 +318,23 @@ def __new__(cls, f):
         # Return zero if expression is trivially constant
         if is_cellwise_constant(f):
             # TODO: Use max topological dimension if there are multiple topological dimensions.
-            dim = extract_unique_domain(f, expand_mesh_sequence=False).topological_dimension()
+            domain_elements = extract_unique_domain(
+                f, expand_mesh_sequence=False
+            ).ufl_coordinate_elements()
+            dim = domain_elements[0].cell.topological_dimension()
+            for e in domain_elements:
+                # TODO: remove this assumption
+                assert e.cell.topological_dimension() == dim
             return Zero(f.ufl_shape + (dim,), f.ufl_free_indices, f.ufl_index_dimensions)
         return CompoundDerivative.__new__(cls)
 
     def __init__(self, f):
         """Initalise."""
         CompoundDerivative.__init__(self, (f,))
-        # TODO: Use max topological dimension if there are multiple topological dimensions.
-        self._dim = extract_unique_domain(f, expand_mesh_sequence=False).topological_dimension()
+        self._dim = max(
+            e.cell.topological_dimension()
+            for e in extract_unique_domain(f, expand_mesh_sequence=False).ufl_coordinate_elements()
+        )
 
     def _ufl_expr_reconstruct_(self, op):
         """Return a new object of the same type with new operands."""