Hi community, I developed a mixed scheme for a nonlinear elasticity problem. The nonlinearity is in the load term and involves a function composition, so it didn't work to use the nonlinear solver. To solve the problem, I implemented the Newton scheme by hand, and it is showing a strange behaviour. As a note, a primal formulation of the problem works correctly, so the issue is in this particular formulation.
I have extracted the portion with the problem, because at first using 'solve' just wouldn't work.
from dolfin import *
M = MixedElement([VectorElement('BDM', mesh.ufl_cell(), 1),
VectorElement('DG', mesh.ufl_cell(), 0),
FiniteElement('DG', mesh.ufl_cell(), 0) ])
W = FunctionSpace(UnitSquareMesh(10,10), M)
tau, v, eta = TestFunctions(W)
w_n = Function(W)
w_n.vector()[:] = 0.
sigma_n, u_n, gamma_n = w_n.split()
#u_n = Function(VectorFunctionSpace(mesh, 'DG', 0))
#sigma_n = Function(VectorFunctionSpace(mesh, 'BDM', 1))
Meta = as_matrix([[0., eta],[-eta, 0.]])
b_below = (dot(v, div(sigma_n)) + inner(Meta, sigma_n))*dx
After running this, for the solver to work, I would need to assemble the form, but executing
assemble(b_below)
throws the following:
Traceback (most recent call last):
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 142, in jit
result = ffc.jit(ufl_object, parameters=p)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 204, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 120, in jit_build
generate=jit_generate)
File "/usr/lib/python2.7/dist-packages/dijitso/jit.py", line 160, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 66, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File "/usr/lib/python2.7/dist-packages/ffc/compiler.py", line 141, in compile_form
prefix, parameters, jit)
File "/usr/lib/python2.7/dist-packages/ffc/compiler.py", line 183, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 94, in analyze_ufl_objects
for form in forms)
File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 94, in
for form in forms)
File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 174, in _analyze_form
do_apply_restrictions=True,
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/compute_form_data.py", line 274, in compute_form_data
form = apply_function_pullbacks(form)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/apply_function_pullbacks.py", line 266, in apply_function_pullbacks
return map_integrand_dags(FunctionPullbackApplier(), expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 58, in map_integrand_dags
form, only_integral_type)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
for itg in form.integrals()]
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 46, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 57, in
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/lib/python2.7/dist-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/usr/lib/python2.7/dist-packages/ufl/corealg/map_dag.py", line 84, in map_expr_dags
r = handlersv._ufl_typecode_
File "/usr/lib/python2.7/dist-packages/ufl/corealg/multifunction.py", line 41, in _memoized_handler
r = handler(self, o)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/apply_function_pullbacks.py", line 255, in form_argument
return apply_single_function_pullbacks(o)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/apply_function_pullbacks.py", line 153, in apply_single_function_pullbacks
assert len(gsh) == 1
AssertionError
For some extra information, this error appears with the following comments above it:
# By placing components in a list and using as_vector at the end,
# we're assuming below that both global function g and its
# reference value r have vector shape, which is the case for most
# elements with the exceptions:
# - TensorElements
# - All cases with scalar subelements and without symmetries
# are covered by the shortcut above
# (ONLY IF REFERENCE VALUE SHAPE PRESERVES TENSOR RANK)
# - All cases with scalar subelements and without symmetries are
# covered by the shortcut above
# - VectorElements of vector-valued basic elements (FIXME)
# - TensorElements with symmetries (FIXME)
But I am not sure about how to treat this.
I also tried using the commented form instead of splitting, but it didn't work. Also, both summands of 'b_below' throw errors, so it isn't some interaction between them. Finally, if I try to compile a form with u_n or gamma_n, it works fine, so there is something with the Tensor elements (or row wise vector) as the comment suggests.
Please let me know of any ideas you have. Best regards!
