Using sympy Spherical Harmonics to customize in-grid evaluation

NGSolve
1

Recently was trying to display real spherical harmonics on a custom coordinate system, and I tried very hard to make it work with ngsolve.bem.SphericalHarmonicsCF, however I gave up and found a more straightforward method using sympy:

import ngsolve as ngs
import sympy as sp
Ynm = sp.functions.special.spherical_harmonics.Ynm

def sympy_to_ngsolve_cf(sp_expr, sub_map):
    “”"Converts a sympy expression into an NGSolve CoefficientFunctionby replacing symbols according to sub_map.“”"
    sp_expr_substituted = sp_expr.expand(func=True).simplify()
    eval_scope = {
    'sqrt': ngs.sqrt,
    'sin': ngs.sin,
    'cos': ngs.cos,
    'exp': ngs.exp,
    'pi': ngs.pi ,
    'I' : 1j
    }

    for sym_obj, ngs_cf_obj in sub_map.items():
        eval_scope[str(sym_obj)] = ngs_cf_obj

    cf_string = sp.printing.sstr(sp_expr_substituted, full_prec=True)

    print(f" cf_string = {cf_string}")
    return eval(cf_string, {}, eval_scope)

def do_the_thing():
    theta_sym, phi_sym = sp.symbols('theta phi')
    real_SH = dict()

    mapping_objects = {
      theta_sym: ngs.CoefficientFunction(....),
      phi_sym: ngs.CoefficientFunction(....)
    }

    
    for L in range(order + 1):
        for M in range(-L, L + 1):
            print(f" building Y[{(L,M)}]:")
            if M == 0:
                real_SH[ (L, 0) ] = sympy_to_ngsolve_cf(Ynm( L, M, theta_sym, phi_sym ), mapping_objects)
            elif M > 0:
                real_SH[ (L, M) ] = sympy_to_ngsolve_cf((1 / sp.sqrt(2)) * ( Ynm( L, M, theta_sym, phi_sym ) + (-1)**M * Ynm( L, -M, theta_sym, phi_sym ) ), mapping_objects)
            elif M < 0:
                real_SH[ (L, M) ] = sympy_to_ngsolve_cf((sp.I / sp.sqrt(2)) * ( Ynm( L, M, theta_sym, phi_sym ) - (-1)**M * Ynm( L, -M, theta_sym, phi_sym ) ), mapping_objects)

Perhaps someone find this technique useful to move easily sympy functions into ngsolve