Imposing initial condition BDM elements

variational formulation
1

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 + dlmbda)*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), vDt)*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!

2

I’m not sure about your specific formulation, but I found the Cahn Hilliard demo to be very helpful for reference when applying initial conditions in a mixed formulation.

3

Thank you so much! I think the problem may arise due to the fact that BDM is in VectorFunctionSpace. How could I write the initial condition of a VectorFunctionSpace or rank 2? It has to be something like:

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

And then interpolate it in the VectorFunctionSpace? Thanks.

4

In order to do that, you must define an InitialCondition Class:

class InitialConditions(Expression):
def eval(self, values, x):
    # For sigma0
    values[0] = Expression(("0.0", "0.0"), degree=1)
    values[1] = Expression(("0.0", "0.0"), degree=1)
    # For v 
    values[2] = 
    # For gamma
   values[3] =
def value_shape(self):
   # Change accordingly 
    return (4,)

You must also define an initial condition for your displacement and rotation. Interpolation works the same as the Cahn-Hilliard Demo from there, except you must define the degree.

# Call the initial class 
init = InitialConditions(degree=?)
Vh.interpolate(init)

Hope that helps!