### latest syntax for F.subs({0.0:gamma_n})

86
views
0
6 weeks ago by
Hi there,

I am trying an old code on nonlinear elasticity which uses the syntax F.subs({0.0:gamma_n}) . F is the deformation gradient i.e. F = I+grad(u)

In this code gamma_n is calculated through a time stepping loop so that gamma_n = 0.5*sin(2*t_n/10*float(pi)) , where t_n is incremented through the loop so that gamma_n first increases and then decreases.

Can any one from the community very kindly suggest latest syntax for F.subs({0.0:gamma_n}) ?
With this syntax I am getting error  : AttributeError: 'Form' object has no attribute 'subs'
Thanks
Community: FEniCS Project
Based on the error, I am guessing that, in the program you are looking at, the variable F has been re-defined as a variational form after being used earlier for the deformation gradient, since the deformation gradient would not be of type Form.  I'm not familiar with the subs() method, but is it possible that you are looking for the replace() function?  See, e.g., the following code snippet:

from dolfin import *
mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh,"Lagrange",1)
u = Function(V)
print(type(F))
import ufl.form
print(isinstance(F,ufl.form.Form))

# Causes a different error from the one given in the question, because F is
# of type Sum, not type Form:
#F.subs({0.0:gamma_n})

# Re-define variable F to be something of type Form, with a terminal object
# named gamma_n to substitute various things for:
gamma_n = Constant(1.0)
F = 3.0*gamma_n*dx(domain=mesh)
print(type(F))
print(isinstance(F,ufl.form.Form))

# Re-produces error from question:
#F.subs({0.0:gamma_n})

# Look at the effect of substituting different things for gamma_n in F:
print(assemble(F))
F_modified = replace(F,{gamma_n:Constant(0.0)})
print(assemble(F_modified))
F_modified = replace(F,{gamma_n:Constant(2.0)})
print(assemble(F_modified))
​
written 6 weeks ago by David Kamensky
# Library imports and settings
from dolfin import*
from fenics import*
#from mshr import *
#import ufl
#ufl.algorithms.preprocess.enable_profiling = True
from numpy import array, arange
#parameters["form_compiler"]["name"] = "sfc"

# Material parameters for Figure 7 in HolzapfelOgden2009
a    =  Constant(0.500)   #kPa
b    =  Constant(8.023)
a_f  =  Constant(16.472)  #kPa
b_f  =  Constant(16.026)
a_s  =  Constant(2.481)   #kPa
b_s  =  Constant(11.120)
a_fs =  Constant(0.356)   #kPa
b_fs =  Constant(11.436)

# Material parameters for compressibility
kappa = Constant(2.0e3)   #kPa
beta  = Constant(9.0)

# Parameters related to time-stepping
T = 10.0
dt = T/100
gamma_max = 0.5

# Parameters related to viscoelasticity
tau = 0.5
beta_inf = 0.25
xi = -dt/(2*tau)

# Strain energy functions for the passive myocardium
# Isochoric part
def psi_iso_inf(I1_bar, I4_f_bar, I4_s_bar, I8_fs_bar, I8_fn_bar):
return(a/(2*b)*exp(b*(I1_bar - 3)) \
+ a_f/(2*b_f)*(exp(b_f*(I4_f_bar - 1)**2) - 1) \
+ a_s/(2*b_s)*(exp(b_s*(I4_s_bar - 1)**2) - 1) \
+ a_fs/(2*b_fs)*(exp(b_fs*I8_fs_bar**2) - 1))

# Volumetric part
def psi_vol_inf(J):
return(kappa*(1/(beta**2)*(beta*ln(J) + 1/(J**beta) - 1)))

# Reference fibre, sheet and sheet-normal directions
f0 = Constant((1, 0, 0))
s0 = Constant((0, 1, 0))
n0 = Constant((0, 0, 1))

# Define kinematic measures in terms of the displacement
def kinematics(u):
I = Identity(3)    # Identity tensor
C = F.T*F                   # Right Cauchy-Green tensor
J = variable(det(F))        # Jacobian
C_bar = J**(-2.0/3.0)*C     # Modified right Cauchy-Green tensor

# Principle isotropic invariants
I1_bar = variable(tr(C_bar))
I2_bar = variable(0.5*(tr(C_bar)**2 - tr(C_bar*C_bar)))

# Anisotropic (quasi) invariants
I4_f_bar = variable(inner(f0, C_bar*f0))
I4_s_bar = variable(inner(s0, C_bar*s0))
I8_fs_bar = variable(inner(f0, C_bar*s0))
I8_fn_bar = variable(inner(f0, C_bar*n0))

return [I, F, C, J, C_bar, I1_bar, I2_bar, \
I4_f_bar, I4_s_bar, I8_fs_bar, I8_fn_bar]

# Define the elastic response of the material
# Isochoric part of the second Piola-Kirchhoff stress
def S_iso_inf(u):
# Define useful kinematic measures
[I, F, C, J, C_bar, I1_bar, I2_bar, \
I4_f_bar, I4_s_bar, I8_fs_bar, I8_fn_bar] = kinematics(u)

# Strain energy functions
psi_iso = psi_iso_inf(I1_bar, I4_f_bar, I4_s_bar, I8_fs_bar, I8_fn_bar)

# Define the second Piola-Kirchhoff stress in terms of the invariants
S_bar =   2*(diff(psi_iso, I1_bar) + diff(psi_iso, I2_bar))*I \
- 2*diff(psi_iso, I2_bar)*C_bar \
+ 2*diff(psi_iso, I4_f_bar)*outer(f0, f0) \
+ 2*diff(psi_iso, I4_s_bar)*outer(s0, s0) \
+ diff(psi_iso, I8_fs_bar)*(outer(f0, s0) + outer(s0, f0)) \
+ diff(psi_iso, I8_fn_bar)*(outer(f0, n0) + outer(n0, f0))
Dev_S_bar = S_bar - (1.0/3.0)*inner(S_bar, C)*inv(C)
S_iso_inf = J**(-2.0/3.0)*Dev_S_bar
return(S_iso_inf)

# Volumetric part of the second Piola-Kirchhoff stress
def S_vol_inf(u):
# Define useful kinematic measures
[I, F, C, J, C_bar, I1_bar, I2_bar, \
I4_f_bar, I4_s_bar, I8_fs_bar, I8_fn_bar] = kinematics(u)
psi_vol = psi_vol_inf(J)
S_vol_inf = J*diff(psi_vol, J)*inv(C)
return(S_vol_inf)

# Cauchy stress
def sigma(u):
[I, F, C, J, C_bar, I1_bar, I2_bar, \
I4_f_bar, I4_s_bar, I8_fs_bar, I8_fn_bar] = kinematics(u)
return(1/J*P(u)*F.T)

# Dimensions and mesh density of the domain
width = 1
depth = 1
height = 1
n = 10
#mesh = Box(0, width, 0, depth, 0, height, n*width, n*depth, n*height)
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(1.0, 1.0, 1.0), 10, 10, 10)

# Function spaces
scalar = FunctionSpace(mesh, "Lagrange", 1)
vector = VectorFunctionSpace(mesh, "Lagrange", 1)
tensor = TensorFunctionSpace(mesh, "Lagrange", 1)

# Functions
du = TrialFunction(vector)            # Incremental displacement
v  = TestFunction(vector)             # Test function
u  = Function(vector)                 # Displacement from previous iteration
S_iso_inf_p = Function(tensor)
S_vol_inf_p = Function(tensor)
Q_p = Function(tensor)

# Boundary conditions
#back_condition   = "x[0] == 0.0 && on_boundary"
#front_condition  = "x[0] == %g && on_boundary" % depth
#left_condition   = "x[1] == 0.0 && on_boundary"
#right_condition  = "x[1] == %g && on_boundary" % width
#bottom_condition = "x[2] == 0.0 && on_boundary"
#top_condition    = "x[2] == %g && on_boundary" % height

#back, front = compile_subdomains([back_condition, front_condition])
#left, right = compile_subdomains([left_condition, right_condition])
#bottom, top = compile_subdomains([bottom_condition, top_condition])
# Boundary conditions
back_condition   = "x[0] == 0.0 && on_boundary"
front_condition  = "x[0] == %g && on_boundary" % depth
left_condition   = "x[1] == 0.0 && on_boundary"
right_condition  = "x[1] == %g && on_boundary" % width
bottom_condition = "x[2] == 0.0 && on_boundary"
top_condition    = "x[2] == %g && on_boundary" % height

back = CompiledSubDomain(back_condition)
front = CompiledSubDomain(front_condition)
left = CompiledSubDomain(left_condition)
right = CompiledSubDomain(right_condition)
bottom = CompiledSubDomain(bottom_condition)
top = CompiledSubDomain(top_condition)
hold = Expression(("0.0", "0.0", "0.0"), degree=0)

# Simple shear along the fs plane
shear = Expression(("0.0", "gamma*depth", "0.0"), gamma=0.0, depth=depth,degree=0)
hold_back = DirichletBC(vector, hold, back)
shear_front = DirichletBC(vector, shear, front)
bcs = [hold_back, shear_front]

# Create files to store output
u_store = TimeSeries('visco/u')
S_vol_inf_store = TimeSeries("visco/S_vol_inf")
S_iso_inf_store = TimeSeries("visco/S_iso_inf")
Q_store = TimeSeries("visco/Q")

# And store initial values
u_store.store(u.vector(), 0.0)
S_iso_inf_store.store(S_iso_inf_p.vector(), 0.0)
S_vol_inf_store.store(S_vol_inf_p.vector(), 0.0)
Q_store.store(Q_p.vector(), 0.0)
import numpy as np
# Define the time range
times = arange(dt, T + dt, dt)
maxshear = np.zeros(1000)
# Subject the body to a known strain protocol and record the stresses
for t_n in times:

# Load previous elastic stress states
t_p = t_n - dt
S_vol_inf_store.retrieve(S_vol_inf_p.vector(), t_p)
S_iso_inf_store.retrieve(S_iso_inf_p.vector(), t_p)
Q_store.retrieve(Q_p.vector(), t_p)

# Compute current shear strain and update the boundary condition
gamma_n = gamma_max*sin(2*t_n/T*float(pi))
shear.gamma = gamma_n

# Update stress state
S_vol_inf_n = S_vol_inf(u)
S_iso_inf_n = S_iso_inf(u)
H_p = exp(xi)*(exp(xi)*Q_p - beta_inf*S_iso_inf_p)
Q_n = beta_inf*exp(xi)*S_iso_inf_n + H_p
S_n = S_vol_inf_n + S_iso_inf_n + Q_n

# Define the variational form for the problem
J = derivative(F, u, du)

# Solve the boundary value problem
solve(F == 0, u, bcs, J=J)

# Project the stress states to the tensor function space
S_iso_inf_n = project(S_iso_inf(u), tensor)
S_vol_inf_n = project(S_vol_inf(u), tensor)
Q_n = project(Q_n, tensor)
S_nn = project(S_n[0,1],scalar)
maxshear = max(S_nn.compute_vertex_values(mesh))
#np.savetxt('visco_shear.txt',maxshear)
print(maxshear)
# Store the displacement and stress state at current time
u_store.store(u.vector(), t_n)
S_iso_inf_store.store(S_iso_inf_n.vector(), t_n)
S_vol_inf_store.store(S_vol_inf_n.vector(), t_n)
Q_store.store(Q_n.vector(), t_n)

# Convert to Cauchy stress for comparison with Dokos et al.
S_n = F.subs({gamma:gamma_n})*S_n*F.T.subs({gamma:gamma_n})
written 6 weeks ago by Ovais
Have added the complete code. I know its a long script to grasp but I am looking for some explanation on the last line of t this code. I will try the suggestion that you have made. Thanks for your interest.
written 6 weeks ago by Ovais
1
Okay, Googling the comment over the last line, it appears that this line is commented-out in the original source.  My guess is that it is just a note by the author of how to post-process the results using sympy, in a separate program.  (All results from Googling 'fenics "subs("' are scripts where subs() is coming from sympy, not an old version of UFL.  However, the above script does not even import sympy.)
written 6 weeks ago by David Kamensky
Hi David, thanks for your interest. I contacted the author and he has recommended me to look for latest FEniCS examples, if there are other cases of substituting a term in a form. He is himself not in touch with fenics and its latest syntax.

Can you suggest something or may be some lead which I may follow to reach my goal .. which is to substitute gamma in the form F.
Can I define (manually) another form and may  be use the other one to compute S_n ?

Thanks.
written 5 weeks ago by Ovais