# Dirichlet boundary conditions in NonlinearVariationalProblem

**URL:** <https://fenicsproject.discourse.group/t/dirichlet-boundary-conditions-in-nonlinearvariationalproblem/3507>\
**Category:** Uncategorized\
**Created:** [June 12, 2020, 5:21pm UTC](https://fenicsproject.discourse.group/t/dirichlet-boundary-conditions-in-nonlinearvariationalproblem/3507 "2020-06-12T17:21:51Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![kausalyam](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kausalyam/32/861_2.png) [@kausalyam](https://fenicsproject.discourse.group/u/kausalyam)\
**Post date:** [June 12, 2020, 5:21pm UTC](https://fenicsproject.discourse.group/t/dirichlet-boundary-conditions-in-nonlinearvariationalproblem/3507/1 "2020-06-12T17:21:51Z")

</div>

I am trying to apply dirichlet boundary conditions that vary at each time step to a nonlinear variational problem but keep seeing the following error when I change the boundary condition inside the time stepping loop.

```auto
*** Error: Unable to check whether point is inside subdomain.
*** Reason: Function inside() not implemented by user.
*** Where: This error was encountered inside SubDomain.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  

```

If I remove the boundary condition definition from outside the loop I do not get an error, but the boundary condition is not applied either. I followed the non-linear poisson demo in the boundary condition definitions.  
A simplified version of my code is below:

```python
from dolfin import *
import numpy as np

parameters["form_compiler"]["quadrature_degree"] = 4

length, H = 1, 0.4

t = 2E-3
# ts = interpolate(Expression('t*(1.0 - (x[0]*x[0] + x[1]*x[1])/(a*a))', t=t, a=a, degree=2), FunctionSpace(mesh, 'CG', 2))
ts = Constant((t))

meshlength = int(100*length)
meshheight = int(100*H)

beta = 8e2#5e3 #dissapation constant
Tend = .003
Nsteps=int(Tend*1e3)
dt=Tend/(Nsteps+1)
betadt = beta/dt

phi_s=1615.5
b0=44898
b1=45156.3

#define quartic energy
cc2 = 41333
cc3 = -181690.6
cc4 = 214515.0
mu = 1150000
eps_s = 0.385

R1=Constant((1,0))
R2=Constant((0.5,sqrt(3)/2))
R3=Constant((-0.5,sqrt(3)/2))

nu = Constant(0.3)
D = 1e-1 *as_tensor([[1., nu, 0.],[nu, 1., 0.],[0., 0., (1. - nu)/2]]) #1e-1
S = Constant(3.6e3)

mesh = RectangleMesh(Point(0., 0.), Point(length, H), meshlength, meshheight)

P1 = FiniteElement("Lagrange", triangle, degree = 1)
P2 = FiniteElement("Lagrange", triangle, degree = 2)
bubble = FiniteElement("B", triangle, degree = 3)
enriched = P1 + bubble

element = MixedElement([VectorElement(P2, dim=2), P2, VectorElement(enriched, dim=2)])

''' -------------------------
BOUNDARY CONDITIONS
------------------------- '''
# Sub domain for Dirichlet boundary condition
class LeftDirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]) < DOLFIN_EPS and on_boundary

class RightDirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - length) < DOLFIN_EPS and on_boundary

''' -------------------------
SET UP VARIATIONAL PROBLEM
------------------------- '''

Q = FunctionSpace(mesh, element)

# Then, we define :py:class:`Function`, :py:class:`TrialFunction` and :py:class:`TestFunction` objects to express the variational forms and we split them into each individual component function::
e_v0 = Expression(("a*x[0]", "a*x[1]") ,a = eps_s,degree=1) #INITIAL STRAIN

e_w0 = Expression(("0*x[0]") ,degree=1)
e_theta0 = Expression(("0*x[0]", "0*x[1]") ,degree=1)

v0 = project(e_v0,Q.sub(0).collapse())
w0 = project(e_w0,Q.sub(1).collapse())
theta0 = project(e_theta0,Q.sub(2).collapse())

q_, q, q_t = Function(Q), TrialFunction(Q), TestFunction(Q)
assign(q_,[v0,w0,theta0])
v_, w_, theta_ = split(q_)

q_old = Function(Q)
assign(q_old,[v0,w0,theta0])
v_old, w_old, theta_old = split(q_old)

''' -------------------------
DIRICHLET BOUNDARY CONDITIONS
------------------------- '''

#Left and right boundary conditions
u_right = Expression(("a*x[0]","a*x[1]"), a = eps_s ,degree=1) #INITIAL STRAIN
u_left = Expression(("a*x[0]","a*x[1]"), a = eps_s ,degree=1) #INITIAL STRAIN

bcl = DirichletBC(Q.sub(0), u_left, LeftDirichletBoundary())
bcr = DirichletBC(Q.sub(0), u_right, RightDirichletBoundary())

rightvert = DirichletBC(Q.sub(1),Constant(0),RightDirichletBoundary())
leftvert = DirichletBC(Q.sub(1),Constant(0),LeftDirichletBoundary())

bcs = [bcr,bcl,rightvert,leftvert]
''' -------------------------
FUNCTIONS TO DEFINE MATERIAL MODEL
------------------------- '''
def eps(v,w):
    return sym(grad(v)) + 0.5*outer(grad(w), grad(w))

def deps(disp,disp_old,z,z_old):
    return (eps(disp,z)-eps(disp_old,z_old))

def eps_dev(disp,z):
    return eps(disp,z)-(tr(eps(disp,z))/2)*Identity(2)

def psi(eps_vol):
    eval = cc2*pow(eps_vol,2) + cc3*pow(eps_vol,3) + cc4*pow(eps_vol,4)
    return eval

def bistable_energy(eps_vol):
    return 3*psi(eps_vol)

def tr_cross(a,b):
    return a[0,0]*b[0,0]+a[0,1]*b[1,0]+a[1,0]*b[0,1]+a[1,1]*b[1,1]

def total_energy(eps_vol,eps_dev,eps_d):
    incr_energy = 3/16 * betadt*(2*tr_cross(eps_d,eps_d)+pow(tr(eps_d),2))
    return bistable_energy(eps_vol)+incr_energy

psi_m = total_energy(tr(eps(v_,w_))/2,eps_dev(v_,w_),deps(v_,v_old,w_,w_old))

#Simplified version without bending
Pi = psi_m*dx
dPi = derivative(Pi, q_, q_t)
J = derivative(dPi, q_, q)

''' NEWTON SOLVER'''
problem = NonlinearVariationalProblem(dPi,q_,bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["relative_tolerance"] = 1e-5
solver.parameters["newton_solver"]["absolute_tolerance"] = 1e-5

time = np.linspace(0, Tend, Nsteps+1)

for i, ti in enumerate(time):
    print("Timestep: "+ str(i)+"/"+str(Nsteps)+", time = "+str(ti))

    u_left = Expression(("a*x[0]+i*0.005","a*x[1]"), a = eps_s,i = i ,degree=1) # if I remove these two lines there is no error
    bcl = DirichletBC(Q.sub(0), u_left, LeftDirichletBoundary())

    bcs = [bcr,bcl,rightvert,leftvert]

    problem = NonlinearVariationalProblem(dPi,q_,bcs, J=J)
    q_old.vector()[:] = q_.vector()

    solver.solve()

    v_h, w_h, theta_h = q_.split(deepcopy=True)

```

---

<div class="post-metadata">

**Author:** ![bhaveshshrimali](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bhaveshshrimali/32/124_2.png) [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)\
**Post date:** [June 12, 2020, 7:03pm UTC](https://fenicsproject.discourse.group/t/dirichlet-boundary-conditions-in-nonlinearvariationalproblem/3507/2 "2020-06-12T19:03:56Z")

</div>

> [@kausalyam](#):
>
> ```auto
> u_left = Expression(("a*x[0]+i*0.005","a*x[1]"), a = eps_s,i = i ,degree=1) # if I remove these two lines there is no error
> bcl = DirichletBC(Q.sub(0), u_left, LeftDirichletBoundary())
> 
> bcs = [bcr,bcl,rightvert,leftvert]
> 
> problem = NonlinearVariationalProblem(dPi,q_,bcs, J=J)
> 
> ```

Hi  
For time dependent BCs, you can simply create a _time-dependent_ expression for `u_left` and update it inside the for loop. For instance,

```python
u_left = Expression(("a*x[0]+i*0.005","a*x[1]"), a = eps_s , i=0, degree=1) #INITIAL STRAIN
...
...
for i, ti in enumerate(time):
    print("Timestep: "+ str(i)+"/"+str(Nsteps)+", time = "+str(ti))
    u_left.i = i

```

instead of creating the DirichletBC object again. Note that there could be other issues with the formulation/solver, since running the problem as is with a _LU_ solver doesn’t yield convergence.