My problem involves solving the equation on the form:
$$\frac{\partial^2 \vec{u}}{\partial t^2} = \frac{1}{\rho}\nabla \cdot P$$
Where u is the displacement vector, $\rho$ is the density and P is the stress tensor given as:
$$P = \lambda(\nabla \cdot \vec{u})I + \mu (\nabla \vec{u} + \nabla \vec{u}^T)$$
I solve this in two domains with different material properties $\lambda$, $\mu$ and $\rho$. I use an explicit second order difference method for the time variable, and FEM for the spatial variables. This gives a variational form:
$$\int_{\Omega}\vec{u}^{(n+1)}\cdot \vec{v}d\Omega = \int_{\Omega}(2\vec{u}^{(n)} - \vec{u}^{(n-1)})\cdot \vec{v} d\Omega + \int_{\Gamma}\vec{n} \cdot P^{(n)} \cdot \vec{v} d\Gamma - \
\int_{\Omega}P^{(n)} : \nabla \vec{v} d\Omega $$
where n denotes the value at time n. I have attached a simplified code with given initial conditions.
from dolfin import *
# Constants
ys = 1
t = 0
dt = 0.01
T = 10
rho1 = 1
rho2 = 10
mu1 = 1
mu2 = 0
lamda1 = 1
lamda2 = 1
# Mesh, functionspace, functions
mesh = RectangleMesh(0,0,2,2,20,20)
V = VectorFunctionSpace(mesh, "CG", 1)
D = VectorFunctionSpace(mesh, "DG", 0)
u = TrialFunction(V)
v = TestFunction(V)
# Define Subdomains
cf = CellFunction("size_t", mesh, 0)
domain1 = AutoSubDomain(lambda x: x[1] > ys)
domain1.mark(cf,1)
# Define variable constants
rhof = Expression(("x[1] > ys ? rho1 : rho2",
"x[1] > ys ? rho1 : rho2"),
ys=ys,rho1=rho1,rho2=rho2)
muf = Expression(("x[1] > ys ? mu1 : mu2",
"x[1] > ys ? mu1 : mu2"),
ys=ys,mu1=mu1,mu2=mu2)
lamdaf = Expression(("x[1] > ys ? lamda1 : lamda2",
"x[1] > ys ? lamda1 : lamda2"),
ys=ys,lamda1=lamda1,lamda2=lamda2)
stepf = Expression(("x[1] > ys ? pow(dt,2)/rho1 : pow(dt,2)/rho2",
"x[1] > ys ? pow(dt,2)/rho1 : pow(dt,2)/rho2"),
ys=ys,dt=dt,rho1=rho1,rho2=rho2)
rho = interpolate(rhof, D)
mu = interpolate(muf, D)
lamda = interpolate(lamdaf, D)
step = interpolate(stepf, D)
# Stress tensor
def sigma(u, lamda, mu):
return lamda*div(u)*Identity(2) + mu*(grad(u) + grad(u).T)
# Initial conditions
Ixy = Constant((1,0))
Vxy = Constant((2,0))
u2 = interpolate(Ixy, V)
u1 = interpolate(Vxy, V)
# Variational forms
F = inner(u,v)*dx - 2*inner(u1,v)*dx + inner(u2,v)*dx + \
step*inner(sigma(u1,lamda,mu),grad(v))*dx
t = 2*dt
left = assemble(lhs(F))
u = Function(V)
# Time stepping
while t < T + DOLFIN_EPS:
plot(u2, rescale=False)
right = assemble(rhs(F))
begin("solving at time step t=%g" % t)
solve(left, u.vector(), right)
end()
u2.assign(u1)
u1.assign(u)
I'm guessing fenics does not like the expression objects for the variable constants, but some advice is needed here. Is there a way to get around this with the current definitions, or is there a way to get around this by following the same logic as in the fenics tutorial for vectorized functions as described in the chapter "Handling Domains with different materials"? Any help is greatly appreciated!
The error message is
daniel@satelite:~/masteroppgave/py$ python question.py
Product can only represent products of scalars.
Traceback (most recent call last):
File "question.py", line 59, in
stepinner(sigma(u1,lamda,mu),grad(v))dx
File "question.py", line 49, in sigma
return lamdadiv(u)Identity(2) + mu(grad(u) + grad(u).T)
File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line 181, in _mul
return _mult(self, o)
File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line 154, in _mult
p = Product(a, b[ii])
File "/usr/lib/python2.7/dist-packages/ufl/algebra.py", line 174, in new
error("Product can only represent products of scalars.")
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 154, in error
raise self._exception_type(self._format_raw(message))
ufl.log.UFLException: Product can only represent products of scalars.
daniel@satelite:~/masteroppgave/py$
