I am solving a problem of uniform thermal expansion in a cylinder (boundary condition: displacement in the z-direction at the top is fixed). When using direct solution methods (the displacement results in the z-direction are shown in Figure 1), the boundary condition was successfully implemented. However, due to insufficient constraints (resulting in an overall shift), I increased the null space and used an iterative method to solve it (displacement results in the z-direction are shown in Figure 2). It is evident that the boundary condition is not met; on the contrary, the displacement at the bottom is larger than at the top. Can you explain what might be happening? I hope to receive some help. Thank you!
Here is my code:
#Direct method
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import io, default_scalar_type
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, grad)
from dolfinx.fem import (functionspace, Constant, dirichletbc, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
#Mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh("unitf1.msh", MPI.COMM_WORLD, gdim=3)
#Parameters
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
def eps(v):
return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
return (lmbda*tr(eps(v))-alpha*(3*lmbda + 2*mu)*dT)*Identity(3)+2.0*mu*eps(v)
#Function space
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim, )))
#Top: z = 0
U_D = Constant(mesh, default_scalar_type(0))
fdim = mesh.topology.dim - 1
facets_top = facet_markers.find(5)
dofs_top_z = locate_dofs_topological(VU.sub(2), fdim, facets_top)
bcU = dirichletbc(U_D, dofs_top_z, VU.sub(2))
#Variational formula
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
#Solve
problem = LinearProblem(aM, LM, bcs=[bcU], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
U = problem.solve()
#Save results
vtkFile = io.VTKFile(MPI.COMM_WORLD, "results/U.pvd", "w")
vtkFile.write_function(U)
vtkFile.close()
#Iterative method using null space
import dolfinx
from dolfinx import la
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import io, default_scalar_type
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, grad)
from dolfinx.fem import (functionspace, FunctionSpace, Function, form, Constant, dirichletbc, locate_dofs_topological)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector, apply_lifting, set_bc)
def null_space(V: FunctionSpace):
_x = Function(V)
gdim = V.mesh.geometry.dim
assert gdim == 2 or gdim == 3
dim = 3 if gdim == 2 else 6
nullspace_basis = [
la.vector(V.dofmap.index_map, bs=V.dofmap.index_map_bs, dtype=PETSc.ScalarType)
for i in range(dim)
]
basis = [b.array for b in nullspace_basis]
dofs = [V.sub(i).dofmap.list.reshape(-1) for i in range(gdim)]
for i in range(gdim):
basis[i][dofs[i]] = 1.0
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.reshape(-1)
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
if gdim == 2:
basis[2][dofs[0]] = -x1
basis[2][dofs[1]] = x0
elif gdim == 3:
basis[3][dofs[0]] = -x1
basis[3][dofs[1]] = x0
basis[4][dofs[0]] = x2
basis[4][dofs[2]] = -x0
basis[5][dofs[2]] = x1
basis[5][dofs[1]] = -x2
for b in nullspace_basis:
b.scatter_forward()
_basis = [x._cpp_object for x in nullspace_basis]
dolfinx.cpp.la.orthonormalize(_basis)
assert dolfinx.cpp.la.is_orthonormal(_basis)
local_size = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
basis_petsc = [
PETSc.Vec().createWithArray(x[:local_size], bsize=gdim, comm=V.mesh.comm)
for x in basis
]
return PETSc.NullSpace().create(comm=V.mesh.comm, vectors=basis_petsc)
#Mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
#Parameters
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
def eps(v):
return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
return (lmbda*tr(eps(v))-alpha*(3*lmbda + 2*mu)*dT)*Identity(3)+2.0*mu*eps(v)
#Function space
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim, )))
#Top: z = 0
U_D = Constant(mesh, default_scalar_type(0))
fdim = mesh.topology.dim - 1
facets_top = facet_markers.find(5)
dofs_top_z = locate_dofs_topological(VU.sub(2), fdim, facets_top)
bcU = dirichletbc(U_D, dofs_top_z, VU.sub(2))
# Variational formula
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
A = assemble_matrix(form(aM), bcs = [bcU])
A.assemble()
b = assemble_vector(form(LM))
apply_lifting(b, [form(aM)], bcs=[[bcU]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bcU])
#Set solver options
opts = PETSc.Options()
opts["ksp_type"] = "cg"
opts["ksp_rtol"] = 1.0e-7
opts["pc_type"] = "gamg"
#Use Chebyshev smoothing for multigrid
opts["mg_levels_ksp_type"] = "chebyshev"
opts["mg_levels_pc_type"] = "jacobi"
#Improve estimate of eigenvalues for Chebyshev smoothing
opts["mg_levels_ksp_chebyshev_esteig_steps"] = 10
#Create PETSc Krylov solver and turn convergence monitoring on
solver = PETSc.KSP().create(mesh.comm)
solver.setFromOptions()
solver.setOperators(A)
null_space = null_space(VU)
A.setNullSpace(null_space)
null_space.remove(b)
#Solve
U = Function(VU)
solver.solve(b, U.vector)
#Save results
vtkFile = io.VTKFile(MPI.COMM_WORLD, "results/U.pvd", "w")
vtkFile.write_function(U)
vtkFile.close()
along with the visualization results and my geo file below:
// Gmsh project created on Tue May 21 09:56:22 2024
SetFactory("OpenCASCADE");
r1 = 3;
N = 31;
Circle(1) = {0, 0, 0, r1, 0, 2*Pi};
Transfinite Curve {1} = N Using Progression 1;
Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 10} {Surface{1};Layers {8};}
Physical Surface("cylinder", 4) = {2};
Physical Surface("top", 5) = {3};
Physical Surface("bottom", 6) = {1};
Physical Volume("v", 7) = {1};