In the below MWE, I am working with a unit square mesh, created with Gmsh (Gmsh script given after the MWE), that grows according to a growth tensor. I am solving the problem \eta \frac{\partial \mathbf{u}}{\partial t}=\nabla. \mathbf{T}, where \eta is the friction coefficient acting on the entire square (the square is immersed in a fluid while it grows), \mathbf{u}=(u(x,y), 0.0, 0.0) is the displacement vector in the X direction only, and \mathbf{T} is the first Piola-Kirchoff stress tensor. The growth of the square is governed by the growth tensor \mathbf{F}_{g}=diag(e^{rt}, 1.0, 1.0), where r is the rate of growth and t is time. In the MWE, I am ramping up time using a \texttt{for} loop. Moreover, the square grows according to the Dirichlet boundary conditions that \mathbf{u}=(0.0, 0.0, 0.0) at the node of the mesh nearest to the center (0.5, 0.5, 0.0), the Y and Z displacements are zero at all time points, and a Neumann boundary condition that \mathbf{T}.\mathbf{n}=0 for the entire boundary of the square where \mathbf{n} = FacetNormal(mesh). I am implementing a hyperelastic setup with a strain energy density function.
According to theoretical calculations, the (2,2) and (3,3) elements of the first Piola-Kirchoff stress tensor should be zero. However, when I am visualizing it in Paraview, I am getting non-zero (2,2) and (3,3) components of the first Piola-Kirchoff stress tensor. I tried P1 and DG0 spaces as well. But the problem persists.
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=True):
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("squaremesh.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_custom = 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_custom = Measure("dx", domain=mesh, subdomain_data=cf)
###################### Defining vector and scalar function spaces ##################
V_bact = VectorFunctionSpace(mesh, 'P', 1)
F_bact = FunctionSpace(mesh, 'P', 1)
T_bact = TensorFunctionSpace(mesh, 'P', 1)
#T_bact = TensorFunctionSpace(mesh, 'DG', 0) # Tried DG 0 as well
#################### Calculating Center of Mass for imposing 0 displacement DBC #################
new_V = VectorFunctionSpace(mesh,"P",1)
new_R = VectorFunctionSpace(mesh,"R",0)
new_position = Function(new_V)
new_position.assign(Expression(["x[0]","x[1]","x[2]"], element=new_V.ufl_element()))
new_c = TestFunction(new_R)
new_bulk = assemble(Constant(1.0)*dx(domain=mesh))
new_centroid = assemble(dot(new_c,new_position)*dx)
new_f = new_centroid/new_bulk
new_f_np = new_f.get_local()
center = Point(0.5,0.5,0.0)
# Find which cell point is located in
tree = mesh.bounding_box_tree()
cell, distance = tree.compute_closest_entity(center)
vertices = mesh.cells()[cell]
# Find the closest vertex in the cell
vertex_coordinates = mesh.coordinates()[vertices]
dist = 1e-1
closest = None
for i, v1 in enumerate(vertex_coordinates):
p_ = Point(v1)
dist_ = center.distance(p_)
print(v1)
if dist_ < dist:
dist = dist_
closest_local = i
closest_vertex = vertices[closest_local]
closestv = mesh.coordinates()[closest_vertex]
######################## Defining function to impose DBC at center ##################
# Find dof corresponding to vertex
vtd_p = vertex_to_dof_map(V_bact)
dof_number = vtd_p[closest_vertex]
print(f"DOF number of closest vertex = {dof_number}")
print(f"closest vertex = {V_bact.tabulate_dof_coordinates()[closest_vertex]}")
class cm(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return near(x[0], closestv[0], tol) and near(x[1], closestv[1], tol) and near(x[2], closestv[2], tol)
####################################################################
vtk6 = File("PK1_ref.pvd")
############################ SOLVING FOR div(v) = 0 with v = grad(u)########################################
u = Function(V_bact)
u_n = Function(V_bact)
v = TestFunction(V_bact)
####################################################################
N = FacetNormal(mesh)
####################################################################
NumStepg = 50
Time = 60 # Units: mins
dt = NumStepg/Time
j = 0
eta = 5.0 # Friction coefficient
r = 0.02 # [r] = /min
mu1 = 1.5 # Stiffness parameter
nu1 = 0.499 # Measure of incompressibility
for T in np.linspace(0,Time,NumStepg):
d = mesh.geometry().dim()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
g_bact = as_tensor([[exp(r*T),0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]]) # Growth Tensor
inv_g_bact = as_tensor([[1.0/exp(r*T),0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]]) # Inverse of the growth tensor
Fe = variable(F*inv_g_bact)
C = variable(Fe.T*Fe)
J = det(Fe)
M = cofac(F)
# Invariants of deformation tensors
I1 = tr(C)
W = (mu1/2)*(I1 - 3) + ((nu1)*(mu1)/(1-2*nu1))*(J-1)**2 - mu1*(ln(J))
PK1 = det(g_bact)*diff(W,Fe)*(inv_g_bact.T) # 1st PK Stress for compressible morphoelasticity
PE = eta*dot(u,v)*dx_custom(5) + dt*inner(PK1,grad(v))*dx_custom(5) - eta*dot(u_n, v)*dx_custom(5)
all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
bcs1 = DirichletBC(V_bact, Constant((0.0,0.0,0.0)), cm(), method='pointwise') # Zero displacement Dirichlet BC near center
bcs2 = DirichletBC(V_bact.sub(1), Constant(0.0), all_domains, 1) # making sure that Y- component of displacement zero
bcs3 = DirichletBC(V_bact.sub(2), Constant(0.0), all_domains, 1) # making sure that Z- component of displacement zero
bcs = [bcs1,bcs2,bcs3]
# Solve variational problem
solve(PE == 0, u , bcs,
solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,
"absolute_tolerance": 5e-10,"maximum_iterations": 50}})
PK1_ref = project(PK1, T_bact)
PK1_ref.rename('PK1_ref','PK1_ref')
vtk6<<(PK1_ref,T)
print(f"Number of iterations: {j}", flush=True)
j += 1
u_n.assign(u)
#######################################################################
print('initial bacteria area =', assemble(Constant(1.0)*dx_custom(5)))
ALE.move(mesh, u)
print('final bacteria area =', assemble(Constant(1.0)*dx_custom(5)))
Gmsh Script:
SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.0, 0, 1.0};
//+
Point(2) = {1.0, 0.0, 0, 1.0};
//+
Point(3) = {1.0, 1.0, 0, 1.0};
//+
Point(4) = {0.0, 1.0, 0, 1.0};
//+
Line(1) = {4, 3};
//+
Line(2) = {3, 2};
//+
Line(3) = {2, 1};
//+
Line(4) = {1, 4};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("top", 1) = {1};
//+
Physical Curve("right", 2) = {2};
//+
Physical Curve("bottom", 3) = {3};
//+
Physical Curve("left", 4) = {4};
//+
Physical Surface("body", 5) = {1};
Here is what the mesh looks like:
I am not sure what I am doing wrong here. Any suggestions, ideas, or comments would be extremely helpful.
Please let me know if you need any further information.
Thanks in advance!