Linear Elasticity with Discontinuous Galerkin in C++ (Interior Penalty)

Question

I am trying to implement the linear elasticity problem with DG method with interior penalty term and using the c++ and uff combination.

My bilinear form is given by

$\sum\limits_K\int_{\partial K} \varepsilon (v) : C : \varepsilon (u) \ ds - \sum\limits_{e \in \mathcal{E}} \int_{e} \left( \lbrack v \rbrack : \lbrace C : \varepsilon (u) \rbrace + \lbrack u \rbrack : \lbrace C :\varepsilon (v) \rbrace - \delta_e \lbrack u \rbrack : \lbrack v \rbrack \right) \ ds $

And my bilinear/linear form in the ufl code is

a = inner(stress(u), strain(v))*dx 
   - inner(avg(dot(stress(u), n)), jump(v))*dS \
   - inner(dot(stress(u), n), v)*ds \
   - eta*inner(jump(u), avg(dot(stress(v), n)))*dS \
   - eta*inner(u, dot(stress(v), n))*ds \
   + penalty/((h('+')+h('-'))**beta)*inner(jump(u), jump(v))*dS \
   + penalty/(h**beta)*inner(u, v)*ds)

L = inner(f,v)*dx + inner(g,v)*ds

With the given information:

cell = triangle
element = VectorElement("DG", cell, 3)

# Trial, test and source functions
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
g = Coefficient(element)

# Normal component, mesh size and right hand side
n = cell.n         
h = cell.volume 

# Penalty term eta and penalty value delta
eta = Constant(cell)
beta = Constant(cell)
penalty = Constant(cell)

# Lame constants
lmbda = Constant(cell)
mu = Constant(cell)

# Define strain functions
def strain(u):
  return 0.5*(grad(u) + grad(u).T) 

# Define stress function
def stress(e):
  return 2.0*mu*strain(e) + lmbda*tr(strain(e))*Identity(cell.geometric_dimension())

In my linear form the last term has this $\delta_e$ penalty term which is depending on the length of each edge

$\delta_e = \delta * \operatorname{area}(e) * ( \frac{1}{vol(K^-)} + \frac{1}{vol(K^+)} )$

for K being the finite element (in this case a triangle) and $\delta > 0$ a fixed penalty parameter.

How do I implement this term of area(e) in my code?

Thanks in advance for your advice :)

Note: The code above is running but there is no convergence and if I change the degree of the VectorElements I get different results.

Linear Elasticity with Discontinuous Galerkin in C++ (Interior Penalty)

1 Answer