This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

arity-mismatch-exception

0 votes

Posting
Dear Fenics users, I am trying to assembly a linear system. It should combine linear elasticity and electrostatics to simulate piezoelectric materials.
I've first implemented the equations for the linear elasticity problem, which worked. For the next step i wanted to calculate the effect of a constant electric field.

As you can see i defined it as a 6x1 constant vector of zeros. When i substract that Constant from the defined sigma (as you cann see in the definition of sigma) i recieve the following error.

ArityMismatchTraceback (most recent call last)
in () 86 # Solve the problem 87 result = Function(VV) ---> 88 solve(a == L, result, boundary) 89 90 # Save solution to VTK format
...
ArityMismatch: Adding expressions with non-matching form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 6664), VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3)), 1, None),) vs ().

Does somebody has any idea what i'm doing wrong? Without subtracting E it works.

from dolfin import *

dielectric = 28000. * 8.8541878176e-12 * Identity(3)
# derived lame constants
Ey = 100.0  # youngs modulus
nu = 0.3   # poisson ratio
mu = Ey/(2.0*(1.0 + nu))
lmbda = Ey*nu/((1.0 + nu)*(1.0 - 2.0*nu))       

# define elasticity tensor
e1 = (1-nu)/(1.0-2.0*nu)
e2 = nu/(1.0-2.0*nu)
# 6x6
c = as_matrix([[e1, e2, e2, 0., 0., 0.], 
               [e2, e1, e2, 0., 0., 0.],
               [e2, e2, e1, 0., 0., 0.],
               [0., 0., 0., 1., 0., 0],
               [0., 0., 0., 0., 1., 0.],
               [0., 0., 0., 0., 0., 1.]]) 

# rewrite the tensor into a vector using Voigt notation
def as_voigt_vector(ten):
    # FEniCS does not know anything about Voigt notation, 
    # so one need to access the components directly as eps[0, 0] etc.   
    return as_vector([ten[0,0],ten[1,1],ten[2,2],ten[1,2],ten[0,2],ten[0,1]])

# rewrite a vector into a tensor using Voigt notation
def from_voigt_vector(vec):
    # FEniCS does not know anything about Voigt notation, 
    # so one need to rewrite the vector by hand into a tensor.   
    return as_tensor([[vec[0], vec[5], vec[4]], 
                      [vec[5], vec[1], vec[3]],
                      [vec[4], vec[3], vec[2]]])

# first define the function for the strain tensor in fenics notation
def epsilon(u):
    # note that it would be shorter to use sym(grad(u)) but the way used now
    # represents the formula
    return 1./2.*(grad(u) + transpose(grad(u)))

# now define the stress tensor in fencis notation
def sigma(u,E):     
    # calcuate sigma with C-matrix in Voigt notation
    # sigma = 6x1 vector; E = 6x1 vector
    sig     = (Ey/(1.+nu))*c* (as_voigt_vector(epsilon(u)) - E(u))    

    # sigma for coupled linear elasicity and electrostatics; d = 3x6 tensor with piezo material constants
    #sig = (Ey/(1.+nu))*c* (as_voigt_vector(epsilon(u))-(transpose(d_piezo)*E(u)))

    # return sigma as symmetric tensor
    return from_voigt_vector(sig)

import numpy 

#boundary conditions:
voltage = 10
q = Constant(0)
f = Constant((0., 0., 0.))
t = Constant((0., 0., -1.))

# build the trial and test function spaces over a domain 
# Setup domain
grid = 10
Omega = UnitCubeMesh(grid, grid, 5*grid)
x = Omega.coordinates()
x[:, 0] *= 0.1
x[:, 1] *= 0.1
x[:, 2] *= 0.01
Omega.bounding_box_tree().build(Omega) 

VV = VectorFunctionSpace(Omega, "Lagrange", 1)

# Get the needed boundaries: 
# 1. To apply a predefined stress on it We keep it simple: the top face sees the same stress vector. 
# 2. For zero Displacement we choose the bottom face of the cube
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.01)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], 0.)

top_boundary  = TopBoundary()
bottom_boundary = BottomBoundary()

# Create the boundary condition for displcement
boundary = []
boundary.append(DirichletBC(VV, Constant((0., 0., 0.)), bottom_boundary))

# As the stress is defined over the outer surface we need to create a 
# integral over that surface
boundaries = FacetFunction("size_t", Omega)
boundaries.set_all(0)
top_boundary.mark(boundaries, 1)
ds = Measure("ds", subdomain_data=boundaries)

E = lambda u: as_vector(Constant((0.,0.,0.,0.,0.,0.)))

# Define the variational problem: get all variables and implement function given above
u = TrialFunction(VV)
v = TestFunction(VV)

a = inner(sigma(u,E), grad(v))*dx        # typically called a
L = inner(f, v)*dx + inner(t, v)*ds(1)   # typically called L

# Solve the problem
result = Function(VV)
solve(a == L, result, boundary)

# Save solution to VTK format
File("displacement.pvd") << result
print "Done"
asked Sep 20, 2016 by andifl FEniCS Novice (120 points)
edited Sep 20, 2016 by andifl

1 Answer

0 votes

Hi,
consider to use

F = inner(sigma(u,E), grad(v))*dx - inner(f, v)*dx + inner(t, v)*ds(1)
a, L = lhs(F), rhs(F)

# Solve the problem
result = Function(VV)
solve(a == L, result, boundary)

The error message indicate that your form a isn't bilinear (because the term E(u) is constant).

answered Sep 21, 2016 by hernan_mella FEniCS Expert (19,460 points)
edited Sep 21, 2016 by hernan_mella

Yes, that works. Thank you for the response!

...