for some optimization problem (shape/topology) I want to penalize asymmetries.
Therefore, I’d like to know the best way to implement a penalty term of the form

short answer: such a non-local integral does not fit into a finite element framework (like NGSolve).

longer answer, but I do not want to make too much hope: It depends on what you want to do with such an integral. Just evaluating it, or computing first and second variations ? E.g. for second variation, the stiffness matrix gets additional entries depending on the argument transformation. For just evaluating the integral, one could maybe combine a coordinate transformation with point-search for the target coordinate.