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()
# %%