diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py
index 7cf0cc391..360f3827a 100644
--- a/psydac/api/discretization.py
+++ b/psydac/api/discretization.py
@@ -592,8 +592,8 @@ def discretize(a, *args, **kwargs):
         domain  = domain_h.domain
         dim     = domain.dim
 
-        # The current implementation of the sum factorization algorithm does not support openMP parallelization
-        # It is not properly tested, whether this way of excluding openMP parallelized code from using the sum factorization algorithm works.
+        # The current implementation of the sum factorization algorithm does not support OpenMP parallelization
+        # It is not properly tested, whether this way of excluding OpenMP parallelized code from using the sum factorization algorithm works.
         backend = kwargs.get('backend')# or None
         assembly_backend = kwargs.get('assembly_backend')# or None
         assembly_backend = backend or assembly_backend
@@ -657,7 +657,11 @@ def discretize(a, *args, **kwargs):
     #     return DiscreteSesquilinearForm(a, kernel_expr, *args, **kwargs)
 
     if isinstance(a, sym_BilinearForm):
-        if kwargs.pop('sum_factorization'):
+        # Currently sum factorization can only be used for interior domains
+        from sympde.expr.evaluation import DomainExpression
+        is_interior_expr = isinstance(kernel_expr, DomainExpression)
+
+        if kwargs.pop('sum_factorization') and is_interior_expr:
             return DiscreteBilinearForm_SF(a, kernel_expr, *args, **kwargs)
         else:
             return DiscreteBilinearForm(a, kernel_expr, *args, **kwargs)
diff --git a/psydac/api/fem_bilinear_form.py b/psydac/api/fem_bilinear_form.py
index 4bd65193a..06b50b946 100644
--- a/psydac/api/fem_bilinear_form.py
+++ b/psydac/api/fem_bilinear_form.py
@@ -12,6 +12,7 @@
 from sympy                  import ImmutableDenseMatrix, Matrix, Symbol, sympify
 from sympy.tensor.indexed   import Indexed, IndexedBase
 from sympy.simplify         import cse_main
+from sympy.printing.pycode  import pycode
 
 from pyccel import epyccel
 
@@ -149,6 +150,7 @@ def __init__(self, expr, kernel_expr, domain_h, spaces, *, nquads,
         self._domain = domain_h.domain
         self._spaces = spaces
         self._matrix = matrix
+        self._func   = None    # The assembly function will be generated later
 
         domain = self.domain
         target = self.target
@@ -251,6 +253,11 @@ def __init__(self, expr, kernel_expr, domain_h, spaces, *, nquads,
 
         #...
 
+        # Determine the type of scalar quantities to be managed in the code
+        dtypes = [getattr(V.symbolic_space, 'codomain_type', 'real') for V in (trial_space, test_space)]
+        assert all(t in ['complex', 'real'] for t in dtypes)
+        self._dtype = 'complex128' if 'complex' in dtypes else 'float64'
+
         # Assuming that all vector spaces (and their Cartesian decomposition,
         # if any) are compatible with each other, extract the first available
         # vector space from which (starts, ends, npts) will be read:
@@ -1216,7 +1223,7 @@ def make_file(self, temps, ordered_stmts, field_derivatives, max_logical_derivat
                 g_mat = g_mat_str.format(u_i=v_j, v_j=u_i)
             else:
                 g_mat = g_mat_str.format(u_i=u_i, v_j=v_j)
-            G_MAT += f'                    {g_mat} : "float64[:,:,:,:,:,:]",\n'
+            G_MAT += f'                    {g_mat} : f"{self._dtype}[:,:,:,:,:,:]",\n'
 
             TT1 += tt1_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:,:]", '
             TT2 += tt2_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:,:]", '
@@ -1282,11 +1289,14 @@ def make_file(self, temps, ordered_stmts, field_derivatives, max_logical_derivat
             global_span_v = global_span_v_str.format(v_j=v_j)
             LOCAL_SPAN += f'        {local_span_v_1} = {global_span_v}1[k_1]\n'
 
+        # Print expressions using SymPy's Python code printer
+        pyc = lambda expr: pycode(expr, fully_qualified_modules=False)
+
         for temp in temps:
-            TEMPS += f'                            {temp.lhs} = {temp.rhs}\n'
+            TEMPS += f'                            {temp.lhs} = {pyc(temp.rhs)}\n'
         for block in blocks:
             for stmt in ordered_stmts[block]:
-                COUPLING_TERMS += f'                            {stmt.lhs} = {stmt.rhs}\n'
+                COUPLING_TERMS += f'                            {stmt.lhs} = {pyc(stmt.rhs)}\n'
         
         KEYS = KEYS_2 + KEYS_3
 
diff --git a/psydac/api/fem_sum_form.py b/psydac/api/fem_sum_form.py
index fc8d05fb1..ca2556c4a 100644
--- a/psydac/api/fem_sum_form.py
+++ b/psydac/api/fem_sum_form.py
@@ -59,6 +59,11 @@ def __init__(self, a, kernel_expr, *args, **kwargs):
         self._kernel_expr = kernel_expr
         operator = None
         for e in kernel_expr:
+
+            # Currently sum factorization can only be used for interior domains
+            from sympde.expr.evaluation import DomainExpression
+            is_interior_expr = isinstance(e, DomainExpression)
+
             if isinstance(a, sym_LinearForm):
                 kwargs['update_ghost_regions'] = False
                 ah = DiscreteLinearForm(a, e, *args, backend=backend, **kwargs)
@@ -74,7 +79,7 @@ def __init__(self, a, kernel_expr, *args, **kwargs):
 
             elif isinstance(a, sym_BilinearForm):
                 kwargs['update_ghost_regions'] = False
-                if sum_factorization:
+                if sum_factorization and is_interior_expr:
                     ah = DiscreteBilinearForm_SF(a, e, *args, assembly_backend=backend, **kwargs)
                 else:
                     ah = DiscreteBilinearForm(a, e, *args, assembly_backend=backend, **kwargs)