Applying Neumann BC in parallel

dolfinx
1

I want to solve a linear elasticity problem of a bar under unilateral static tension. The bar is clamped at x=0 and a Neumann BC is applied at the free end at x=2.

import numpy as np
from dolfinx import cpp, Constant, VectorFunctionSpace, FunctionSpace, Function, DirichletBC, BoxMesh
from ufl import Identity, TestFunction, TrialFunction, grad, inner, sym, tr, Measure, dot
from dolfinx.fem import locate_dofs_geometrical, assemble_scalar, assemble_vector, assemble_matrix
from dolfinx.mesh import locate_entities_boundary, MeshTags
from petsc4py import PETSc
from dolfinx.io import XDMFFile
from mpi4py import MPI

solver = PETSc.KSP().create(MPI.COMM_WORLD)
mesh = BoxMesh(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]), np.array([2,1,1])], [10,4,4], cell_type=cpp.mesh.CellType.hexahedron)

V = VectorFunctionSpace(mesh, ("Lagrange",1))
u = TrialFunction(V)
v = TestFunction(V)
u_k = Function(V)

E = 1.0
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
lmbda = nu / ((1.0 - 2 * nu)*(1 + nu))

def eps(v_):
  return sym(grad(v_))

def sigma(v_):
  return 2.0 * mu * eps(v_) + lmbda * tr(eps(v_)) * Identity(3)

def clamped_end(x):
    return np.isclose(x[0], 0.0)

def free_end(x):
  return np.isclose(x[0], 2)

uvw_clamped_end = Function(V)

with uvw_clamped_end.vector.localForm() as uvw_clamped_end_loc:
  uvw_clamped_end_loc.set(0)

facets_clamped_end = locate_entities_boundary(mesh, mesh.topology.dim-1, clamped_end)
facets_free_end = locate_entities_boundary(mesh, mesh.topology.dim-1, free_end)

bc_u = [''] * 1
bc_u[0] = DirichletBC(uvw_clamped_end, locate_dofs_geometrical(V, clamped_end))

dx = Measure("dx")
mt = MeshTags(mesh, mesh.topology.dim-1, facets_free_end, 1)
ds = Measure("ds", subdomain_data=mt)

t = Constant(mesh, (1,0,0))
a = inner(sigma(u), eps(v)) * dx
l = dot(t, v) * ds(1)
elastic_energy = inner(sigma(u_k), eps(u_k)) * dx

A = assemble_matrix(a, bc_u)
A.assemble()
b = assemble_vector(l)
solver.setOperators(A)
solver.solve(b, u_k.vector)
u_k.x.scatter_forward()

local_elastic_energy = assemble_scalar(elastic_energy)
gathered_elastic_energy = MPI.COMM_WORLD.allgather(local_elastic_energy)
global_elastic_energy = sum(gathered_elastic_energy)
print('Rank: ',MPI.COMM_WORLD.rank,'local energy: ', local_elastic_energy, 'global energy:', global_elastic_energy)

filename_u = '/u.xdmf'
with XDMFFile(MPI.COMM_WORLD, filename_u, "w", encoding=XDMFFile.Encoding.HDF5) as file:
  file.write_mesh(mesh)
  file.write_function(u_k)

Using mpirun the global elastic energies are as expected and do not differ for -np 1, -np 2 or -np 3. However, for -np 4 the result is different. Inspecting the displacements in paraview shows that there might be something wrong with the Neumann BC at the free end where two mesh partitions border each other. How can I fix this?

displ_paraview
displ_paraview1524×1068 213 KB

1 Like
2

You have forgotten two things in your code:

  1. Using ghost update on your RHS
  2. Applying bcs through lifting on the RHS

The following additions fixes your issues:

from dolfinx.fem import apply_lifting, set_bc
A = assemble_matrix(a, bc_u)
A.assemble()
b = assemble_vector(l)
apply_lifting(b, [a], [bc_u])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bc_u)
solver.setOperators(A)
solver.solve(b, u_k.vector)
u_k.x.scatter_forward()
2 Likes