Hello,
I am having mixed form problem.
from dolfin import *
import numpy as np
from numpy import linalg as LA
mesh = UnitCubeMesh(5, 5, 5)
# bc
def left_boundary(x):
return x[0] < DOLFIN_EPS
def right_boundary(x):
return (x[1] < (1.0 - DOLFIN_EPS))
V = VectorFunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "CG", 1)
Z = MixedFunctionSpace([V,W,Q])
z = TrialFunction(Z)
(u,c,b) = split(z)
zn = TestFunction(Z)
(du, dc, db) = split(zn)
z0 = Function(Z)
(u0, c0, b0) = split(z0)
f = Constant((1.0, 0.0, 0.0))
# params
E, nu = 10.0, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
l = Constant(0.7, name="l")
gc = Constant(1e-3, name="gc")
k = Constant(1.0e-6, name="eps")
kappa = Constant(180.0, name="kappa")
eta = 1e-3
ffp = lambda d: (d+abs(d))/2.0
ffm = lambda d: (d-abs(d))/2.0
def g(c):
return (1-c)**2
def eps(v):
return sym(grad(v))
def sigma_0(v):
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/(1.0 - nu**2)
return 2.0*mu*(eps(v)) + lmbda*tr(eps(v))*Identity(3)
def sigma(u,c):
return (g(c)+k)*sigma_0(u)
def elastic_energy(u):
return 0.5*inner(sigma(u,c),eps(u))
dt = 1e-6
F1 = inner(sigma(u, c), grad(du))*dx
F2 = inner(l*grad(c), grad(dc))*dx
F3 = c*db*dx - c0*db*dx + dt*inner(-1/eta * ffp(b + 2 * (1-c) *(elastic_energy(u))/gc - c/l), db) * dx
F = F1 + F2 + F3
zz = Function(Z)
dF = derivative(F,z, zn)
tension = Constant((0.0, -1.0, 0.0))
fixed = Constant((0.0, 0.0, 0.0))
bct = DirichletBC(V, tension, right_boundary)
bcf = DirichletBC(V, fixed, left_boundary)
bcs = [bct, bcf]
zz0 = Function(Z)
# c = Function(W)
# c0 = Function(W)
iter = 0
maxiter = 2
# Iterations
# while iter<maxiter:
# c0.vector()[:] = c.vector()[:]
solve(F == 0, zz0, bcs, J=dF)
# iter=iter+1
But I am running to following error, does anybody have suggestion, what could be possibly wrong?
Invalid coefficient type for Argument(MixedElement(*[VectorElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, dim=3, quad_scheme=None), FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, quad_scheme=None), FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, quad_scheme=None)], **{'value_shape': (5,) }), 1, None)
Traceback (most recent call last):
File "linear_phase_field_without_decomposition.py", line 78, in <module>
dF = derivative(F,z, zn)
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/fem/formmanipulations.py", line 80, in derivative
return ufl.derivative(form, u, du, coefficient_derivatives)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/formoperators.py", line 244, in derivative
coefficients, arguments = _handle_derivative_arguments(form, coefficient, argument)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/formoperators.py", line 197, in _handle_derivative_arguments
ufl_assert(isinstance(c, Indexed), "Invalid coefficient type for %s" % repr(c))
File "/home/fenics/build/lib/python2.7/site-packages/ufl/assertions.py", line 37, in ufl_assert
if not condition: error(*message)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Invalid coefficient type for Argument(MixedElement(*[VectorElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, dim=3, quad_scheme=None), FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, quad_scheme=None), FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, quad_scheme=None)], **{'value_shape': (5,) }), 1, None)
