# Imposing initial condition BDM elements

**URL:** <https://fenicsproject.discourse.group/t/imposing-initial-condition-bdm-elements/1677>\
**Category:** variational formulation\
**Created:** [October 6, 2019, 2:42am UTC](https://fenicsproject.discourse.group/t/imposing-initial-condition-bdm-elements/1677 "2019-10-06T02:42:53Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Simone\_Puel](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/simone_puel/32/520_2.png) [@Simone\_Puel](https://fenicsproject.discourse.group/u/Simone_Puel)\
**Post date:** [October 6, 2019, 2:42am UTC](https://fenicsproject.discourse.group/t/imposing-initial-condition-bdm-elements/1677/1 "2019-10-06T02:42:53Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![ia267](https://avatars.discourse-cdn.com/v4/letter/i/ed8c4c/32.png) [@ia267](https://fenicsproject.discourse.group/u/ia267)\
**Post date:** [October 6, 2019, 9:27pm UTC](https://fenicsproject.discourse.group/t/imposing-initial-condition-bdm-elements/1677/2 "2019-10-06T21:27:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Simone\_Puel](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/simone_puel/32/520_2.png) [@Simone\_Puel](https://fenicsproject.discourse.group/u/Simone_Puel)\
**Post date:** [October 7, 2019, 5:23am UTC](https://fenicsproject.discourse.group/t/imposing-initial-condition-bdm-elements/1677/3 "2019-10-07T05:23:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ia267](https://avatars.discourse-cdn.com/v4/letter/i/ed8c4c/32.png) [@ia267](https://fenicsproject.discourse.group/u/ia267)\
**Post date:** [October 7, 2019, 12:36pm UTC](https://fenicsproject.discourse.group/t/imposing-initial-condition-bdm-elements/1677/4 "2019-10-07T12:36:01Z")

</div>

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!