Green's function, "double mesh", and weak singularities.


I am trying to compute a first order corrector for a perturbed resonance problem. To do this I need to calculate \int_B' \int_B G(x',x)u(x')u*(x)dx dx'

Where B and B’ are the unit squares
x = [x_1, x_2]
x’ = [x_1’,x_2’]
G is the Hankel function of the first kind of order 0 (it has a log singularity when x=x’)
u and u* are found using an Arnoldi solver (similar to 1.7.1 in the tutorials) and are a resonance function and its conjugate.

I can do this using base python but the estimated runtime is very large and I know NGSolve is much smarter about doing these types of problems than my basic python for loops.

My questions are as follows:
First, to evaluate the integral I need u, and u* to be defined on duplicate mesh (same mesh twice once for u and one for u*) but with different variables and I’m unsure how to do such a thing. Do I need to define a new mesh (different from the one used in the Arnoldi) or can I just define functions such that u(x_1, x_2) and u*(x_1’, x_2’) and then use the same mesh twice somehow.

Second, can the integral (or assemble or apply) handle the weak singularity since G=H(k sqrt(x_1-x_1')^2+(x_2-x_2')^2) or is there some built-in feature that deals with this issue. (k is the resonance produced from the Arnoldi solver). If not how would I go about modifying integrate, or assemble, or apply such that if NaN occurs do this instead while still retaining the speed of NGSolve.

Any help or ideas would be greatly appreciated