working with complex valued functions

[gcdiwan]

Hi David,
Do you know how I could get the library working with complex numbers? For eg. what is the correct way of installing tIGAr with dolfinx using its docker image (see https://github.com/FEniCS/dolfinx)?
Many thanks,
Ganesh

[david-kamensky]

Getting tIGAr working with DOLFIN-X would involve significant modifications of the library itself. For now, the only way to implement a complex PDE is to rewrite it as a coupled system of two real PDEs, where the real and imaginary parts of the complex solution are separate unknown fields.

[gcdiwan]

Yes, I am trying the following but the solution does not make sense to me.

` normal = spline.n #normals?
ur = TrialFunction(spline.V)
ui = TrialFunction(spline.V)
vr = TestFunction(spline.V)
vi = TestFunction(spline.V)
omega = math.pi
# real part domain integral
ar1robin = (inner(nabla_grad(ur), nabla_grad(vr)) - inner(nabla_grad(ui), nabla_grad(vi))) * spline.dx
- (omega2)*inner(ur, vr)*spline.dx + (omega2)*inner(ui, vi)*spline.dx
# imag part domain integral
ai1robin = (inner(nabla_grad(ur), nabla_grad(vi)) + inner(nabla_grad(ui), nabla_grad(vr))) * spline.dx
- ( omega2)*inner(ur, vi)*spline.dx - (omega2)*inner(ui, vr)*spline.dx
# boundary integral on robin bc contributing to weak form: real
IrRobin = -omega * inner(ui, vr)*spline.ds - omega * inner(ur, vi)*spline.ds
# boundary integral on robin bc contributing to weak form: imag
IiRobin = omega * inner(ur, vr)*spline.ds - omega * inner( ui, vi)*spline.ds

x = spline.spatialCoordinates()
## rhs boundary integrals:
gr = - omega*sin(omega*x[0]) - normal[0]*omega*sin(omega*x[0])
gi = omega*cos(omega*x[0]) + normal[0]*omega*cos(omega*x[0])

uh = Function(spline.V)
arobin = ar1robin + ai1robin + (IrRobin + IiRobin) 
rhsRealInt = - vi*(omega*cos(omega*x[0]) + normal[0]*omega*cos(omega*x[0]))*spline.ds - vr*(omega*sin(omega*x[0]) + normal[0]*omega*sin(omega*x[0]))*spline.ds    
rhsImagInt = vr*(omega*cos(omega*x[0]) + normal[0]*omega*cos(omega*x[0]))*spline.ds - vi*(omega*sin(omega*x[0]) + normal[0]*omega*sin(omega*x[0]))*spline.ds
Lrobin = (rhsRealInt + rhsImagInt)   
spline.solveLinearVariationalProblem(arobin==Lrobin,uh)
# uhr, uhi = uh.split(True) ## fails
File("soln.pvd") << uh ## check what fenics is writing`

[david-kamensky]

I don't see the spline setup, but, from the definition of the trial and test functions, it looks like maybe spline.V is a scalar space. What you want to do is set up a spline space with two fields (e.g., splineGenerator = EqualOrderSpline(2,controlMesh), so spline.V is vector-valued with two components, not scalar-valued). Then, instead of defining ur and ui as two independent TrialFunctions in a scalar space, you would define

u = TrialFunction(spline.V)
ur,ui = split(u)

and likewise for the TestFunction.

Another note: I've previously run into bugs with complicated boundary integrands on quad meshes in FEniCS (independent of tIGAr), where the ds integrals are incorrect. As another debugging step, you might try extracting splines to a triangle mesh, by passing the optional useRect=False keyword argument to the constructor for the control mesh. (This is less efficient, but the triangle element implementation in FEniCS is more mature than the relatively-recent quad element implementation, so not all functionality is available for quad elements yet.)

[gcdiwan]

I don't see the spline setup, but, from the definition of the trial and test functions, it looks like maybe spline.V is a scalar space.

Yes, it indeed was! I simply copied the Poisson example code and hoped it would work with complex-valued functions. After setting up the correct space as you suggest above, I now get the correct solution to the problem.

've previously run into bugs with complicated boundary integrands on quad meshes in FEniCS (independent of tIGAr), where the ds integrals are incorrect

Thanks for the tip!

Best,
Ganesh