I am in the process of learning two-field methods with the purpose of avoiding the issue of volumetric locking. I haven’t implemented a two-field method in FEniCS yet. This is my first try. I read a lot of posts and books and developed the following MWE with information collected from various posts.
In the following MWE, I have an annulus with inner radius 1 and outer radius 2. I assigned the physical tag 1 to the inner boundary and 2 to the outer boundary. I have assigned the physical tag 5 to the annular region. I created my geometries in Gmsh.
The MWE has two parts: (1) solves a hyperelastic problem using a single vector function space (works fine), (2) tries to solve the same problem using a mixed function space (spits errors).
The main feature of both the parts is that I am imposing a DirichletBC to make the third component of the displacement vector 0 to avoid any axial stretch. Thanks to @Kei_Yamamoto, @conpierce8, and @dokken, I am able to do it successfully for problem (1) with their suggestions in How to keep the third component of the displacement vector 0 while solving a hyperelasticity problem? - #11 by avishmj
However, for problem (2), I am not able to apply the same strategy.
Here is the MWE:
from dolfin import *
import meshio
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from ufl import cofac, sqrt
from scipy import linalg as la
from ufl import *
from ufl.classes import *
import pandas as pd
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
######################################################################################################################
msh = meshio.read("singleannulus.msh")
######################################################################################################################
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds = Measure("dS", domain=mesh, subdomain_data=mf) # Defining global boundary measure, for global surface measure, plug in "dS" instead of "ds", which is used for subdomain boundary measures
dx = Measure("dx", domain=mesh, subdomain_data=cf)
###################################### DEFING A SINGLE VECTOR FUNCTION SPACE ###########################################
V_ima = VectorFunctionSpace(mesh, 'P', 1)
############################### DEFINING MATERIAL PARAMETERS ########################
nu1 = 0.4
mu1 = 27.9
############################### (1): SOLVING FOR KINEMATICS USING THE SIMPLE VECTOR FUNCTION SPACE #####################
u1 = Function(V_ima)
v1 = TestFunction(V_ima)
d = mesh.geometry().dim()
I = Identity(d) # Identity tensor
F = I + grad(u1) # Deformation gradient u = (u1, u2, 0)
Fe = variable(F)
Ce = Fe.T*Fe
Je = det(Fe)
I1 = tr(Ce)
M = cofac(F)
# Stored strain energy density (incompressible neo-Hookean model)
psi1 = (mu1/2)*(I1 - 3) + ((nu1)*(mu1)/(1-2*nu1))*(Je-1)**2 - mu1*(ln(Je))
TT_i = diff(psi1,Fe)
Pi = inner(TT_i,grad(v1))*dx(5)
all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
bcs = DirichletBC(V_ima.sub(2), Constant(0.0), all_domains, 1) # Making the Z-component of displacement 0
solve(Pi == 0, u1 , bcs,
solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,
"absolute_tolerance": 5e-9,"maximum_iterations": 10}})
print(f"Solved IMA displacement using a single vector function space")
############################# SETTING UP THE SAME KINEMATICS WITH MIXED FUNCTION SPACE #############################
############################### DEFINING A MIXED FUNCTION SPACE ######################
P22 = VectorElement("P", mesh.ufl_cell(), 4)
P11 = FiniteElement("P", mesh.ufl_cell(), 1)
TH = MixedElement([P22,P11])
W = FunctionSpace(mesh, TH)
(u2, p2) = TrialFunctions(W)
(v2, q2) = TestFunctions(W)
############################## (2): SOLVING FOR KINEMATICS USING MIXED FUNCTION SPACE ##########################
d = mesh.geometry().dim()
I = Identity(d) # Identity tensor
F = I + grad(u2) # Deformation gradient u = (u1, u2, 0)
Fe = variable(F)
Ce = Fe.T*Fe
Je = det(Fe)
I1 = tr(Ce)
M = cofac(F)
# Stored strain energy density (incompressible neo-Hookean model)
psi1 = (mu1/2)*(I1 - 3) + p2*(Je-1)**2 - mu1*(ln(Je))
TT_i = diff(psi1,Fe)
Pi = inner(TT_i,grad(v2))*dx(5)
U = Function(W)
u2, p2 = U.split()
all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
bcs = DirichletBC(W.sub(2), Constant(0.0), all_domains, 1) # The error says it cannot extract the subspace
solve(Pi == 0, U , bcs,
solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,
"absolute_tolerance": 5e-9,"maximum_iterations": 10}})
print(f"Solved IMA displacement using a mixed function space")
u2_, p2_ = U.split()
#############################################################################################################
The second problem (2) spits the following error:
bcs = DirichletBC(W.sub(2), Constant(0.0), all_domains, 1)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/functionspace.py", line 119, in sub
raise ValueError("Can only extract SubSpaces with i = 0 ... %d" %
ValueError: Can only extract SubSpaces with i = 0 ... 1
I tried W.sub(0), which gives the error
Error: Unable to create Dirichlet boundary condition. *** Reason: Expecting a vector-valued boundary value but given function is scalar.
I also tried W.sub(1), which gives the error
raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b))) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_1',).
Since I am trying mixed function spaces for the first time, I am not sure why these errors happen and how to fix them.
Any suggestions would be greatly appreciated.
Thank you.
`