Support point/vertex cells (#30)

* Return zero() for determinant of Zero matrix

* don't scale integral over vertex

* allow empty elements in MixedElement

* update AUTHORS

AUTHORS

@@ -28,3 +28,4 @@ Contributors:
   | Corrado Maurini       <corrado.maurini@upmc.fr>
   | Jack S. Hale          <mail@jackhale.co.uk>
   | Tuomas Airaksinen     <tuomas.airaksinen@gmail.com>
+  | Reuben W. Hill        <reuben.hill10@imperial.ac.uk>

ufl/algorithms/apply_integral_scaling.py

@@ -28,11 +28,15 @@ def compute_integrand_scaling_factor(integral):
     # Polynomial degree of integrand scaling
     degree = 0
     if integral_type == "cell":
-        detJ = JacobianDeterminant(domain)
-        degree = estimate_total_polynomial_degree(apply_geometry_lowering(detJ))
-        # Despite the abs, |detJ| is polynomial except for
-        # self-intersecting cells, where we have other problems.
-        scale = abs(detJ) * weight
+        if tdim > 0:
+            detJ = JacobianDeterminant(domain)
+            degree = estimate_total_polynomial_degree(apply_geometry_lowering(detJ))
+            # Despite the abs, |detJ| is polynomial except for
+            # self-intersecting cells, where we have other problems.
+            scale = abs(detJ) * weight
+        else:
+            # No need to scale 'integral' over a vertex
+            scale = 1
 
     elif integral_type.startswith("exterior_facet"):
         if tdim > 1:

ufl/compound_expressions.py

@@ -13,6 +13,7 @@
 from ufl.core.multiindex import indices, Index
 from ufl.tensors import as_tensor, as_matrix, as_vector
 from ufl.operators import sqrt
+from ufl.constantvalue import Zero, zero
 
 
 # Note: To avoid typing errors, the expressions for cofactor and
@@ -82,7 +83,9 @@ def pseudo_inverse_expr(A):
 def determinant_expr(A):
     "Compute the (pseudo-)determinant of A."
     sh = A.ufl_shape
-    if sh == ():
+    if isinstance(A, Zero):
+        return zero()
+    elif sh == ():
         return A
     elif sh[0] == sh[1]:
         if sh[0] == 1:

ufl/finiteelement/mixedelement.py

@@ -56,9 +56,12 @@ def __init__(self, *elements, **kwargs):
 
         # Check that all elements use the same quadrature scheme TODO:
         # We can allow the scheme not to be defined.
-        quad_scheme = elements[0].quadrature_scheme()
-        if not all(e.quadrature_scheme() == quad_scheme for e in elements):
-            error("Quadrature scheme mismatch for sub elements of mixed element.")
+        if len(elements) == 0:
+            quad_scheme = None
+        else:
+            quad_scheme = elements[0].quadrature_scheme()
+            if not all(e.quadrature_scheme() == quad_scheme for e in elements):
+                error("Quadrature scheme mismatch for sub elements of mixed element.")
 
         # Compute value sizes in global and reference configurations
         value_size_sum = sum(product(s.value_shape()) for s in self._sub_elements)
@@ -292,10 +295,8 @@ def __init__(self, family, cell=None, degree=None, dim=None,
         # Initialize element data
         MixedElement.__init__(self, sub_elements, value_shape=value_shape,
                               reference_value_shape=reference_value_shape)
-        # FIXME: Storing this here is strange, isn't that handled by
-        # subclass?
-        self._family = sub_element.family()
-        self._degree = sub_element.degree()
+        FiniteElementBase.__init__(self, sub_element.family(), cell, sub_element.degree(), quad_scheme,
+                                   value_shape, reference_value_shape)
         self._sub_element = sub_element
 
         # Cache repr string