Create a custom initial condition for a field

General
1

I would like to to define a custom initial condition for a density field rho_n
I tried several things, including :

domain = Circle(Point(0, 0), R) # Define geometry for background
mesh = generate_mesh(domain, 50)
K_CG1 = FunctionSpace(mesh, dFE_CG1,constrained_domain=pbc)

class MyExpression0(UserExpression):
    def eval(self, value, r):
        m = (1 - numpy.arctan((R - R0) / dr)) / 2.0
        value[0] =  (1 - numpy.arctan((r - R0) / dr)) / 2.0 - m
    def value_shape(self):
        return (1,)

f0 = MyExpression0()
x = SpatialCoordinate(mesh)
r = project(sqrt(x[0]**2+x[1]**2))
rho_n = project(f0(r),K_CG1)

This yields an error : NotImplementedError: Cannot take length of non-vector expression.

How should I do this ?

2

Here is an even more minimal exemple :


domain = Circle(Point(0, 0), R) # Define geometry for background
mesh = generate_mesh(domain, 50) 
K_CG1 = FunctionSpace(mesh, dFE_CG1,constrained_domain=pbc)

f0 = Expression('1.0', degree=2)
x = SpatialCoordinate(mesh)
r = project(sqrt(x[0]**2+x[1]**2))
rho_n = project(f0(x),K_CG1)

Yielding :
TypeError: expected scalar arguments for the coordinates

3

I’d use ufl directly, something like

x = SpatialCoordinate(mesh)
r = ufl.sqrt(x[0]**2+x[1]**2)
R = whatever
R0 = whatever
dr = whatever
m = (1 - ufl.atan((R - R0) / dr)) / 2.0
f0 =  (1 - ufl.atan((r - R0) / dr)) / 2.0 - m

For future posts, please notice that a minimal example must be reproducible, i.e. something that I could copy in my own terminal and run. What you produced are snippets which are very minimal (thanks!), but not reproducible.