Merge pull request #25 from mscroggs/cross_normal_nedelec

Correct higher order Nedelec

Author: mscroggs

.github/workflows/run-tests.yml

@@ -21,7 +21,7 @@ jobs:
             sudo apt-get update
             sudo apt-get install -y libeigen3-dev
             pip install --upgrade pip
-            pip install scikit-build pytest numpy
+            pip install scikit-build pytest numpy scipy
             pip install pytest-xdist
         name: Install dependencies
       - uses: actions/checkout@v2

symfem/__init__.py

@@ -31,7 +31,8 @@
                 and issubclass(_element, _FiniteElement)
                 and _element != _FiniteElement
             ):
-                _elementlist.append(_element)
+                if _element not in _elementlist:
+                    _elementlist.append(_element)
                 for _n in _element.names:
                     if _n in _elementmap:
                         assert _element == _elementmap[_n]

symfem/core/finite_element.py

@@ -1,7 +1,8 @@
 """Abstract finite element classes and functions."""
 
 import sympy
-from .symbolic import x, zero, subs
+from .symbolic import x, zero, subs, sym_sum
+from . import mappings
 
 
 class FiniteElement:
@@ -20,6 +21,10 @@ def __init__(self, reference, basis, dofs, domain_dim, range_dim,
         self._basis_functions = None
         self._reshaped_basis_functions = None
 
+    def entity_dofs(self, entity_dim, entity_number):
+        """Get the numbers of the DOFs associated with the given entity."""
+        return [i for i, j in enumerate(self.dofs) if j.entity == (entity_dim, entity_number)]
+
     def get_polynomial_basis(self, reshape=True):
         """Get the polynomial basis for the element."""
         if reshape and self.range_shape is not None:
@@ -32,24 +37,26 @@ def get_polynomial_basis(self, reshape=True):
 
         return self.basis
 
+    def get_dual_matrix(self):
+        """Get the dual matrix."""
+        mat = []
+        for b in self.basis:
+            row = []
+            for d in self.dofs:
+                row.append(d.eval(b))
+            mat.append(row)
+        return sympy.Matrix(mat)
+
     def get_basis_functions(self, reshape=True):
         """Get the basis functions of the element."""
         if self._basis_functions is None:
-            mat = []
-            for b in self.basis:
-                row = []
-                for d in self.dofs:
-                    row.append(d.eval(b))
-                mat.append(row)
-            minv = sympy.Matrix(mat).inv("LU")
+            minv = self.get_dual_matrix().inv("LU")
             self._basis_functions = []
             if self.range_dim == 1:
                 # Scalar space
                 for i, dof in enumerate(self.dofs):
-                    b = zero
-                    for c, d in zip(minv.row(i), self.basis):
-                        b += c * d
-                    self._basis_functions.append(b)
+                    self._basis_functions.append(
+                        sym_sum(c * d for c, d in zip(minv.row(i), self.basis)))
             else:
                 # Vector space
                 for i, dof in enumerate(self.dofs):
@@ -108,80 +115,15 @@ def tabulate_basis(self, points, order="xyzxyz"):
             return output
         raise ValueError(f"Unknown order: {order}")
 
+    def map_to_cell(self, f, vertices):
+        """Map a function onto a cell using the appropriate mapping for the element."""
+        map = self.reference.get_map_to(vertices)
+        inverse_map = self.reference.get_inverse_map_to(vertices)
+        return getattr(mappings, self.mapping)(f, map, inverse_map, self.reference.tdim)
+
     @property
     def name(self):
         """Get the name of the element."""
         return self.names[0]
 
     names = []
-
-
-def make_integral_moment_dofs(
-    reference,
-    vertices=None, edges=None, faces=None, volumes=None,
-    cells=None, facets=None, ridges=None, peaks=None
-):
-    """Generate DOFs due to integral moments on sub entities.
-
-    Parameters
-    ----------
-    reference: symfem.references.Reference
-        The reference cell.
-    vertices: tuple
-        DOFs on dimension 0 entities.
-    edges: tuple
-        DOFs on dimension 1 entities.
-    faces: tuple
-        DOFs on dimension 2 entities.
-    volumes: tuple
-        DOFs on dimension 3 entities.
-    cells: tuple
-        DOFs on codimension 0 entities.
-    facets: tuple
-        DOFs on codimension 1 entities.
-    ridges: tuple
-        DOFs on codimension 2 entities.
-    peaks: tuple
-        DOFs on codimension 3 entities.
-    """
-    from symfem import create_reference
-    dofs = []
-
-    # DOFs per dimension
-    for dim, moment_data in enumerate([vertices, edges, faces, volumes]):
-        if moment_data is not None:
-            IntegralMoment, SubElement, order = moment_data
-            if order >= SubElement.min_order:
-                sub_type = reference.sub_entity_types[dim]
-                if sub_type is not None:
-                    assert dim > 0
-                    for i, vs in enumerate(reference.sub_entities(dim)):
-                        sub_ref = create_reference(
-                            sub_type,
-                            vertices=[reference.reference_vertices[v] for v in vs])
-                        sub_element = SubElement(sub_ref, order)
-                        for f, d in zip(
-                            sub_element.get_basis_functions(), sub_element.dofs
-                        ):
-                            dofs.append(IntegralMoment(sub_ref, f, d))
-
-    # DOFs per codimension
-    for _dim, moment_data in enumerate([peaks, ridges, facets, cells]):
-        dim = reference.tdim - 3 + _dim
-        if moment_data is not None:
-            IntegralMoment, SubElement, order = moment_data
-            if order >= SubElement.min_order:
-                sub_type = reference.sub_entity_types[dim]
-                if sub_type is not None:
-                    assert dim > 0
-                    for i, vs in enumerate(reference.sub_entities(dim)):
-                        sub_ref = create_reference(
-                            sub_type,
-                            vertices=[reference.reference_vertices[v] for v in vs])
-                        sub_element = SubElement(sub_ref, order)
-                        for f, d in zip(
-                            sub_element.get_basis_functions(), sub_element.dofs
-                        ):
-                            dofs.append(IntegralMoment(sub_ref, f, d))
-
-    return dofs

symfem/core/functionals.py

@@ -7,6 +7,9 @@
 class BaseFunctional:
     """A functional."""
 
+    def __init__(self, entity=(None, None)):
+        self.entity = entity
+
     def eval(self, fun):
         """Apply to the functional to a function."""
         raise NotImplementedError
@@ -21,17 +24,17 @@ def dof_direction(self):
 
     def entity_dim(self):
         """Get the dimension of the entitiy this DOF is associated with."""
-        return None
+        return self.entity[0]
 
     name = None
 
 
 class PointEvaluation(BaseFunctional):
     """A point evaluation."""
 
-    def __init__(self, point, entity_dim=None):
+    def __init__(self, point, entity=(None, None)):
+        super().__init__(entity)
         self.point = point
-        self._entity_dim = entity_dim
 
     def eval(self, function):
         """Apply to the functional to a function."""
@@ -41,20 +44,16 @@ def dof_point(self):
         """Get the location of the DOF in the cell."""
         return self.point
 
-    def entity_dim(self):
-        """Get the dimension of the entitiy this DOF is associated with."""
-        return self._entity_dim
-
     name = "Point evaluation"
 
 
 class PointDirectionalDerivativeEvaluation(BaseFunctional):
     """A point evaluation of a derivative in a fixed direction."""
 
-    def __init__(self, point, direction, entity_dim=None):
+    def __init__(self, point, direction, entity=(None, None)):
+        super().__init__(entity)
         self.point = point
         self.dir = direction
-        self._entity_dim = entity_dim
 
     def eval(self, function):
         """Apply to the functional to a function."""
@@ -68,18 +67,14 @@ def dof_direction(self):
         """Get the direction of the DOF."""
         return self.dir
 
-    def entity_dim(self):
-        """Get the dimension of the entitiy this DOF is associated with."""
-        return self._entity_dim
-
     name = "Point evaluation of directional derivative"
 
 
 class PointNormalDerivativeEvaluation(PointDirectionalDerivativeEvaluation):
     """A point evaluation of a normal derivative."""
 
-    def __init__(self, point, edge):
-        super().__init__(point, edge.normal(), entity_dim=edge.tdim)
+    def __init__(self, point, edge, entity=(None, None)):
+        super().__init__(point, edge.normal(), entity=entity)
         self.reference = edge
 
     name = "Point evaluation of normal derivative"
@@ -88,10 +83,10 @@ def __init__(self, point, edge):
 class PointComponentSecondDerivativeEvaluation(BaseFunctional):
     """A point evaluation of a component of a second derivative."""
 
-    def __init__(self, point, component, entity_dim=None):
+    def __init__(self, point, component, entity=(None, None)):
+        super().__init__(entity)
         self.point = point
         self.component = component
-        self._entity_dim = entity_dim
 
     def eval(self, function):
         """Apply to the functional to a function."""
@@ -101,20 +96,16 @@ def dof_point(self):
         """Get the location of the DOF in the cell."""
         return self.point
 
-    def entity_dim(self):
-        """Get the dimension of the entitiy this DOF is associated with."""
-        return self._entity_dim
-
     name = "Point evaluation of Jacobian component"
 
 
 class PointInnerProduct(BaseFunctional):
     """An evaluation of an inner product at a point."""
 
-    def __init__(self, point, vec, entity_dim=None):
+    def __init__(self, point, vec, entity=(None, None)):
+        super().__init__(entity)
         self.point = point
         self.vec = vec
-        self._entity_dim = entity_dim
 
     def eval(self, function):
         """Apply to the functional to a function."""
@@ -132,20 +123,16 @@ def dof_direction(self):
         """Get the location of the DOF in the cell."""
         return self.vec
 
-    def entity_dim(self):
-        """Get the dimension of the entitiy this DOF is associated with."""
-        return self._entity_dim
-
     name = "Point inner product"
 
 
 class DotPointEvaluation(BaseFunctional):
     """A point evaluation in a given direction."""
 
-    def __init__(self, point, vector, entity_dim=None):
+    def __init__(self, point, vector, entity=(None, None)):
+        super().__init__(entity)
         self.point = point
         self.vector = vector
-        self._entity_dim = entity_dim
 
     def eval(self, function):
         """Apply to the functional to a function."""
@@ -159,23 +146,22 @@ def dof_direction(self):
         """Get the direction of the DOF."""
         return self.vector
 
-    def entity_dim(self):
-        """Get the dimension of the entitiy this DOF is associated with."""
-        return self._entity_dim
-
     name = "Dot point evaluation"
 
 
 class IntegralMoment(BaseFunctional):
     """An integral moment."""
 
-    def __init__(self, reference, f, dof):
+    def __init__(self, reference, f, dof, entity=(None, None)):
+        super().__init__(entity)
         self.reference = reference
         self.dof = dof
         self.f = subs(f, x, t)
         if isinstance(self.f, tuple):
+            # TODO: is this one of the mappings?
             self.f = tuple(
-                sum(self.reference.scaled_axes()[j][i] * c for j, c in enumerate(self.f))
+                sum(self.reference.axes[j][i] * c / self.reference.jacobian()
+                    for j, c in enumerate(self.f))
                 for i, o in enumerate(self.reference.origin)
             )
 
@@ -210,18 +196,14 @@ def dof_direction(self):
             for i in range(self.reference.gdim)
         )
 
-    def entity_dim(self):
-        """Get the dimension of the entitiy this DOF is associated with."""
-        return self.reference.tdim
-
     name = "Integral moment"
 
 
 class VecIntegralMoment(IntegralMoment):
     """An integral moment applied to a component of a vector."""
 
-    def __init__(self, reference, f, dot_with, dof):
-        super().__init__(reference, f, dof)
+    def __init__(self, reference, f, dot_with, dof, entity=(None, None)):
+        super().__init__(reference, f, dof, entity=entity)
         self.dot_with = dot_with
 
     def dot(self, function):
@@ -238,16 +220,16 @@ def dof_direction(self):
 class TangentIntegralMoment(VecIntegralMoment):
     """An integral moment in the tangential direction."""
 
-    def __init__(self, reference, f, dof):
-        super().__init__(reference, f, reference.tangent(), dof)
+    def __init__(self, reference, f, dof, entity=(None, None)):
+        super().__init__(reference, f, reference.tangent(), dof, entity=entity)
 
     name = "Tangential integral moment"
 
 
 class NormalIntegralMoment(VecIntegralMoment):
     """An integral moment in the normal direction."""
 
-    def __init__(self, reference, f, dof):
-        super().__init__(reference, f, reference.normal(), dof)
+    def __init__(self, reference, f, dof, entity=(None, None)):
+        super().__init__(reference, f, reference.normal(), dof, entity=entity)
 
     name = "Normal integral moment"