On 08/28/2018 12:43 AM, Irwin Zaid wrote:

Before answering your question, let me point out this. When doing electromagnetic simulations, AFAIK, vector finite elements from the H(curl) space are needed (e.g. Nedelec elements). I have no experience with such elements, and those are not implemented in SfePy. So unless you are sure that H1 elements are OK for your problem, SfePy might not be the right tool to use.

This is very clear, thanks. To be clear on my end, because this is a scattering problem, this actually reduces to solving the vector Helmholtz equation for a single field, either E or H. It's not the full Maxwell's equations, this is a standard approach for light scattering. The other field is easily determined from the first, as I mentioned above for H.

My understanding right now is that I can solve the vector Helmholtz equation within SfePy, and I believe I've implemented it correctly, but could be wrong.

Yes, the vector Helmholtz equation should be possible.

It is not currently possible, but it could be rather easily implemented, assuming your term evaluation function would be pure Python (or Cython). The reason of the current implementation is more or less historical - a number of term evaluation functions are implemented in C, so they expect just simple float64 arrays.

I'm happy to do pure Python or Cython. How would you suggest I proceed to get the complex numbers? Is there any way I can make this work with SfePy as it is now? I'd be willing to do PRs for SfePy in the near future if I am sure I can solve my problem within it -- which is why I'd like to work with it as it is now.

You would need to modify eval_complex() in sfepy/terms/terms.py so that it does not split the term arguments into the real and imaginary parts. For this, try using the term_mode argument. The call would then be something like:

val = problem.evaluate('your_term.i.Omega(E)', mode='eval', term_mode='complex')

and in eval_complex() you just check for term_mode == 'complex' and in that case just call

status = self.call_function(out, fargs)

I'll think of a cleaner way in the meantime (some other small changes in the term evaluation code might be needed - maybe adding 'ceval' to possible values of mode argument.

If this is meant just for post-processing, then the curl operator could be added. However, for solving the full coupled Maxwell equations, H(curl) elements would be required, as mentioned above. This is currently out of scope for me to add.

Yes, it is meant for just post-processing. Basically I solve for E, and get from H = const * curl(E). I'd also be interested in using the curl operator on parameter fields in the problem, but definitely not for test or trial fields. Again, is there a way to do this now with custom terms, or does it require changes to the SfePy source?

A new term similar to ev_grad would do it for the moment. It might be also a good idea to add curl as a possible mode to FieldVariable.evaluate().

I think we're on the same page -- I

*believe*I can use SfePy for this because I'm solving the equations for light scattering (vector Helmholtz equation), rather than the full Maxwell equations. But please correct me if I'm wrong. :)

Yes, the Helmholtz equation should work :)

r.