PDE complex numbers
Can anyone confirm that complex numbers are still not supported on fenics?
Is there a good example on how to solve the wave equation in 2D in the frequency domain (or something similar) somewhere?
I've tried to follow this https://fenicsproject.org/qa/6775/gmres-convergence-stronly-depends-on-frequency, but I get problems (TypeError: unsupported operand type(s) for *: 'FunctionSpace' and 'FunctionSpace') with the following:
V = FunctionSpace(mesh, 'Lagrange', 1)
V2 = V * V
To solve the Helmholtz equation with 1st order absorbing boundary condition, for example, you can use something like:
Vr = FiniteElement('P', mesh.ufl_cell(), degree) Vi = FiniteElement('P', mesh.ufl_cell(), degree) Vcomplex = Vr * Vi V = FunctionSpace(mesh, Vcomplex) #Real and Imaginary parts of the trial and test functions ur, ui = TrialFunctions(V) vr, vi = TestFunctions(V) a = - inner(grad(ur), grad(vr))*dx + (k**2)*inner(ur, vr)*dx\ - inner(grad(ui), grad(vi))*dx + (k**2)*inner(ui, vi)*dx \ + k*(ur*vi*ds) - k*(ui*vr*ds) # Source fr - real fi - imaginary L = (fr*vr + fi*vi)*dx A = assemble(a) b = assemble(L) u = Function(V) solve(A,u.array(),b) # Split u in real and imaginary parts ur , ui = u.split()
However, If you want to get the complex matrix, it's just a matter of finding the real and imaginary dofs:
# dofs Real Part dofs_real = V.sub(0).dofmap().dofs() # dofs Imaginary Part dofs_imag = V.sub(1).dofmap().dofs()
To correctly extract the real and imaginary matrices, in one of scipy sparse matrix classes:
# Slice Matrix - Real Part A_real = A[dofs_real,:] A_real = A_real[:,dofs_real] b_real = b [dofs_real] # Slice Matrix - Real Part A_imag = A[dofs_imag,:] A_imag = A_imag[:,dofs_real] b_imag = b [dofs_imag]
mixed_element = MixedElement(V.ufl_element(),V.ufl_element()) V2 = FunctionSpace(mesh,mixed_element)
I believe fenics does not support complex numbers so you will have to solve for the real and imaginary part equations separately.
V2 = VectorFunctionSpace(mesh, 'CG', pdim, dim=2) (u_r, u_i) = TrialFunction(V2) (phi_r, phi_i) = TestFunction(V2) ... u = Function(V2) A_tot = A_block(u_r,phi_r,omega_) - A_block(u_i,phi_i,omega_) + A_block(u_r,phi_i,omega_) + A_block(u_i,phi_r,omega_) b_tot = b_block(source_termR,phi_r) + b_block(source_termI,phi_i) A = assemble(A_tot) b = assemble(b_tot) [bc.apply(A,b) for bc in bcs] solve(A, u.vector(), b) (u_real, u_imag) = u.split(True)
but I don't get the expected result.. and I can't tell if it is a bug somewhere else or if this "split" doesn't work properly.
Ideally I'd like to do something like in this paper http://pe.org.pl/articles/2017/3/40.pdf
Does anyone has an example ???