poroelasticmodelforacutemyocarditis's Introduction

PoroelasticModelForAcuteMyocarditis

What

Finite element simulation of immune system dynamics coupled with large-strain poromechanics. It relies on the FEniCS package (tested with v.2019.1).

Context

The immune system dynamics are as follows

And the simulations produce the following type of solutions

How to cite

Further details on the model and discretisation can be found in the following two references [1] [2]

@rticle{Barnafi21,
  url     = {http://arxiv.org/abs/2111.04206},
  doi     = {2111.04206xxx},
  author  = {Barnafi, Nicolas and G\'omez-Vargas, Bryan  and Louren\c{c}o, Wesley de Jesus and 
              Reis, Ruy Freitas and Rocha, Bernardo Martins and Lobosco, Marcelo and 
              Ruiz-Baier, Ricardo and Weber dos Santos, Rodrigo},
  title   = {Mixed methods for large-strain poroelasticity/chemotaxis models 
              simulating the formation of myocardial oedema},
  year    = {2021},
  journal = {arXiv preprint}
}

@article{Lourenco22,
  doi = {10.3389/fphys.2022.888515},
  year    = {2022},
  volume  = {13}, 
  pages ={e888515(1--14)},
  author  = {Louren\c{c}o, Wesley de Jesus and Reis, Ruy Freitas and 
              Ruiz-Baier, Ricardo and Rocha, Bernardo Martins and 
              Weber dos Santos, Rodrigo and Lobosco, Marcelo},
  title   = {A poroelastic approach for modelling myocardial oedema 
              in acute myocarditis},
  journal = {Frontiers in Physiology}
}

poroelasticmodelforacutemyocarditis's People

Contributors

Stargazers

Watchers

poroelasticmodelforacutemyocarditis's Issues

The results do not coincide with the paper.

I have made modifications to the code in ex04 to replicate the example from Fig. 4 in the Barnafi21 paper, but unfortunately, it did not succeed. Numerically, our results are significantly smaller even though we used the same parameters.

Here is my updated code:

from fenics import *

parameters["form_compiler"]["representation"] = "uflacs"
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 2

# parameter-dependent functions

normalize_vec = lambda u : u/sqrt(dot(u,u))
kappa  = lambda J,phi: kappa0*(1 + (1-phi0)**2/phi0**3*J*phi**3*(J-phi)**2)  # isotropic Kozeny-Carman

# time constants
t = Expression("t", t=0,degree=1)
dt = 1
Tfinal = 100
cont = 0

# ********* model constants ******* #

E = Constant(3.e4)
nu = Constant(0.2)
ell = Constant(0.0)

rhos = Constant(2.e-3)
rhof = Constant(1.e-3)
kappa0 = Constant(1.e-8)
alpha = Constant(1)
muf = Constant(1.e-3)
bb = Constant((0, 0))
I = Identity(2)

# neo-Hookean
mu_s = Constant(E/(2.*(1. + nu)))
lmbdas = Constant(E*nu/((1. + nu)*(1. - 2.*nu)))

# ********* Mesh and I/O ********* #
mesh = RectangleMesh(Point(0.0, 0.0), Point(8.0, 5.0), 100, 100)
bdry = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)


class Boundary_1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)
        

class Boundary_2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 8.0)


class Boundary_3(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0) 


class Boundary_4(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 5.0)


class Boundary_5(SubDomain):
    def inside(self, x, on_boundary):
       return on_boundary and near(x[1], 5.0) and (0.0 <= x[0] <= 1.0)  


boundary_1 = Boundary_1()
boundary_2 = Boundary_2()
boundary_3 = Boundary_3()
boundary_4 = Boundary_4()
boundary_5 = Boundary_5()


boundary_1.mark(bdry, 1)
boundary_2.mark(bdry, 2)
boundary_3.mark(bdry, 3)
boundary_4.mark(bdry, 4)
boundary_5.mark(bdry, 5)


fileO = File("outputs/boundary.pvd")
fileO << bdry

ds = Measure("ds", subdomain_data= bdry)  
nn = FacetNormal(mesh)

fileO_p = File("outputs/drainage1/p.pvd")
fileO_u = File("outputs/drainage1/u.pvd")
fileO_phi = File("outputs/drainage1/phi.pvd")

# ********* Finite dimensional spaces ********* #

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Bub = FiniteElement("Bubble", mesh.ufl_cell(), 4)
P1b = VectorElement(P1 + Bub)
P2vec = VectorElement("CG", mesh.ufl_cell(), 1)
Hh = FunctionSpace(mesh, MixedElement([P2vec, P1, P1]))

print("**************** Total Dofs = ", Hh.dim())

Sol = Function(Hh)
dSol = TrialFunction(Hh)
u, phi, p = split(Sol)
v, psi, q = TestFunctions(Hh)

# ******* Boundary conditions ********** #
u_D = project(Constant((0, 0)), Hh.sub(0).collapse())
bcU = DirichletBC(Hh.sub(0), u_D, bdry, 3)
bcU1 = DirichletBC(Hh.sub(0).sub(0), Constant(0), bdry, 1)
bcU2 = DirichletBC(Hh.sub(0).sub(0), Constant(0), bdry, 2)
bcP1 = DirichletBC(Hh.sub(2), 0, bdry, 4)
bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 5)
# bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 1)
# bcP3 = DirichletBC(Hh.sub(2), 0, bdry, 2)
# bcP4 = DirichletBC(Hh.sub(2), 0, bdry, 3)
bcH = [bcU, bcP1,bcP2, bcU1, bcU2]

# ******* Initial conditions ********** #
phi0 = Constant(0.33)
uold = project(Constant((0, 0)), Hh.sub(0).collapse())
phiold = project(phi0, Hh.sub(1).collapse())
pold = interpolate(Constant(0), Hh.sub(2).collapse())


# ********  Weak form ********** #

ndim = u.geometric_dimension()

F = I + grad(u)
J = det(F)
invF = inv(F)

Peff = mu_s*(F - invF.T) + lmbdas*ln(J)*invF.T
t_1 = 0.2*(lmbdas+2.0*mu_s)*sin(2*pi*t/15.0)

# Poromechanics
F1 = inner(Peff - alpha*p*J*invF.T, grad(v)) * dx \
    + psi*(J-1-phi+phi0)*dx \
    + rhof*J/dt*(phi-phiold)*q * dx \
    + rhof*dot(phi*J*invF*kappa(J, phi)/muf*invF.T*grad(p), grad(q)) * dx \
    - rhos*dot(bb, v) * dx \
    - dot(-J*t_1*invF.T*nn, v) * ds(5) \
    - rhof*J*ell*q*dx
    
    

Tang = derivative(F1, Sol, dSol)

# ********* Time loop ************* #

while (t.t < Tfinal):

    t.t += dt

    print("t=%.2f" % t.t)

    # ********* Solving ************* #
    
    solve(F1 == 0, Sol, J=Tang, bcs=bcH,
          solver_parameters={'newton_solver': {'linear_solver': 'mumps',
                                               'absolute_tolerance': 1.0e-7,
                                               'relative_tolerance': 1.0e-7}})
    uh, phih, ph = Sol.split()

    # ********* Writing ************* #

    fileO_p << ph
    fileO_u << uh
    fileO_phi << phih

    cont += 1

    # ********* Updating ************* #
    assign(uold, uh)
    assign(phiold, phih)
    assign(pold, ph)

# ************* End **************** #

Recommend Projects