Implement variational form with conditional statements in C++

1

Dear all,

I am working on solving a non-linear PDE with a coefficient depending on the solution in a complex way. So far, I have used UserExpression to define the coefficient. For example, I have a coefficient q like below with UserExpression in Python (the detail of the variational form does not matter).

N = 100
X = 1.0
mesh = dl.IntervalMesh(N, 0, X)
degree = 1 # linear basis functions
Vh = dl.FunctionSpace(mesh, "Lagrange", degree)
u = dl.interpolate(dl.Expression("x[0]", degree = 1), Vh)
alpha = 2.0

class user_expr(dl.UserExpression):
    def set_alpha(self, u, alpha): # in the user_expr class
        self.u = u
        self.alpha = alpha

    def eval(self, value, x): # in the user_expr class
        if self.u(x) > 0.5:
            value[0] = 0.5
        else:  
            value[0] = self.alpha * self.u(x)* self.u(x)
        
    def value_shape(self):
        return ()

user_expr2 = user_expr(degree = 0)
user_expr2.set_alpha(u, alpha)

q = dl.interpolate(user_expr2, dl.FunctionSpace(mesh, "Lagrange", 1))
print(assemble((q*dl.dx))

However, to solve the PDE using the Newton method, it might be better to implement the variational form without using UserExpression (well, the Newton can be implemented with UserExpression, but mathematically a bit different from what I expected). If the coeffcient q is simple, such as

def q(u):
    return 1 + u**2

, then there is no problem in implementing it in Python like above and evaluate the variational form

print(assemble((q(u)*dl.dx))

However, my coefficient is very long and complex (includes conditional statements), so implementing it in Python leaded to unacceptable computational time. So, I want to implement the variational form in C++ to speed up the calculation. Is there any good ways?

Thank you for all the help from the community, I really appreciate that.