Discretizing a norm with user defined mapping

[alex-m-h]

When I tried to compute the norm of a vector valued FemField, the discretize function generates a file, that is not executable. Here an example where I define a spline approximated torus mapping. The TorusCallableMapping class is defined via the BasicCallableMapping interface as follows:

class TorusCallableMapping(BasicCallableMapping):
    def __init__(self, R0: float) -> None:
        super().__init__()
        self._R0 = R0
    
    def __call__(self, *eta):
        """ Evaluate mapping at location eta. """
        x = self._R0 + eta[0] * np.cos(eta[1]) * np.cos(eta[2])
        y = self._R0 + eta[0] * np.cos(eta[1]) * np.sin(eta[2])
        z = self._R0 + eta[0] * np.sin(eta[1])

        return (x,y,z)


    def jacobian(self, *eta):
        """ Compute Jacobian matrix at location eta. """
        raise NotImplementedError

    def jacobian_inv(self, *eta):
        """ Compute inverse Jacobian matrix at location eta.
            An exception should be raised if the matrix is singular.
        """
        raise NotImplementedError

    def metric(self, *eta):
        """ Compute components of metric tensor at location eta. """
        raise NotImplementedError

    def metric_det(self, *eta):
        """ Compute determinant of metric tensor at location eta. """
        raise NotImplementedError

    @property
    def ldim(self):
        """ Number of logical/parametric dimensions in mapping
            (= number of eta components).
        """
        return 3
    
    @property
    def pdim(self):
        """ Number of physical dimensions in mapping
            (= number of x components).
        """
        return 3

Then in the following code, we define the domain, spaces and try to compute the norm of B_h:

cube = Cube(name='cube', bounds1=(0.1, 0.3), bounds2=(0,2*np.pi), bounds3=(0,0.25))
cube_discrete = discretize(cube, ncells=[8,8,4], periodic=[False, True, False])

V0_for_domain = ScalarFunctionSpace(name='V0_for_domain', domain=cube, kind='H1')

V0h_for_domain = discretize(V0_for_domain, cube_discrete, degree=[3,3,3])
torus_callable_mapping = TorusCallableMapping(R0=1.0)
torus_spline_mapping = SplineMapping.from_mapping(V0h_for_domain, torus_callable_mapping)

torus_mapping = Mapping(name='torus_mapping', dim=3)
torus_mapping.set_callable_mapping(torus_spline_mapping)
domain = torus_mapping(cube)

domain_h = discretize(domain, ncells=[8,8,4], periodic = [False, True, False])
derham = Derham(domain, ['H1', 'Hcurl', 'Hdiv', 'L2'])
degree = [2, 2, 2]
derham_h : DiscreteDerham = discretize(derham, domain_h, degree=degree)

projector_hcurl = Projector_Hcurl(derham_h.V1)
B1_R = lambda r, th, ze: r
B1_th = lambda r, th, ze: np.sin(th)
B1_ze = lambda r, th, ze: np.cos(ze)
B_h = projector_hcurl((B1_R, B1_th, B1_ze))

B, v = elements_of(derham.V1, 'B, v')
l2_norm_B_sym = Norm(ImmutableDenseMatrix([B[0],B[1], B[2]]), domain, kind='l2')
l2_norm_B_discrete = discretize(l2_norm_B_sym, domain_h, derham_h.V1)
l2_norm_B = l2_norm_B_discrete.assemble(B=B_h)
print(l2_norm_B)

I get the following traceback:

Traceback (most recent call last):
  File "isolatebug_normB.py", line 137, in <module>
    isolate_bug()
  File "isolatebug_normB.py", line 133, in isolate_bug
    l2_norm_B = l2_norm_B_discrete.assemble(B=B_h)
  File "/home/alex/repos/psydac/psydac/api/fem.py", line 1558, in assemble
    v = self._func(*args)
  File "/home/alex/repos/psydac/psydac/api/tests/__psydac__/dependencies_qxx3hnq7.py", line 152, in assemble_scalar_qxx3hnq7
    x = arr_x[i_quad_1,i_quad_2,i_quad_3]
NameError: name 'arr_x' is not defined