Dirichlet boundary conditions in NonlinearVariationalProblem

1

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.

*** 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:

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)
2

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,

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.

1 Like