Write Expression in c++ code

mathematics
1

Hi all. I have a problem in which I need to project a Fourier series onto a finite element mesh, and for that, I’ve been using a subclass of dolfin.UserExpression to write the series. However, due to the complexity and required number of calls to the eval() method, it simply takes too long to calculate. Could you please help me translate the FunctionEvolution() class in the provided MWE into C++ language so that I could use a JIT-compiled Expression and reduce the computational time?

Thanks in advance for your help, and sorry for the relatively long MWE.

import dolfin as dlf
import ufl
import numpy as np

n = 4
mesh = dlf.BoxMesh(dlf.Point(0., 0., 0.), dlf.Point(1., 1., 1.), n, n, 1)

V_scalar   = dlf.FunctionSpace(mesh, "CG", 1)
point_source = dlf.Function(V_scalar)
func = dlf.Function(V_scalar, name="Function")

class Face(dlf.SubDomain):
    def inside(self, x, on_boundary):
        return dlf.near(x[2], 0.) and on_boundary
mesh_face = Face()
facets = dlf.MeshFunction(
    "size_t", mesh, mesh.topology().dim()-1
)
facets.set_all(0)
mesh_face.mark(facets, 1)
ds_face = dlf.Measure('ds', domain=mesh, subdomain_data=facets)

exp_params = {
    "ds"                     : ds_face,
    "func_space"             : V_scalar,
    "point_source"           : point_source,
}

class FunctionEvolution(dlf.UserExpression):
    def __init__(self, params, **kwargs):
        super().__init__(**kwargs)
        self._ds       = params["ds"]
        self._V        = params["func_space"]
        self._l        = params["point_source"]


    def Bmn(self, m, n, x1, x2):
        self._l.vector()[:] = 0.
        p_source = dlf.PointSource(
            self._V,
            dlf.Point(
                0.5, 
                0.5 + 0.016,
                0.
            ),  
            1. 
        )
        p_source.apply(self._l.vector())

        return (
            dlf.assemble(
                ufl.sin(m * np.pi * x1)
                * ufl.sin(n * np.pi * x2)
                * self._l(x1, x2, 0.) * self._ds(1)
            )
        )

    def eval(self, values, x):
        values[0] = sum([
            sum([
                self.Bmn(m+1, n+1, x[0], x[1]) 
                * ufl.sin(m+1 * np.pi * x[0])
                * ufl.exp(((m+1)**2 * np.pi**2 + (n+1)**2 * np.pi**2))
                for n in range(4)
            ])
        for m in range(4)
        ])

temp_expr = FunctionEvolution(
    params=exp_params, element=V_scalar.ufl_element()
)

func.assign(dlf.interpolate(temp_expr, V_scalar))

file_results = dlf.XDMFFile("outfile.xdmf")
file_results.write(func, 0.)
2

Your python code very inefficient because you’re using so many convenience functions.

    values[0] = sum([
        sum([
            self.Bmn(m+1, n+1, x[0], x[1]) 
            * ufl.sin(m+1 * np.pi * x[0])
            * ufl.exp(((m+1)**2 * np.pi**2 + (n+1)**2 * np.pi**2))
            for n in range(4)
        ])
    for m in range(4)
    ])

UFL is designed for symbolic formulations and not numeric evaluation. Furthermore, generating forms for each eval call is very expensive.

self.Bmn(m+1, n+1, x[0], x[1])

It’s not obvious if you use dolfin.Constants to reduce JIT compile time for parameters like m and n either.

        self._l.vector()[:] = 0.
        p_source = dlf.PointSource(
            self._V,
            dlf.Point(
                0.5, 
                0.5 + 0.016,
                0.
            ),  
            1. 
        )
        p_source.apply(self._l.vector())

Your use of PointSource also doesn’t make much sense.

                * self._l(x1, x2, 0.) * self._ds(1)

Calling dolfin.Function.__call__ with only coordinates is extremely expensive since you have to traverse the mesh’s BoundingBoxTree to find the required cell each time.

            dlf.assemble(
                ufl.sin(m * np.pi * x1)
                * ufl.sin(n * np.pi * x2)
                * self._l(x1, x2, 0.) * self._ds(1)
            )

If your domain is just a box are you sure you can’t just compute this integral analytically?

I’d recommend refactoring so you only compute what’s required in eval rather than recomputing for every single call. You would have to do this for a C++ implementation anyway since you won’t have the convenience of JIT compiling.

2 Likes