Dear all,

I am dealing with a 3D nonlinear PDE which looks like, with u - scalar solution, f - vector solution:

div(f) = 0;

(a+b*norm2(f)**c)*f/norm2(f)=-grad(u)

where norm2(f) is the L2 norm. and a b c are predefined constant.

As you can see dividing f by its norm can cause dividing-by-zero, as a can be nonzero. In simple numpy scripts I can just do a if statement to make such case output a 0 vector. i.e. I tried

sele = FiniteElement(“CG”, mesh.ufl_cell(), 1)

vele = VectorElement(“CG”, mesh.ufl_cell(), 1)

mixed_space = FunctionSpace(mesh, MixedElement([sele, vele]))

U = Function(mixed_space)

u,f = split(U)

V = TestFunction(mixed_space)

v,tau = split(V)

def cfunc(f):

return conditional(le(norma(f),Constant(1e-6)),0*f,(a+b*norma(f)**c)/norma(f))

However this woudn’t compile with error

raise self._exception_type(self._format_raw(*message))

UFLException: Shape mismatch between conditional branches.

How should I get around with this or does fenics have its all function to deal with it?

Thanks a lot in advance!