Hyperelastic deformation with applied pressue on internal surface

Chris_S · May 9, 2024, 4:45am

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

dokken · May 9, 2024, 4:50am

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:

github.com/jorgensd/dolfinx-tutorial

Add one-sided integration example as example of custom integration

opened 08:04PM - 28 Nov 23 UTC

jorgensd

```python # Computing a onesided integral over an interior facet with consisten…t orientation # Enabled by: https://github.com/FEniCS/dolfinx/pull/2269 # Copyright 2023 Jørgen S. Dokken # SPDX: MIT import dolfinx import numpy as np import ufl from mpi4py import MPI mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet) # Create connectivties required for defining integration entities tdim = mesh.topology.dim fdim = tdim - 1 mesh.topology.create_entities(fdim) mesh.topology.create_connectivity(fdim, tdim) mesh.topology.create_connectivity(tdim, fdim) # Get number of cells on process cell_map = mesh.topology.index_map(tdim) num_cells = cell_map.size_local + cell_map.num_ghosts # Create markers for each size of the interface cell_values = np.ones(num_cells, dtype=np.int32) cell_values[dolfinx.mesh.locate_entities( mesh, tdim, lambda x: x[0] <= 0.5+1e-13)] = 2 ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange( num_cells, dtype=np.int32), cell_values) facet_map = mesh.topology.index_map(fdim) num_facets = facet_map.size_local + facet_map.num_ghosts # Create facet markers facet_values = np.ones(num_facets, dtype=np.int32) facet_values[dolfinx.mesh.locate_entities( mesh, fdim, lambda x: np.isclose(x[0], 0.5))] = 2 ft = dolfinx.mesh.meshtags(mesh, fdim, np.arange( num_facets, dtype=np.int32), facet_values) # Give a set of facets marked with a value (in this case 2), get a consistent orientation for an interior integral facets_to_integrate = ft.find(2) f_to_c = mesh.topology.connectivity(fdim, tdim) c_to_f = mesh.topology.connectivity(tdim, fdim) # Compute integration entities for a single facet of a cell. # Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process # local_facet_index is the local indexing of a facet for a given cell integration_entities = [] for i, facet in enumerate(facets_to_integrate): # Only loop over facets owned by the process to avoid duplicate integration if facet >= facet_map.size_local: continue # Find cells connected to facet cells = f_to_c.links(facet) # Get value of cells marked_cells = ct.values[cells] # Get the cell marked with 2 correct_cell = np.flatnonzero(marked_cells == 2) assert len(correct_cell) == 1 # Get local index of facet local_facets = c_to_f.links(cells[correct_cell[0]]) local_index = np.flatnonzero(local_facets == facet) assert len(local_index) == 1 # Append integration entities integration_entities.append(cells[correct_cell[0]]) integration_entities.append(local_index[0]) # Create custom integration measure for one-sided integrals breakpoint() ds = ufl.Measure("ds", domain=mesh, subdomain_data=[ (8, np.asarray(integration_entities, dtype=np.int32))]) n = ufl.FacetNormal(mesh) x = ufl.SpatialCoordinate(mesh) # Exact integral is [y/2**2]_0^1= 1/2 L = ufl.dot(ufl.as_vector((x[1], 0)), n)*ds(8) L_compiled = dolfinx.fem.form(L) print( f"Correct integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L_compiled), op=MPI.SUM)}") # Create reference implementation where we use a restricted two-sided integral with no notion of orientation dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft, subdomain_id=2) n = ufl.FacetNormal(mesh) x = ufl.SpatialCoordinate(mesh) L2 = ufl.dot(ufl.as_vector((x[1], 0)), n("+"))*dS L2_compiled = dolfinx.fem.form(L2) print( f"Wrong integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_compiled), op=MPI.SUM)}") ``` Source: https://fenicsproject.discourse.group/t/wrong-facetnormal-vector-on-internal-boundaries/12887/2?u=dokken

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

Chris_S · May 9, 2024, 5:19am

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

dokken · May 9, 2024, 5:43am

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

Chris_S · May 12, 2024, 6:10am

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

dokken · May 13, 2024, 7:29pm

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

github.com/jorgensd/dolfinx-tutorial

Add one-sided integration example as example of custom integration

opened 08:04PM - 28 Nov 23 UTC

jorgensd

```python # Computing a onesided integral over an interior facet with consisten…t orientation # Enabled by: https://github.com/FEniCS/dolfinx/pull/2269 # Copyright 2023 Jørgen S. Dokken # SPDX: MIT import dolfinx import numpy as np import ufl from mpi4py import MPI mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet) # Create connectivties required for defining integration entities tdim = mesh.topology.dim fdim = tdim - 1 mesh.topology.create_entities(fdim) mesh.topology.create_connectivity(fdim, tdim) mesh.topology.create_connectivity(tdim, fdim) # Get number of cells on process cell_map = mesh.topology.index_map(tdim) num_cells = cell_map.size_local + cell_map.num_ghosts # Create markers for each size of the interface cell_values = np.ones(num_cells, dtype=np.int32) cell_values[dolfinx.mesh.locate_entities( mesh, tdim, lambda x: x[0] <= 0.5+1e-13)] = 2 ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange( num_cells, dtype=np.int32), cell_values) facet_map = mesh.topology.index_map(fdim) num_facets = facet_map.size_local + facet_map.num_ghosts # Create facet markers facet_values = np.ones(num_facets, dtype=np.int32) facet_values[dolfinx.mesh.locate_entities( mesh, fdim, lambda x: np.isclose(x[0], 0.5))] = 2 ft = dolfinx.mesh.meshtags(mesh, fdim, np.arange( num_facets, dtype=np.int32), facet_values) # Give a set of facets marked with a value (in this case 2), get a consistent orientation for an interior integral facets_to_integrate = ft.find(2) f_to_c = mesh.topology.connectivity(fdim, tdim) c_to_f = mesh.topology.connectivity(tdim, fdim) # Compute integration entities for a single facet of a cell. # Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process # local_facet_index is the local indexing of a facet for a given cell integration_entities = [] for i, facet in enumerate(facets_to_integrate): # Only loop over facets owned by the process to avoid duplicate integration if facet >= facet_map.size_local: continue # Find cells connected to facet cells = f_to_c.links(facet) # Get value of cells marked_cells = ct.values[cells] # Get the cell marked with 2 correct_cell = np.flatnonzero(marked_cells == 2) assert len(correct_cell) == 1 # Get local index of facet local_facets = c_to_f.links(cells[correct_cell[0]]) local_index = np.flatnonzero(local_facets == facet) assert len(local_index) == 1 # Append integration entities integration_entities.append(cells[correct_cell[0]]) integration_entities.append(local_index[0]) # Create custom integration measure for one-sided integrals breakpoint() ds = ufl.Measure("ds", domain=mesh, subdomain_data=[ (8, np.asarray(integration_entities, dtype=np.int32))]) n = ufl.FacetNormal(mesh) x = ufl.SpatialCoordinate(mesh) # Exact integral is [y/2**2]_0^1= 1/2 L = ufl.dot(ufl.as_vector((x[1], 0)), n)*ds(8) L_compiled = dolfinx.fem.form(L) print( f"Correct integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L_compiled), op=MPI.SUM)}") # Create reference implementation where we use a restricted two-sided integral with no notion of orientation dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft, subdomain_id=2) n = ufl.FacetNormal(mesh) x = ufl.SpatialCoordinate(mesh) L2 = ufl.dot(ufl.as_vector((x[1], 0)), n("+"))*dS L2_compiled = dolfinx.fem.form(L2) print( f"Wrong integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_compiled), op=MPI.SUM)}") ``` Source: https://fenicsproject.discourse.group/t/wrong-facetnormal-vector-on-internal-boundaries/12887/2?u=dokken

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.

Chris_S · May 23, 2024, 1:08am

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:

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).
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])

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

Define DG-0 space:

V_c = fem.functionspace(domain, ("DG", 0))

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

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)

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

Topic		Replies	Views
Nonlinear incompressible inclusion problem with SNES solver dolfinx	0	409	May 30, 2023
Hyperelasticity example issues General	9	189	May 3, 2024
Nonlinear hyperelasticity in dolfin-x: Minimal example, and a few issues dolfinx	15	4653	May 29, 2022
Penalty Contact with circular shaped boundary Condition variational formulation	0	544	April 21, 2023
Hyperelastic model with two different materials General	1	187	February 4, 2024

Hyperelastic deformation with applied pressue on internal surface

Related Topics