Problem with grad and interpolate

Errors
1

Hello everybody,
I’m working on the Cahn-Hilliard demo (Cahn-Hilliard equation — DOLFIN documentation), where I’m trying to change the initial conditions. In particular I am editing

class InitialConditions(UserExpression):
    #def __init__(self, **kwargs):
        #random.seed(2 + MPI.rank(MPI.comm_world))
        #super().__init__(**kwargs)
    def eval(self, values, x):
        #values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[0] = exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2)))
        values[1] = 200*(exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))*(1-exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))*(1-2*exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))-lmbda*div(grad(exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2)))))
        #values[1] = 1
    def value_shape(self):
        return (2)

But it seems to me that the grad function causes problems as to follow is present:

u.interpolate(u_init)

Which gives me the following error:

Cannot determine geometric dimension from expression.
Traceback (most recent call last):
  File "/Users/lorenzomarta/Desktop/Simulazioni_demo_fenics/Cahn-Hilliard/demo_cahn-hilliard.py", line 256, in <module>
    u.interpolate(u_init)
  File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/function/function.py", line 363, in interpolate
    self._cpp_object.interpolate(u._cpp_object)
  File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/function/expression.py", line 53, in wrapped_eval
    self.user_expression.eval(values, x)
  File "/Users/lorenzomarta/Desktop/Simulazioni_demo_fenics/Cahn-Hilliard/demo_cahn-hilliard.py", line 136, in eval
    values[1] = 200*(exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))*(1-exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))*(1-2*exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))-lmbda*div(grad(exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2)))))
  File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/operators.py", line 381, in grad
    return Grad(f)
  File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/differentiation.py", line 152, in __new__
    dim = find_geometric_dimension(f)
  File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/domain.py", line 384, in find_geometric_dimension
    error("Cannot determine geometric dimension from expression.")
  File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cannot determine geometric dimension from expression.
2

The div and grad operators are ufl.operators, and work on ufl objects. You are trying to use them inside a function that takes in numpy arrays (values and x), and therefore you see the error message.

I would suggest computing the div(grad(your_expression)) by hand, or use sympy to do it.

Please note that it is easier to do these things in DOLFINx, which is our new recommended software to use if you are just starting with DOLFINx,
see: The new DOLFINx solver is now recommended over DOLFIN

There I would create the expression with ufl:

from ufl import pow, exp, div, grad, as_vector

# Define mesh, V and u
# …
x = ufl.SpatialCoordinate(mesh)
v0 = exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2)))
v1 = 200*(exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))*(1-exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))*(1-2*exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2))))-lmbda*div(grad(exp(-(pow((x[0]-0.5),2)+pow((x[1]-0.5),2)))))
v = as_vector([v0,v1])
u_init = dolfinx.fem.Expression(v, V.element.interpolation_points)
u.interpolate(u_init)