Error in implementing weak form for a hyperelastic time-dependent PDE

General
1

Hi all,

I have a rectangular domain, immersed in a fluid, with the corner points at (0,0), (5,0), (5,0.1), and (0,0.1) that is exponentially growing in the X direction under the growth tensor diag(e^{rt}, 1, 1) over a period of 1 minute. On this domain, I need to solve \eta\frac{du}{dt}=-\nabla \Pi[u], where \mathbf{u}(x,y) is the displacement vector, \eta is the friction coefficient from the fluid-substrate interaction, and \Pi is the stored potential energy in the system due to growth. A zero displacement Dirichlet BC is imposed at the center of the beam (2.5, 0.05, 0.0).
Discretizing the time derivative and multiplying with a test function \mathbf{v}, the equation becomes \eta\int_{\Omega}\mathbf{v}\mathbf{u}dA+\Delta t\int_{\Omega}\mathbf{v}\nabla\Pi[\mathbf{u}]dA-\eta\int_{\Omega}\mathbf{v}\mathbf{u}_{(n)}dA, where \mathbf{u}_{n} is the solution from the previous time step. In the below MWE, I am calculating \Pi[\mathbf{u}] with the first Piola-Kirchoff stress tensor.
I am getting error when calculating the integral \int_{\Omega}\mathbf{v}\nabla\Pi[\mathbf{u}]dA

MWE:

msh = meshio.read("mesh.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)

####################  Finding the closest vertex to (2.5, 0.05, 0.0)  ############

center = Point(2.5,0.05,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)

############################  Defining functions  #########################

u = Function(V_bact) 
u_n = Function(V_bact) 
v = TestFunction(V_bact)

####################################################################

N = FacetNormal(mesh)

####################################################################

NumStepg = 10
Time = 60  # Units: mins
dt = NumStepg/Time
j = 0
eta = 0.5 # Friction coefficient
r = 0.02  # [r] = /min
mu1 = 10.0             

for T in np.linspace(0,60,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]])  
 
 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)  # Strain energy density function for Neo Hookean model
 
 TT = det(g_bact)*diff(W,Fe)*(inv(g_bact).T)  #   1st PK Stress
 
 Pi1 = inner(TT,grad(v))*dx_custom(5)   # Stored potential energy, Pi, of the system
 
 grad_Pi = grad(Pi1)
 
 VF = eta*dot(u,v)*dx_custom(5) + dt*inner(v,grad_Pi)*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 Y- component of displacement zero
 bcs = [bcs1,bcs2]
  
 # Solve variational problem
 solve(PE == 0, u , bcs,
      solver_parameters={"newton_solver": {"relative_tolerance": 5e-3,   
      "absolute_tolerance": 5e-4,"maximum_iterations": 50}})

 j += 1
 
 u_n.assign(u)

When executed, it generates the error

grad_Pi = grad(Pi1)
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 369, in grad
    f = as_ufl(f)
  File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 437, in as_ufl
    raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: { ({ A | A_{i_{16}, i_{17}} = sum_{i_{18}} ({ A | A_{i_{14}, i_{15}} = (d/d[var0({ A | A_{i_8, i_9} = sum_{i_{10}} (I + (grad(f_37)))[i_8, i_{10}] * (([
  [1.0, 0, 0],
  [0, 1.0, 0],
  [0, 0, 1.0]
])^-1)[i_{10}, i_9]  })] (5.0 * (-3 + (tr(var1({ A | A_{i_{11}, i_{12}} = sum_{i_{13}} (var0({ A | A_{i_8, i_9} = sum_{i_{10}} (I + (grad(f_37)))[i_8, i_{10}] * (([
  [1.0, 0, 0],
  [0, 1.0, 0],
  [0, 0, 1.0]
])^-1)[i_{10}, i_9]  }))[i_{13}, i_{12}] * ((var0({ A | A_{i_8, i_9} = sum_{i_{10}} (I + (grad(f_37)))[i_8, i_{10}] * (([
  [1.0, 0, 0],
  [0, 1.0, 0],
  [0, 0, 1.0]
])^-1)[i_{10}, i_9]  }))^T)[i_{11}, i_{13}]  }))))))[i_{14}, i_{15}] * (det([
  [1.0, 0, 0],
  [0, 1.0, 0],
  [0, 0, 1.0]
])) })[i_{16}, i_{18}] * ((([
  [1.0, 0, 0],
  [0, 1.0, 0],
  [0, 0, 1.0]
])^-1)^T)[i_{18}, i_{17}]  }) : (grad(v_0)) } * dx(<Mesh #0>[5], {}) can not be converted to any UFL type.

Also, I am confused about the weak form of the integral \Delta t\int_{\Omega}\mathbf{v}\nabla\Pi[\mathbf{u}]dA

Any suggestions would be extremely helpful.

Thanks in advance.

2

It’s unlikely anyone will be able to help, because your MWE does not contain imports needed to run the code, and you have not provided the mesh to be used there. Consider replacing your custom mesh with a built in (like a cube) for ease of reading.

1 Like
3

Thanks @francesco-ballarin for the suggestion. It was my mistake not to post the entire code along with the Gmsh script. However, I have solved the problem. The equation was not well defined.