How to add a constant stress tensor to the stress expression? (linear elasticity)

1

Hi, a new user here, so a very basic question.
I have started to modify an elasticity example code. The initial code describes a steady-state deformation of a rectangular 3D beam, clamped on one end and deforming under its own gravity.
I have changed it to 2D and then wanted to include pre-stress, basically adding a constant (2,2) tensor to the stress function (see the code below).
I could not make the shapes of the addition arguments to match. I have also tried adding a constant scalar to a scalar field (Poisson example) in a similar fashion (taking the trial function u, and defining a new one) and it worked well, however, adding a Constant((1,0)) vector to the gradient of the test function didn’t work.
What am I doing wrong here? Thanks in advance :slight_smile:

MWE code:

from fenics import *
from ufl import nabla_div, nabla_grad

# Scaled variables
L = 1; W = 0.2 # length and width
mu, lambda_ = 1, 1.25 #Lame coeffs 
rho = 1
delta = W/L
g = 0.4*delta**2
# Create mesh and define function space
mesh = RectangleMesh(Point(0, 0), Point(L, W), 50, 15, diagonal = "crossed")
V = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)
# Define strain and stress
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
def sigma(u):
    # The problem is in the line below. It works perfectly fine if the last term is not there
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u) + as_tensor(((1,2),(3,4)))

##### Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, -rho*g))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f,v)*dx 

u = Function(V)
solve(a == L, u, bc)

The error is:

Traceback (most recent call last):
  File "MWE_elasticity.py", line 46, in <module>
    solve(a == L, u, bc)
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 241, in _solve_varproblem
    problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,
  File "/usr/lib/python3/dist-packages/dolfin/fem/problem.py", line 56, in __init__
    a = Form(a, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
  File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_1',).
1 Like
2

Hello
Your bilinear part (a) should contain both trial and test functions. As you are adding a constant value to your stress term, then definition of the bilinear part is not correct in your implementation. So instead of:

a = inner(sigma(u), epsilon(v))*dx
L = dot(f,v)*dx

You can replace:

F = inner(sigma(u), epsilon(v))*dx - dot(f,v)*dx
a, L = lhs(F), rhs(F)
1 Like
3

It works, thanks a lot! You have saved me a lot of hours of frustration :slight_smile:
If I understand correctly, my added term becomes linear
inner(my_constant_tensor,epsilon(v))*dx
and probably moves to the right hand side?

4

That is right. As it only includes the test function, then this term should appear in the linear part (e.g. L).

5

Hi Leo

I have a similar issue, but i need to add a depth dependent stress tensor like:

t=np.array([[22.26,4.29],[4.29,13.8]])
T=as_tensor(t)
y_insitu=849.9 #depth stress in situ medition
y_max=1179.735 #upper limit
stress_in = Expression(("(y_max-x[1])/(y_max-y_insitu)",0),degree=1,y_max=y_max,y_insitu=y_insitu)
#Functions

def eps(v):
  return sym(grad(v))
def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) + T*stress_in

When i add the tensor as a constant its works, but when apply the depht dependent function y have the error: UFLException: Can’t add expressions with different shapes.

I hope you can help me pls

6

From the error it is clear that you are not allowed to multiply an expression of rank 1 (vector) to a tensor (rank 2). In other words, the shapes must be consistent. You may want to take a look at Expression class in FEniCS

Here is a working example:

from fenics import *
from ufl import nabla_div, nabla_grad

import numpy as np


# Scaled variables
L = 1; W = 0.2 # length and width
mu, lmbda = 1, 1.25 #Lame coeffs 
rho = 1
delta = W/L
g = 0.4*delta**2
# Create mesh and define function space
mesh = RectangleMesh(Point(0, 0), Point(L, W), 50, 15, diagonal = "crossed")
V = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)



x = SpatialCoordinate(mesh)


t=np.array([[22.26,4.29],[4.29,13.8]])
T=as_tensor(t)
y_insitu=849.9 #depth stress in situ medition
y_max=1179.735 #upper limit

S = as_tensor([[(y_max-x[1])/(y_max-y_insitu), 0],[0, (y_max-x[1])/(y_max-y_insitu)]])

def eps(v):
  return sym(grad(v))

def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) + T*S


##### Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, -rho*g))
F = inner(sigma(u), eps(v))*dx - dot(f,v)*dx
a, L = lhs(F), rhs(F)

u = Function(V)
solve(a == L, u, bc)
2 Likes
7

Thanka a lot, it works perfect!