diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index 17a1701346b3e199ca80135a8fc490d27d429131..3a297fd1f320e28d3f2308ed46e79c878e97bb90 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -909,6 +909,13 @@ def dd_axis(axis, ambient_dim, operand):
     """Return the derivative along (XYZ) axis *axis*
     (in *ambient_dim*-dimensional space) of *operand*.
     """
+    from pytools.obj_array import is_obj_array, with_object_array_or_scalar
+    if is_obj_array(operand):
+        def dd_axis_comp(operand_i):
+            return dd_axis(axis, ambient_dim, operand_i)
+
+        return with_object_array_or_scalar(dd_axis_comp, operand)
+
     d = Derivative()
 
     unit_vector = np.zeros(ambient_dim)
@@ -1409,7 +1416,9 @@ def n_cross(vec, where=None):
 
 def div(vec):
     ambient_dim = len(vec)
-    return sum(dd_axis(iaxis, ambient_dim, vec) for iaxis in range(ambient_dim))
+    return sum(
+            dd_axis(iaxis, ambient_dim, vec[iaxis])
+            for iaxis in range(ambient_dim))
 
 
 def curl(vec):
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index f81a5e550d06e879dd2c9309813d97381bef97b6..3e50aa94a36d8e5710a3bd3b6b4fbaa95ba855fd 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -389,16 +389,20 @@ def test_perf_data_gathering(ctx_getter, n_arms=5):
 
 # {{{ test 3D jump relations
 
-@pytest.mark.parametrize("relation", ["sp", "nxcurls"])
+@pytest.mark.parametrize("relation", ["sp", "nxcurls", "div_s"])
 def test_3d_jump_relations(ctx_factory, relation, visualize=False):
-    #logging.basicConfig(level=logging.INFO)
+    # logging.basicConfig(level=logging.INFO)
 
     pytest.importorskip("pyfmmlib")
 
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
-    target_order = 4
+    if relation == "div_s":
+        target_order = 3
+    else:
+        target_order = 4
+
     qbx_order = target_order
 
     from pytools.convergence import EOCRecorder
@@ -469,6 +473,16 @@ def test_3d_jump_relations(ctx_factory, relation, visualize=False):
                     - (sym.Sp(knl, density_sym, qbx_forced_limit="avg")
                         - 0.5*density_sym))
 
+        elif relation == "div_s":
+
+            density = m.cos(2*x) * m.cos(2*y) * m.cos(z)
+            density_sym = sym.var("density")
+
+            jump_identity_sym = (
+                    sym.div(sym.S(knl, sym.normal(3).as_vector()*density_sym,
+                        qbx_forced_limit="avg"))
+                    + sym.D(knl, density_sym, qbx_forced_limit="avg"))
+
         else:
             raise ValueError("unexpected value of 'relation': %s" % relation)