diff --git a/elastodynamicsx/pde/materials/elasticmaterial.py b/elastodynamicsx/pde/materials/elasticmaterial.py
index 0fa2a44..26a48dd 100644
--- a/elastodynamicsx/pde/materials/elasticmaterial.py
+++ b/elastodynamicsx/pde/materials/elasticmaterial.py
@@ -4,7 +4,7 @@
 #
 # SPDX-License-Identifier: MIT
 
-from dolfinx import fem
+from dolfinx import fem, default_scalar_type
 from petsc4py import PETSc
 import ufl
 import numpy as np
@@ -222,8 +222,8 @@ def DG_numerical_flux(self) -> 'function':
     def DG_SIPG_regularization_parameter(self) -> 'dolfinx.fem.Constant':
         """Regularization parameter for the Symmetric Interior Penalty Galerkin methods (SIPG)"""
         degree = self._function_space.ufl_element().degree()
-        gamma  = fem.Constant(self._function_space.mesh, PETSc.ScalarType(degree*(degree+1) + 1)) #+1 otherwise blows with elements of degree 1
-        #gamma  = fem.Constant(self._function_space.mesh, PETSc.ScalarType(160)) #+1 otherwise blows with elements of degree 1
+        gamma  = fem.Constant(self._function_space.mesh, default_scalar_type(degree*(degree+1) + 1)) #+1 otherwise blows with elements of degree 1
+        #gamma  = fem.Constant(self._function_space.mesh, default_scalar_type(160)) #+1 otherwise blows with elements of degree 1
         P_mod  = self.P_modulus
         R_     = gamma*P_mod
         return R_
@@ -476,23 +476,23 @@ def k1_CG(self) -> 'function':
 
     @property
     def k2_CG(self) -> 'function':
-        return lambda u,v: fem.Constant(u.ufl_function_space(), PETSc.ScalarType(0)) * ufl.inner(u,v) * self._dx
+        return lambda u,v: fem.Constant(u.ufl_function_space().mesh, default_scalar_type(0)) * ufl.inner(u,v) * self._dx
 
     @property
     def k3_CG(self) -> 'function':
         return lambda u,v: ufl.inner(self.mu*self._L_onaxis(u), self._L_onaxis(v)) * self._dx
-        
+
     @property
     def DG_numerical_flux_SIPG(self) -> 'function':
         inner, avg, jump = ufl.inner, ufl.avg, ufl.jump
         V = self._function_space
         n = ufl.FacetNormal(V)
-        h = ufl.MinCellEdgeLength(V) #works!
+        h = ufl.MinCellEdgeLength(V)  # works!
         h_avg  = (h('+') + h('-'))/2.0
         R_ = self.DG_SIPG_regularization_parameter()
         dS = self._dS
         sigma = self.sigma
-        
+
         k_int_facets = lambda u,v: \
                        -           inner(avg(sigma(u)), jump(v,n)    ) * dS \
                        -           inner(jump(u,n)    , avg(sigma(v))) * dS \
@@ -510,7 +510,7 @@ def DG_numerical_flux_NIPG(self) -> 'function':
         R_ = self.DG_SIPG_regularization_parameter()
         dS = self._dS
         sigma = self.sigma
-        
+
         k_int_facets = lambda u,v: \
                        -           inner(avg(sigma(u)), jump(v,n)    ) * dS \
                        +           inner(jump(u,n)    , avg(sigma(v))) * dS \
diff --git a/test/test_materials.py b/test/test_materials.py
index 33d10c4..f1ca8ca 100644
--- a/test/test_materials.py
+++ b/test/test_materials.py
@@ -6,7 +6,7 @@
 
 # TODO: currently only testing for build & compile, not for the validity of a result.
 
-from dolfinx import mesh, fem
+from dolfinx import mesh, fem, default_scalar_type
 from mpi4py import MPI
 from petsc4py import PETSc
 
@@ -27,29 +27,29 @@ def create_mesh(dim):
 def tst_scalar_material(dim, eltname="Lagrange"):
     # FE domain
     V = fem.FunctionSpace(create_mesh(dim), (eltname, 1))
-    
+
     # Material
-    const = lambda x: fem.Constant(V.mesh, PETSc.ScalarType(x))
+    const = lambda x: fem.Constant(V.mesh, default_scalar_type(x))
     mat  = material(V, 'scalar', rho=const(1), mu=const(1))
     mats = [mat]
-    
+
     # PDE
     pde = PDE(V, materials=mats)
-    
+
     # Compile some matrices
     _, _, _ = pde.M() ,  pde.C() ,  pde.K()
     _, _, _ = pde.K1(),  pde.K2(),  pde.K3()
-    
+
     # The end
 
 
 
 def tst_vector_materials(dim, nbcomps, eltname="Lagrange"):
     # FE domain
-    V = fem.VectorFunctionSpace(create_mesh(dim), (eltname, 1), dim=nbcomps)
+    V = fem.FunctionSpace(create_mesh(dim), (eltname, 1, (nbcomps,)))
     
     # Material
-    const = lambda x: fem.Constant(V.mesh, PETSc.ScalarType(x))
+    const = lambda x: fem.Constant(V.mesh, default_scalar_type(x))
     rho  = const(1)
     Coo  = const(1) #a dummy Cij
     types= ('isotropic', 'cubic', 'hexagonal', 'trigonal', 'tetragonal', 'orthotropic', 'monoclinic', 'triclinic')