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 ?
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
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.