I am trying to implement a model for nonlinear elasticity where I need to calculate the Frobenius norm of the deviatoric part of the strain tensor, as well as quantities inside sqrt and exp functions. As my model gets more and more complex, I also get more of the ArityMismatch errors, which has become a kind of a barrier that keeps me from understanding the logic behind forms, UFL and general implementation using Fenics.
Am I getting these errors because I am violating a rule that I do not know yet? I am guessing that there must be some limitations to any algebra system that seeks to be elegant, general and efficient at the same time, so is it the case; or am I just using it wrong? When I try to construct the equation using a == L, the shapes match and it succeeds, but I get an error when I try to solve the equation. Below is the code,
from __future__ import print_function
from dolfin import *
import time
import numpy
# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
# Sub domain for clamp at left end
def left(x, on_boundary):
return x[0] < 0.001 and on_boundary
# Sub domain for rotation at right end
def right(x, on_boundary):
return x[0] > 0.99 and on_boundary
# Load mesh and define function space
mesh = UnitCubeMesh(20, 20, 20)
# Define function space
V = VectorFunctionSpace(mesh, "CG", 1)
# Test and trial functions
u = TrialFunction(V)
delta_u = TestFunction(V)
lmbda = 10.
mu_0 = 10.
mu_inf = 2.
eta = 2.
# External forces (body and applied tractions
f = Constant((10.0, 0.0, 0.0))
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
force_boundary = AutoSubDomain(right)
force_boundary.mark(boundary_subdomains, 3)
# Define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)
# Set up boundary condition at left end
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, left)
def eps(u):
return sym(grad(u))
def tensor_norm(T):
# return 1 # errors disappear when the norm terms are set as constant
return sqrt(inner(T,T))
# Stress tensor
def sigma(u):
return lmbda * tr(eps(u)) * Identity(len(u)) \
+ (2 * mu_hat(u) + d_mu_hat_d_dev_eps(u) * tensor_norm(dev(eps(u)))) * dev(eps(u))
# The terms in the nonlinear deviatoric part
def mu_hat(u):
return mu_0 + (mu_inf - mu_0) * (1 - exp(-1 * tensor_norm(dev(eps(u))) / eta))
def d_mu_hat_d_dev_eps(u):
return (mu_inf - mu_0) * exp(-1 * tensor_norm(dev(eps(u))) / eta)
# Forms
a = inner(sigma(u), eps(delta_u))*dx
L = inner(f, delta_u)*dss(3)
# Time-stepping
u = Function(V)
xdmf_file = XDMFFile(mesh.mpi_comm(), "out.xdmf")
start = time.time()
print("Started solution")
solve(a == L, u,
bcs=bc,
solver_parameters={"linear_solver": "cg"},)
And I get such an error when I try to run it:
fenics@90d9768619d6:~/shared/linear_elastic$ python nonlinear_elasticity.py
Started solution
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
File "nonlinear_elasticity.py", line 88, in <module>
solver_parameters={"linear_solver": "cg"},)
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/fem/solving.py", line 300, in solve
_solve_varproblem(*args, **kwargs)
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/fem/solving.py", line 325, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/fem/solving.py", line 82, in __init__
a = Form(a, form_compiler_parameters=form_compiler_parameters)
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/fem/form.py", line 94, in __init__
mpi_comm=mesh.mpi_comm())
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 65, in mpi_jit
return local_jit(*args, **kwargs)
File "/home/fenics/build/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 124, in jit
result = ffc.jit(ufl_object, parameters=p)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/jitcompiler.py", line 201, in jit
module = jit_build_with_instant(ufl_object, module_name, parameters)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/jitcompiler.py", line 98, in jit_build_with_instant
code_h, code_c = jit_generate(ufl_object, module_name, parameters)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/jitcompiler.py", line 62, in jit_generate
parameters=parameters, jit=True)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/compiler.py", line 154, in compile_form
analysis = analyze_forms(forms, parameters)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/analysis.py", line 60, in analyze_forms
parameters) for form in forms)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/analysis.py", line 60, in <genexpr>
parameters) for form in forms)
File "/home/fenics/build/lib/python2.7/site-packages/ffc/analysis.py", line 142, in _analyze_form
form_data = compute_form_data(form)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/algorithms/compute_form_data.py", line 362, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is
File "/home/fenics/build/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 143, in check_form_arity
check_integrand_arity(itg.integrand(), arguments)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 137, in check_integrand_arity
args = map_expr_dag(rules, expr, compress=False)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/corealg/map_dag.py", line 82, in map_expr_dags
r = handlers[v._ufl_typecode_](v)
File "/home/fenics/build/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 34, in nonlinear_operator
raise ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.__name__, t))
ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator Exp to expression depending on form argument v_1.
Aborted
I tried to read the UFL documentation in order to get some clues, but it did not help. Do you have any other material that can explain?
Thanks in advance.
