Hi! This is my very first project using FEniCS and I’m having some issues when running my code. The error I receive is the following:
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
Traceback (most recent call last):
File “partc1.py”, line 38, in
u0 = interpolate(indata,V)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/interpolation.py”, line 71, in interpolate
Pv.interpolate(v._cpp_object)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py”, line 383, in interpolate
self._cpp_object.interpolate(u)
RuntimeError:
*** -------------------------------------------------------------------------
*** Error: Unable to interpolate function into function space.
*** Reason: Rank of function (0) does not match rank of function space (1).
*** Where: This error was encountered inside FunctionSpace.cpp.
This is what my code looks like:
from dolfin import *
from fenics import *
# Create mesh and define function space
mesh = Mesh("circle.xml")
# Construct the finite element space
V = VectorFunctionSpace (mesh, 'P', 1)
# Define parameters :
T = 150
dt = 0.5
alpha = 0.4
beta = 2
gamma = 0.8
delta = 1
# Class representing the intial conditions
class InitialConditions(UserExpression):
def eval(self,values,x):
#values[0] = Expression(("(4/25)-2*pow(10,-7)*(x[0]-0.1*x[1]-225)*(x[0]-0.1*x[1]-675)"),degree=2)
#values[1] = Expression(("(22/45)-3*pow(10,-5)*(x[0]-450)-1.2*pow(10,-4)*(x[1]-150)"),degree=2)
values = Expression(("(4/25)-2*pow(10,-7)*(x[0]-0.1*x[1]-225)*(x[0]-0.1*x[1]-675)","(22/45)-3*pow(10,-5)*(x[0]-450)-1.2*pow(10,-4)*(x[1]-150)"), degree=2)
def value_shape(self):
return (2 ,)
# Define initial condition
indata = InitialConditions(degree=2)
u0 = Function(V)
u0 = interpolate(indata,V)
# Test and trial functions
u,v = TrialFunction(V), TestFunction(V)
# Create bilinear and linear forms
s0 = (pow(u0[0],2) + (u0[0]*u[1])/(u0[0]+alpha))*dx
s1 = (u0[0]*u[1])/(u0[0]+alpha)*dx
a0 = (u[0]*v[0]*dx) - (dt*u[0]*v[0]*dx) - (0.5*delta*dt*inner(grad(u[0]),grad(v[0]))*dx)
a1 = (u[1]*v[1]*dx) + (dt*gamma*u[1]*v[1]*dx)(0.5*delta*dt*inner(grad(u[1]),grad(v[1]))*dx)
L0 = (u0[0]*v[0]*dx) + (dt*u0[0]*v[0]*dx) + (0.5*delta*dt*inner(grad(u0[0]),grad(v[0]))*dx) - (dt*u0[0]*v[0]*dx*s0)
L1 = (u0[1]*v[1]*dx) - (dt*gamma*u0[1]*v[1]*dx) + (0.5*delta*dt*inner(grad(u0[1]),grad(v[1]))*dx) + (dt*beta*u0[1]*v[1]*dx*s1)
a = a0 + a1
L = L0 + L1
#Set up boundary condition
g = Constant([0.0,0.0])
bc = DirichletBC(V,u_initial,DirichletBoundary())
bc = [] #NEUMANN
#Assemble matrix
A = assemble(a)
# Set an output file
out_file = File("Solution_partc1.pvd","compressed")
# Set initial condition
u = Function(V)
u.assign(u0)
t = 0.0
out_file << (u,t)
u_initial = Function(V)
u_initial.assign(u0)
t_save = 0
num_samples = 20
# Time - stepping
while t < T:
# assign u0
u0.assign(u)
#Assemble vector and apply boundary conditions
A = assemble(a)
b = asseble(L)
t_save += dt
if t_save > T/ num_samples or t >= T-dt:
print("Saving!")
#Save the solution to file
out_file << (uv,t)
t_save = 0
#Move to next interval and adjust boundary condition
t += dt
I think that my error is located where I set my initial conditions and define values[0] and values[1], however I am not completely sure.
Thanks in advance for any tips you might have!
UPDATE!
I noticed that def values_shape was not inside the class which triggered the first error, however I’m having a new problem now and this the error I get:
Expecting a tuple of Index objects to A^indices := as_tensor(A, indices).
Traceback (most recent call last):
File “partc1.py”, line 52, in
L0 = (u0[0]v[0]dx) + (dtu0[0]v[0]dx) + (0.5deltadtinner(grad(u0[0]),grad(v[0]))dx) - (dtu0[0]*v[0]dx)((u0[0]^2) + (u0[0]*u[1])/(u0[0]+alpha)*dx)
File “/usr/lib/python3/dist-packages/ufl/exproperators.py”, line 91, in _as_tensor
error(“Expecting a tuple of Index objects to A^indices := as_tensor(A, indices).”)
File “/usr/lib/python3/dist-packages/ufl/log.py”, line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Expecting a tuple of Index objects to A^indices := as_tensor(A, indices).
This is my updated code:
from dolfin import *
from fenics import *
# Create mesh and define function space
nx = ny = nz = 30
mesh = UnitSquareMesh(nx,ny)
# Construct the finite element space
V = VectorFunctionSpace (mesh, 'P', 1)
# Define parameters :
T = 150
dt = 0.5
alpha = 0.4
beta = 2
gamma = 0.8
delta = 1
# Class representing the intial conditions
class InitialConditions ( UserExpression ):
def eval (self , values , x):
values[0] =(4/25)-2*pow(10,-7)*(x[0]-0.1*x[1]-225)*(x[0]-0.1*x[1]-675)
values[1] =(22/45)-3*pow(10,-5)*(x[0]-450)-1.2*pow(10,-4)*(x[1]-150)
def value_shape ( self ):
return (2 ,)
# Define initial condition
indata = InitialConditions(degree=2)
u0 = Function(V)
u0 = interpolate(indata,V)
# Test and trial functions
u,v = TrialFunction(V), TestFunction(V)
# Create bilinear and linear forms
a0 = (u[0]*v[0]*dx) - (dt*u[0]*v[0]*dx) - (0.5*delta*dt*inner(grad(u[0]),grad(v[0]))*dx)
a1 = (u[1]*v[1]*dx) + (dt*gamma*u[1]*v[1]*dx) - (0.5*delta*dt*inner(grad(u[1]),grad(v[1]))*dx)
L0 = (u0[0]*v[0]*dx) + (dt*u0[0]*v[0]*dx) + (0.5*delta*dt*inner(grad(u0[0]),grad(v[0]))*dx) - (dt*u0[0]*v[0]*dx)*((u0[0]^2) + (u0[0]*u[1])/(u0[0]+alpha)*dx)
L1 = (u0[1]*v[1]*dx) - (dt*gamma*u0[1]*v[1]*dx) + (0.5*delta*dt*inner(grad(u0[1]),grad(v[1]))*dx) + (dt*beta*u0[1]*v[1]*dx*(u0[0]*u[1])/(u0[0]+alpha)*dx)
a = a0 + a1
L = L0 + L1
If I change the u0[0]^2 to pow(u0[0],2) in L0 I get the following error:
Traceback (most recent call last):
File “partc1.py”, line 52, in
L0 = (u0[0]v[0]dx) + (dtu0[0]v[0]dx) + (0.5deltadtinner(grad(u0[0]),grad(v[0]))dx) - (dtu0[0]*v[0]dx)(pow(u0[0],2) + (u0[0]*u[1])/(u0[0]+alpha)*dx)
TypeError: unsupported operand type(s) for +: ‘Power’ and ‘Form’