Error while updating a matrix using a for loop

1

I created a deformed 3-layer annulus and tried to update a matrix, which is involved in the dynamics of the problem, through a “for loop”. The problem is set up in 3-D. But at the same time, while experimenting with an MWE, I didn’t encounter any error while updating the matrix through a “for loop”. Below is the part of my code that produces the error (mentioned after the code) followed by the MWE that runs perfectly without producing an error.

My code:

i = 0.0
for j in range(0,4):

  bcs=[]

  # Kinematics
  d = mesh.geometry().dim()
  print(d) 
  I = Identity(d)             
  F = I + grad(u)             
  g_int = as_matrix([[1.0+i,0.0,0.0],[0.0,1.0+i,0.0],[0.0,0.0,1.0]]) 
  g_med = Identity(d)  # Growth tensor for the Media
  g_ad = Identity(d)  # Growth tensor for the Adventitia
  C_int = (F*inv(g_int)).T*(F*inv(g_int))                   
  C_med = (F*inv(g_med)).T*(F*inv(g_med))                  
  C_ad = (F*inv(g_ad)).T*(F*inv(g_ad))                   
  M = cofac(F)
  # Invariants of deformation tensors
  Ic_int = tr(C_int)     # Ic for the Intima
  Ic_med = tr(C_med)     # Ic for the Media
  Ic_ad = tr(C_ad)     # Ic for the Adventitia
  J_int  = det(F*inv(g_int))
  J_med = det(F*inv(g_med))
  J_ad = det(F*inv(g_ad))
  

  # Elasticity parameters 
  nu = 0.3   
  mu1, mu2, mu3 = 27.9, 1.27, 7.56 
  eta_int, eta_med, eta_ad = 263.66, 21.6, 38.57
  beta_int, beta_med, beta_ad = 170.88, 8.21, 85.03
  rho_int, rho_med, rho_ad = 0.51, 0.25, 0.55

# Stored strain energy density 
  psi1 = (mu1/2)*(Ic_int - 3) + ((nu)*(mu1)/(1-2*nu))*(J_int-1)**2 - mu1*(ln(J_int)) + (eta_int/beta_int)*(exp(beta_int*((1-rho_int)*(Ic_int - 3)**2))-1)
  psi2 = (mu2/2)*(Ic_med - 3) + ((nu)*(mu2)/(1-2*nu))*(J_med-1)**2 - mu2*(ln(J_med)) + (eta_med/beta_med)*(exp(beta_med*((1-rho_med)*(Ic_med - 3)**2))-1)
  psi3 = (mu3/2)*(Ic_ad - 3) + ((nu)*(mu3)/(1-2*nu))*(J_ad-1)**2 - mu3*(ln(J_ad)) + (eta_ad/beta_ad)*(exp(beta_ad*((1-rho_ad)*(Ic_ad - 3)**2))-1)
  # Total potential energy
  Pi = (det(g_int))*(psi1*dx_custom(5)) + (det(g_med))*(psi2*dx_custom(6)) + (det(g_ad))* 
  (psi3*dx_custom(7)) - (0.5)*dot(M*P*(-n), u)*ds_custom(1)

# Compute first variation of Pi (directional derivative about u in the direction of v)
  F1 = derivative(Pi, u, v)

# Compute Jacobian of F
  Ju = derivative(F1, u, du)

# Solve variational problem
  solve(F1 == 0, u, bcs,
    form_compiler_parameters=ffc_options)

  i +=0.01

This produces the following error after the first iteration (but the error doesn’t show up when I manually update the matrix without a “for loop”):

Pi = (det(g_int))*(psi1*dx_custom(5)) + (det(g_med))*(psi2*dx_custom(6)) + (det(g_ad))*(psi3*dx_custom(7)) - (0.5)*dot(M*P*(-n), u)*ds_custom(1)
  File "/usr/lib/python3/dist-packages/ufl/measure.py", line 406, in __rmul__
    error("Can only integrate scalar expressions. The integrand is a "
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (3,) and free indices with labels ().

Now here is an MWE that is similar and runs perfectly with a similar “for loop”:

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C3 = Circle(Point(0, 0),3)
C4 = Circle(Point(0, 0),4)
C_1 = C2 - C1
C_2 = C3 - C2
C_3 = C4 - C3


class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class b2(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class b3(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 9, eps=1e-3) and on_boundary
class b4(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 16, eps=1e-3) and on_boundary
class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4
class s2(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 4 and x[0]**2 + x[1]**2 < 9
class s3(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 9 and x[0]**2 + x[1]**2 < 16

mesh = generate_mesh(C4-C1, 5)


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
b1().mark(boundary_markers, 1)
b2().mark(boundary_markers, 2)
b3().mark(boundary_markers, 3)
b4().mark(boundary_markers, 4)
s1().mark(surface_markers, 5)
s2().mark(surface_markers, 6)
s3().mark(surface_markers, 7)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)


V = VectorFunctionSpace(mesh, "Lagrange", 1)
n=FacetNormal(mesh)

du = TrialFunction(V)
v  = TestFunction(V)             
u  = Function(V)                 
P = Constant(5.0)

i = 0.0
for k in range(0,4):

  bcs=[]

  # Kinematics
  d = mesh.geometry().dim()
  print(d) 
  I = Identity(d)             
  F = I + grad(u)             
  g_int = as_matrix([[1.0+i,0.0],[0.0,1.0+i]])
  J_g_i = det(g_int)
  g_med = Identity(d)  
  g_ad = Identity(d)  
  C_int = (F*inv(g_int)).T*(F*inv(g_int))                   
  C_med = (F*inv(g_med)).T*(F*inv(g_med))                   
  C_ad = (F*inv(g_ad)).T*(F*inv(g_ad))                   
  M = cofac(F)
  
  # Invariants of deformation tensors
  Ic_int = tr(C_int)     
  Ic_med = tr(C_med)     
  Ic_ad = tr(C_ad)     
  J_int  = det(F*inv(g_int))
  J_med = det(F*inv(g_med))
  J_ad = det(F*inv(g_ad))
  #J_i = det(F)

  # Elasticity parameters 
  nu = 0.3
  mu1, mu2, mu3 = 27.9, 1.27, 7.56
  eta_int, eta_med, eta_ad = 263.66, 21.6, 38.57
  beta_int, beta_med, beta_ad = 170.88, 8.21, 85.03
  rho_int, rho_med, rho_ad = 0.51, 0.25, 0.55

# Stored strain energy density
  psi1 = (mu1/2)*(Ic_int - 3) + ((nu)*(mu1)/(1-2*nu))*(J_int-1)**2 - mu1*(ln(J_int)) + (eta_int/beta_int)*(exp(beta_int*((1-rho_int)*(Ic_int - 3)**2))-1)
  psi2 = (mu2/2)*(Ic_med - 3) + ((nu)*(mu2)/(1-2*nu))*(J_med-1)**2 - mu2*(ln(J_med)) + (eta_med/beta_med)*(exp(beta_med*((1-rho_med)*(Ic_med - 3)**2))-1)
  psi3 = (mu3/2)*(Ic_ad - 3) + ((nu)*(mu3)/(1-2*nu))*(J_ad-1)**2 - mu3*(ln(J_ad)) + (eta_ad/beta_ad)*(exp(beta_ad*((1-rho_ad)*(Ic_ad - 3)**2))-1)
  Coords = Expression(('x[0]','x[1]'),degree=0)
  coords = project(Coords, V)
  phi = coords+u
  Pi = (det(g_int))*(psi1*dx(5)) + (det(g_med))*(psi2*dx(6)) + (det(g_ad))*(psi3*dx(7)) - (0.5)*dot(M*P*(-n), phi)*ds(1)

# Compute first variation of Pi (directional derivative about u in the direction of v)
  F1 = derivative(Pi, u, v)

# Compute Jacobian of F
  Ju = derivative(F1, u, du)

# Solve variational problem
  solve(F1 == 0, u, bcs,
    form_compiler_parameters=ffc_options)

  i +=0.01

I am confused as to what is causing the error in my code while the MWE works perfectly fine!

2

Please note that it is quite hard to help you, as you for the code that is failing have not supplied a mesh, or a fully working code (i.e. variables such as dx_custom is not defined, and you are missing ffc_options, and several imports).
Additionally, your usage of brackets is highly confusing. Try to reduce the number of them, i.e.

can be written as
det(g_int)*psi1*dx_custom(5).

Other than that, I would suggest you remove a term from Pi one by one, until you are left with one term causing the issue, as it is then alot easier to debug (and you might be able to further reduce the lenght of the failing code by removing unused variables).