Incompressible Mooney-Rivlin with Component Wise BCs

Hey All,

I am trying to setup unaxial extension (z-axis) on an annulus generated with gmsh. I want to use incompressible Mooney-Rivlin or at least nearly-incompressible. I have been scouring the discussion board and found quite a few examples of people doing similar things (i.e. Neo-Hookean, Component Wise BCs, Mooney-Rivlin).

However, my current setup is causing errors on running:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0

I imagine my errors are occurring in either my FunctionSpace setup:

n.ufl_cell(), degree=2)  
FE1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), degree=1)  
state_space = fem.FunctionSpace(domain, VE2 * FE1)
V, _ = state_space.sub(0).collapse()

Or BCs:

base_dofs_z = fem.locate_dofs_topological((state_space.sub(0).sub(2), V), facet_tag.dim, facet_tag.find(1))
u0_bc = fem.Function(V)
u0 = lambda x: np.zeros_like(x, dtype=default_scalar_type)
u0_bc.interpolate(u0)
bc_base_z = fem.dirichletbc(u0_bc, base_dofs_z, state_space.sub(0).sub(2))

But I am unsure, thank you for your help. The full code is below. Sorry for the superfluous commenting, I am trying to learn how it all works.

# +==+==+
# Setup
# += Imports
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 = 1.5

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
    # += Read .msh into domain for FEniCSx
    #    (1): File name .msh
    #    (2): Multiprocessing assignment
    #    (3): Rank of multiprocessing
    #    (4): Dimension of mesh
    domain, _, facet_markers = io.gmshio.read_from_msh("gmsh_msh/" + test_name + ".msh", MPI.COMM_WORLD, 0, gdim=MESH_DIM)
    # += Create Vector Element
    #    (1): Interpolation style
    #    (2): Cell from mesh domain
    #    (3): Degree of interpolation style
    VE2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), degree=2)  
    FE1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), degree=1)  
    # += Create Function Space
    #    (1): Mesh domain space
    #    (2): Finite element setup
    # += Find total Function Space which contains both interpolation schemes
    state_space = fem.FunctionSpace(domain, VE2 * FE1)
    # Collapse into just Geometric Space
    V, _ = state_space.sub(0).collapse()

    # +==+==+
    # Determine Boundaries
    fdim = MESH_DIM-1
    #    (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_x = fem.locate_dofs_topological((state_space.sub(0), V.sub(0)), facet_tag.dim, facet_tag.find(1))
    # base_dofs_y = fem.locate_dofs_topological((state_space.sub(0), V.sub(1)), facet_tag.dim, facet_tag.find(1))
    base_dofs_z = fem.locate_dofs_topological((state_space.sub(0).sub(2), V), facet_tag.dim, facet_tag.find(1))
    # += Set Dirichlet BCs of (1) on (2)
    # bc_base_x = fem.dirichletbc(default_scalar_type(0), base_dofs_x, (state_space.sub(0), V.sub(0)))
    # bc_base_y = fem.dirichletbc(default_scalar_type(0), base_dofs_y, (state_space.sub(0), V.sub(1)))
    u0_bc = fem.Function(V)
    u0 = lambda x: np.zeros_like(x, dtype=default_scalar_type)
    u0_bc.interpolate(u0)
    bc_base_z = fem.dirichletbc(u0_bc, base_dofs_z, state_space.sub(0).sub(2))
    # += Top Dirichlet BCs
    #    (1): Interpolation space
    #    (2): Internal function to determine if on z = 0
    top_dofs_z = fem.locate_dofs_topological((state_space.sub(0).sub(2), V), facet_tag.dim, facet_tag.find(2))
    # += Set Dirichlet BCs of (1) on (2)
    u1_bc = fem.Function(V)
    u1 = lambda x: np.zeros_like(x, dtype=default_scalar_type)
    u1_bc.interpolate(u1)
    bc_top_z = fem.dirichletbc(u1_bc, top_dofs_z, state_space.sub(0).sub(2))
    # bc_top_z = fem.dirichletbc(default_scalar_type(0.05), top_dofs_z, (state_space.sub(0), V.sub(2)))
    # += Concatenate boundaries
    bc = [bc_base_z, bc_top_z]

    # # +==+==+ 
    ## FOR LATER
    # # Pressure Setup
    # # += Pressure expression, contribution of pressure at internal radius
    # p = ufl.exp(('p*x[0]/R', 'p*x[1]/R'), R = R0)
    # # += Pressure Space (decreasing order of interpolation by 1) 
    # pre_interp = fem.FunctionSpace(domain, ("Lagrange", test_order-1))
    # pressure = fem.Function(pre_interp)
    # # +=  Interpolate expression over points
    # pre_expr = fem.Expression(p, pre_interp.element.interpolation_points())
    # pressure.interpolate(pre_expr)

    # +==+==+
    # Setup Parameteres for Variational Equation
    # += Test and Trial Functions
    state = fem.Function(state_space, name="state")
    test_state = ufl.TestFunctions(state_space)
    v, q = test_state
    u, p = ufl.split(state)
    # += 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)
    #    (3): λ1^2*λ2^2*λ3^2; det(C) = J^2
    Ic = ufl.variable(ufl.tr(C))
    # uniMod_Ic = ufl.variable(J**(-2/3) * Ic)
    #    (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)
    # uniMod_IIc = ufl.variable(J**(-4/3) * IIc)
    J = ufl.variable(ufl.det(F))
    # IIIc = ufl.variable(ufl.det(C))
    # += Material Parameters
    c1 = 2
    c2 = 6
    # kappa = 1
    # += Mooney-Rivlin Strain Energy Density Function
    # psi = c1 * (Ic - 3) + c2 *(IIc - 3) + kappa*(IIIc-1)**2
    psi = c1 * (Ic - 3) + c2 *(IIc - 3) 
    # Terms
    gamma1 = ufl.diff(psi, Ic) + Ic * ufl.diff(psi, IIc)
    gamma2 = -ufl.diff(psi, IIc)
    # sPK_term1 = ufl.diff(psi, Ic) + ufl.diff(Ic, E) + ufl.diff(psi, IIc) + ufl.diff(IIc, E)
    # sPK_term2 = ufl.diff(psi, Ic) + ufl.diff(Ic, E.T) + ufl.diff(psi, IIc) + ufl.diff(IIc, E.T)
    # sPK = 0.5 (sPK_term1 + sPK_term2)
    firstPK = F * 2 * ufl.diff(psi, C)
    # += First Piola Stress
    # piola = ufl.diff(psi, F)
    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=facet_tag, 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 
    # += Problem Setup and Solver
    problem = NonlinearProblem(R, state, 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(state)
    u_sol, p_sol = state.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", [u], engine="BP4") as vtx:
        vtx.write(0.0)

# +==+==+
# Main check for script operation.
#   Will operate argparse to take values from terminal. 
#   Then runs main() program for script execution 
if __name__ == '__main__':
    
    # +==+ Intake Arguments
    argparser = argparse.ArgumentParser("FEniCSz Program for Passive Annulus Inflation")
    # += Name for file writing
    argparser.add_argument("test_ID")
    # += Type of generation, i.e. Tetrahedral or Hexagonal
    argparser.add_argument("gen_type")
    # += Order of generated mesh
    argparser.add_argument("gen_order", type=int)
    # += Refinement level
    argparser.add_argument("refinement", type=bool)
    # += Capture arguments and store accordingly
    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)

Going to answer my own question here. It seems to be a matter of setting up the problem as a Mixed Element:

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)
V, _ = W.sub(0).collapse()
Vz, _ = V.sub(2).collapse()

And improving the description of the BCs:

base_dofs_z = fem.locate_dofs_topological((W.sub(0), V), facet_tag.dim, facet_tag.find(1))
u0_bc = fem.Function(V)
u0 = lambda x: np.zeros_like(x, dtype=default_scalar_type)
u0_bc.interpolate(u0)
bc_base_z = fem.dirichletbc(u0_bc, base_dofs_z, W.sub(0))
top_dofs_z = fem.locate_dofs_topological((W.sub(0).sub(2), Vz), facet_tag.dim, facet_tag.find(2))
u1_bc = fem.Function(Vz)
u1 = lambda x: np.full(x.shape[1], default_scalar_type(0.2))
u1_bc.interpolate(u1)
bc_top_z = fem.dirichletbc(u1_bc, top_dofs_z, W.sub(0).sub(2))
bc = [bc_base_z, bc_top_z]

Not entirely sure why, but I’ll take it. :partying_face: