Right boundary conditions for linear elastic compression of a cube

dolfinx
1

What are the right boundary conditions to add in order to obtain the right results on uniaxial plate compression of a cube. I need to kill some translations and rotations. So here some BCs are missing.

I tried everything I was thinking about, but the results are never the ones intended.

Here is the code:

from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
W = 1.0
mu = 1.0
rho = 1.0
d = W/20
beta = 1.25
lambda_ = beta

############# Domain

domain = mesh.create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([W, W, W])],
    [20, 20, 20],
    cell_type=mesh.CellType.hexahedron,
)

V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))


######## scalar BC for compression

# Function subspace
Vx = V.sub(0)
Vy = V.sub(1)
Vz = V.sub(2)

# Boundary markers
def boundary_left(x):
    return np.isclose(x[0], 0)

def boundary_right(x):
    return np.isclose(x[0], W)



fdim = domain.topology.dim - 1
facets_l = mesh.locate_entities_boundary(domain, fdim, boundary_left)
facets_r = mesh.locate_entities_boundary(domain, fdim, boundary_right)

# Locate dofs (IMPORTANT tuple syntax)
dofs_l_x = fem.locate_dofs_topological((Vx, V), fdim, facets_l)[1]
dofs_r_x = fem.locate_dofs_topological((Vx, V), fdim, facets_r)[1]

# Dirichlet BCs on ux only
bc_l = fem.dirichletbc(
    fem.Constant(domain, default_scalar_type(d)),
    dofs_l_x,
    Vx
)

bc_r = fem.dirichletbc(
    fem.Constant(domain, default_scalar_type(0.0)),
    dofs_r_x,
    Vx
)


# --- Combine with left/right compression BCs ---
bcs = [
    bc_l,  # ux on x=0
    bc_r
]

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

ds = ufl.Measure("ds", domain=domain)


########## Variational formulation

def epsilon(u):
    return ufl.sym(
        ufl.grad(u)
    )  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(T, v) * ds



########### Solve the linear variational problem

problem = LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"}#,
    #petsc_options_prefix="linear_elasticity", ########### -> error: not recognized: TypeError: LinearProblem.__init__() got an unexpected keyword argument 'petsc_options_prefix'
)
uh = problem.solve()

The core of the code comes from this tutorial.

(Note that petsc_options_prefix=“linear_elasticity” is not working for me for print(dolfinx.version)
→0.9.0).

The code gives this kind of results (or variations of it following the other BC I tried)

image
image601×551 131 KB

Thx in advance

2

Simplest is to consider only a quarter of the cube and enfore uy=0 on the y=0 plane and uz=0 on the z=0 plane

3

This is not how I would apply the boundary conditions.
I would do:

V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))


######## scalar BC for compression

# Function subspace
Vx = V.sub(0)
Vy = V.sub(1)
Vz = V.sub(2)

# Boundary markers
def boundary_left(x):
    return np.isclose(x[0], 0)

def boundary_right(x):
    return np.isclose(x[0], W)



fdim = domain.topology.dim - 1
facets_l = mesh.locate_entities_boundary(domain, fdim, boundary_left)
facets_r = mesh.locate_entities_boundary(domain, fdim, boundary_right)

# Locate dofs (IMPORTANT tuple syntax)
dofs_l_x = fem.locate_dofs_topological(Vx, fdim, facets_l)
dofs_r_x = fem.locate_dofs_topological(Vx, fdim, facets_r)
# Dirichlet BCs on ux only
bc_l = fem.dirichletbc(
    fem.Constant(domain, default_scalar_type(d)),
    dofs_l_x,
    Vx
)

bc_r = fem.dirichletbc(
    fem.Constant(domain, default_scalar_type(0)),
        dofs_r_x,
    Vx
)

You can also try to set the null space to PETSc, i.e.

from dolfinx import mesh, fem, plot, io, default_scalar_type
from petsc4py import PETSc
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
W = 1.0
mu = 1.0
rho = 1.0
d = W/20
beta = 1.25
lambda_ = beta

############# Domain

domain = mesh.create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([W, W, W])],
    [20, 20, 20],
    cell_type=mesh.CellType.hexahedron,
)

V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))


######## scalar BC for compression

# Function subspace
Vx = V.sub(0)
Vy = V.sub(1)
Vz = V.sub(2)

# Boundary markers
def boundary_left(x):
    return np.isclose(x[0], 0)

def boundary_right(x):
    return np.isclose(x[0], W)



fdim = domain.topology.dim - 1
facets_l = mesh.locate_entities_boundary(domain, fdim, boundary_left)
facets_r = mesh.locate_entities_boundary(domain, fdim, boundary_right)

# Locate dofs (IMPORTANT tuple syntax)
dofs_l_x = fem.locate_dofs_topological(Vx, fdim, facets_l)
dofs_r_x = fem.locate_dofs_topological(Vx, fdim, facets_r)
# Dirichlet BCs on ux only
bc_l = fem.dirichletbc(
    fem.Constant(domain, default_scalar_type(d)),
    dofs_l_x,
    Vx
)

bc_r = fem.dirichletbc(
    fem.Constant(domain, default_scalar_type(0)),
        dofs_r_x,
    Vx
)


# --- Combine with left/right compression BCs ---
bcs = [
    bc_l,  # ux on x=0
    bc_r
]

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

ds = ufl.Measure("ds", domain=domain)


########## Variational formulation

def epsilon(u):
    return ufl.sym(
        ufl.grad(u)
    )  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(T, v) * ds



########### Solve the linear variational problem
from dolfinx import la
from dolfinx import default_scalar_type as dtype
def build_nullspace(V):
    """Build PETSc nullspace for 3D elasticity."""
    # Create vectors that will span the nullspace
    bs = V.dofmap.index_map_bs
    length0 = V.dofmap.index_map.size_local
    # Remove translation in x-direction as it is not in nullspace
    basis = [la.vector(V.dofmap.index_map, bs=bs, dtype=dtype) for i in range(3)]
    b = [b.array for b in basis]

    # Get dof indices for each subspace (x, y and z dofs)
    dofs = [V.sub(i).dofmap.list.flatten() for i in range(3)]

    # Set the three translational rigid body modes
    for i in range(2):
        b[i][dofs[i+1]] = 1.0

    # Set the three rotational rigid body modes
    x = V.tabulate_dof_coordinates()
    dofs_block = V.dofmap.list.flatten()
    x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
    b[2][dofs[2]] = x1
    b[2][dofs[1]] = -x2

    la.orthonormalize(basis)

    basis_petsc = [
        PETSc.Vec().createWithArray(x[: bs * length0], bsize=3, comm=V.mesh.comm) for x in b
    ]
    return PETSc.NullSpace().create(vectors=basis_petsc)

nsp = build_nullspace(V)

problem = LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                   "ksp_error_if_not_converged": True,
                   "ksp_monitor": None,
                   "pc_factor_mat_solver_type": "mumps"},
    petsc_options_prefix="a"
)
problem.A.setNullSpace(nsp)
uh = problem.solve()

assert nsp.test(problem.A)

import dolfinx.fem.petsc
with io.VTXWriter(domain.comm, "u.bp", [uh]) as bp:
    bp.write(0.0)

which would remove the y-z translation + and rotation in the y-z plane.