Hyperelastic deformation with applied pressue on internal surface

Hi Everyone, I posted this on slack already but here may be a better place. I am simulating an applied pressure onto a boundary of a hyperelastic material, which works well enough. That boundary is a hole in a rectangle like the fluid dynamics tutorial structure from Jsdokken.com . Now if I want to also simulate the deformation within the hole by adding a mesh there (using gmsh fragment instead of cut), the simulation converges but there is no applied force and just a stationary mesh. I am not sure if I am missing anything. This is the fragmentation code:

# %%

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import matplotlib.pyplot as plt
import pyvista
from dolfinx import nls
import numpy as np
import ufl
import matplotlib.pyplot as plt
from dolfinx.plot import vtk_mesh
import meshio

from mpi4py import MPI
from dolfinx import fem, mesh, plot
import gmsh
import json
import sys
import pandas as pd
from dolfinx.io.gmshio import model_to_mesh

from dolfinx.io import XDMFFile
from dolfinx.plot import vtk_mesh

from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)


# %%
#============================== making geometry ==============================
gmsh.initialize()
gmsh.clear()
gmsh.model.add("lung")


L = 2.2
H = 0.41
c_x = c_y = 0.205
r = 0.15
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

    obstacle_mass = gmsh.model.occ.getMass(2,obstacle)

if mesh_comm.rank == model_rank:
    fluid = gmsh.model.occ.fragment([(gdim, rectangle)], [(gdim, obstacle)],removeObject=True,removeTool=True)
    gmsh.model.occ.synchronize()

# %%
# removing dups and healing
gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()
gmsh.model.occ.healShapes()
gmsh.model.occ.synchronize()
gmsh.write("trial1.brep")

# %%
# ====================== physical groups ======================
fluid_marker = 1
obstacle_srf_marker = 6
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=gdim)
    assert (len(volumes) == 2)
    for vol in volumes:
        if np.isclose(gmsh.model.occ.getMass(2,vol[1]),obstacle_mass):
            gmsh.model.addPhysicalGroup(vol[0], [vol[1]], obstacle_srf_marker)
            gmsh.model.setPhysicalName(vol[0], obstacle_srf_marker, "Obstacle_srf")            
        else:
            gmsh.model.addPhysicalGroup(vol[0], [vol[1]], fluid_marker)
            gmsh.model.setPhysicalName(vol[0], fluid_marker, "Fluid")


inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = [gmsh.model.getBoundary([x], oriented=False) for x in volumes]
    boundaries = list(set([x for xs in boundaries for x in xs]))
    #boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H / 2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H / 2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

# %%
# looking at output

gmsh.model.occ.synchronize()
gmsh.write("trial2.brep")

# %%
# ====================== Generating Mesh ======================
# Create distance field from obstacle.
# Add threshold of mesh sizes based on the distance field
# LcMax -                  /--------
#                      /
# LcMin -o---------/
#        |         |       |
#       Point    DistMin DistMax
res_min = r / 3
if mesh_comm.rank == model_rank:
    distance_field = gmsh.model.mesh.field.add("Distance")
    gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
    threshold_field = gmsh.model.mesh.field.add("Threshold")
    gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
    gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
    gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
    gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
    gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
    min_field = gmsh.model.mesh.field.add("Min")
    gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
    gmsh.model.mesh.field.setAsBackgroundMesh(min_field)


# %%

if mesh_comm.rank == model_rank:
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")


# %%
# finalizing mesh
domain, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
ft.name = "Facet_markers"
gmsh.finalize()

V = fem.VectorFunctionSpace(domain, ("Lagrange", 2))

# %%
# plotting
pyvista.start_xvfb()
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(domain, domain.topology.dim))
num_local_cells = domain.topology.index_map(domain.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices < num_local_cells]
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
cell_tag_fig = plotter.screenshot("cell_tags.png")
if not pyvista.OFF_SCREEN:
    plotter.show()


# %%
# =========================== simulation ===========================

fdim = domain.topology.dim - 1
lung_facets = ft.indices[ft.values == obstacle_marker] # obstacle marker


marked_facets = np.hstack([lung_facets])
marked_values = np.hstack([np.full_like(lung_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

# %%
# adding bc
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

rib_dofs = fem.locate_dofs_topological(V, fdim, ft.find(outlet_marker))
bcs = [fem.dirichletbc(u_bc, rib_dofs, V)]


# %%

B = fem.Constant(domain, default_scalar_type((0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

# %%
# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))


# %%
#E: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2814213/#:~:text=They%20found%20that%20the%20transversal,offer%20better%20resistance%20to%20deformation.
#nu: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4644690/#:~:text=The%20Poisson%27s%20ratio%20(v)%20of,et%20al.%2C%202012).
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2693244/
# https://clinmedjournals.org/articles/jor/jor-4-048-table1.html
# Elasticity parameters
E = default_scalar_type(24.7e3) # default = 1.0e4, muscle = 2.47e4
nu = default_scalar_type(0.47) # default = 0.3, muscle = 0.47
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

# %%
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata) # external traction
#ds = ufl.Measure('ds', domain=domain, subdomain_data=lung_facet_tag, metadata=metadata) # internal pressure
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
#F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(20)


##adding pressure
n = ufl.FacetNormal(domain)
P_intr = fem.Constant(domain, default_scalar_type(-0))  # Internal pressure value, adjust as needed
T_intr = P_intr * n  # Internal traction due to pressure
F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T_intr) * ds(2)

# %%
problem = NonlinearProblem(F, u, bcs)

# %%
solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

# %%
pyvista.start_xvfb()
plotter = pyvista.Plotter()
plotter.open_gif("pipe_obstacle.gif", fps=3)
plotter.show_axes()

#2D camera position
camera_position = (1.1, 0.21, 4.75)  # This positions the camera above the XY plane
camera_focal_point = (1.1, 0.21, 0)  # This focuses the camera on the center of your mesh on the XY plane
camera_view_up = (0, 1, 0)  # This sets the 'up' direction of the camera to be aligned with the Y-axis

plotter.camera_position = [camera_position, camera_focal_point, camera_view_up]

topology, cells, geometry = plot.vtk_mesh(u.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)

#--> Assume 'ct' has correct cell tags from FEniCSx
function_grid.cell_data['PhysicalGroup'] = ct.values

values = np.zeros((geometry.shape[0], 3))
values[:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
function_grid["u"] = values
function_grid.set_active_vectors("u")

# Warp mesh by deformation
warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])

# Compute magnitude of displacement to visualize in GIF
Vs = fem.FunctionSpace(domain, ("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array

# %%
# Simulation time parameters

#def pressure(t, amplitude, frequency, offset = 0):
#    """Calculate sinusoidal pressure variation."""
#    return amplitude * np.cos(2 * np.pi * frequency * t) + offset

def pressure(t, amplitude, frequency, offset = 0):
    """Calculate sinusoidal pressure variation."""
    return amplitude * np.sin(2 * np.pi * frequency * t) + offset

# Sinusoidal pressure parameters
amp = - 1000  # Amplitude of pressure (Pa) variation
freq = 18/60  # Frequency of breathing cycles 
#offset = (-1 * amp) + amp*0.2 #assuming patient is inhaled, hence offset sine wave with small 10% extra inhale
offset = 0

start_time = 0.0
end_time = 1/freq
num_steps = 12
dt = (end_time - start_time) / num_steps
time = start_time

# %%
log.set_log_level(log.LogLevel.INFO)


#tval0 = -1.5
# = 10
for n in range(num_steps):
    

    current_pressure = pressure(time,amp,freq,offset)
    P_intr.value = current_pressure

    num_its, converged = solver.solve(u)

    assert (converged)
    u.x.scatter_forward()
    print(f"step {n}, time {time}, Number of iterations {num_its}, Load {P_intr.value}")
    function_grid["u"][:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
    magnitude.interpolate(us)
    warped.set_active_scalars("mag")
    warped_n = function_grid.warp_by_vector(factor=1)
    plotter.update_coordinates(warped_n.points.copy(), render=False)
    plotter.update_scalar_bar_range([0, 10])
    plotter.update_scalars(magnitude.x.array)
    plotter.write_frame()


    #plt.scatter(warped_n.points.copy()[:,0]/SI_conv,warped_n.points.copy()[:,1]/SI_conv)
    #plt.show()


    # Update time
    time += dt

plotter.close()
# %%

ds is an integration measure over exterior facets, thus those that are connected to only one cell. By adding a new domain in the circle this becomes and interior facet.

I would follow the notes in:

on how to do a one-sided interior facet integral.

This is great! Is there an equivalent “dolfinx.mesh.GhostMode.shared_facet” for gmsh?

Model to mesh takes in a partitioner:
dolfinx.io.gmshio — DOLFINx 0.7.3 documentation which by default is GhostMode.none. But you can create a shared facet partitioner, as described in
Mesh creation in serial and parallel — FEniCSx Documentation

1 Like

I am still working trying to work it out but am currently getting “Discontinuous type Jacobian must be restricted.” I have changed ds → dS and specified v(‘+’). Anyone know whatI need to change?

# %%

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import matplotlib.pyplot as plt
import pyvista
from dolfinx import nls
import numpy as np
import ufl
import matplotlib.pyplot as plt
from dolfinx.plot import vtk_mesh
import meshio

from mpi4py import MPI
from dolfinx import fem, mesh, plot
import gmsh
import json
import sys
import pandas as pd
from dolfinx.io.gmshio import model_to_mesh

from dolfinx.io import XDMFFile
from dolfinx.plot import vtk_mesh

from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)


# %%
#============================== making geometry ==============================
gmsh.initialize()
gmsh.clear()
gmsh.model.add("lung")


L = 2.2
H = 0.41
c_x = c_y = 0.205
r = 0.15
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

    obstacle_mass = gmsh.model.occ.getMass(2,obstacle)

if mesh_comm.rank == model_rank:
    fluid = gmsh.model.occ.fragment([(gdim, rectangle)], [(gdim, obstacle)],removeObject=True,removeTool=True)
    gmsh.model.occ.synchronize()

# %%
# removing dups and healing
gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()
gmsh.model.occ.healShapes()
gmsh.model.occ.synchronize()
gmsh.write("trial1.brep")

# %%
# ====================== physical groups ======================
fluid_marker = 1
obstacle_srf_marker = 6
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=gdim)
    assert (len(volumes) == 2)
    for vol in volumes:
        if np.isclose(gmsh.model.occ.getMass(2,vol[1]),obstacle_mass):
            gmsh.model.addPhysicalGroup(vol[0], [vol[1]], obstacle_srf_marker)
            gmsh.model.setPhysicalName(vol[0], obstacle_srf_marker, "Obstacle_srf")            
        else:
            gmsh.model.addPhysicalGroup(vol[0], [vol[1]], fluid_marker)
            gmsh.model.setPhysicalName(vol[0], fluid_marker, "Fluid")


inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = [gmsh.model.getBoundary([x], oriented=False) for x in volumes]
    boundaries = list(set([x for xs in boundaries for x in xs]))
    #boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H / 2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H / 2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

# %%
# looking at output

gmsh.model.occ.synchronize()
gmsh.write("trial2.brep")

# %%
# ====================== Generating Mesh ======================
# Create distance field from obstacle.
# Add threshold of mesh sizes based on the distance field
# LcMax -                  /--------
#                      /
# LcMin -o---------/
#        |         |       |
#       Point    DistMin DistMax
res_min = r / 3
if mesh_comm.rank == model_rank:
    distance_field = gmsh.model.mesh.field.add("Distance")
    gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
    threshold_field = gmsh.model.mesh.field.add("Threshold")
    gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
    gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
    gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
    gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
    gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
    min_field = gmsh.model.mesh.field.add("Min")
    gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
    gmsh.model.mesh.field.setAsBackgroundMesh(min_field)


# %%

if mesh_comm.rank == model_rank:
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")


# %%
# finalizing mesh
domain, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
ft.name = "Facet_markers"
gmsh.finalize()

V = fem.VectorFunctionSpace(domain, ("Lagrange", 2))

# %%
# plotting
pyvista.start_xvfb()
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(domain, domain.topology.dim))
num_local_cells = domain.topology.index_map(domain.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices < num_local_cells]
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
cell_tag_fig = plotter.screenshot("cell_tags.png")
if not pyvista.OFF_SCREEN:
    plotter.show()


# %%
# =========================== simulation ===========================

fdim = domain.topology.dim - 1
lung_facets = ft.indices[ft.values == obstacle_marker] # obstacle marker


marked_facets = np.hstack([lung_facets])
marked_values = np.hstack([np.full_like(lung_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

# %%
# adding bc
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

rib_dofs = fem.locate_dofs_topological(V, fdim, ft.find(outlet_marker))
bcs = [fem.dirichletbc(u_bc, rib_dofs, V)]


# %%

B = fem.Constant(domain, default_scalar_type((0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

# %%
# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))


# %%
#E: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2814213/#:~:text=They%20found%20that%20the%20transversal,offer%20better%20resistance%20to%20deformation.
#nu: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4644690/#:~:text=The%20Poisson%27s%20ratio%20(v)%20of,et%20al.%2C%202012).
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2693244/
# https://clinmedjournals.org/articles/jor/jor-4-048-table1.html
# Elasticity parameters
E = default_scalar_type(24.7e3) # default = 1.0e4, muscle = 2.47e4
nu = default_scalar_type(0.47) # default = 0.3, muscle = 0.47
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

# %%
metadata = {"quadrature_degree": 4}
dS = ufl.Measure('dS', domain=domain, subdomain_data=facet_tag, metadata=metadata,subdomain_id=5) # external traction
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata) # external traction
#ds = ufl.Measure('ds', domain=domain, subdomain_data=lung_facet_tag, metadata=metadata) # internal pressure
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
#F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(20)


##adding pressure
n = ufl.FacetNormal(domain)
P_intr = fem.Constant(domain, default_scalar_type(-0))  # Internal pressure value, adjust as needed
T_intr = P_intr * n  # Internal traction due to pressure
F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v('+'), T_intr) * dS

# %%
problem = NonlinearProblem(F, u, bcs)

# %%
solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

# %%
pyvista.start_xvfb()
plotter = pyvista.Plotter()
plotter.open_gif("pipe_obstacle.gif", fps=3)
plotter.show_axes()

#2D camera position
camera_position = (1.1, 0.21, 4.75)  # This positions the camera above the XY plane
camera_focal_point = (1.1, 0.21, 0)  # This focuses the camera on the center of your mesh on the XY plane
camera_view_up = (0, 1, 0)  # This sets the 'up' direction of the camera to be aligned with the Y-axis

plotter.camera_position = [camera_position, camera_focal_point, camera_view_up]

topology, cells, geometry = plot.vtk_mesh(u.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)

#--> Assume 'ct' has correct cell tags from FEniCSx
function_grid.cell_data['PhysicalGroup'] = ct.values

values = np.zeros((geometry.shape[0], 3))
values[:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
function_grid["u"] = values
function_grid.set_active_vectors("u")

# Warp mesh by deformation
warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])

# Compute magnitude of displacement to visualize in GIF
Vs = fem.FunctionSpace(domain, ("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array

# %%
# Simulation time parameters

#def pressure(t, amplitude, frequency, offset = 0):
#    """Calculate sinusoidal pressure variation."""
#    return amplitude * np.cos(2 * np.pi * frequency * t) + offset

def pressure(t, amplitude, frequency, offset = 0):
    """Calculate sinusoidal pressure variation."""
    return amplitude * np.sin(2 * np.pi * frequency * t) + offset

# Sinusoidal pressure parameters
amp = - 1000  # Amplitude of pressure (Pa) variation
freq = 18/60  # Frequency of breathing cycles 
#offset = (-1 * amp) + amp*0.2 #assuming patient is inhaled, hence offset sine wave with small 10% extra inhale
offset = 0

start_time = 0.0
end_time = 1/freq
num_steps = 12
dt = (end_time - start_time) / num_steps
time = start_time

# %%
log.set_log_level(log.LogLevel.INFO)


#tval0 = -1.5
# = 10
for n in range(num_steps):
    

    current_pressure = pressure(time,amp,freq,offset)
    P_intr.value = current_pressure

    num_its, converged = solver.solve(u)

    assert (converged)
    u.x.scatter_forward()
    print(f"step {n}, time {time}, Number of iterations {num_its}, Load {P_intr.value}")
    function_grid["u"][:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
    magnitude.interpolate(us)
    warped.set_active_scalars("mag")
    warped_n = function_grid.warp_by_vector(factor=1)
    plotter.update_coordinates(warped_n.points.copy(), render=False)
    plotter.update_scalar_bar_range([0, 10])
    plotter.update_scalars(magnitude.x.array)
    plotter.write_frame()


    #plt.scatter(warped_n.points.copy()[:,0]/SI_conv,warped_n.points.copy()[:,1]/SI_conv)
    #plt.show()


    # Update time
    time += dt

plotter.close()
# %%

You need to restrict n to the correct side.
There are two ways of doing this, either as shown in:

or by creating a DG-0 marker function, as done in TEAM30/utils.py at v0.7.0 · Wells-Group/TEAM30 · GitHub
where you create a DG-0 space with 0 for all cells on the “wrong side” of the interface, and 1 on the “right side” and assemble the sum of the (“+”) and (“-”) restriction of T_intr.

1 Like

I got it to work, putting code snippets here for anyone else in the community trying to get it to work. Dokken’s proposed solution works great, where you create a DG-0 space. In particular form the TEAM30 link. I did the following in more detail:

  1. Trial with TEAM30 mesh “generate_team30_meshes.py”; The mesh has an external square boundary ( facet tag = 1), internal circular boundary (facet tag = 2) and inner and outer cells with respect to the internal circular boundary already labeled with cell tags 3 (inner) & 2 (outer).
  2. Apply traction to internal boundary by identifying internal boundary facets:
obstacle_marker = internal_boundary_marker = 2

fdim = domain.topology.dim - 1
lung_facets = ft.indices[ft.values == obstacle_marker]


marked_facets = np.hstack([lung_facets])
marked_values = np.hstack([np.full_like(lung_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
  1. Define internal integral measure dS:
metadata = {"quadrature_degree": 4}
dS = ufl.Measure('dS', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dS_inner = dS(2) # 2 is for the inner boundary facet tag
  1. Define DG-0 space:
V_c = fem.functionspace(domain, ("DG", 0))
  1. Restrict to inner or outer side using the labeled cells (this is from Team30 utils.py):
gap_markers = (2,3) # 2,3 are the inner and outer cell tags with respect to inner boundary
_restriction = fem.Function(V_c)
_restriction.interpolate(
    lambda x: np.ones(x.shape[1], dtype=default_scalar_type), cells=ct.find(gap_markers[1])
)
_restriction.interpolate(
    lambda x: np.zeros(x.shape[1], dtype=default_scalar_type), cells=ct.find(gap_markers[0])
)
_restriction.x.scatter_forward()
  1. Apply traction inner surface, here normal to surface
##adding pressure
n = ufl.FacetNormal(domain)
P_intr = fem.Constant(domain, default_scalar_type(+1500))  # Internal pressure value, adjust as needed
#T_intr = P_intr * n  # Internal traction due to pressure

#traction_integral = fem.form(((g('+') * ufl.dot(T, n('+')) + g('-') * ufl.dot(T, n('-'))) * ufl.dS))

traction_vector_plus = P_intr * n('+')
traction_vector_minus = P_intr * n('-')

# Dot product of traction vector with test function to form a scalar
traction_integral = ((_restriction('+') * ufl.inner(traction_vector_plus, v('+')) + 
                              _restriction('-') * ufl.inner(traction_vector_minus, v('-'))) * dS_air)
  1. Note I applied a regular fixed boundary condition on the outer square boundary, if you do that the deformation will be heavily restricted and it is hard to see the internal deformation but it is there.

Hope this will help others

1 Like