JohnsonMercier on manifolds

FIAT/expansions.py

@@ -92,7 +92,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
     dX = pad_jacobian(Jinv, pad_dim)
 
     phi0 = numpy.array([sum((ref_pts[i] - ref_pts[i] for i in range(dim)), 0.0)])
-    results = [numpy.zeros((num_members,) + (dim,)*k + phi0.shape[1:], dtype=phi0.dtype)
+    results = [numpy.zeros((num_members,) + (len(dX[0]),)*k + phi0.shape[1:], dtype=phi0.dtype)
                for k in range(order+1)]
 
     phi, dphi, ddphi = results + [None] * (2-order)
@@ -263,22 +263,23 @@ def __new__(cls, *args, **kwargs):
         reference element."""
         if cls is not ExpansionSet:
             return super().__new__(cls)
+        ref_el = args[0]
+        shape = ref_el.get_shape()
         try:
-            ref_el = args[0]
             expansion_set = {
                 reference_element.POINT: PointExpansionSet,
                 reference_element.LINE: LineExpansionSet,
                 reference_element.TRIANGLE: TriangleExpansionSet,
                 reference_element.TETRAHEDRON: TetrahedronExpansionSet,
-            }[ref_el.get_shape()]
-            return expansion_set(*args, **kwargs)
+            }[shape]
         except KeyError:
-            raise ValueError("Invalid reference element type.")
+            raise ValueError(f"Invalid reference element type {type(ref_el).__name__}.")
+        return expansion_set(*args, **kwargs)
 
     def __init__(self, ref_el, scale=None, variant=None):
         self.ref_el = ref_el
         self.variant = variant
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         top = ref_el.get_topology()
         base_ref_el = reference_element.default_simplex(sd)
         base_verts = base_ref_el.get_vertices()
@@ -303,7 +304,7 @@ def reconstruct(self, ref_el=None, scale=None, variant=None):
 
     def get_scale(self, n, cell=0):
         scale = self.scale
-        sd = self.ref_el.get_spatial_dimension()
+        sd = self.ref_el.get_topological_dimension()
         if isinstance(scale, str):
             vol = self.ref_el.volume_of_subcomplex(sd, cell)
             scale = scale.lower()
@@ -335,11 +336,12 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         ref_pts = numpy.add(numpy.dot(pts, A.T), b).T
         Jinv = A if direction is None else numpy.dot(A, direction)[:, None]
         sd = self.ref_el.get_spatial_dimension()
+        tdim = self.ref_el.get_topological_dimension()
         scale = self.get_scale(n, cell=cell)
-        phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
+        phi = dubiner_recurrence(tdim, n, lorder, ref_pts, Jinv,
                                  scale, variant=self.variant)
         if self.continuity == "C0":
-            phi = C0_basis(sd, n, phi)
+            phi = C0_basis(tdim, n, phi)
 
         # Pack linearly independent components into a dictionary
         result = {(0,) * sd: numpy.asarray(phi[0])}
@@ -418,7 +420,8 @@ def tabulate_normal_jumps(self, n, ref_pts, facet, order=0):
         :returns: a numpy array of tabulations of normal derivative jumps.
         """
         sd = self.ref_el.get_spatial_dimension()
-        transform = self.ref_el.get_entity_transform(sd-1, facet)
+        tdim = self.ref_el.get_topological_dimension()
+        transform = self.ref_el.get_entity_transform(tdim-1, facet)
         pts = transform(ref_pts)
         cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=False)
         cell_node_map = self.get_cell_node_map(n)
@@ -507,7 +510,7 @@ def get_dmats(self, degree, cell=0):
         if degree == 0:
             return cache.setdefault(key, numpy.zeros((self.ref_el.get_spatial_dimension(), 1, 1), "d"))
 
-        D = self.ref_el.get_dimension()
+        D = self.ref_el.get_topological_dimension()
         top = self.ref_el.get_topology()
         verts = self.ref_el.get_vertices_of_subcomplex(top[D][cell])
         pts = reference_element.make_lattice(verts, degree, variant="gl")
@@ -556,7 +559,7 @@ def __eq__(self, other):
 class PointExpansionSet(ExpansionSet):
     """Evaluates the point basis on a point reference element."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 0:
+        if ref_el.get_topological_dimension() != 0:
             raise ValueError("Must have a point")
         super().__init__(ref_el, **kwargs)
 
@@ -570,8 +573,8 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
 class LineExpansionSet(ExpansionSet):
     """Evaluates the Legendre basis on a line reference element."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 1:
-            raise Exception("Must have a line")
+        if ref_el.get_topological_dimension() != 1:
+            raise ValueError("Must have a line")
         super().__init__(ref_el, **kwargs)
 
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
@@ -600,16 +603,16 @@ class TriangleExpansionSet(ExpansionSet):
     """Evaluates the orthonormal Dubiner basis on a triangular
     reference element."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 2:
-            raise Exception("Must have a triangle")
+        if ref_el.get_topological_dimension() != 2:
+            raise ValueError("Must have a triangle")
         super().__init__(ref_el, **kwargs)
 
 
 class TetrahedronExpansionSet(ExpansionSet):
     """Collapsed orthonormal polynomial expansion on a tetrahedron."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 3:
-            raise Exception("Must be a tetrahedron")
+        if ref_el.get_topological_dimension() != 3:
+            raise ValueError("Must be a tetrahedron")
         super().__init__(ref_el, **kwargs)
 
 
@@ -627,7 +630,7 @@ def polynomial_dimension(ref_el, n, continuity=None):
     elif continuity == "C0":
         space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top)
     else:
-        dim = ref_el.get_spatial_dimension()
+        dim = ref_el.get_topological_dimension()
         space_dimension = math.comb(n + dim, dim) * len(top[dim])
     return space_dimension
 
@@ -641,7 +644,7 @@ def polynomial_entity_ids(ref_el, n, continuity=None):
     :returns: a dict of dicts mapping dimension and entity id to basis functions.
     """
     top = ref_el.get_topology()
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
     entity_ids = {}
     cur = 0
     for dim in sorted(top):
@@ -668,7 +671,7 @@ def polynomial_cell_node_map(ref_el, n, continuity=None):
     :returns: a numpy array mapping cell id to basis functions supported on that cell.
     """
     top = ref_el.get_topology()
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
 
     entity_ids = polynomial_entity_ids(ref_el, n, continuity)
     ref_entity_ids = polynomial_entity_ids(ref_el.construct_subelement(sd), n, continuity)
@@ -697,7 +700,7 @@ def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12):
     :returns: a dict mapping cell id to the point ids nearest to that cell.
     """
     top = ref_el.get_topology()
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
     if len(top[sd]) == 1:
         return {0: Ellipsis}
 

FIAT/finite_element.py

@@ -190,7 +190,7 @@ def tabulate(self, order, points, entity=None):
                      default cell-wise tabulation is performed.
         """
         if entity is None:
-            entity = (self.ref_el.get_spatial_dimension(), 0)
+            entity = (self.ref_el.get_topological_dimension(), 0)
 
         entity_dim, entity_id = entity
         transform = self.ref_el.get_entity_transform(entity_dim, entity_id)

FIAT/functional.py

@@ -361,7 +361,7 @@ class IntegralMomentOfNormalDerivative(IntegralMomentOfDerivative):
     def __init__(self, ref_el, facet_no, Q_face, f_at_qpts):
         n = ref_el.compute_normal(facet_no)
         # map points onto facet
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         Q = quadrature.FacetQuadratureRule(ref_el, sd-1, facet_no, Q_face, avg=True)
         super().__init__(ref_el, Q, f_at_qpts, n, nm="IntegralMomentOfNormalDerivative")
 

FIAT/johnson_mercier.py

@@ -13,13 +13,13 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
             raise ValueError(f"Johnson-Mercier does not have the {variant} variant")
         ref_el = ref_complex.get_parent()
         top = ref_el.get_topology()
-        sd = ref_el.get_spatial_dimension()
         entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
         nodes = []
 
         # Face dofs: bidirectional (nn and nt) moments
-        dim = sd - 1
-        R = numpy.array([[0, 1], [-1, 0]])
+        tdim = ref_el.get_topological_dimension()
+        dim = tdim - 1
+        n = list(map(ref_el.compute_scaled_normal, sorted(top[dim])))
         ref_facet = ref_el.construct_subelement(dim)
         Qref = parse_quadrature_scheme(ref_facet, 2*degree, quad_scheme)
         P = polynomial_set.ONPolynomialSet(ref_facet, degree)
@@ -28,9 +28,9 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
             cur = len(nodes)
             Q = FacetQuadratureRule(ref_el, dim, f, Qref, avg=True)
             thats = ref_el.compute_tangents(dim, f)
-            if sd == 2:
+            if tdim == 2:
                 # Face moments of sigma.n against n P1 and t P1
-                nhat = numpy.dot(R, *thats)
+                nhat = n[f]
                 components = (nhat, *thats)
             else:
                 # Face moments of sigma.n against n P1 and (n x t_j) P1
@@ -44,14 +44,14 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
 
         cur = len(nodes)
         # Interior dofs: moments for each independent component
-        n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
         Q = parse_quadrature_scheme(ref_complex, 2*degree-1, quad_scheme)
         P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
+        sd = ref_el.get_spatial_dimension()
         phis = P.tabulate(Q.get_points())[(0,) * sd]
         nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
-                     for phi in phis for i in range(sd) for j in range(i, sd))
+                     for phi in phis for i in range(tdim) for j in range(i, tdim))
 
-        entity_ids[sd][0].extend(range(cur, len(nodes)))
+        entity_ids[tdim][0].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 

FIAT/macro.py

@@ -30,15 +30,18 @@ def xy_to_bary(verts, pts, result=None):
     verts = numpy.asarray(verts)
     pts = numpy.asarray(pts)
     npts = pts.shape[0]
-    sdim = verts.shape[1]
+    nverts = verts.shape[0]
 
-    mat = numpy.vstack((verts.T, numpy.ones((1, sdim+1))))
+    A = numpy.vstack((verts.T, numpy.ones((1, nverts))))
     b = numpy.vstack((pts.T, numpy.ones((1, npts))))
-    foo = numpy.linalg.solve(mat, b)
+    if A.shape[0] == A.shape[1]:
+        bary = numpy.linalg.solve(A, b)
+    else:
+        bary, *_ = numpy.linalg.lstsq(A, b)
     if result is None:
-        return numpy.copy(foo.T)
+        result = bary.T.copy()
     else:
-        result[:, :] = foo[:, :].T
+        numpy.copyto(result, bary.T)
     return result
 
 
@@ -129,7 +132,7 @@ def __init__(self, parent, vertices, topology):
         self._child_to_parent = child_to_parent
         self._parent_to_children = parent_to_children
 
-        sd = parent.get_spatial_dimension()
+        sd = parent.get_topological_dimension()
         inv_top = invert_cell_topology(topology)
 
         # dict mapping cells to their boundary facets for each dimension,
@@ -186,7 +189,7 @@ def construct_subelement(self, dimension):
         return self.get_parent().construct_subelement(dimension)
 
     def get_facet_element(self):
-        dimension = self.get_spatial_dimension()
+        dimension = self.get_topological_dimension()
         return self.construct_subelement(dimension - 1)
 
     def is_macrocell(self):
@@ -256,7 +259,7 @@ class PowellSabinSplit(SplitSimplicialComplex):
     """
     def __init__(self, ref_el, dimension=1):
         self.split_dimension = dimension
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         top = ref_el.get_topology()
         connectivity = ref_el.get_connectivity()
         new_verts = list(ref_el.get_vertices())
@@ -293,7 +296,7 @@ def construct_subcomplex(self, dimension):
         """Constructs the reference subcomplex of the parent complex
         specified by subcomplex dimension.
         """
-        if dimension == self.get_dimension():
+        if dimension == self.get_topological_dimension():
             return self
         parent = self.get_parent_complex()
         subcomplex = parent.construct_subcomplex(dimension)
@@ -316,7 +319,7 @@ def __new__(cls, ref_el):
             return ref_el._split_cache.setdefault(cls, self)
 
     def __init__(self, ref_el):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         super().__init__(ref_el, dimension=sd)
 
 
@@ -332,7 +335,7 @@ def __new__(cls, ref_el):
             return ref_el._split_cache.setdefault(cls, self)
 
     def __init__(self, ref_el):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         super().__init__(ref_el, dimension=sd-1)
 
 
@@ -389,7 +392,7 @@ class MacroQuadratureRule(QuadratureRule):
     :returns: A quadrature rule on the sub entities of the simplicial complex.
     """
     def __init__(self, ref_el, Q_ref, parent_facets=None):
-        parent_dim = Q_ref.ref_el.get_spatial_dimension()
+        parent_dim = Q_ref.ref_el.get_topological_dimension()
         if parent_facets is not None:
             parent_to_children = ref_el.get_parent_to_children()
             facets = []
@@ -409,7 +412,7 @@ def __init__(self, ref_el, Q_ref, parent_facets=None):
 
         # Collapse repeated points if any of them lie on facets
         atol = 1E-10
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         top = ref_el.get_topology()
         for cell in top[sd]:
             bary = ref_el.compute_barycentric_coordinates(pts, entity=(sd, cell))
@@ -450,7 +453,7 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs):
         if not isinstance(order, (int, dict)):
             raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}")
 
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         if isinstance(order, int):
             order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)}
         if vorder is not None:
@@ -524,10 +527,17 @@ def hdiv_conforming_coefficients(U, order=0):
     ref_el = U.get_reference_element()
     coeffs = U.get_coeffs()
     shape = U.get_shape()
+
+    sd = ref_el.get_topological_dimension()
+    if sd != ref_el.get_spatial_dimension():
+        parent = ref_el.get_parent() or ref_el
+        thats = parent.compute_tangents(sd, 0)
+        mapping = "double contravariant piola" if len(shape) == 2 else "contravariant piola"
+        coeffs = pullback(coeffs, mapping, J=thats.T)
+        shape = coeffs.shape[1:-1]
+
     expansion_set = U.get_expansion_set()
     k = 1 if expansion_set.continuity == "C0" else 0
-
-    sd = ref_el.get_spatial_dimension()
     facet_el = ref_el.construct_subelement(sd-1)
 
     phi_deg = 0 if sd == 1 else degree - k
@@ -565,7 +575,7 @@ class HDivPolynomialSet(polynomial_set.PolynomialSet):
     :kwarg scale: The scale for the underlying ExpansionSet.
     """
     def __init__(self, ref_el, degree, order=0, **kwargs):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         U = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), **kwargs)
         coeffs = hdiv_conforming_coefficients(U, order=order)
         super().__init__(ref_el, degree, degree, U.expansion_set, coeffs)
@@ -613,9 +623,15 @@ def pullback(phi, mapping, J=None, Jinv=None, Jdet=None):
     if J is None:
         J = numpy.linalg.pinv(Jinv)
     if Jinv is None:
-        Jinv = numpy.linalg.pinv(J)
+        if J.shape[0] == J.shape[1]:
+            Jinv = numpy.linalg.inv(J)
+        else:
+            Jinv = numpy.linalg.pinv(J)
     if Jdet is None:
-        Jdet = numpy.linalg.det(J)
+        if J.shape[0] == J.shape[1]:
+            Jdet = numpy.linalg.det(J)
+        else:
+            Jdet = numpy.linalg.det(J.T @ J)**0.5
     F1 = Jinv.T
     F2 = J / Jdet
 
@@ -643,7 +659,7 @@ class MacroPolynomialSet(polynomial_set.PolynomialSet):
     """
     def __init__(self, ref_el, element):
         top = ref_el.get_topology()
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
 
         mapping, = set(element.mapping())
         base_ref_el = element.get_reference_element()