Hey again,
I’ve been running into a problem regarding the formulation of a linear elasticity problem. The issue may result from the way I am imposing the transformation conditions of the normal vector from the undeformed to deformed configuration.
I am currently trying to couple my 2D nonlinear reaction-diffusion PDE system with a pseudo steady-state linear elasticity problem using the following set of equations:
PDE Problem: \nabla \cdot \boldsymbol{\sigma} + \gamma \boldsymbol{\dot{u}} = 0
Boundary Condition: \boldsymbol{\sigma}^T \textbf{n} = \boldsymbol{\sigma} \cdot \textbf{n} = F_{\pm}(C_A) \textbf{n}
Coupling Term: F_{\pm}(C_A) = \frac{f_R}{1 + e^{-2 s_1 (C_A - s_0)}} = \frac{0.3}{1 + e^{-2 \cdot 10 \cdot (C_A - 1)}}
Piola Transform: \textbf{n} = \frac{J \textbf{F}^{-T} \textbf{N}}{\lVert J \textbf{F}^{-T} \textbf{N} \rVert}
\boldsymbol{\sigma} = \boldsymbol{\sigma}^e + \boldsymbol{\sigma}^v
\boldsymbol{\sigma}^e = \lambda \text{tr}(\varepsilon) \textbf{I} + 2 \mu \varepsilon
\boldsymbol{\sigma}^v = \beta (\lambda \text{tr}(\dot{\varepsilon}) \textbf{I} + 2 \mu \dot{\varepsilon})
\boldsymbol{\varepsilon} = \frac{1}{2} (\nabla \textbf{u} + (\nabla \textbf{u})^T)
The code for the linear elasticity portion MWE is as follows:
#--------
# IMPORTS
#--------
import ufl
import numpy as np
import gmsh
from petsc4py import PETSc
from mpi4py import MPI
import dolfinx
from dolfinx import fem, mesh, plot, log, default_scalar_type, geometry
from dolfinx.io import gmshio
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
import basix
#------------
# CREATE MESH
#------------
nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-1, -1]), np.array([1, 1])], [nx, ny], mesh.CellType.triangle)
#--------------------
# CREAT FEM FUNCTIONS
#--------------------
# Global Spatial Coordinates as a Vector-Valued Expression i.e. x[0] = x, x[1] = y
x_spatial = ufl.SpatialCoordinate(domain)
metadata = {"quadrature_degree": 4}
# Define Measure Over Cells of Mesh
dx = ufl.Measure('dx', domain=domain, metadata=metadata)
# Define Measure Over Boundary of Mesh
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
# Time Discretisation Process
t = 0 # Initial Time
T = 60 # Final Time
dt_original = 0.1
dt = fem.Constant(domain, dt_original) # Time Step Size
Nsteps = int(T/dt_original) # Number of Time Steps
#----------------------------
# LINEAR ELASTICITY FUNCTIONS
#----------------------------
# Vector-Valued Function Space for Linear Elasticity
v_cg1 = basix.ufl.element(family="Lagrange", cell=domain.topology.cell_name(), degree=1, shape=(domain.geometry.dim, ))
V = fem.functionspace(mesh=domain, element=v_cg1)
# Trial and Test Functions for Linear Elasticity Problem (Mechanics Vectors)
# Initial Conditions: Mechanical System is at Rest
v = ufl.TestFunction(V)
u = ufl.TrialFunction(V)
u_n = fem.Function(V)
u_n.name = "u_n"
#-----------------------------
# REACTION-DIFFUSION FUNCTIONS
#-----------------------------
# Mixed-Element Function Space for Reaction Diffusion
P1 = basix.ufl.element(family="Lagrange", cell=domain.topology.cell_name(), degree=1)
element = basix.ufl.mixed_element([P1, P1])
C = fem.functionspace(mesh=domain, element=element)
# Trial and Test Functions for Reaction-Diffusion Problem (Concentrations)
c = fem.Function(C)
c.name = "c"
c_A, c_I = ufl.split(c)
# Initial Species Concentration Conditions of Species I=2 and Species A=0
c_A_0 = fem.Constant(domain, 0.0)
c_I_0 = fem.Constant(domain, 2.0)
c.sub(0).interpolate(lambda x: np.full((x.shape[1],), c_A_0))
c.sub(1).interpolate(lambda x: np.full((x.shape[1],), c_I_0))
c.sub(0).collapse()
c.sub(1).collapse()
#---------------------------
# LINEAR ELASTICITY PROBLEM
#---------------------------
# Spatial Dimension
spa_dim = len(x_spatial)
# Identity Tensor
I = ufl.Identity(spa_dim)
# Deformation Gradient Tensor
F = I + ufl.grad(u)
F_inv = ufl.inv(F) # F^-1
F_inv_tra = ufl.inv(F).T # F^-T
# Jacobian Function
J = ufl.det(F)
J_n = J
# Calculate Normal Vectors along Boundary in Reference/Undeformed Configuration (might need minus to point outwards)
n0 = ufl.FacetNormal(domain) # Stationary Domain Normal Vector
piola_transform = J*F_inv_tra*n0
piola_transform_norm = ufl.inner(piola_transform, piola_transform)
n = piola_transform/ufl.sqrt(piola_transform_norm)
# Mechanics Constants
youngs_modulus_E = 0.30 # kPa
poisson_ratio_nu = 0.45
stoke_drag_gamma = 20.0 # kPa*s/um^2 for final velocity of 0.17 um/s
lambda_ = fem.Constant(domain, youngs_modulus_E*poisson_ratio_nu/((1+poisson_ratio_nu)*(1-2*poisson_ratio_nu))) # ~0.93 kPa
mu_ = fem.Constant(domain, youngs_modulus_E/(2*(1+poisson_ratio_nu))) # ~0.1 kPa
beta_ = 10.0 # s
# Polymerisation Force Constants
s0 = fem.Constant(domain, 1.0)
s1 = fem.Constant(domain, 10.0)
f_R = fem.Constant(domain, 0.3) # nN/um
# Cauchy Stress Tensor
def sigma(u_current, u_previous, del_t):
identity_ = ufl.Identity(len(u_current))
eps_current = epsilon(u_current)
eps_previous = epsilon(u_previous)
eps_dot = (eps_current - eps_previous)/del_t
stress = (lambda_ * ufl.tr(eps_current) * identity_ + 2 * mu_ * eps_current) \
+ (beta_ * (lambda_ * ufl.tr(eps_dot) * identity_ + 2 * mu_ * eps_dot))
return stress
# Small-Scale Strain
def epsilon(u):
small_strain = 0.5 * (ufl.nabla_grad(u) + ufl.transpose(ufl.nabla_grad(u)))
return small_strain
# Heaviside Polymerisation Force from Actin
def poly_force(cA):
# C_A is active species
polymerisation_force = f_R / (1 + ufl.exp(-2 * s1 * (cA - s0)))
return polymerisation_force
#-----------------------------
# WEAK FORM PROBLEM DERIVATION
#-----------------------------
term1 = ufl.inner(sigma(u, u_n, dt), epsilon(v)) * dx
term2 = poly_force(c_A) * ufl.inner(n, v) * ds
term3 = - stoke_drag_gamma * ufl.dot(u-u_n, v) * dx
Func = term1 + term2 + term3
a = ufl.lhs(Func)
L = ufl.rhs(Func)
problem = LinearProblem(a, L, bcs=[])
In implementing and running the code, the error produced is:
line 165, in <module>
a = ufl.lhs(Func)
^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/formoperators.py", line 78, in lhs
return compute_form_lhs(form)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 365, in compute_form_lhs
return compute_form_with_arity(form, 2)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 334, in compute_form_with_arity
return map_integrands(_transform, form)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 30, in map_integrands
map_integrands(function, itg, only_integral_type) for itg in form.integrals()
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
return itg.reconstruct(function(itg.integrand()))
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 329, in _transform
e, provides = pe.visit(e)
^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 109, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 260, in linear_indexed_type
part, provides = self.visit(expression)
^^^^^^^^^^^^^^^^^^^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/transformer.py", line 113, in visit
r = h(o)
^^^^
File "/Users/vozcoban/micromamba/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/formtransformations.py", line 204, in division
raise ValueError(
ValueError: Found Argument in denominator of <Division id=5119406928>, this is an invalid expression.
I’ve honed in on where the problem is occurring and has something to do with the deformation gradient tensor F and the Jacobian J (which are used in term2 towards the end of the code) in the following section of code:
# Spatial Dimension
spa_dim = len(x_spatial)
# Identity Tensor
I = ufl.Identity(spa_dim)
# Deformation Gradient Tensor
F = I + ufl.grad(u)
F_inv = ufl.inv(F) # F^-1
F_inv_tra = ufl.inv(F).T # F^-T
# Jacobian Function
J = ufl.det(F)
J_n = J
# Calculate Normal Vectors along Boundary in Reference/Undeformed Configuration (might need minus to point outwards)
n0 = ufl.FacetNormal(domain) # Stationary Domain Normal Vector
piola_transform = J*F_inv_tra*n0
piola_transform_norm = ufl.inner(piola_transform, piola_transform)
n = piola_transform/ufl.sqrt(piola_transform_norm)
If you remove both J and F_inv_tra from piola_transform the error dissappears. Placing just F_inv_tra back into the multiplication results in the division error. Additionally, if you formulate the problem as piola_transform=J*n0 another error results.
I had a few additional questions on top of the error:
- Is the piola transform introducing some sort of non-linearity?
- Is it fine to create two function spaces: a vector-valued function space for my linear elasticity problem and a mixed-element function space for my reaction-diffusion problem? I ask this as I’m not sure whether there are any conflicts when trying to pass through values from the mixed-element function space (C_A in this case) into a function which will be used in the linear elasticity weak formulation (see poly_force(cA) function).
- Should this problem be solved as a linear or nonlinear problem in FEniCSx? The report from which I’m trying to replicate this from states this problem was solved using the in built solver in FEniCS (unfortunately, I am unable to acquire this code).
Cheers,
Volkan