In the below MWE, I am compiling the outer three annuli into a single domain. So I am left out with two meshes, the mesh on the innermost disc and the single combined mesh of the outer three annuli. Next, I am applying a pressure on the innermost circle (for the pressure, I am considering a for loop to slowly increase the pressure, otherwise the code diverges) and solving for the deformation (u1) of the combined mesh by minimizing a variational form that consists of identity matrices and a pressure term. The deformation (u2), due to deformation u1, of the disc mesh is computed by minimizing another variational form and passing the displacement of the nodes on the innermost circle as a DirichletBC. This code works perfectly fine and gives me accurate results.
Next, I am applying ALE.move() on both meshes by their corresponding solutions.
Next, there is another similar loop. In this loop, the pressure doesn’t increase. Instead the entries of a diagonal matrix changes with iterations which is incorporated into the strain energy density function.
MWE:
I have considered the following geometry:
GMSH script (in the MWE, the msh file is named MWE.msh):
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 1, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 1.5, 0, 2*Pi};
//+
Circle(3) = {0, 0, 0, 2.5, 0, 2*Pi};
//+
Circle(4) = {0, 0, 0, 3.5, 0, 2*Pi};
//+
Hide "*";
//+
Show {
Curve{1}; Curve{2}; Curve{3}; Curve{4};
}
//+
Physical Curve("C1", 1) = {1};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {3};
//+
Physical Curve("C4", 4) = {4};
//+
Curve Loop(1) = {1};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {2};
//+
Curve Loop(3) = {1};
//+
Plane Surface(2) = {2, 3};
//+
Curve Loop(4) = {3};
//+
Curve Loop(5) = {2};
//+
Plane Surface(3) = {4, 5};
//+
Curve Loop(6) = {4};
//+
Curve Loop(7) = {3};
//+
Plane Surface(4) = {6, 7};
//+
Physical Surface("S5", 5) = {2};
//+
Physical Surface("S6", 6) = {3};
//+
Physical Surface("S7", 7) = {4};
//+
Physical Surface("S8", 8) = {1};
The mesh looks like
Here is the MWE:
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
#***********IMPORTING GEOMETRY FROM GMSH*****************************
msh = meshio.read("MWE.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()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
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)
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)
#***********************EXTRACTING DISC MESH AND DEFINING MEASURES******************
lumen=SubMesh(mesh,cf,8)
boundary_marker_lum = MeshFunction("size_t", lumen, lumen.topology().dim()-1, 0)
ncells2 = MeshFunction("size_t", lumen, lumen.topology().dim())
surface_marker_lum = MeshFunction("size_t", lumen, lumen.topology().dim(), 0)
vmap2 = lumen.data().array("parent_vertex_indices", 0)
cmap2 = lumen.data().array("parent_cell_indices", lumen.topology().dim())
n = 0
for c in cells(lumen):
parent_cell = Cell(mesh, cmap2[c.index()])
surface_marker_lum.array()[c.index()] = cf.array()[parent_cell.index()]
for f in facets(parent_cell):
for g in facets(c):
g_vertices = vmap2[g.entities(0)]
if set(f.entities(0)) == set(g_vertices):
boundary_marker_lum.array()[g.index()] = mf.array()[f.index()]
n=n+1
ds_lum = Measure('ds', domain=lumen, subdomain_data=boundary_marker_lum)
dx_lum = Measure('dx', domain=lumen, subdomain_data=surface_marker_lum)
#************COMPILING OUTER THREE ANNULI TO SOLVE THE HYPERELASTICITY PROBLEM***************
combined_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
cs_2 = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
combined_subdomains.array()[cf.array()==5] = 1
combined_subdomains.array()[cf.array()==6] = 1
combined_subdomains.array()[cf.array()==7] = 1
mesh_ima = SubMesh(mesh, combined_subdomains, 1)
boundary_marker_ima = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim()-1, 0)
ncells0 = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim())
surface_marker_ima = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim(), 0)
vmap0 = mesh_ima.data().array("parent_vertex_indices", 0)
cmap0 = mesh_ima.data().array("parent_cell_indices", mesh_ima.topology().dim())
n = 0
for c in cells(mesh_ima):
parent_cell = Cell(mesh, cmap0[c.index()])
surface_marker_ima.array()[c.index()] = cf.array()[parent_cell.index()]
for f in facets(parent_cell):
for g in facets(c):
g_vertices = vmap0[g.entities(0)]
if set(f.entities(0)) == set(g_vertices):
boundary_marker_ima.array()[g.index()] = mf.array()[f.index()]
n=n+1
ds_ima = Measure("ds", domain=mesh_ima, subdomain_data=boundary_marker_ima)
dx_ima= Measure("dx", domain=mesh_ima, subdomain_data=surface_marker_ima)
n = FacetNormal(mesh_ima)
#*****************SOLVING HYPERELASTICITY PROBLEM ON THE OUTER THREE ANNULI COMBINED**************
V1 = VectorFunctionSpace(mesh_ima, "Lagrange", 1) # Introducing a vector function space on intima, media and adventitia combined, i.e., ima submesh
V2 = VectorFunctionSpace(lumen, "Lagrange", 1) # Introducing a vector function space on lumen
#----------DEFINING FUNCTIONS NEEDED FOR THE VARIATIONAL FORM--------------------------
u1 = Function(V1)
u2 = Function(V2)
v1 = TestFunction(V1)
v2 = TestFunction(V2)
u_ima = Function(V1)
u_lum = Function(V2)
P = 0.0
Pmax = 0.0
#---------------------DEFINING FUNCTIONS NEEDED FOR THE VARIATIONAL FORM FOR THE DISC---------------------
dp = 0.1 # Small pressure increment
NumSteps = 10
for i in range(0,NumSteps): # For loop for pressure increment
dl = lumen.geometry().dim()
Il = Identity(dl) # Identity tensor
Fl = Il + grad(u2)
J_e_l = det(Fl)
C_lumen = Fl.T*Fl
Ic_lumen = tr(C_lumen)
alpha, mu_lumen = 0.16666,1.0
psi_lumen = (mu_lumen/2)*(Ic_lumen - 3) + (alpha)*(J_e_l-1)**2 - mu_lumen*ln(J_e_l)
Pi_lumen = (psi_lumen*dx(domain=lumen))
#-----------------------------------------------------SETTING THE KINEMATICS FOR IMA-----------------------------------------------------------------------------------------------------
d = mesh_ima.geometry().dim()
I = Identity(d) # Identity tensor
F = I + grad(u1) # Deformation gradient
g_int = as_matrix([[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])
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_e = det(F)
C_lumen = F.T*F
Ic_lumen = tr(C_lumen)
# Elasticity parameters
nu = 0.497
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 (compressible neo-Hookean model)
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)
###########################################################################################
Pi = (det(g_int))*(psi1*dx_ima(5)) + (det(g_med))*(psi2*dx_ima(6)) + (det(g_ad))*(psi3*dx_ima(7)) - (0.5)*inner(M*(P+(i*dp))*(-n), u1)*ds_ima(1)
Pi2 = psi_lumen*dx_lum(8)
# Compute first variation of Pi (directional derivative about u in the direction of v)
F1 = derivative(Pi, u1, v1)
F2 = derivative(Pi2, u2, v2)
bcs = []
# Solve variational problem
solve(F1 == 0, u1, bcs,
form_compiler_parameters=ffc_options)
u1.set_allow_extrapolation(True)
bc2 = DirichletBC(V2, u1, boundary_marker_lum, 1)
solve(F2 == 0, u2, bc2,
form_compiler_parameters=ffc_options)
u_ima.assign(u1)
u1 = interpolate(u_ima, V1)
u_lum.assign(u2)
u2 = interpolate(u_lum, V2)
print('Pressure at the end of each loop', P+(i*dp))
print('END OF PRESSURE FOR LOOP')
Pmax = P + (NumSteps-1)*dp
print('Pmax =', Pmax)
#************************MOVING MESHES BY THE SOLUTIONS**********************
ALE.move(mesh_ima,u1) # Without ALE.move, the code works fine
ALE.move(lumen,u2)
#*****************************FOR LOOP FOR MATRIX****************************
V_lum = VectorFunctionSpace(lumen, 'P', 1)
V_ima = VectorFunctionSpace(mesh_ima, 'P', 1)
ug1 = Function(V_ima)
ug2 = Function(V_lum)
vg1 = TestFunction(V_ima)
vg2 = TestFunction(V_lum)
ug_ima = Function(V_ima)
ug_lum = Function(V_lum)
ng = FacetNormal(mesh_ima)
dt = 0.05
NumStepg = 10
for j in range(0,NumStepg): # For loop for matrix increment
dlg = lumen.geometry().dim()
Ilg = Identity(dlg) # Identity tensor
Flg = Ilg + grad(ug2)
J_e_lg = det(Flg)
Cg_lumen = Flg.T*Flg
Icg_lumen = tr(Cg_lumen)
alpha, mu_lumen = 0.16666,1.0
psig_lumen = (mu_lumen/2)*(Icg_lumen - 3) + (alpha)*(J_e_lg-1)**2 - mu_lumen*ln(J_e_lg) # Accoding to Gou_Pence
#-----------------------------------------------------SETTING THE KINEMATICS FOR IMA-----------------------------------------------------------------------------------------------------
dg = mesh_ima.geometry().dim()
Ig = Identity(dg) # Identity tensor
Fg = Ig + grad(ug1) # Deformation gradient
g_i = as_matrix([[1.0+(j*dt),0.0,0.0],[0.0,1.0+(j*dt),0.0],[0.0,0.0,1.0]])
g_med = Identity(dg)
g_ad = Identity(dg)
Cg_int = (Fg*inv(g_i)).T*(Fg*inv(g_i))
Cg_med = (Fg*inv(g_med)).T*(Fg*inv(g_med))
Cg_ad = (Fg*inv(g_ad)).T*(Fg*inv(g_ad))
Mg = cofac(Fg)
# Invariants of deformation tensors
Icg_int = tr(Cg_int)
Icg_med = tr(Cg_med)
Icg_ad = tr(Cg_ad)
Jg_int = det(Fg*inv(g_int))
Jg_med = det(Fg*inv(g_med))
Jg_ad = det(Fg*inv(g_ad))
Jg_e = det(Fg)
# Stored strain energy density (compressible neo-Hookean model)
psig1 = (mu1/2)*(Icg_int - 3) + ((nu)*(mu1)/(1-2*nu))*(Jg_int-1)**2 - mu1*(ln(Jg_int)) + (eta_int/beta_int)*(exp(beta_int*((1-rho_int)*(Icg_int - 3)**2))-1)
psig2 = (mu2/2)*(Icg_med - 3) + ((nu)*(mu2)/(1-2*nu))*(Jg_med-1)**2 - mu2*(ln(Jg_med)) + (eta_med/beta_med)*(exp(beta_med*((1-rho_med)*(Icg_med - 3)**2))-1)
psig3 = (mu3/2)*(Icg_ad - 3) + ((nu)*(mu3)/(1-2*nu))*(Jg_ad-1)**2 - mu3*(ln(Jg_ad)) + (eta_ad/beta_ad)*(exp(beta_ad*((1-rho_ad)*(Icg_ad - 3)**2))-1)
Pig = det(g_i)*(psig1*dx_ima(5)) + (det(g_med))*(psig2*dx_ima(6)) + (det(g_ad))*(psig3*dx_ima(7)) - (0.5)*inner(Mg*(Pmax)*(-ng), ug1)*ds_ima(1)
Pig2 = psig_lumen*dx_lum(8)
# Compute first variation of Pi (directional derivative about u in the direction of v)
Fg1 = derivative(Pig, ug1, vg1)
Fg2 = derivative(Pig2, ug2, vg2)
bcgs_ima = []
# Solve variational problem
solve(Fg1 == 0, ug1, bcgs_ima,
form_compiler_parameters=ffc_options)
ug1.set_allow_extrapolation(True)
bcgs_lum = DirichletBC(V_lum, ug1, boundary_marker_lum, 1)
solve(Fg2 == 0, ug2, bcgs_lum,
form_compiler_parameters=ffc_options)
ug_ima.assign(ug1)
ug1 = interpolate(ug_ima, V_ima)
ug_lum.assign(ug2)
ug2 = interpolate(ug_lum, V_lum)
print('g_int(11) =', 1.0+j*dt)
print('end of a growth loop')
ALE.move(mesh_ima,ug1)
ALE.move(lumen,ug2)
This provides the following error:
File "MWE.py", line 317, in <module>
solve(Fg2 == 0, ug2, bcgs_lum,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 314, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to create mesh entity.
*** Reason: Mesh entity index -1 out of range [0, 378] for entity of dimension 2.
*** Where: This error was encountered inside MeshEntity.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
The code works fine if ALE.move() is not applied in between the loops.
I am unable to figure out the problem here. Any help would be greatly appreciated. Thanks.