Hi, I am trying to solve a time-dependent non-linear elasticity problem. The boundary conditions are implemented using a Nitsche method which results in a nonlinear problem. Discretization results in a mixed problem, which has to be solved in each iteration.
I am not sure, if derivative() can also handle the part
F2 = gamma/h**2 * dot(ppos(dot(u, n)), dot(phi2, n)) * ds(1)
My idea was to use derivative on the form F1 and provide the derivative of F2 by hand and sum up both to get the Jacobian for the NonLinearVariational solver.
Unfortunately I get the error:
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).
Is this the right idea to handle the derivative?
I provided a MWE below.
Thank you
from vtkplotter.dolfin import plot as dolfin_plot
import numpy as np
import matplotlib.pyplot as plt
from ufl import nabla_div
import sys
from dolfin import *
import ufl
from mshr.cpp import Rectangle, generate_mesh
if __name__ == '__main__':
#TODO: https://fenicsproject.org/qa/12129/computing-jacobian-manually-coupled-problem-hilliard-equation/
k = 10
# Define geometry
l_end = -3.0
r_end = 3.0
t_end = 2.0
b_end = 0.0
m_t_end = 1.0
m_l_end = -1.5
m_r_end = 1.5
l_b_end = 1.0
r_b_end = 1.0
t_l_r_end = -0.75
t_r_l_end = 0.75
rect1 = Rectangle(Point(l_end, l_b_end), Point(t_l_r_end, t_end))
rect2 = Rectangle(Point(t_r_l_end, r_b_end), Point(r_end, t_end))
rect3 = Rectangle(Point(m_l_end, b_end), Point(m_r_end, m_t_end))
rect = rect1 + rect2 + rect3
mesh = generate_mesh(rect, k)
# Plot mesh
#dolfin_plot(mesh, title="Mesh")
# Degree of interpolation polynomials
degree = 2
# parameters
rho = 1
# Lamé constants
mu = 1
beta = 1.25
lambda_ = beta
# constants
gamma_nitsche = Constant(100.0)
# Time stepping parameters
T_start = 0
T_final = 10 # final time
t = 0
num_steps = 10 # number of time steps
dt = (T_final - T_start)/num_steps # time step size
# Define function space for system of PDEs
P1 = VectorElement('Lagrange', cell=triangle, degree=degree)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
tol = 1E-14
def clamped_boundary(x, on_boundary): # beam is only fixed on the left end
return on_boundary and near(x[0], 0.0, tol) #x[0] < tol
u_D1 = Constant((0, 0))
bc1 = DirichletBC(V.sub(0), u_D1, clamped_boundary)
bc2 = DirichletBC(V.sub(1), u_D1, clamped_boundary)
def force_on_rhs_boundary(x, on_boundary):
return on_boundary and near(x[0], 60.0, tol)
class topbot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and 0.0 + tol < x[0] < 60.0 - tol
# force on rhs boundary
# f1 = 1
# u_D2 = Constant((f1, 0))
# time dependent force: f(t) = 1-e^{-\kappa t}
kappa = 0.5
tt = np.linspace(0, 10, 100)
S = 1.0
SW = 0
y = S-(S-0)*np.exp(-kappa*tt)
fct = str(S) + '-(' + str(S) + '-' + str(SW) + ')*exp(-kappa*t)'
u_D22 = Expression((fct, '0'), degree=1, kappa=kappa, t=0)
u_D21 = Constant((0, 0))
bc3 = DirichletBC(V.sub(0), u_D21, force_on_rhs_boundary)
bc4 = DirichletBC(V.sub(1), u_D22, force_on_rhs_boundary)
# Mark boundaries as "ds"
boundaries = MeshFunction("size_t", mesh, dim=1)
boundaries.set_all(0)
# Mark boundaries where u*n=0 as "1"
topbot = topbot()
topbot.mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Boundary conditions
bcs = [bc3, bc4]
# Define strain and stress
def epsilon(z): # strain
return 0.5*(nabla_grad(z) + nabla_grad(z).T)
# return sym(nabla_grad(z))
def sigma(z): # stress
return lambda_*nabla_div(z) * Identity(d) + 2 * mu * epsilon(z)
# Define test functions
vu = Function(V)
v, u = split(vu)
print(type(u))
vnun = Function(V)
v_n, u_n = split(vnun)
(phi1, phi2) = TestFunctions(V)
d = v.geometric_dimension()
# Define source terms
f_1 = Constant((0, 0))
f_2 = Constant((0, 0))
# Define variational problem
gamma = Constant(gamma_nitsche)
h = CellDiameter(mesh)
n = FacetNormal(mesh)
def ppos(x):
return (x+abs(x))/2
# Weak formulation
F1 = (dot(v-v_n, phi2) * dx
- dt * inner(sigma(u), epsilon(phi2)) * dx
- dt * dot(v, phi1) * dx
+ dot(u-u_n, phi1) * dx)
# Add Nitsche term:
F1 -= dot(dot(n, dot(sigma(u), n)), dot(phi2, n)) * ds(1)
# Add penalty term
F2 = gamma/h**2 * dot(ppos(dot(u, n)), dot(phi2, n)) * ds(1)
F = F1 + F2
vuh = Function(V)
dvu = Function(V)
dv, du = split(dvu)
for n in range(num_steps):
# Update current time
t += dt
print('t: ', t)
# Update bc
u_D22.t = t
#J = derivative(F, vuh, du)
# Compute derivative of F for Newton method
R = action(F1, vuh)
DR = derivative(R, vuh, dvu)
print(type(DR))
def heavyside(u):
return conditional(ufl.operators.And(ge(u[0], 0), ge(u[1], 0)), 1, 0)
DR2 = heavyside(u) * dot(du, phi1) * dx
#DR2 = action(DR2, vuh)
# Sum over directional derivatives:
DR = DR + DR2
problem = NonlinearVariationalProblem(R, vuh, bcs, DR)
solver = NonlinearVariationalSolver(problem)
solver.solve()
# Split solutions
vh, uh = vuh.split(True)
U = VectorFunctionSpace(mesh, 'P', degree)
print('Lösung u: Dolfin Plot')
uh_plot = project(uh, U)
dolfin_plot(uh_plot, style='paraview', lw=0)
U = FunctionSpace(mesh, 'P', degree)
u_magnitude = sqrt(dot(uh_plot, uh_plot))
u_magnitude = project(u_magnitude, U)
print('Displacement magnitude')
print('Max/Min Displacement')
print('min/max u:',
u_magnitude.vector().min(),
u_magnitude.vector().max())
# Update previous solution
v_n = vh
u_n = uh