exact cut of a Function

Question

I want to do the the cut on a Function, namely
rho = 0 if rho<0, rho if 0=<rho=<1, 1 if rho>1

I have three variants to realize

def rho_cut_1(rho):
    func_sp = rho.function_space()
    rho_cut = Max(Constant(0.), Min(Constant(1.), rho))
    return project(rho_cut, func_sp)

def rho_cut_2(rho):
    func_sp = rho.function_space()
    upper_bound = conditional(gt(rho, 1.), 1., rho)
    rho_cut = conditional(lt(upper_bound, 0.), 0., upper_bound)
    return project(rho_cut, func_sp)

class rho_cut_3(Expression):
    def __init__(self, rho, rho_range=(0., 1.), val_range=(0., 1.)):
        self.rho = rho
        self.rho_range = rho_range
        self.val_range = val_range

    def eval(self, values, x):
        if self.rho(x) >= self.rho_range[1]:
            values[0] = self.val_range[1]
        elif self.rho_range[0] < self.rho(x) < self.rho_range[1] :
            values[0] = self.rho(x)
        else:
            values[0] = self.val_range[0]

To check them I simply did the following,

mesh = UnitSquareMesh(50, 50)
V = FunctionSpace(mesh, 'Lagrange', 1)

rho = interpolate(Expression('x[0]*x[0]+x[1]*x[1]'), V)
rho_1 = rho_cut_1(rho)
rho_2 = rho_cut_2(rho)
rho_3 = interpolate(rho_cut_3(rho), V)

Here are the results, rho_cut_1 and rho_cut_2 use project, and not end in an exact answer, but faster. rho_cut_3 gives the best output but very slow.

The more exact solution in rho_cut_3 may be due to interpolate. However interpolate only works for Expression.

Can I get an accurate answer like interpolate by using project?
Or do I need to go for JIT by using a cpp code Expression? How can I input my rho and share it in a cpp Expression?

Thx!

Answer 1

You can do

expr = "rho < 0. ? 0 : (1. < rho ? 1. : rho)"
rho_4 = interpolate(Expression(expr,rho=rho),V)

If you do not need parallel performance, you can work directly on the arrays.

Comments

cool! it's faster than the subclass one and in the same accuracy