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

252

views

0

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

Never mind, it works now! Thanks for the help

written
7 months ago by
Mats Brun

1

Please, add comments using "Add comment" button rather than a new answer.

written
7 months ago by
Jan Blechta

### 1 Answer

3

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.
# Define external load
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

gives

`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

Please login to add an answer/comment or follow this question.

and added the following terms to my variational form

where skw() is defined as in your code. However, the term

now gives me the error "This integral is missing an integration domain.". Any idea what is wrong?

Thanks again.