diff --git a/hog/forms.py b/hog/forms.py
index 482dcf86dcd2785f99d081265cd52b68a57fd7c8..44b9500d4dcc123efd9b293896610ba917f8d66d 100644
--- a/hog/forms.py
+++ b/hog/forms.py
@@ -1631,10 +1631,10 @@ Weak formulation
 
     if trial != test:
         raise HOGException(
-            "Trial space must be equal to test space to assemble diffusion matrix."
+            "Trial space must be equal to test space to assemble second derivative matrix."
         )
 
-    with TimedLogger("assembling advection matrix", level=logging.DEBUG):
+    with TimedLogger("assembling second derivative matrix", level=logging.DEBUG):
         tabulation = Tabulation(symbolizer)
 
         jac_affine = symbolizer.jac_ref_to_affine(geometry.dimensions)
@@ -1642,89 +1642,90 @@ Weak formulation
         jac_affine_det = symbolizer.abs_det_jac_ref_to_affine()
 
         if isinstance(blending, ExternalMap):
+            TimedLogger("ExternalMap is not tested well").log()
             jac_blending = jac_affine_to_physical(geometry, symbolizer)
+            jacinvjac_blending = jacinvjac_affine_to_physical(geometry, symbolizer)
         else:
             affine_coords = trafo_ref_to_affine(geometry, symbolizer)
             jac_blending = blending.jacobian(affine_coords)
+            if not isinstance(blending, IdentityMap):
+                hessian_blending_map = blending.hessian(affine_coords)
 
         jac_blending_det = abs(det(jac_blending))
         with TimedLogger("inverting blending Jacobian", level=logging.DEBUG):
             jac_blending_inv = inv(jac_blending)
 
-        phi_eval_symbols = tabulation.register_phi_evals(
-            velocity_function_space.shape(geometry)
-        )
-        
-        ux, _ = fem_function_on_element(
-            velocity_function_space,
-            geometry,
-            symbolizer,
-            domain="reference",
-            function_id="ux",
-            basis_eval=phi_eval_symbols,
-        )
+        mat = create_empty_element_matrix(trial, test, geometry)
+        it = element_matrix_iterator(trial, test, geometry)
 
-        uy, _ = fem_function_on_element(
-            velocity_function_space,
-            geometry,
-            symbolizer,
-            domain="reference",
-            function_id="uy",
-            basis_eval=phi_eval_symbols,
-        )
+        if velocity_function_space != None:
+            u_eval_symbols = tabulation.register_phi_evals(
+                velocity_function_space.shape(geometry)
+            )
 
-        uz, _ = fem_function_on_element(
-            velocity_function_space,
-            geometry,
-            symbolizer,
-            domain="reference",
-            function_id="uz",
-            basis_eval=phi_eval_symbols,
-        )
+            ux, _ = fem_function_on_element(
+                velocity_function_space,
+                geometry,
+                symbolizer,
+                domain="reference",
+                function_id="ux",
+                basis_eval=u_eval_symbols,
+            )
 
-        if geometry.dimensions == 2:
-            u = sp.Matrix([[ux], [uy]])
-        elif geometry.dimensions == 3:
-            u = sp.Matrix([[ux], [uy], [uz]])
+            uy, _ = fem_function_on_element(
+                velocity_function_space,
+                geometry,
+                symbolizer,
+                domain="reference",
+                function_id="uy",
+                basis_eval=u_eval_symbols,
+            )
 
-        mat = create_empty_element_matrix(trial, test, geometry)
-        it = element_matrix_iterator(trial, test, geometry)
+            if isinstance(geometry, TetrahedronElement):
+                uz, _ = fem_function_on_element(
+                    velocity_function_space,
+                    geometry,
+                    symbolizer,
+                    domain="reference",
+                    function_id="uz",
+                    basis_eval=u_eval_symbols,
+                )
+                u = sp.Matrix([[ux], [uy], [uz]])
+            else:
+                u = sp.Matrix([[ux], [uy]])
 
         for data in it:
-            grad_phi = data.trial_shape_grad
             psi = data.test_shape
-            grad_psi = data.test_shape_grad
-            if blending != IdentityMap():
-                jac_affine_inv_T_grad_phi_symbols = tabulation.register_factor(
-                    "jac_affine_inv_T_grad_phi",
-                    jac_affine_inv.T * grad_phi,
-                )
-                jac_affine_inv_T_grad_psi_symbols = tabulation.register_factor(
-                    "jac_affine_inv_T_grad_psi",
-                    jac_affine_inv.T * grad_psi,
+            grad_phi = data.trial_shape_grad
+
+            jac_affine_inv_T_grad_phi_symbols = tabulation.register_factor(
+                "jac_affine_inv_T_grad_phi",
+                jac_affine_inv.T * grad_phi,
+            )
+
+            # jac_blending_inv_T_jac_affine_inv_T_grad_phi_symbols = tabulation.register_factor(
+            #     "jac_affine_inv_T_grad_psi",
+            #     jac_blending_inv.T * jac_affine_inv_T_grad_phi_symbols,
+            # )
+
+            if isinstance(blending, IdentityMap):
+                form = (
+                    dot(u, jac_affine_inv_T_grad_phi_symbols)
+                    * psi
+                    * jac_affine_det
                 )
+            else:
                 form = (
-                    dot(
-                        jac_blending_inv.T
-                        * sp.Matrix(jac_affine_inv_T_grad_phi_symbols),
-                        u,
-                    )
+                    dot(u, jac_blending_inv.T * jac_affine_inv_T_grad_phi_symbols)
+                    * psi
                     * jac_affine_det
                     * jac_blending_det
                 )
-            else:
-                jac_affine_inv_T_grad_phi = jac_affine_inv.T * grad_phi
-
-                jac_affine_inv_grad_phi_jac_affine_inv_grad_psi_det = (
-                    dot(jac_affine_inv_T_grad_phi, u) 
-                    * psi 
-                    * jac_affine_det,
-                )[0]
-                form = jac_affine_inv_grad_phi_jac_affine_inv_grad_psi_det
+                # HOGException("Only for testing with Blending map")
 
             mat[data.row, data.col] = form
 
-    return Form(mat, tabulation, symmetric=True, docstring=docstring)
+    return Form(mat, tabulation, symmetric=False, docstring=docstring)
 
 def supg_diffusion(
     trial: FunctionSpace,