JohnsonMercier on manifolds

Author: pbrubeck

FIAT/expansions.py

@@ -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()
@@ -334,7 +335,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         A, b = self.affine_mappings[cell]
         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()
+        sd = self.ref_el.get_topological_dimension()
         scale = self.get_scale(n, cell=cell)
         phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
                                  scale, variant=self.variant)
@@ -417,7 +418,7 @@ 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()
+        sd = self.ref_el.get_topological_dimension()
         transform = self.ref_el.get_entity_transform(sd-1, facet)
         pts = transform(ref_pts)
         cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=False)
@@ -507,7 +508,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")
@@ -519,14 +520,14 @@ def get_dmats(self, degree, cell=0):
     def tabulate(self, n, pts):
         if len(pts) == 0:
             return numpy.array([])
-        sd = self.ref_el.get_spatial_dimension()
+        sd = self.ref_el.get_topological_dimension()
         return self._tabulate(n, pts)[(0,) * sd]
 
     def tabulate_derivatives(self, n, pts):
         from FIAT.polynomial_set import mis
         vals = self._tabulate(n, pts, order=1)
         # Create the ordinary data structure.
-        sd = self.ref_el.get_spatial_dimension()
+        sd = self.ref_el.get_topological_dimension()
         v = vals[(0,) * sd]
         dv = [vals[alpha] for alpha in mis(sd, 1)]
         data = [[(v[i, j], [vi[i, j] for vi in dv])
@@ -537,7 +538,7 @@ def tabulate_derivatives(self, n, pts):
     def tabulate_jet(self, n, pts, order=1):
         vals = self._tabulate(n, pts, order=order)
         # Create the ordinary data structure.
-        sd = self.ref_el.get_spatial_dimension()
+        sd = self.ref_el.get_topological_dimension()
         v0 = vals[(0,) * sd]
         data = [v0]
         for r in range(1, order+1):
@@ -556,7 +557,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 +571,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,15 +601,15 @@ 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:
+        if ref_el.get_topological_dimension() != 2:
             raise Exception("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:
+        if ref_el.get_topological_dimension() != 3:
             raise Exception("Must be a tetrahedron")
         super().__init__(ref_el, **kwargs)
 
@@ -627,7 +628,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 +642,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 +669,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 +698,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/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()
+        sd = ref_el.get_topological_dimension()
         entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
         nodes = []
 
         # Face dofs: bidirectional (nn and nt) moments
+        n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
         dim = sd - 1
-        R = numpy.array([[0, 1], [-1, 0]])
         ref_facet = ref_el.construct_subelement(dim)
         Qref = parse_quadrature_scheme(ref_facet, 2*degree, quad_scheme)
         P = polynomial_set.ONPolynomialSet(ref_facet, degree)
@@ -30,7 +30,7 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
             thats = ref_el.compute_tangents(dim, f)
             if sd == 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,7 +44,6 @@ 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")
         phis = P.tabulate(Q.get_points())[(0,) * sd]

FIAT/macro.py

@@ -30,11 +30,15 @@ def xy_to_bary(verts, pts, result=None):
     verts = numpy.asarray(verts)
     pts = numpy.asarray(pts)
     npts = pts.shape[0]
-    sdim = verts.shape[1]
+    sdim = verts.shape[0]
 
-    mat = numpy.vstack((verts.T, numpy.ones((1, sdim+1))))
+    mat = numpy.vstack((verts.T, numpy.ones((1, sdim))))
     b = numpy.vstack((pts.T, numpy.ones((1, npts))))
-    foo = numpy.linalg.solve(mat, b)
+
+    if mat.shape[0] == mat.shape[1]:
+        foo = numpy.linalg.solve(mat, b)
+    else:
+        foo, *_ = numpy.linalg.lstsq(mat, b)
     if result is None:
         return numpy.copy(foo.T)
     else:
@@ -129,7 +133,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 +190,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 +260,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 +297,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 +320,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 +336,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 +393,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 +413,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 +454,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 +528,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 +576,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 +624,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 +660,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()

FIAT/polynomial_set.py

@@ -224,7 +224,7 @@ class ONSymTensorPolynomialSet(PolynomialSet):
     def __init__(self, ref_el, degree, size=None, **kwargs):
         expansion_set = expansions.ExpansionSet(ref_el, **kwargs)
 
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         if size is None:
             size = sd
 

FIAT/quadrature.py

@@ -174,7 +174,7 @@ class CollapsedQuadratureSimplexRule(QuadratureRule):
     from the hypercube to the simplex."""
 
     def __init__(self, ref_el, m):
-        dim = ref_el.get_spatial_dimension()
+        dim = ref_el.get_topological_dimension()
         Ref1 = reference_element.default_simplex(dim)
         pts_ref, wts_ref = simplexgausslegendre(dim, m)
         pts, wts = map_quadrature(pts_ref, wts_ref, Ref1, ref_el)