The code below is based on the hyperelasticity demo but breaks at J = Form(J) with
/usr/local/lib/python2.7/dist-packages/ffc/analysis.pyc in _attach_integral_metadata(form_data, parameters)
258 # Update scheme for QuadratureElements
259 if all_equal(quad_schemes):
--> 260 scheme = quad_schemes[0]
261 else:
262 scheme = "canonical"
IndexError: list index out of range
Something to do with quad_schemes that I don't understant. My hyperelasticity demo works fine. Any help?
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": False, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
class Primitive(Expression):
def __init__(self,f,mesh,**kwargs):
self._f = f
self._mesh = mesh
subdomains = CellFunction('size_t',mesh)
self._subdomains = subdomains
dx = Measure('dx',domain=mesh)[subdomains]
self._dx=dx
self._F = Form(f*dx(1))
def eval(self, values, x):
inside_function = lambda y: y<x[0]
domain = AutoSubDomain(inside_function=inside_function)
self._subdomains.set_all(0)
domain.mark(self._subdomains,1)
F = self._F
values[0] = assemble(F)
mesh = IntervalMesh(100,0,1.0)
V = FunctionSpace(mesh, "CG", 1)
right = CompiledSubDomain("(std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary")
class LinearMismatchProfile(Expression):
def __init__(self, G, W, W_b, f_0,**kwargs):
self._W = W
self._G = G
self._W_b = W_b
self._f_0 = f_0
def eval(self, values, x):
W = self._W
W_b = self._W_b
f_0 = self._f_0
G = self._G
if x[0] <= W-W_b:
values[0] = G*x[0] + f_0
else:
values[0] = G*(W-W_b)
Ed = Constant(8.138931143006406e-09) #energy per unit length of dislocation
c = Constant(92017082108.12033) #elastic constant for biaxial strain
b = Constant(2.161453136448972e-10) #part of Burger's vector responsible for relaxation (in-plane component)
W = 1e-6
W_b = 100e-9
G = 1e5 #PQ(0.1,'1/mum').inBaseUnits().value
f_0 = 0.0
f = LinearMismatchProfile(G=G,W=W,W_b=W_b,f_0=f_0)
u = Function(V) # dislocation density profile
v = TestFunction(V)
du = TrialFunction(V)
epsilon = f - b*Primitive(f=u,mesh=mesh,) #in-plane strain
E = 2*Ed*u*dx(mesh) + epsilon**2*dx(mesh)
F = derivative(E,u,v)
J = derivative(F, u, du)
F = Form(F)
J = Form(J)
bcr = DirichletBC(V,0.0,right)
bcs = [bcr]
# Solve non-linear variational problem
solve(F == 0, u, bcs, J=J,
form_compiler_parameters=ffc_options,
solver_parameters={'newton_solver':
{"linear_solver": "lu",
"maximum_iterations": 100,
'relative_tolerance': 1e-30
}})
# Plot and hold solution
plot(u, interactive = False)
