Projector_Hcurl, ProjectorH_div not working with complex

[e-moral-sanchez]

Projectors of Hcurl and Hdiv do not work when using a complex underlying field. A work-around is to declare a deRham sequence with underlying real field, project lambda function and save coefficients of the real and imaginary part separately to a file and rerun simulation with loaded coefficients.

Minimal example (analogous for Projector_Hdiv):

from sympy import cos
from psydac.feec.global_projectors  import Projector_Hcurl, Projector_Hdiv

domain = Cube('domain', bounds1=(0, 1), bounds2=(0,1), bounds3=(0,1))
derham = Derham(domain)
derham.V0.codomain_type='complex'
derham.V1.codomain_type='complex'
derham.V2.codomain_type='complex'
derham.V3.codomain_type='complex'
domain_h = discretize(domain, ncells=[60,2,2], periodic=[False,True,True])
derham_h = discretize(derham, domain_h, degree=[3,1,1])
P1 = Projector_Hcurl(derham_h.V1, [5,2,2])
e0 = P1([lambda x,y,z: 0*z,lambda x,y,z: 0*z,lambda x,y,z: cos(x)])

Error in line 199 of psydac/linalg/direct_solvers.py :
out[:] = self._splu.solve(rhs.T, trans='T' if transposed else 'N').T
Message:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/elena-hp/psydac/psydac/feec/global_projectors.py", line 437, in __call__
    return super().__call__(fun)
  File "/home/elena-hp/psydac/psydac/feec/global_projectors.py", line 312, in __call__
    coeffs = self._solver.dot(self._rhs)
  File "/home/elena-hp/psydac/psydac/linalg/block.py", line 655, in dot
    self._func(self._blocks_as_args, v, out, **self._args)
  File "/home/elena-hp/psydac/psydac/linalg/block.py", line 672, in _dot
    out[i] += Lij.dot(v[j], out=inc[i])
  File "/home/elena-hp/psydac/psydac/linalg/kron.py", line 539, in dot
    return self.solve(v, out=out)
  File "/home/elena-hp/psydac/psydac/linalg/kron.py", line 567, in solve
    self._solve_nd(inslice, outslice)
  File "/home/elena-hp/psydac/psydac/linalg/kron.py", line 585, in _solve_nd
    self._solver_passes[i].solve_pass(temp1, temp2)
  File "/home/elena-hp/psydac/psydac/linalg/kron.py", line 681, in solve_pass
    self._solver.solve(view, out=view)
  File "/home/elena-hp/psydac/psydac/linalg/direct_solvers.py", line 199, in solve
    out[:] = self._splu.solve(rhs.T, trans='T' if transposed else 'N').T
TypeError: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe'

out and rhs both have dtype 'complex128', but the L, U matrices from the LU decomposition have the dtype 'float64' (because the matrix of which the decomposition is computed is also dtype 'float64' ).