How to define the cross section at the center of a 3D beam and impose a zero displacement Dirichlet BC on the vertices lying on the cross section defined at the center?

General
1

Hi! In the below MWE, I have a boxmesh defined by mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(10.0, 0.1, 0.1), 100, 10, 10). On this mesh, I am solving the PDE \eta\frac{\partial\mathbf u}{\partial t} = (\text{CSA})\nabla.\mathbf{T}, where \mathbf{u}, t, \eta, CSA, and \mathbf{T} represent displacement, time, friction coefficient, cross section area, and the 1st Piola–Kirchoff stress tensor. I am applying a hyperelastic framework with a growth tensor and a strain energy density function. I solve for the displacement \mathbf{u} that minimizes the stress.

I need to impose a zero displacement DirichletBC on the vertices lying on cross section at the center of the beam (the cross section at the center is a rectangle with the corners located at (5.0, 0.0, 0.0), (5.0, 0.0, 0.1), (5.0, 0.1, 0.1), and (5.0, 0.1, 0.0)). I tried to locate the vertices of the tetrahedra that lies close to the center center = Point(5.0, 0.05, 0.05). But it returns only a point [1.9 0.04 0.08] which is far away from the defined center and also returns only one point. I need all vertices lying on the cross section at the center of the beam. I am not sure what I am doing wrong in the below MWE

from dolfin import *
import ufl as ufl
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# Creating Rectangular mesh #

L = 10.0
H = 0.1
W = 0.1

mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(10.0, W, H), 100, 10, 10)

# Defining the volume measure of the rectangle #
    
volume_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
 
dv_custom = Measure("dx", domain=mesh, subdomain_data=volume_markers)

# Defining vector and scalar function spaces  #

F_bact = VectorFunctionSpace(mesh, 'P', 1)
V_bact = VectorFunctionSpace(mesh, 'P', 1)
T_bact = TensorFunctionSpace(mesh, 'P', 1)

###############

vtk1 = File("solutionbacteriagrowth.pvd")
vtk2 = File("sigmareference.pvd")
vtk3 = File("sigmadeformed.pvd")
vtk4 = File("velocity_time.pvd")

##############

u = Function(V_bact) 
u_n = Function(V_bact) 
v = TestFunction(V_bact)

vel = Function(V_bact)

############

csa_second_moment = (H**4)/12
csa = H**2  

NumStepg = 100 
Time = 60  
dt = Time/NumStepg
j = 0
eta = 5.0  
r = 0.020  
mu1 = 2500.0/csa_second_moment 
nu1 = 0.4999

#  Defining function to impose DBC on the cross section at center 

center = Point(L/2, H/2, W/2)

# Find which tetrahedra point is located in

tree = mesh.bounding_box_tree()
cell, distance = tree.compute_closest_entity(center)
vertices = mesh.cells()[cell]

# Find the closest vertex in the cell
vertex_coordinates = mesh.coordinates()[vertices]
dist = 1e-1
closest = None
for i, v1 in enumerate(vertex_coordinates):
    p_ = Point(v1)
    dist_ = center.distance(p_)
    print(v1)
    if dist_ < dist:
        dist = dist_
        closest_local = i
closest_vertex = vertices[closest_local]


closestv = mesh.coordinates()[closest_vertex]

# Find dof corresponding to vertex

vtd_p = vertex_to_dof_map(F_bact)
dof_number = vtd_p[closest_vertex]

print(f"DOF number of closest vertex = {dof_number}")
print(f"closest vertex = {V_bact.tabulate_dof_coordinates()[closest_vertex]}")

class cm(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return near(x[0], closestv[0], tol)
        
#############

for T in np.linspace(0,Time,NumStepg): 

 d = mesh.geometry().dim()
 I = Identity(d)             # Identity tensor
 F = I + grad(u)             # Deformation gradient
 Fg_theta = as_tensor([[exp(r*T),0.0, 0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])  
 
 Fe = variable(F*inv(Fg_theta))
 C = variable(Fe.T*Fe)
 J  = det(Fe)
 
 # Invariants of deformation tensors
 
 I1 = tr(C)  
 
 W = (mu1/2)*(I1 - 3) + ((nu1)*(mu1)/(1-(2*nu1)))*(J-1)**2 - mu1*(ln(J)) 

 TT = det(Fg_theta)*diff(W,Fe)*(inv(Fg_theta).T) 

 sigma = (2/J)*((mu1/2)*(Fe*Fe.T - I) + (nu1*mu1*(J-1)*J*I)/(1-2*nu1))
 
 PE = (eta)*dot(u,v)*dx_custom() + csa*dt*inner(sigma,grad(v))*dx_custom() - (eta)*dot(u_n, v)*dx_custom()
 
 all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
 bcs1 = DirichletBC(V_bact.sub(0), Constant(0.0), cm(), method='pointwise') 
 bcs2 = DirichletBC(V_bact.sub(1), Constant(0.0), cm(), method='pointwise')
 bcs3 = DirichletBC(V_bact.sub(2), Constant(0.0), cm(), method='pointwise')

 bcs = [bcs1, bcs2, bcs3]
 
  
 # Solve variational problem
 solve(PE == 0, u , bcs,
      solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,   
      "absolute_tolerance": 5e-10,"maximum_iterations": 50}})
 
 
 print(f"Number of iterations: {j}", flush=True)
 print(f"Time = {T} mins", flush=True)
 
 j += 1
 
 vel.assign((u - u_n)/dt)
 
 u_n.assign(u)
 
 disp = project(u, V_bact)
 disp.rename('disp','disp')

 velocity = project(vel, V_bact)
 velocity.rename("velocity","velocity")

 vtk1<<(disp, T)
 vtk4<<(velocity, T)  

################  END  ###############

In addition to the above problem, I also need to keep the cross-section constant with time across the length of the beam as it grows and deforms with time. I am not sure how to implement that.

Any suggestions would be greatly appreciated!
Thanks in advance!

2

Consider the following MWE:

from dolfin import *
mesh = UnitCubeMesh(6, 2,2)
V = VectorFunctionSpace(mesh, "Lagrange", 1)



class CrossSection(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.5)
    

mf = MeshFunction("size_t", mesh, 0, 0)
CrossSection().mark(mf, 1)

vertices = mf.where_equal(1)

vertex_to_dofmap = vertex_to_dof_map(V)
# As each vertex goes to three dofs, we reshape the array such that each row corresponds to a vertex
vertex_to_dofmap = vertex_to_dofmap.reshape((-1, mesh.geometry().dim()))

dofs = vertex_to_dofmap[vertices]
print(V.tabulate_dof_coordinates()[dofs])
1 Like
3

@dokken Thanks for the suggestion. That solves the problem. I have a follow-up question that I mentioned at the end of my original post.

In the MWE, as the 3D beam deforms, due to elastic stretch, the cross section area of the beam along its length does not remain constant. How do I impose a condition such that the cross section areas along the length of the beam remain constant as the beam deforms?