Hi,
I am trying to convert the following legacy FEniCS code to FEniCSx:
from dolfin import *
# Define the mesh
N = 32
mesh = RectangleMesh(Point(0, 0), Point(2.5, 2.5), N, N)
# Define the mixed finite elements
Ue = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Us = FunctionSpace(mesh, Ue)
element = MixedElement([Ue, Ue])
ME = FunctionSpace(mesh, element)
# Lower the log level to get more verbose output
set_log_level(40)
def BC_U_and_stim(t):
## returns tuple of BCs to apply (to V) and stimulus term at time t
## What happens when you play with this parameters?
start_stim_end_time = 6
stimstarttime = 150
stimduration = 6.
stimamp = 0.8
if t < start_stim_end_time:
# set a Dirichlet boundary condition of u=1 at left edge of the domain, but no stimulus term S
subdomain = CompiledSubDomain("near(x[0], 0) && on_boundary")
return [DirichletBC(ME.sub(0), Constant(1), subdomain)], Constant(0)
elif t > stimstarttime and t < stimstarttime + stimduration:
# have no Dirichlet boundary conditions, but a stimulus in the bottom left quarter of the domain
class Stimulus(UserExpression):
def eval(self, value, x):
if x[0] < 1.25 and x[1] < 1.25:
value[0] = stimamp
else:
value[0] = 0
def value_shape(self):
return ()
return [], Stimulus(degree=0)
else:
# no stimulus and no Dirichlet boundrary conditions
return [], Constant(0)
My tentative so far has been
# Define the mesh
N = 32
length = 2.5
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([length,length])], [N,N], mesh.CellType.triangle)
# Define the mixed finite elements
Ue = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
element = ufl.MixedElement([Ue, Ue])
ME = fem.FunctionSpace(domain, element)
ME_0 = ME.sub(0).collapse()[0]
# Lower the log level to get more verbose output
log.set_log_level(log.LogLevel.WARNING)
def BC_U_and_stim(domain,stim0,t):
## returns tuple of BCs to apply (to V) and stimulus term at time t
## What happens when you play with this parameters?
start_stim_end_time = 6
stimstarttime = 150
stimduration = 6.
stimamp = 0.8
def left(x):
return np.isclose(x[0],0)
left_dofs = fem.locate_dofs_geometrical((ME.sub(0),ME_0), left)
if t < start_stim_end_time:
# set a Dirichlet boundary condition of u=1 at left edge of the domain, but no stimulus term S
return [fem.dirichletbc(1.0, left_dofs[0], ME.sub(0))], 0.0
elif t > stimstarttime and t < stimstarttime + stimduration:
def stimulus(x):
return stimamp*np.all([(x[0] < 1.25),(x[1] < 1.25)],axis=0)
return [], stim0.interpolate(stimulus)
else:
# no stimulus and no Dirichlet boundrary conditions
return [], 0.0
My problem is that, in the main function solving for each timestep I have the following:
u_p, v_p = Up.sub(0), Up.sub(1)
U = fem.Function(ME)
dU = ufl.TrialFunction(ME)
(psiu, psiv) = ufl.TestFunctions(ME)
u, v = ufl.split(U)
a = fem.Constant(domain, 0.1)
eps = fem.Constant(domain, 0.01)
beta = fem.Constant(domain, 0.5)
gamma = fem.Constant(domain, 1.0)
D = fem.Constant(domain, 1e-4)
dtc = fem.Constant(domain, dt)
stimf = fem.Function(ME)
stim0, stim1 = stimf.split() # .split() returns a function object, split(stimf) return an indexable (non dolfin!) object
bcs, stim = BC_U_and_stim(domain,stim0,t)
def RHS_E1(v, w):
return - D*ufl.inner(ufl.grad(u), ufl.grad(psiu)) + psiu *u*(1-u)*(u-a) - w*psiu + stim*psiu
def RHS_E2(v, w):
return eps*(beta*v - gamma*w)*psiv
F1 = (u - u_p)/dtc * psiu * ufl.dx - (theta * RHS_E1(u, v_p) + (1-theta) * RHS_E1(u_p, v_p)) * ufl.dx
F2 = (v - v_p)/dtc * psiv * ufl.dx - (theta * RHS_E2(u, v) + (1-theta) * RHS_E2(u_p, v_p)) * ufl.dx
But I keep getting the error
return - D*ufl.inner(ufl.grad(u), ufl.grad(psiu)) + psiu *u*(1-u)*(u-a) - w*psiu + stim*psiu
~~~~^~~~~
TypeError: unsupported operand type(s) for *: 'NoneType' and 'Indexed'
Since I am new both to FEniCS and FEniCSx and this is my first porting from the legacy to the maintained library, I do not have any clue on how to solve it and if this is the best practice to port such a code. Thank you in advance