diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py
index ac59ef4..cb3e176 100644
--- a/sympde/calculus/core.py
+++ b/sympde/calculus/core.py
@@ -93,10 +93,11 @@
 
 from operator  import mul, add
 from functools import reduce
+from math      import sqrt
 
 from sympy                    import Indexed, sympify
 from sympy                    import Matrix, ImmutableDenseMatrix
-from sympy                    import cacheit
+from sympy                    import cacheit, conjugate
 from sympy.core               import Basic
 from sympy.core               import Add, Mul, Pow
 from sympy.core.containers    import Tuple
@@ -182,6 +183,17 @@ def is_zero(x):
     else:
         return x == 0
 
+# def is_Function_complex(list):
+#     res=False
+#     for expr in list:
+#         if isinstance(expr, (ScalarFunction, VectorFunction)):
+#             res = res or expr.is_complex
+#         else:
+#             res = res or is_Function_complex(expr.args)
+#         if res:
+#             break
+#     return res
+
 #==============================================================================
 class BasicOperator(CalculusFunction):
     """
@@ -266,6 +278,7 @@ def is_scalar(atom):
     """
     return is_constant(atom) or isinstance(atom, ScalarFunction)
 
+
 #==============================================================================
 # TODO add dot(u,u) +2*dot(u,v) + dot(v,v) = dot(u+v,u+v)
 # now we only have dot(u,u) + dot(u,v)+ dot(v,u) + dot(v,v) = dot(u+v,u+v)
@@ -351,9 +364,14 @@ def __new__(cls, arg1, arg2, **options):
         args_2 = [i for i in b if not i.is_commutative]
         c2     = [i for i in b if not i in args_2]
 
+
+        c = Mul(*c1)*Mul(*c2)
+        # 1D case where everything is commutative
+        if args_1==[] and args_2==[]:
+            return(c)
+
         a = reduce(mul, args_1)
         b = reduce(mul, args_2)
-        c = Mul(*c1)*Mul(*c2)
 
         if str(a) > str(b):
             a,b = b,a
@@ -505,12 +523,19 @@ def __new__(cls, arg1, arg2, **options):
         args_2 = [i for i in b if not i.is_commutative]
         c2     = [i for i in b if not i in args_2]
 
+        # if is_Function_complex(b):
+        #     args_2 = [i if isinstance(i, conjugate) else conjugate(i) for i in args_2]
+        #     c2     = [i if isinstance(i, conjugate) else conjugate(i) for i in c2]
+
+        c = Mul(*c1)*Mul(*c2)
+        if args_1==[] and args_2==[]:
+            return(c)
+
         a = reduce(mul, args_1)
         b = reduce(mul, args_2)
-        c = Mul(*c1)*Mul(*c2)
 
-        if str(a) > str(b):
-            a,b = b,a
+        # if str(a) > str(b):
+        #     a,b = b,a
 
         obj = Basic.__new__(cls, a, b)
 
@@ -780,7 +805,6 @@ def eval(cls, expr):
                 raise ArgumentTypeError(msg)
 
         return cls(expr, evaluate=False)
-
 #==============================================================================
 class Curl(DiffOperator):
     """
diff --git a/sympde/expr/basic.py b/sympde/expr/basic.py
index 9666096..4db54cf 100644
--- a/sympde/expr/basic.py
+++ b/sympde/expr/basic.py
@@ -62,14 +62,16 @@ class BasicExpr(Expr):
     is_Function   = True
     is_linear     = False
     is_bilinear   = False
+    # TODO put is_sesquilinear=False here but it don't work
+    #is_sesquilinear   = False
     is_functional = False
 
     @property
     def fields(self):
         atoms  = self.expr.atoms(ScalarFunction, VectorFunction)
-        if self.is_bilinear or self.is_linear:
+        if self.is_bilinear or self.is_sesquilinear or self.is_linear:
             args = self.variables
-            if self.is_bilinear:
+            if self.is_bilinear or self.is_sesquilinear :
                 args = args[0]+args[1]
             fields = tuple(atoms.difference(args))
         else:
diff --git a/sympde/expr/equation.py b/sympde/expr/equation.py
index 5b4ed1a..66cb759 100644
--- a/sympde/expr/equation.py
+++ b/sympde/expr/equation.py
@@ -11,10 +11,10 @@
 from sympde.topology.space import VectorFunction, IndexedVectorFunction
 from sympde.topology import Boundary, NormalVector, TangentVector
 from sympde.topology import Trace, trace_0, trace_1
-from sympde.calculus import grad, dot
+from sympde.calculus import grad, dot, inner
 from sympde.core.utils import random_string
 
-from .expr import BilinearForm, LinearForm
+from .expr import BilinearForm, LinearForm, SesquilinearForm
 from .expr import linearize
 from .errors import ( UnconsistentLhsError, UnconsistentRhsError,
                       UnconsistentArgumentsError, UnconsistentBCError )
@@ -192,8 +192,8 @@ class Equation(Basic):
 
     def __new__(cls, lhs, rhs, trials, tests, bc=None, constraint=None):
         # ...
-        if not isinstance(lhs, BilinearForm):
-            raise UnconsistentLhsError('> lhs must be a bilinear')
+        if not isinstance(lhs, (BilinearForm, SesquilinearForm)):
+            raise UnconsistentLhsError('> lhs must be a bilinear or a Sesquilinear')
 
         if not isinstance(rhs, LinearForm):
             raise UnconsistentRhsError('> rhs must be a linear')
@@ -380,7 +380,10 @@ def __new__(cls, form, fields, bc=None, trials=None):
 def find(trials, *, forall, lhs, rhs, bc=None, constraint=None):
 
     tests = forall
-    lhs = BilinearForm((trials, tests), lhs)
+    if tests.space.codomain_complex:
+        lhs = SesquilinearForm((trials, tests), lhs)
+    else:
+        lhs = BilinearForm((trials, tests), lhs)
     rhs =   LinearForm(         tests , rhs)
 
     return Equation(lhs, rhs, trials, tests, bc=bc, constraint=constraint)
diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py
index e75189b..89a9674 100644
--- a/sympde/expr/evaluation.py
+++ b/sympde/expr/evaluation.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 from sympy import Abs, S, cacheit
-from sympy import Indexed, Matrix, ImmutableDenseMatrix
+from sympy import Indexed, Matrix, ImmutableDenseMatrix, conjugate
 from sympy import expand
 from sympy.core import Basic, Symbol
 from sympy.core import Add, Mul, Pow
@@ -139,7 +139,7 @@ def _get_trials_tests(expr, *, flatten=False):
     if not isinstance(expr, (BasicForm, BasicExpr)):
         raise TypeError("Expression must be of type BasicForm or BasicExpr, got '{}' instead".format(type(expr)))
 
-    if expr.is_bilinear:
+    if expr.is_bilinear or expr.is_sesquilinear:
         trials = _unpack_functions(expr.variables[0]) if flatten else expr.variables[0]
         tests  = _unpack_functions(expr.variables[1]) if flatten else expr.variables[1]
 
@@ -828,6 +828,10 @@ def eval(cls, expr, domain):
             expr      = cls(expr.expr, domain=domain)
             return LogicalExpr(expr, domain=domain)
 
+        elif isinstance(expr, conjugate):
+            expr      = cls(expr.args[0], domain=domain)
+            return conjugate(expr)
+
         return expr
 
 
diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py
index 1061b75..777e63b 100644
--- a/sympde/expr/expr.py
+++ b/sympde/expr/expr.py
@@ -5,6 +5,7 @@
 
 from sympy import Dummy
 from sympy import Matrix, ImmutableDenseMatrix
+from sympy import conjugate
 from sympy.core import Basic, S
 from sympy.core import Expr, Add, Mul
 from sympy.core.numbers import Zero as sy_Zero
@@ -30,6 +31,7 @@
 
 __all__ = (
     'BilinearForm',
+    'SesquilinearForm',
     'Functional',
     'IntAdd',
     'Integral',
@@ -287,7 +289,8 @@ def _sympystr(self, printer):
 
 #==============================================================================
 class Functional(BasicForm):
-    is_functional = True
+    is_functional   = True
+    is_sesquilinear = False
 
     def __new__(cls, expr, domain, evaluate=True, **options):
 
@@ -332,6 +335,7 @@ def space(self):
 #==============================================================================
 class LinearForm(BasicForm):
     is_linear = True
+    is_sesquilinear   = False
 
     def __new__(cls, arguments, expr, **options):
 
@@ -347,11 +351,16 @@ def __new__(cls, arguments, expr, **options):
         args = _sanitize_arguments(arguments, is_linear=True)
 
         # Check linearity with respect to the given arguments
-        if not is_linear_expression(expr, args):
-            msg = 'Expression is not linear w.r.t [{}]'.format(args)
-            raise UnconsistentLinearExpressionError(msg)
+        if args[0].is_complex:
+            if not is_antilinear_expression(expr, args):
+                msg = 'Complex case : Expression is not antilinear w.r.t [{}]'.format(args)
+                raise UnconsistentLinearExpressionError(msg)
+        else:
+            if not is_linear_expression(expr, args):
+                msg = 'Real case :Expression is not linear w.r.t [{}]'.format(args)
+                raise UnconsistentLinearExpressionError(msg)
 
-        # Create new object of type LinearForm
+            # Create new object of type LinearForm
         obj = Basic.__new__(cls, args, expr)
 
         # Compute 'domain' property (scalar or tuple)
@@ -411,6 +420,7 @@ def __call__(self, *tests, **kwargs):
 #==============================================================================
 class BilinearForm(BasicForm):
     is_bilinear = True
+    is_sesquilinear   = False
     _is_symmetric = None
 
     def __new__(cls, arguments, expr, **options):
@@ -430,6 +440,12 @@ def __new__(cls, arguments, expr, **options):
         # Distinguish between trial and test functions
         trial_functions, test_functions = args
 
+        if trial_functions[0].space.codomain_complex != test_functions[0].space.codomain_complex:
+            raise TypeError('Trial space and Test space should both be real.')
+        if trial_functions[0].space.codomain_complex:
+            raise TypeError('Trial space and Test space should be real. In the complex case, a SesquilinearForm should be called ')
+
+
         # Check linearity with respect to trial functions
         if not is_linear_expression(expr, trial_functions):
             msg = ' Expression is not linear w.r.t trial functions {}'\
@@ -515,9 +531,122 @@ def __call__(self, trials, tests, **kwargs):
 
         return expr
 
+#==============================================================================
+class SesquilinearForm(BasicForm):
+    is_sesquilinear = True
+    _is_hermitian   = None
+
+    def __new__(cls, arguments, expr, **options):
+
+        # Trivial case: null expression
+        if expr == 0:
+            return sy_Zero()
+
+        # Check that integral expression is given
+
+        if not expr.atoms(Integral):
+            raise ValueError('Expecting integral Expression')
+
+        # TODO: why do we 'sanitize' here?
+        args = _sanitize_arguments(arguments, is_bilinear=True)
+
+        # Distinguish between trial and test functions
+        trial_functions, test_functions = args
+        if trial_functions[0].space.codomain_complex != test_functions[0].space.codomain_complex:
+            raise TypeError('Trial space and Test space should both be complex.')
+        if not trial_functions[0].space.codomain_complex:
+            raise TypeError('Trial space and Test space should be complex. In the real case, a BilinearForm should be called ')
+
+
+        # Check linearity with respect to trial functions
+        if not is_linear_expression(expr, trial_functions):
+            msg = ' Real Case : Expression is not linear w.r.t trial functions {}'\
+                    .format(trial_functions)
+            raise UnconsistentLinearExpressionError(msg)
+
+        # Check linearity with respect to test functions
+        if not is_antilinear_expression(expr, test_functions):
+            msg = ' Complex Case : Expression is not linear w.r.t test functions {}'\
+                    .format(test_functions)
+            raise UnconsistentLinearExpressionError(msg)
+
+        # Create new object of type SesquilinearForm
+        obj = Basic.__new__(cls, args, expr)
+
+        # Compute 'domain' property (scalar or tuple)
+        # TODO: is this is useful?
+        obj._domain = _get_domain(expr)
+
+        return obj
+
+    @property
+    def variables(self):
+        return self._args[0]
+
+    @property
+    def expr(self):
+        return self._args[1]
+
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def test_functions(self):
+        return self.variables[1]
+
+    @property
+    def trial_functions(self):
+        return self.variables[0]
+
+    @property
+    def test_spaces(self):
+        return [u.space for u in self.test_functions]
+
+    @property
+    def trial_spaces(self):
+        return [u.space for u in self.trial_functions]
+
+    @property
+    def coordinates(self):
+        return self.test_spaces[0].coordinates
+
+    @property
+    def ldim(self):
+        return self.test_spaces[0].ldim
+
+    @property
+    def is_hermitian(self):
+        if self._is_hermitian is None:
+            left, right = self.variables
+            a1 = self(left, right)
+            a2 = self(right, left)
+            self._is_hermitian = (conjugate(a1) == a2)
+        return self._is_hermitian
+
+    def __call__(self, trials, tests, **kwargs):
+
+        # Use free variables if given and available
+        expr = self._update_free_variables(**kwargs)
+
+        # If needed, convert positional arguments to lists
+        if not is_sequence(trials): trials = [trials]
+        if not is_sequence(tests ): tests  = [tests ]
+
+        # Concatenate input values into single list
+        values = [*trials, *tests]
+
+        # Substitute variables with given values in bilinear expression
+        variables = [*self.variables[0], *self.variables[1]]
+        subs      = dict(zip(variables, values))
+        expr, _   = expr._xreplace(subs)
+
+        return expr
+
 #==============================================================================
 class Norm(Functional):
     is_norm = True
+    is_sesquilinear   = False
 
     def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
 #        # ...
@@ -549,36 +678,40 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
             exponent = 2
 
             if not is_vector:
-                expr = expr*expr
+                expr = expr * conjugate(expr)
 
             else:
                 if not( expr.shape[1] == 1 ):
                     raise ValueError('Wrong expression for Matrix. must be a row')
 
-                v = Tuple(*expr[:,0])
-                expr = Dot(v, v)
+                args = expr[:, 0]
+                args_conj=[conjugate(i) for i in args]
+                v = Tuple(*args)
+                v_conj = Tuple(*args_conj)
+                expr = Dot(v, v_conj)
 
         elif kind == 'h1'and evaluate :
             exponent = 2
 
             if not is_vector:
                 a    = Grad(expr)
-                expr = Dot(a, a)
+                expr = Dot(a, conjugate(a))
 
             else:
                 if not( expr.shape[1] == 1 ):
                     raise ValueError('Wrong expression for Matrix. must be a row')
-
-                v = Tuple(*expr[:,0])
-                a = Grad(v)
-                expr = Inner(a, a)
+                args = expr[:, 0]
+                args_conj=[conjugate(i) for i in args]
+                v = Tuple(*args)
+                v_conj = Tuple(*args_conj)
+                expr = Dot(Grad(v), Grad(v_conj))
 
         elif kind == 'h2'and evaluate :
             exponent = 2
 
             if not is_vector:
                 a    = Hessian(expr)
-                expr = Dot(a, a)
+                expr = Dot(a, conjugate(a))
 
             else:
                 raise NotImplementedError('TODO')
@@ -724,6 +857,78 @@ def is_linear_expression(expr, args, integral=True, debug=True):
 
     left_expr = expr.subs(list(zip(args, left_args)))
     left_expr = coeff * left_expr
+    if not( (newexpr-left_expr).expand() == 0 or newexpr.expand()==left_expr.expand()):
+        # TODO use a warning or exception?
+        if debug:
+            a=(newexpr-left_expr).expand()
+            b=newexpr.expand()
+            c=left_expr.expand()
+            print('Failed to assert multiplication property')
+            print('{} != {}'.format(newexpr, left_expr))
+        return False
+    # ...
+
+    return True
+
+#==============================================================================
+def is_antilinear_expression(expr, args, debug=True):
+    """checks if an expression is linear with respect to the given arguments."""
+    # ...
+    left_args  = []
+    right_args = []
+
+    for arg in args:
+        tag    = random_string( 4 )
+
+        if isinstance(arg, ScalarFunction):
+            left  = ScalarFunction(arg.space, name='l_' + tag)
+            right = ScalarFunction(arg.space, name='r_' + tag)
+
+        elif isinstance(arg, VectorFunction):
+            left  = VectorFunction(arg.space, name='l_' + tag)
+            right = VectorFunction(arg.space, name='r_' + tag)
+        else:
+            raise TypeError('argument must be a {Scalar|Vector}Function')
+
+        left_args  += [left]
+        right_args += [right]
+    # ...
+
+    # ... check addition
+    newargs = [left + right for left, right in zip(left_args, right_args)]
+
+    newexpr    = expr.subs(zip(args, newargs))
+    left_expr  = expr.subs(zip(args, left_args))
+    right_expr = expr.subs(zip(args, right_args))
+
+    a = newexpr
+    b = left_expr + right_expr
+
+    if not( (a-b).expand() == 0 or a.expand() == b.expand()):
+        # TODO use a warning or exception?
+        if debug:
+            print('Failed to assert addition property')
+            print('{} != {}'.format(a.expand(), b.expand()))
+        return False
+
+    # ...
+
+    # ... check multiplication
+    tag   = random_string( 4 )
+    coeff = Constant('alpha_' + tag)
+
+    newexpr = expr
+    for arg, left in zip(args, left_args):
+        newarg  = coeff * left
+        newexpr = newexpr.subs(arg, newarg)
+
+    atoms     = list(newexpr.atoms(BasicOperator))
+    subs      = [e.func(*e.args, evaluate=True) for e in atoms]
+    newexpr   = newexpr.subs(zip(atoms, subs))
+
+
+    left_expr = expr.subs(list(zip(args, left_args)))
+    left_expr = conjugate(coeff) * left_expr
     if not( (newexpr-left_expr).expand() == 0 or newexpr.expand()==left_expr.expand()):
         # TODO use a warning or exception?
         if debug:
diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py
index f5d3537..cad1d53 100644
--- a/sympde/topology/derivatives.py
+++ b/sympde/topology/derivatives.py
@@ -15,7 +15,7 @@
 from sympy import S
 from sympy import Indexed
 from sympy import diff
-from sympy import log
+from sympy import log, conjugate, Abs
 from sympy import preorder_traversal
 from sympy import cacheit
 from sympde.old_sympy_utilities import is_sequence
@@ -162,6 +162,12 @@ def eval(cls, expr):
             v = (log(b)*cls(e, evaluate=True) + e*cls(b, evaluate=True)/b) * b**e
             return v
 
+        elif isinstance(expr, conjugate):
+            return conjugate(cls(expr.args[0], evaluate=True))
+
+        elif isinstance(expr, Abs):
+            return Abs(cls(expr.args[0], evaluate=True))
+
         else:
             msg = '{expr} of type {type}'.format(expr=expr, type=type(expr))
             raise NotImplementedError(msg)
diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py
index 5c8b658..2bd41b9 100644
--- a/sympde/topology/mapping.py
+++ b/sympde/topology/mapping.py
@@ -80,7 +80,7 @@ def get_logical_test_function(u):
     kind            = space.kind
     dim             = space.ldim
     logical_domain  = space.domain.logical_domain
-    l_space         = type(space)(space.name, logical_domain, kind=kind)
+    l_space         = type(space)(space.name, logical_domain, kind=kind, codomain_complex=space.codomain_complex)
     el              = l_space.element(u.name)
     return el
 
@@ -895,7 +895,7 @@ def eval(cls, expr, domain, **options):
         """."""
 
         from sympde.expr.evaluation import TerminalExpr, DomainExpression
-        from sympde.expr.expr import BilinearForm, LinearForm, BasicForm, Norm
+        from sympde.expr.expr import BilinearForm, LinearForm, BasicForm, Norm, SesquilinearForm
         from sympde.expr.expr import Integral
 
         types = (ScalarFunction, VectorFunction, DifferentialOperator, Trace, Integral)
@@ -1221,6 +1221,12 @@ def eval(cls, expr, domain, **options):
             body    = cls.eval(expr.expr, domain)
             return BilinearForm((trials, tests), body)
 
+        elif isinstance(expr, SesquilinearForm):
+            tests   = [get_logical_test_function(a) for a in expr.test_functions]
+            trials  = [get_logical_test_function(a) for a in expr.trial_functions]
+            body    = cls.eval(expr.expr, domain)
+            return SesquilinearForm((trials, tests), body)
+
         elif isinstance(expr, LinearForm):
             tests   = [get_logical_test_function(a) for a in expr.test_functions]
             body    = cls.eval(expr.expr, domain)
diff --git a/sympde/topology/space.py b/sympde/topology/space.py
index cce1f55..2b88570 100644
--- a/sympde/topology/space.py
+++ b/sympde/topology/space.py
@@ -135,7 +135,7 @@ class BasicFunctionSpace(Basic):
     _kind       = None
     _regularity = None # TODO pass it as arg to __new__
     _is_broken  = None
-    def __new__(cls, name, domain, shape, kind):
+    def __new__(cls, name, domain, shape, kind, *, codomain_complex=False):
 
         if not isinstance(domain, BasicDomain):
             raise TypeError('> Expecting a BasicDomain object for domain')
@@ -144,7 +144,7 @@ def __new__(cls, name, domain, shape, kind):
         obj._name   = name
         obj._domain = domain
         obj._shape  = shape
-
+        obj._codomain_complex=codomain_complex
         # ...
         if kind is None:
             kind = 'undefined'
@@ -206,6 +206,10 @@ def is_broken(self):
     def regularity(self):
         return self._regularity
 
+    @property
+    def codomain_complex(self):
+        return self._codomain_complex
+
     @property
     def coordinates(self):
         return self.domain.coordinates
@@ -224,9 +228,9 @@ class ScalarFunctionSpace(BasicFunctionSpace):
     """
     Represents a basic continuous scalar Function space.
     """
-    def __new__(cls, name, domain, kind=None):
+    def __new__(cls, name, domain, kind=None, *, codomain_complex=False):
         shape = 1
-        return BasicFunctionSpace.__new__(cls, name, domain, shape, kind)
+        return BasicFunctionSpace.__new__(cls, name, domain, shape, kind, codomain_complex=codomain_complex)
 
     def element(self, name):
         return ScalarFunction(self, name)
@@ -236,9 +240,9 @@ class VectorFunctionSpace(BasicFunctionSpace):
     """
     Represents a basic continuous vector Function space.
     """
-    def __new__(cls, name, domain, kind=None):
+    def __new__(cls, name, domain, kind=None, *, codomain_complex=False):
         shape = domain.dim
-        return BasicFunctionSpace.__new__(cls, name, domain, shape, kind)
+        return BasicFunctionSpace.__new__(cls, name, domain, shape, kind, codomain_complex=codomain_complex)
 
     def element(self, name):
         return VectorFunction(self, name)
@@ -246,7 +250,7 @@ def element(self, name):
 #=============================================================================
 class Derham:
     """."""
-    def __init__(self, domain, sequence=None):
+    def __init__(self, domain, sequence=None, *, codomain_complex=False):
         shape = domain.dim
         self._V0  = None
         self._V1  = None
@@ -254,8 +258,8 @@ def __init__(self, domain, sequence=None):
         self._V3  = None
 
         if shape == 1:
-            spaces = [ScalarFunctionSpace('H1', domain, kind='H1'),
-                      ScalarFunctionSpace('L2', domain, kind='L2')]
+            spaces = [ScalarFunctionSpace('H1', domain, kind='H1', codomain_complex=codomain_complex),
+                      ScalarFunctionSpace('L2', domain, kind='L2', codomain_complex=codomain_complex)]
 
             self._V0  = spaces[0]
             self._V1  = spaces[1]
@@ -264,19 +268,19 @@ def __init__(self, domain, sequence=None):
             assert sequence is not None
 
             space = sequence[1]
-            spaces = [ScalarFunctionSpace('H1', domain, kind='H1'),
-                      VectorFunctionSpace(space, domain, kind=space),
-                      ScalarFunctionSpace('L2', domain, kind='L2')]
+            spaces = [ScalarFunctionSpace('H1', domain, kind='H1' ,codomain_complex=codomain_complex),
+                      VectorFunctionSpace(space, domain, kind=space, codomain_complex=codomain_complex),
+                      ScalarFunctionSpace('L2', domain, kind='L2', codomain_complex=codomain_complex)]
 
             self._V0  = spaces[0]
             self._V1  = spaces[1]
             self._V2  = spaces[2]
 
         elif shape == 3:
-            spaces = [ScalarFunctionSpace('H1', domain, kind='H1'),
-                      VectorFunctionSpace('Hcurl', domain, kind='Hcurl'),
-                      VectorFunctionSpace('Hdiv', domain, kind='Hdiv'),
-                      ScalarFunctionSpace('L2', domain, kind='L2')]
+            spaces = [ScalarFunctionSpace('H1', domain, kind='H1'      , codomain_complex=codomain_complex),
+                      VectorFunctionSpace('Hcurl', domain, kind='Hcurl', codomain_complex=codomain_complex),
+                      VectorFunctionSpace('Hdiv', domain, kind='Hdiv'  , codomain_complex=codomain_complex),
+                      ScalarFunctionSpace('L2', domain, kind='L2'      , codomain_complex=codomain_complex)]
 
             self._V0  = spaces[0]
             self._V1  = spaces[1]
@@ -434,13 +438,15 @@ class ScalarFunction(Symbol):
     is_commutative = True
     _space         = None
     _projection_of = None
+    is_complex     = False
 
     def __new__(cls, space, name):
         if not isinstance(space, ScalarFunctionSpace):
             raise ValueError('Expecting a ScalarFunctionSpace')
-        obj = Expr.__new__(cls)
-        obj._space = space
-        obj._name  = name
+        obj            = Expr.__new__(cls)
+        obj._space     = space
+        obj._name      = name
+        obj.is_complex = space.codomain_complex
         return obj
 
     @property
@@ -543,12 +549,14 @@ class VectorFunction(Symbol, IndexedBase):
     is_commutative = False
     _space         = None
     _projection_of = None
+    is_complex     = False
 
     def __new__(cls, space, name):
         if not isinstance(space, VectorFunctionSpace):
             raise ValueError('Expecting a VectorFunctionSpace')
         obj        = Expr.__new__(cls)
         obj._space = space
+        obj.is_complex=space.codomain_complex
         obj._name = name
         return obj