Setting vector component-wise Dirichlet BCs

dolfinx
1

Hi

I have been successfully running a nonlinear mixed vector field problem (simple uni-axial tension of a rod) using a non-homogeneous Dirichlet boundary condition incrementally applied as a vector, i.e [ u_x, u_y, u_z] = [ t*1.0, 0.0, 0.0 ] at the right end of the rod (the left end is clamped):

msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([L, H, H])], [20, 1, 1],
                      mesh.CellType.hexahedron)
P1 = ufl.VectorElement("CG", msh.ufl_cell(), 2, dim=3)  # Lagrange family, 2nd order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))

v, yf = ufl.TestFunctions(VY)     # test functions
uwf = fem.Function(VY)            # current displacement and change of director
u, wf = ufl.split(uwf)

def right(x):
    return np.isclose(x[0], L)

class incremented_displacement_expression:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        # linearly incremented displacement
        return np.stack((np.full(x.shape[1], self.t * 1.0), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

V0, u_dofs = VY.sub(0).collapse()
incremented_displacement_expr = incremented_displacement_expression()
incremented_displacement_expr.t = 0
incremented_displacement = fem.Function(V0)
incremented_displacement.interpolate(incremented_displacement_expr.eval)
right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), right)
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, VY.sub(0))

skipping some not relevant lines

for n in range(1, 10):

    # update the inhomogeneous displacement boundary condition
    incremented_displacement_expr.t = n
    incremented_displacement.interpolate(incremented_displacement_expr.eval)

    # solve the equation system
    num_its, converged = solver.solve(uwf)

I realized now that the actual problem only prescribes the u_x component but u_y and u_z are both kept free, i.e. I only need to apply the Dirichlet BC component-wise for u_x.

I am following the example showing vector component-wise https://jorgensd.github.io/dolfinx-tutorial/chapter3/component_bc.html but wanted to stick to geometrical DOF localization and I also had to transfer the example code into the mixed environment. I got it running without errors but I don’t get a converged solution anymore, although the load (the non-homogeneous Dirichlet BC) has effectively not changed:

class incremented_x_displacement_expression:
    def __init__(self):
        self.t = 0.0

    def eval(self, x):
        # linearly incremented displacement
        return np.full(x.shape[1], self.t * 1.0)

V0, u_dofs = VY.sub(0).collapse()
V0x, u_x_dofs = V0.sub(0).collapse()

incremented_x_displacement_expr = incremented_x_displacement_expression()
incremented_x_displacement_expr.t = 0
incremented_x_displacement = fem.Function(V0x)
incremented_x_displacement.interpolate(incremented_x_displacement_expr.eval)
right_ux_dofs = fem.locate_dofs_geometrical((V0.sub(0),V0x), right)
bc2 = fem.dirichletbc(incremented_x_displacement, right_ux_dofs, V0.sub(0))

skipping some not relevant lines

for n in range(1, 10):

    # update the inhomogeneous displacement boundary condition
    incremented_x_displacement_expr.t = n
    incremented_x_displacement.interpolate(incremented_x_displacement_expr.eval)

    # solve the equation system
    num_its, converged = solver.solve(uwf)

Please advise. Many thanks in advance.

2

You can apply the .sub() function repeatedly without collapsing the subspaces, i.e.:

right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0).sub(0),V0x), right)
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, VY.sub(0).sub(0))

When you call sub(), a subspace is returned that retains a link to its parent subspace. When you collapse a subspace, a new FunctionSpace is created that has no knowledge of its parent function space. When you pass the collapsed subspace to fem.dirichletbc, the Dirichlet bc has no knowledge of the parent function space and thus the bc cannot affect the dofs of the function you are trying to solve (since it is in the parent subspace).

Here’s an adapted code based on my previous answer that applies the incremented bc only to u_x:

from contextlib import ExitStack

import numpy as np

import ufl
from dolfinx import fem, nls, mesh, plot, log, io
from mpi4py import MPI
from petsc4py import PETSc

import pdb


# ----------------------------------------------------------------------------------------------------------------------
# Create a box mesh of a beam
L = 20.0
H = 1.0
msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([L, H, H])], [10, 1, 1],
                      mesh.CellType.hexahedron)
P1 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 1, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))


# ----------------------------------------------------------------------------------------------------------------------
# Boundary conditions and body loads
def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], L)

def top(x):
    return np.isclose(x[1], 1)

def bottom(x):
    return np.isclose(x[1], 0)


def fixed_displacement_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def incremented_displacement_expression(x):
    return np.full(x.shape[1], 1.0e-02)

def fixed_director_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))

# ----------------------------------------------
# Dirichlet boundary application to affected displacement dofs

# Fix displacement DOF on left face and apply axial displacement on right face
V0, submap = VY.sub(0).collapse()
V0x, subsubmap = V0.sub(0).collapse()
fixed_displacement = fem.Function(V0)
fixed_displacement.interpolate(fixed_displacement_expression)
incremented_displacement = fem.Function(V0x)
incremented_displacement.interpolate(incremented_displacement_expression)
left_u_dofs = fem.locate_dofs_geometrical((VY.sub(0),V0), left)
right_u_dofs = fem.locate_dofs_geometrical((VY.sub(0).sub(0),V0x), right)
bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, VY.sub(0))
bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, VY.sub(0).sub(0))

# Fix ALL change of director degrees of freedom in the domain
V1, submap = VY.sub(1).collapse()
fixed_director = fem.Function(V1)
fixed_director.interpolate(fixed_director_expression)
top_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), top)
bottom_wf_dofs = fem.locate_dofs_geometrical((VY.sub(1),V1), bottom)
bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, VY.sub(1))
bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,VY.sub(1))

bcs = [bc1,  bc2, bc3, bc4]
pdb.set_trace()

# ----------------------------------------------------------------------------------------------------------------------
# Function space
v, yf = ufl.TestFunctions(VY)     # test functions
uwf = fem.Function(VY)            # current displacement and change of director
u, wf = ufl.split(uwf)

# ----------------------------------------------------------------------------------------------------------------------
# Kinematics

# 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
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

# ----------------------------------------------------------------------------------------------------------------------
# Constitutive law

# Set the elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(msh, E/(2*(1 + nu)))
lmbda = fem.Constant(msh, E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi_m = (mu / 2.0) * (Ic - 3.0) - mu * ufl.ln(J) + (lmbda / 2.0) * (ufl.ln(J))**2

# 1st Piola Kirchhoff stress tensor
P = ufl.diff(psi_m, F)                                     # 1st Piola Kirchhoff stress tensor
df = fem.Constant(msh, PETSc.ScalarType((0.0, 0.0, 0.0)))  # dummy director stress

# ----------------------------------------------------------------------------------------------------------------------
# Variational formulation
metadata = {"quadrature_degree": 3}
dx = ufl.Measure("dx", domain=msh, metadata=metadata)

# Define form Pi (we want to find u such that Pi(u) = 0)
Pi = ufl.inner(P, ufl.grad(v))*dx + ufl.inner(df, yf)*dx

# Initialize nonlinear problem
problem = fem.petsc.NonlinearProblem(Pi, uwf, bcs)

# ----------------------------------------------------------------------------------------------------------------------
# Initialise Newton Rhapson solver
solver = nls.petsc.NewtonSolver(msh.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

# ----------------------------------------------------------------------------------------------------------------------
# Incremental application of the traction loading
log.set_log_level(log.LogLevel.INFO)
for n in range(1, 2):
    num_its, converged = solver.solve(uwf)
    assert(converged)

with io.XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    xdmf.write_function(uwf.split()[0])
6 Likes
3

Thanks a lot, on the dot!