Hi,
I am trying to select the right facets (on an r=0.25 circle on top of a cube). To do this I am using locate_entities_boundary() and meshtags():
T = fem.Constant(domain, default_scalar_type((0.0, 0.0, -1.0)))
f = fem.Constant(domain, default_scalar_type((0.0, 0.0, 0.0)))
traction_facets = mesh.locate_entities_boundary(domain, 1, lambda x : np.sqrt((x[0]-grasp_node[0])**2 + (x[1]-grasp_node[1])**2 + (x[2]-grasp_node[2])**2) < force_radius)
mt = mesh.meshtags(domain, domain.geometry.dim-1, traction_facets, 1)
ds = ufl.Measure("ds", subdomain_data=mt)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(T, v) * ufl.dx + ufl.dot(T, v) * ds(1)
When plotting in ParaView it seems the traction is being applied possibly as a body force but not as a surface traction on top of the cube. Is there anything wrong with the code? A minimal example code is shown below:
from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, default_scalar_type, io
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, set_bc #,create_vector
import basix
import numpy as np
def distToCenter(tri_vertex, ctr):
return np.sqrt((tri_vertex[0]-ctr[0])**2 + (tri_vertex[1]-ctr[1])**2 + (tri_vertex[2]-ctr[2])**2)
# Import nodes and connectivity matrix into Fenics domain and mesh objects
points = np.array([[0.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 1.0, 1.0],
[0.0, 0.0, 1.0],
[1.0, 0.0, 0.0],
[1.0, 1.0, 0.0],
[1.0, 1.0, 1.0],
[1.0, 0.0, 1.0]])
cells = np.array([[0, 1, 3, 4],
[7, 3, 4 ,6],
[1, 2, 3, 6],
[3, 1, 6, 4],
[5, 1, 6, 4]])
finiteElement = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, ))
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, points, finiteElement)
fdim = domain.topology.dim - 1
for n in np.arange(4):
domain.topology.create_entities(1)
domain = mesh.refine(domain)
domain.topology.create_connectivity(fdim, fdim+1)
# Variables
fix_radius = 0.25
force_radius = 0.25
E = 3.5e9
nu = 0.35
lambda_ = E * nu / ((1.0 + nu) * (1 - 2.0*nu))
mu = 0.5 * E / (1.0 + nu)
# FEA functions
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)
# Calculate grasp node
grasp_node = np.array([0.5, 0.5, 1.0])
grasp_facet_node = np.array([0.5, 0.5, 0.0])
# Get facets that are within sphere with center at grasp_facet_node
free_end_facets = mesh.exterior_facet_indices(domain.topology)
free_end_connectivity2 = mesh.entities_to_geometry(domain, fdim, free_end_facets)
fixed_facets = np.empty((0), dtype=int)
for p in range(free_end_connectivity2.shape[0]):
dist1 = distToCenter(domain.geometry.x[free_end_connectivity2[p,0],:].squeeze(), grasp_facet_node)
dist2 = distToCenter(domain.geometry.x[free_end_connectivity2[p,1],:].squeeze(), grasp_facet_node)
dist3 = distToCenter(domain.geometry.x[free_end_connectivity2[p,2],:].squeeze(), grasp_facet_node)
if (dist1 < fix_radius) and (dist2 < fix_radius) and (dist3 < fix_radius):
fixed_facets = np.append(fixed_facets, p)
## FEA SETUP DEPENDING ON DOMAIN2
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
T = fem.Constant(domain, default_scalar_type((0.0, 0.0, -1.0)))
f = fem.Constant(domain, default_scalar_type((0.0, 0.0, 0.0)))
traction_facets = mesh.locate_entities_boundary(domain, 1, lambda x : np.sqrt((x[0]-grasp_node[0])**2 + (x[1]-grasp_node[1])**2 + (x[2]-grasp_node[2])**2) < force_radius)
mt = mesh.meshtags(domain, domain.geometry.dim-1, traction_facets, 1)
ds = ufl.Measure("ds", subdomain_data=mt)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(T, v) * ufl.dx + ufl.dot(T, v) * ds(1)
## Solution setup
a_compiled = fem.form(a)
L_compiled = fem.form(L)
# Create solution function
uh = fem.Function(V)
# Solver
solver = PETSc.KSP().create(domain.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
## SOLVE FEA
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, free_end_facets[fixed_facets]), V)
# Assemble system, applying boundary conditions
A = assemble_matrix(a_compiled, bcs=[bc])
A.assemble()
b = assemble_vector(L_compiled)
apply_lifting(b, [a_compiled], bcs=[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# Solver
solver.setOperators(A)
# Compute solution
solver.solve(b, uh.x.petsc_vec)
with io.XDMFFile(domain.comm, "deformation" + str(0) + ".xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
uh.name = "Deformation"
xdmf.write_function(uh)
