add volume scaling for integral dots

.github/workflows/test.yml

@@ -24,7 +24,7 @@ jobs:
         arch: [default]
     runs-on: [self-hosted, Linux]
     container:
-      image: ubuntu:latest
+      image: ubuntu:25.10
     env:
       OMPI_ALLOW_RUN_AS_ROOT: 1
       OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1

fuse/cells.py

@@ -839,6 +839,34 @@ def quadrature(self, degree):
         pts, wts = Q.get_points(), Q.get_weights()
         return pts, wts
 
+    def volume(self):
+        vertices = np.asarray(self.ordered_vertex_coords())
+        if self.get_spatial_dimension() == 0:
+            return 1
+        elif self.get_spatial_dimension() == 1:
+            return abs(vertices[1] - vertices[0])[0]
+        elif self.get_spatial_dimension() == 2:
+            x = vertices[:, 0]
+            y = vertices[:, 1]
+            return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
+        elif self.get_spatial_dimension() == 3:
+            vertices = np.asarray(vertices)
+            V = 0.0
+            for face in self.d_entities(2):
+                pts = np.array([self.get_node(v, return_coords=True) for v in face.ordered_vertices()])
+                c = pts.mean(axis=0)
+
+                n = np.zeros(3)
+                for i in range(len(pts)):
+                    v0 = pts[i]
+                    v1 = pts[(i + 1) % len(pts)]
+                    n += np.cross(v0, v1)
+                V += np.dot(c, n)
+            V = abs(V) / 3.0
+            return V
+        else:
+            raise NotImplementedError("Dimension not accounted for")
+
     def cartesian_to_barycentric(self, pts):
         verts = np.array(self.ordered_vertex_coords())
         v_0 = self.ordered_vertex_coords()[0]

fuse/dof.py

@@ -35,7 +35,7 @@ def __call__(self, kernel, v, cell):
         return v(*kernel.pt)
 
     def tabulate(self):
-        return 1
+        return np.eye(self.entity.dim())
 
     def add_entity(self, entity):
         res = DeltaPairing()
@@ -80,8 +80,8 @@ def tabulate(self):
         if self.orientation:
             new_bvs = np.array(self.entity.orient(self.orientation).basis_vectors())
             basis_change = np.matmul(np.linalg.inv(new_bvs), bvs)
-            return basis_change
-        return np.eye(bvs.shape[0])
+            return (1/self.entity.volume())*basis_change
+        return (1/self.entity.volume())*np.eye(bvs.shape[0])
 
     def add_entity(self, entity):
         res = L2Pairing()
@@ -423,13 +423,7 @@ def convert_to_fiat(self, ref_el, interpolant_degree, value_shape=tuple()):
     def to_quadrature(self, arg_degree, value_shape):
         Qpts, Qwts = self.cell_defined_on.quadrature(self.kernel.degree(arg_degree))
         Qwts = Qwts.reshape(Qwts.shape + (1,))
-        dim = self.cell_defined_on.get_spatial_dimension()
-        if dim > 0:
-            bvs = np.array(self.cell_defined_on.basis_vectors())
-            new_bvs = np.array(self.cell_defined_on.orient(self.pairing.orientation).basis_vectors())
-            basis_change = np.matmul(np.linalg.inv(new_bvs), bvs)
-        else:
-            basis_change = np.eye(dim)
+        basis_change = self.pairing.tabulate()
 
         if self.immersed and (isinstance(self.kernel, VectorKernel) or isinstance(self.kernel, BarycentricPolynomialKernel) or isinstance(self.kernel, PolynomialKernel)):
             def immersed(pt):

test/test_2d_examples_docs.py

@@ -344,6 +344,7 @@ def test_nd_example():
     for dof in ned.generate():
         assert [np.allclose(1, dof.eval(basis_func).flatten()) for basis_func in basis_funcs].count(True) == 1
         assert [np.allclose(0, dof.eval(basis_func).flatten()) for basis_func in basis_funcs].count(True) == 2
+    ned.to_fiat()
 
 
 def construct_rt(tri=None):

test/test_cells.py

@@ -7,7 +7,7 @@
 from test_convert_to_fiat import helmholtz_solve
 
 
-@pytest.fixture(scope='module', params=[0, 1, 2])
+@pytest.fixture(scope='module', params=[0, 1, 2, 3])
 def C(request):
     dim = request.param
     if dim == 0:
@@ -16,6 +16,8 @@ def C(request):
         return Point(1, [Point(0), Point(0)], vertex_num=2)
     elif dim == 2:
         return polygon(3)
+    elif dim == 3:
+        return make_tetrahedron()
 
 
 def test_vertices(C):