Hey All,
Was just reading the following: Wrong FacetNormal vector on internal boundaries
I am trying to implement the solution to see if it can fix my inflation test on the inner boundary of an annulus. The code uses a gmsh generated annulus with full incompressibility and Mooney-Rivlin. It worked fine for axial extension but trying to implement pressure has been an issue which brought me to learn about the facetNormals
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from numpy import linalg as LNG
from mpi4py import MPI
import numpy as np
import argparse
import meshio
import gmsh
import ufl
# += Parameters
MESH_DIM = 3
LEN_Z = 1
R0 = 1
R1 = 1.5
P_INNER = 0.15
P_OUTER = 0
def gen_annulus(test_name, test_type, test_order, refine_check):
# +==+==+
# Initialise and begin geometry
gmsh.initialize()
gmsh.model.add(test_name)
# += Setup base of annulus, surface
innerCircle = gmsh.model.occ.addCircle(x=0, y=0, z=0, r=R0, tag=1)
outerCircle = gmsh.model.occ.addCircle(x=0, y=0, z=0, r=R1, tag=2)
innerCircleCurve = gmsh.model.occ.addCurveLoop([innerCircle], 1)
outerCircleCurve = gmsh.model.occ.addCurveLoop([outerCircle], 2)
baseSurface = gmsh.model.occ.addPlaneSurface(wireTags=[outerCircle, innerCircle], tag=1)
# += Synchronize and add physical group
basePhysicalGroup = gmsh.model.addPhysicalGroup(dim=2, tags=[baseSurface], tag=100, name="Annulus Base Surface")
gmsh.model.occ.synchronize()
# += Extrude from geometry
baseExtrusion = gmsh.model.occ.extrude(dimTags=[(2, baseSurface)], dx=0, dy=0, dz=LEN_Z)
gmsh.model.occ.synchronize()
# += Create Physical Group on volume
volumePhysicalGroup = gmsh.model.addPhysicalGroup(3, [baseExtrusion[1][1]], tag=1000, name="Internal Volume")
innerPhysicalGroup = gmsh.model.addPhysicalGroup(2, [baseExtrusion[3][1]], tag =101, name="Inner Surface")
outerPhysicalGroup = gmsh.model.addPhysicalGroup(2, [baseExtrusion[2][1]], tag =102, name="Outer Surface")
topPhysicalGroup = gmsh.model.addPhysicalGroup(2, [baseExtrusion[0][1]], tag =103, name="Annulus Top Surface")
# += Create Mesh
gmsh.model.mesh.generate(MESH_DIM)
if refine_check == "True":
gmsh.model.mesh.refine()
gmsh.model.mesh.setOrder(test_order)
# += Write File
gmsh.write("gmsh_msh/" + test_name + ".msh")
gmsh.finalize()
def main(test_name, test_type, test_order, refine_check):
# +==+==+
# Mesh Generation
gen_annulus(test_name, test_type, test_order, refine_check)
# +==+==+
# Load Domain & Interpolation
domain, _, facet_markers = io.gmshio.read_from_msh("gmsh_msh/" + test_name + ".msh", MPI.COMM_WORLD, 0, gdim=MESH_DIM)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_entities(fdim)
domain.topology.create_connectivity(fdim, tdim)
domain.topology.create_connectivity(tdim, fdim)
VE2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), degree=2)
FE1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), degree=1)
W = fem.FunctionSpace(domain, ufl.MixedElement([VE2, FE1]))
w = fem.Function(W)
# += Collapse into geometric space
V, _ = W.sub(0).collapse()
Vx, _ = V.sub(0).collapse()
Vy, _ = V.sub(1).collapse()
Vz, _ = V.sub(2).collapse()
# +==+==+
# Determine Boundaries
# (1): Base
base_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[2], 0))
# (2): Top
top_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[2], LEN_Z))
# (3): Inner Cylinder Surface
def inner_bound(x):
return np.isclose(np.sqrt(x[0]**2+x[1]**2), R0)
inner_facets = mesh.locate_entities_boundary(domain, fdim, inner_bound)
# (4): Outer Cylinder Surface
def outer_bound(x):
return np.isclose(np.sqrt(x[0]**2+x[1]**2), R1)
outer_facets = mesh.locate_entities_boundary(domain, fdim, outer_bound)
# += Collate Marked boundaries
marked_facets = np.hstack([base_facets, top_facets, inner_facets, outer_facets])
# += Assign boundaries IDs
marked_values = np.hstack([
np.full_like(base_facets, 1),
np.full_like(top_facets, 2),
np.full_like(inner_facets, 3),
np.full_like(outer_facets, 4)
])
# += Sort and assign
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
# +==+==+
# Boundary Conditions
# += Base Dirichlet BCs but broken per subspace
# (1): Interpolation space
# (2): Internal function to determine if on z = 0
base_dofs_z = fem.locate_dofs_topological((W.sub(0).sub(2), Vz), facet_tag.dim, facet_tag.find(1))
# += Set Dirichlet BCs of (1) on (2)
u0_bc = fem.Function(Vz)
u0 = lambda x: np.full(x.shape[1], default_scalar_type(0.0))
u0_bc.interpolate(u0)
bc_base_z = fem.dirichletbc(u0_bc, base_dofs_z, W.sub(0).sub(2))
# += Top Dirichlet BCs
top_dofs_z = fem.locate_dofs_topological((W.sub(0).sub(2), Vz), facet_tag.dim, facet_tag.find(2))
# # += Set Dirichlet BCs of (1) on (2)
u1_bc = fem.Function(Vz)
u1 = lambda x: np.full(x.shape[1], default_scalar_type(0.0))
u1_bc.interpolate(u1)
bc_top_z = fem.dirichletbc(u1_bc, top_dofs_z, W.sub(0).sub(2))
# += Concatenate boundaries
bc = [bc_base_z, bc_top_z]
# # +==+==+
# # Pressure Setup
p_inner = fem.Constant(domain, default_scalar_type(P_INNER))
p_outer = fem.Constant(domain, default_scalar_type(P_OUTER))
n = ufl.FacetNormal(domain)
## START OF SOLUTION [Wrong FacetNormal vector on internal boundaries]
cell_map = domain.topology.index_map(tdim)
num_cells = cell_map.size_local + cell_map.num_ghosts
cell_values = np.ones(num_cells, dtype=np.int32)
cell_values[mesh.locate_entities(
domain, tdim, inner_bound
)] = 3
cell_tag = mesh.meshtags(
domain, tdim,
np.arange(num_cells, dtype=np.int32), cell_values
)
facet_map = domain.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts
facets_to_integrate = facet_tag.find(3)
f_to_c = domain.topology.connectivity(fdim, tdim)
c_to_f = domain.topology.connectivity(tdim, fdim)
integration_entities = []
for i, facet in enumerate(facets_to_integrate):
if facet >= facet_map.size_local:
continue
cells = f_to_c.links(facet)
marked_cells = cell_tag.values[cells]
correct_cell = np.flatnonzero(marked_cells == 2)
assert len(correct_cell) == 1
local_facets = c_to_f.links(cells[correct_cell[0]])
local_index = np.flatnonzero(local_facets == facet)
assert len(local_index) == 1
integration_entities.append(cells[correct_cell[0]])
integration_entities.append(local_index[0])
## END OF SOLUTION [Wrong FacetNormal vector on internal boundaries]
# +==+==+
# Setup Parameteres for Variational Equation
# += Test and Trial Functions
v, q = ufl.TestFunctions(W)
u, p = ufl.split(w)
# # += Identity Tensor
I = ufl.variable(ufl.Identity(MESH_DIM))
# += Deformation Gradient Tensor
# F = ∂u/∂X + I
F = ufl.variable(I + ufl.grad(u))
# += Right Cauchy-Green Tensor
C = ufl.variable(F.T * F)
# += Invariants
# (1): λ1^2 + λ2^2 + λ3^2; tr(C)
Ic = ufl.variable(ufl.tr(C))
# (2): λ1^2*λ2^2 + λ2^2*λ3^2 + λ3^2*λ1^2; 0.5*[(tr(C)^2 - tr(C^2)]
IIc = ufl.variable((Ic**2 - ufl.inner(C,C))/2)
# (3): λ1^2*λ2^2*λ3^2; det(C) = J^2
J = ufl.variable(ufl.det(F))
# IIIc = ufl.variable(ufl.det(C))
# += Material Parameters
c1 = 2
c2 = 6
# += Mooney-Rivlin Strain Energy Density Function
psi = c1 * (Ic - 3) + c2 *(IIc - 3)
# Terms
gamma1 = ufl.diff(psi, Ic) + Ic * ufl.diff(psi, IIc)
gamma2 = -ufl.diff(psi, IIc)
# += First Piola Stress
firstPK = 2 * F * (gamma1*I + gamma2*C) + p * J * ufl.inv(F).T
# +==+==+
# Setup Variational Problem Solver
# += Gaussian Quadrature
metadata = {"quadrature_degree": 4}
# += Domains of integration
ds = ufl.Measure('ds', domain=domain, subdomain_data=[
(3, np.asarray(integration_entities, dtype=np.int32))], metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# += Residual Equation (Variational, for solving)
R = ufl.inner(ufl.grad(v), firstPK) * dx + q * (J - 1) * dx - p_inner * J * ufl.dot(ufl.inv(F).T * n, v) * ds(3)
problem = NonlinearProblem(R, w, bc)
solver = NewtonSolver(domain.comm, problem)
# += Tolerances for convergence
solver.atol = 1e-8
solver.rtol = 1e-8
# += Convergence criteria
solver.convergence_criterion = "incremental"
num_its, converged = solver.solve(w)
u_sol, p_sol = w.split()
if converged:
print(f"Converged in {num_its} iterations.")
else:
print(f"Not converged after {num_its} iterations.")
# +==+==+
# ParaView export
with io.VTXWriter(MPI.COMM_WORLD, test_name + ".bp", w.sub(0).collapse(), engine="BP4") as vtx:
vtx.write(0.0)
vtx.close()
if __name__ == '__main__':
# +==+ Intake Arguments
argparser = argparse.ArgumentParser("FEniCSz Program for Passive Annulus Inflation")
argparser.add_argument("test_ID")
argparser.add_argument("gen_type")
argparser.add_argument("gen_order", type=int)
argparser.add_argument("refinement", type=bool)
args = argparser.parse_args()
test_name = args.test_ID
test_type = args.gen_type
test_order = args.gen_order
refine_check = args.refinement
# +==+ Feed Main()
main(test_name, test_type, test_order, refine_check)
# python annulus_inflation.py quick_check Tetra 2 False
(this can be run as $python file_name.py TEST Tetra 2 False)
If anyone could just explain how to properly implement pressure on the internal boundary that would be super.
Thank you 