Hey All,
I want to conduct a simple deformation over a 2D square with subdomains defined by constitutive equations. An example geometry is shown here:
The constitutive equations are a Guccione for the straight thin section, and simple Mooney-Rivlin for the remainder.
I have these regions defined via GMSH and have been following the example setup in the relevant Subdomain Tutorial however because my material property setup is quite complex it doesn’t seem like I can apply the same steps of applying a constant:
Tutorial:
kappa = Function(Q)
bottom_cells = ct.find(bottom_marker)
kappa.x.array[bottom_cells] = np.full_like(bottom_cells, 1, dtype=default_scalar_type)
My material setup:
# += Identify the subdomains of interest
str_myo = ct.find(MYO_STRAIGHT)
cytosol = ct.find(CYTOSOL)
# += Test and Trial Parameters
u, p = ufl.split(w)
v, q = ufl.TestFunctions(W)
# += Curvilinear Coordinates
x = Function(V)
x.interpolate(lambda x: (x[0], x[1]))
# += Tensor Indices
i, j, k, a, b = ufl.indices(5)
# += Curvilinear Mapping
Push = ufl.as_matrix([
[ufl.cos(ROT), -ufl.sin(ROT)],
[ufl.sin(ROT), ufl.cos(ROT)]
])
x_nu = ufl.inv(Push) * x
u_nu = ufl.inv(Push) * u
nu = ufl.inv(Push) * (x + u_nu)
# += Metric Tensors
Z_un = ufl.grad(x_nu).T * ufl.grad(x_nu)
Z_co = ufl.grad(nu).T * ufl.grad(nu)
Z_ct = ufl.inv(Z_co)
# += Christoffel Symbol | Γ^{i}_{j, a}
gamma = ufl.as_tensor((
0.5 * Z_ct[k, a] * (
ufl.grad(Z_co)[a, i, j] + ufl.grad(Z_co)[a, j, i] - ufl.grad(Z_co)[i, j, a]
)
), [k, i, j])
# += Covariant Derivative
covDev = ufl.grad(v) - ufl.as_tensor(v[k]*gamma[k, i, j], [i, j])
# += Kinematics
I = ufl.variable(ufl.Identity(MESH_DIM))
F = ufl.as_tensor(I[i, j] + ufl.grad(u)[i, j], [i, j]) * Push
C = ufl.variable(ufl.as_tensor(F[k, i]*F[k, j], [i, j]))
E = ufl.variable(ufl.as_tensor((0.5*(Z_co[i,j] - Z_un[i,j])), [i, j]))
J_myo = ufl.variable(ufl.det(F))
# += Material Setup | Guccione
Q = (
GCC_CONS[1] * E[0,0]**2 +
GCC_CONS[2] * (E[1,1]**2) +
GCC_CONS[3] * (2*E[0,1]*E[1,0])
)
piola = GCC_CONS[0]/4 * ufl.exp(Q) * ufl.as_matrix([
[4*GCC_CONS[1]*E[0,0], 2*GCC_CONS[3]*(E[1,0] + E[0,1])],
[2*GCC_CONS[3]*(E[0,1] + E[1,0]), 4*GCC_CONS[2]*E[1,1]],
]) - p * Z_un
myo_term = ufl.as_tensor(piola[a, b] * F[j, b] * covDev[j, a])
Which is repeated but for a Mooney-Rivlin.
I = ufl.variable(ufl.Identity(MESH_DIM))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
IIc = ufl.variable((Ic**2 - ufl.inner(C,C))/2)
J_cyt = ufl.variable(ufl.det(F))
c1 = 2
c2 = 6
# += Mooney-Rivlin Strain Energy Density Function
psi = c1 * (Ic - 3) + c2 *(IIc - 3)
gamma1 = ufl.diff(psi, Ic) + Ic * ufl.diff(psi, IIc)
gamma2 = -ufl.diff(psi, IIc)
piola = 2 * F * (gamma1*I + gamma2*C) + p * J_cyt * ufl.inv(F).T
cyt_term = ufl.inner(ufl.grad(v), piola)
My question is:
Can I define my variational terms “myo_term” and “cyt_term” in a similar way to the tutorial (i.e. kappa.x.array[bottom_cells]) as myo_term.x.array[bottom_cells] does not work or myoterm[bottom_cells] like:
metadata = {"quadrature_degree": quad_order}
ds = ufl.Measure('ds', domain=domain, subdomain_data=ft, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
term.x.array[str_myo] = myo_term
term.x.array[cytosol] = cyt_term
R = term * dx + q * (J - 1) * dx
Or
Is there a way to achieve this with measures like below:
dx_myo = ufl.Measure("dx", domain=domain, subdomain_data=MYO_STRAIGHT, metadata=metadata)
dx_cyt = ufl.Measure("dx", domain=domain, subdomain_data=CYTOSOL, metadata=metadata)
R = myo_term * dx_myo + q * (J_myo - 1) * dx_myo + cyt_term * dx_cyt + q * (J_cyt - 1) * dx_cyt
Neither approaches are working thus far.
I apologise for any obvious ignorance or issues with my question, happy to amend.
Thank you 