### Symmetric tensor space, alternative to Arnold-Winther?

252
views
0
7 months ago by
I have a question regarding symmetric tensor valued space in a mixed setting. I would like to construct the following mixed space (the S-element should be a tensor element, but i think the "AW" option means it is tensor valued already?)
mesh = UnitSquareMesh(32, 32)

U = VectorElement('DG', mesh.ufl_cell(), 0)
S = FiniteElement('AW', mesh.ufl_cell(), None)
P = FiniteElement('DG', mesh.ufl_cell(), 0)
W = FiniteElement('RT', mesh.ufl_cell(), 1)
T = FiniteElement('DG', mesh.ufl_cell(), 0)
R = FiniteElement('RT', mesh.ufl_cell(), 1)
elem = MixedElement([U, S, P, W, T, R])   # mixed element
MS = FunctionSpace(mesh, elem)            # mixed space​

def boundary(x, on_boundary):
return on_boundary

# Dirichlet BC
bc_su = DirichletBC(MS.sub(0), Constant((0.0, 0.0)), boundary)
bc_wp = DirichletBC(MS.sub(2), Constant(0.0), boundary)
bc_rt = DirichletBC(MS.sub(4), Constant(0.0), boundary)
bcs = [bc_su, bc_wp, bc_rt]

# trial and test functions
v, tau, q, z, S, y = TestFunctions(MS)
u, sig, p, w, T, r = TrialFunctions(MS)

but apparently the "AW" option is not available. In my variational form I have a term which involves div(sig), so my question is:
is there another way I can construct the mixed space I need without using the "AW"-option? Or perhaps the "AW"-option will be included in the near future?

Many thanks,
Mats

Community: FEniCS Project
Thanks for the reply. I have now edited my code as follows
mesh = UnitSquareMesh(32, 32)
mesh_size = mesh.hmax()

# define function spaces
U_el = VectorElement('DG', mesh.ufl_cell(), 0)
S_el = VectorElement('BDM', mesh.ufl_cell(), 1)
P_el = FiniteElement('DG', mesh.ufl_cell(), 0)
W_el = FiniteElement('RT', mesh.ufl_cell(), 1)
T_el = FiniteElement('DG', mesh.ufl_cell(), 0)
R_el = FiniteElement('RT', mesh.ufl_cell(), 1)
x_el = FiniteElement('DG', mesh.ufl_cell(), 0)
el = MixedElement([U_el, S_el, P_el, W_el, T_el, R_el, x_el])         # mixed element
MS = FunctionSpace(mesh, el)                                          # mixed space

def boundary(x, on_boundary):
return on_boundary

# Dirichlet BC for displacement and pressure
bc_su = DirichletBC(MS.sub(0), Constant((0.0, 0.0)), boundary)
bc_wp = DirichletBC(MS.sub(2), Constant(0.0), boundary)
bc_rt = DirichletBC(MS.sub(4), Constant(0.0), boundary)
bcs = [bc_su, bc_wp, bc_rt]

# trial and test functions
v, tau, q, z, S, y, e = TestFunctions(MS)
u, sig, p, w, T, r, x = TrialFunctions(MS)​

and added the following terms to my variational form
inner(skw(x), tau)*dx + inner(sig, skw(e))*dx​

where skw() is defined as in your code. However, the term
inner(sig,tau)*dx​
now gives me the error "This integral is missing an integration domain.". Any idea what is wrong?

Thanks again.
written 7 months ago by Mats Brun
Never mind, it works now! Thanks for the help
written 7 months ago by Mats Brun
1
written 7 months ago by Jan Blechta

3
7 months ago by
I don't think the AW elements are fully supported in FEniCS.  But you can use the mixed finite elements for elasticity with weakly imposed symmetry (often called the AFW elements).  Here's a simple program that uses them.  Note that they work fine in 2D as well.

"""
Solve the 2D elasticity problem

A sigma = eps(u),  div(sigma) = f

using the mixed formulation with weak symmetry and AFW elements.
Here, only homogeneous displacement boundary conditions are given.
"""

from dolfin import *

# skew matrix determined by scalar
def skw(r):
return as_matrix([[0, r], [-r, 0]])

# Create a mesh
mesh = UnitSquareMesh(40, 40);

# set material parameters
mu = 4.
lam = 1.

f = Expression(("sin(2*pi*x[0])", "0."), degree=5)

# Define function spaces
deg = 1
BDMxBDM = VectorElement('BDM', mesh.ufl_cell(), deg)  # BDM matrices
DGxDG = VectorElement('DG', mesh.ufl_cell(), deg-1) # DG vectors
DG = FiniteElement('DG', mesh.ufl_cell(), deg-1)
X = FunctionSpace(mesh, MixedElement(BDMxBDM, DGxDG, DG))

# Create trial and test functions
(sigma, u, r) = TrialFunctions(X)
(tau, v, q) = TestFunctions(X)

# Define variational formulation
a = .5/mu * ( inner(sigma, tau) - lam/(2*mu + 2*lam) * tr(sigma) * tr(tau) ) * dx \
+ dot(u, div(tau))*dx + inner(skw(r), tau)*dx \
+ inner(div(sigma), v)*dx + inner(sigma, skw(q))*dx
L = inner(f, v)*dx
x = Function(X)

# Solve the system
solve(a == L, x)
(sigma, u, r) = x.split()
plot(u[0])
interactive()
​
Indeed

import ffc
import pprint
pprint.pprint(ffc.supported_elements)

gives

['Brezzi-Douglas-Fortin-Marini',
 'Brezzi-Douglas-Marini',
 'BrokenElement',
 'Bubble',
 'Crouzeix-Raviart',
 'Discontinuous Lagrange',
 'Discontinuous Raviart-Thomas',
 'Discontinuous Taylor',
 'EnrichedElement',
 'Gauss-Legendre',
 'Gauss-Lobatto-Legendre',
 'HDiv Trace',
 'Hellan-Herrmann-Johnson',
 'Lagrange',
 'Nedelec 1st kind H(curl)',
 'Nedelec 2nd kind H(curl)',
 'NodalEnrichedElement',
 'Raviart-Thomas',
 'Regge',
 'TensorProductElement']
written 7 months ago by Jan Blechta
A small typo.  I meant to say that they work fine in 3D as well.
written 7 months ago by Douglas N Arnold