Hi, I’ve read quite a few posts on here about properly setting up a nonlinear variational problem in a mixed finite element space but to no avail I cannot get my code to run.
I believe I’ve set up Function(), TrialFunctions(), TestFunctions() correctly for the mixed element spaces and have used the derivative() properly for my residual and Jacobian. I’ve read somewhere that it’s tricky using derivative() for a mixed function space, maybe that’s it?
The current error thrown is ufl.log.UFLException: Coefficient and argument shapes do not match! when trying to get my residual. Any help is greatly appreciated.
# ----------------------------------------------------------------
# FEniCS implementation: Cosserat fiber-reinforced elasticity
# Written by: Ryan McAvoy F24
# ----------------------------------------------------------------
from dolfin import *
from dolfin import mesh
import ufl
import sys
import time
import math
import numpy as np
from mpi4py import MPI
tic = time.time()
# ----------------------------------------------------------------
# Input parameters
# ----------------------------------------------------------------
# Material parameters
mu = 20.0e6 # Matrix shear modulus [Pa]
nu = 0.33 # Matrix Poisson's ratio
kappa = 2*mu*(1+nu)/(3*(1-2*nu)) # Matrix bulk modulus
lamda = 2*mu*nu/(1 - 2*nu) # Matrix first Lame Parameter
r = 1e-3 # Fiber cross-section radius [m]
E = 23e9 # Fiber Young's modulus [Pa]
A = math.pi*r**2 # Fiber cross-section area [m^2]
I = math.pi*(r**4)/2 # Fiber moment of inertia [m^4]
k = E*A # Fiber axial stiffness [N]
F = E*I # Fiber flexural stiffness [Nm^2]
# Define the domain corners
x_min = 0.0
x_max = 1.0
y_min = 0.0
y_max = 1.0
# Define the number of divisions in the x and y directions
nx = 100
ny = 100
alpha = 0.2
beta = 2*math.pi*20
D1 = Expression("1 / (pow((a*sin(b*x[0])),2) + 1)", degree=1, a = alpha, b = beta)
D2 = Expression("a*sin(b*x[0]) / (pow((a*sin(b*x[0])),2) + 1)", degree=1, a = alpha, b = beta)
# ----------------------------------------------------------------
# Read mesh
# ----------------------------------------------------------------
mesh = RectangleMesh(Point(x_min, y_min), Point(x_max, y_max), nx, ny)
dim = mesh.geometric_dimension()
mesh_coord = mesh.coordinates()
# ----------------------------------------------------------------
# Define function spaces
# ----------------------------------------------------------------
u_elem = VectorElement('CG', mesh.ufl_cell(), 2) # displacement
theta_elem = FiniteElement('CG', mesh.ufl_cell(), 1) # rotation
lambda_elem = FiniteElement('CG', mesh.ufl_cell(), 1) # Lagrange multiplier
mixedUT = MixedElement([u_elem, theta_elem, lambda_elem])
V = FunctionSpace(mesh, mixedUT)
U, T, Lam = V.split()
U_0, U_1 = U.split() # 2 dofs for disp. (U_0, U_1) / 1 dof for rotation (T)
# ----------------------------------------------------------------
# Define boundary conditions
# ----------------------------------------------------------------
left = CompiledSubDomain("near(x[0], mesh_xmin) && on_boundary", mesh_xmin = x_min)
right = CompiledSubDomain("near(x[0], mesh_xmax) && on_boundary", mesh_xmax = x_max)
# Dirichlet boundary
BC_left_x0 = DirichletBC(U_0, Constant(0.0), left)
BC_left_y0 = DirichletBC(U_1, Constant(0.0), left)
BC_left_theta = DirichletBC(T, Constant(0.0), left)
BC_right_x0 = DirichletBC(U_0,Constant(1.0), right)
BC_right_y0 = DirichletBC(U_1,Constant(0.0), right)
BC_right_theta = DirichletBC(T, Constant(0.0), right)
BC_right_x0 = DirichletBC(U_0,Constant(1.0), right)
BC = [BC_left_x0, BC_left_theta, BC_left_y0, BC_right_x0, BC_right_theta, BC_right_y0]
# ----------------------------------------------------------------
# Define energy
# ----------------------------------------------------------------
uh = Function(V)
u, theta, lam = split(uh)
x_test = TestFunctions(V)
delx = TrialFunctions(V)
rot = as_tensor([[cos(theta), sin(theta)],[-sin(theta), cos(theta)]])
Ften = Identity(2) + grad(u)
Eten = (rot.T)*Ften
C = Eten.T * Eten
IC = det(C)
J = det(Eten)
phi = as_vector([ theta.dx(0), theta.dx(1) ])
phiprime = phi[0]*D1 + phi[1]*D2
ED = as_vector([ D1*Eten[0,0] + D2*Eten[0,1], D1*Eten[1,0] + D2*Eten[1,1] ])
Lambda2 = ED[0]**2 + ED[1]**2
fibep = (Lambda2-1)/2
DpED = -ED[0]*D2 + ED[1]*D1
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(IC - 3) - mu*ln(J) + (lamda/2)*(ln(J))**2 + (k/2)*(fibep)**2 + (F/2)*(phiprime**2) + lam*DpED
Pi = psi*dx
# ----------------------------------------------------------------
# Variational problem
# ----------------------------------------------------------------
F = derivative(Pi, uh, x_test)
Jac = derivative(F, uh, delx)
solve(F == 0, uh, BC, J=Jac)
u_new, theta_new, lam_new = uh.split()
vtkfile_u << u_new
vtkfile_theta << theta_new
vtkfile_lamda << lam_new
toc = time.time() - tic
print('Elapsed CPU time: ', toc, '[sec]')
