Hi everybody!

I am using FEniCS 2019.1.0 running on MacOS. I am trying to solve a linear Maxwell viscoelasticity problem using a mixed formulation and following the approach of Rognes and Winther (2008), where stress is calculated using BDM elements, velocity and rotation DG.

I think I wrote the variational formulation correctly, but when I write the initial condition of zero stress, I have a problem in the assemble of the rhs. This is part of my code:

# define functions

def asym(z):

return z[0,1] - z[1,0]

def AEsigma(s):

return 1./(2.*mu)*(s - lmbda/(2.*mu + d*lmbda)*tr(s)*Identity(d))

def AVsigma(s):

return (1./(2.*eta))*s

# define function spaces

k = 1

BDM = VectorFunctionSpace(mesh, “BDM”, k) # stress

DGv = VectorFunctionSpace(mesh, “DG”, k - 1) # displacement

DGs = FunctionSpace(mesh, “DG”, k - 1) # rotation

# mixed Function Space

Vh_element = MixedElement([BDM.ufl_element(), DGv.ufl_element(), DGs.ufl_element()])

Vh = FunctionSpace(mesh, Vh_element)

n = FacetNormal(mesh)

# zero stress Dirichlet BCs

zero_tensor = Expression(((“0.0”, “0.0”), (“0.0”, “0.0”)), degree=1)

bc1 = DirichletBC(Vh.sub(0), zero_tensor, boundaries, topleft_id)

bc2 = DirichletBC(Vh.sub(0), zero_tensor, boundaries, topright_id)

bcs = [bc1, bc2]

# initial condition

sigma0 = Function(BDM)

v0 = Function(DGv)

Dt = Constant(dt)

# define trial and test functions

(sigma, v, gamma) = TrialFunctions(Vh)

(tau, w, q) = TestFunctions(Vh)

# write the weak form: lhs and rhs

lhs = inner(AVsigma(sigma)*Dt, tau)*dx + inner(AEsigma(sigma), tau)*dx + inner(div(tau), v*Dt)*dx

+ inner(asym(tau)*Dt, gamma)*dx

+ inner(div(sigma), w)*dx + inner(asym(sigma), q)*dx

rhs = inner(v0, tau*n)*ds(zerodispl_id)

+ inner(AEsigma(sigma0), tau)*dx

t = 0.

sol = Function(Vh)

while t < t_end:

```
t += dt
```

# apply prescribed displacement at t = dt (after first iteration)

```
if t == dt:
slip = -1
rhs += inner(slip, dot(tau('+')*n('+'), tangent(n)('+')))*dS(fault_id)
# Assemble the linear system
A = assemble(lhs)
b = assemble(rhs)
[bc.apply(A) for bc in bcs]
[bc.apply(b) for bc in bcs]
solve(A, sol.vector(), b)
(sigma, v, gamma) = sol.split(deepcopy=True)
# Update previous solutions
sigma0.assign(sigma)
```

If could anyone help me, I would appreciate it. Thank you so much in advance!