diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py
index a229b48efec686c48ebc1c53efdc553737463b0f..5f67e8080a0922cb7549f8bcb4c401ccacee8ae5 100644
--- a/pytential/symbolic/execution.py
+++ b/pytential/symbolic/execution.py
@@ -127,24 +127,24 @@ class EvaluationMapper(EvaluationMapperBase):
         assert not is_obj_array(result)  # numpy bug with obj_array.imag
         return result.imag
 
-    def apply_sqrt(self, args):
+    def apply_conj(self, args):
         arg, = args
-        return cl.clmath.sqrt(self.rec(arg))
+        return self.rec(arg).conj()
 
     def apply_abs(self, args):
         arg, = args
         return abs(self.rec(arg))
 
-    def apply_conj(self, args):
-        arg, = args
-        return self.rec(arg).conj()
-
     # }}}
 
     def map_call(self, expr):
-        from pytential.symbolic.primitives import Function
-        if isinstance(expr.function, Function):
+        from pytential.symbolic.primitives import EvalMapperFunction, CLMathFunction
+
+        if isinstance(expr.function, EvalMapperFunction):
             return getattr(self, "apply_"+expr.function.name)(expr.parameters)
+        elif isinstance(expr.function, CLMathFunction):
+            return getattr(cl.clmath, expr.function.name)(
+                    *(self.rec(arg) for arg in expr.parameters), queue=self.queue)
         else:
             return EvaluationMapperBase.map_call(self, expr)
 
diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index 27c2e24c950627a13cc363303c809b9ff200cc5f..57b3ff97cc0dcdb9469deb46b59a5b8ae154f14a 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -153,14 +153,43 @@ class Function(var):
 
             return with_object_array_or_scalar(make_op, operand)
         else:
-            return var.__call__(self, operand)
+            return var.__call__(self, operand, *args, **kwargs)
 
 
-real = Function("real")
-imag = Function("imag")
-conj = Function("conj")
-sqrt = Function("sqrt")
-abs = Function("abs")
+class EvalMapperFunction(Function):
+    pass
+
+
+class CLMathFunction(Function):
+    pass
+
+
+real = EvalMapperFunction("real")
+imag = EvalMapperFunction("imag")
+conj = EvalMapperFunction("conj")
+abs = EvalMapperFunction("abs")
+
+sqrt = CLMathFunction("sqrt")
+
+sin = CLMathFunction("sin")
+cos = CLMathFunction("cos")
+tan = CLMathFunction("tan")
+
+asin = CLMathFunction("asin")
+acos = CLMathFunction("acos")
+atan = CLMathFunction("atan")
+atan2 = CLMathFunction("atan2")
+
+sinh = CLMathFunction("sinh")
+cosh = CLMathFunction("cosh")
+tanh = CLMathFunction("tanh")
+
+asinh = CLMathFunction("asinh")
+acosh = CLMathFunction("acosh")
+atanh = CLMathFunction("atanh")
+
+exp = CLMathFunction("exp")
+log = CLMathFunction("log")
 
 
 class DiscretizationProperty(Expression):
@@ -267,6 +296,9 @@ def parametrization_derivative(ambient_dim, dim, where=None):
 
 
 def pseudoscalar(ambient_dim, dim=None, where=None):
+    """
+    Same as the outer product of all parametrization derivative columns.
+    """
     if dim is None:
         dim = ambient_dim - 1
 
@@ -299,7 +331,9 @@ def normal(ambient_dim, dim=None, where=None):
     pder = (
             pseudoscalar(ambient_dim, dim, where)
             / area_element(ambient_dim, dim, where))
-    return cse(pder.I | pder, "normal",
+    return cse(
+            # Dorst Section 3.7.2
+            pder << pder.I.inv(),
             cse_scope.DISCRETIZATION)