How to impose a DirchletBC at the center of mass?

General
1

I have a unit square (constructed with Gmsh). I am solving a hyperelastic problem where the displacement of the domain is governed by a growth tensor growth tensor. I want to impose a zero displacement Dirichlet BC at the center of mass. However, the center of mass might not lie at a node if the mesh is not heavily refined. I want the zero Dirichlet BC to be imposed at a node that is closest to the center of mass. How do I achieve that? This is what I tried:

msh = meshio.read("square.msh") # Simple unit square mesh created in Gmsh

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

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_custom = Measure("dS", domain=mesh, subdomain_data=mf) # Defining global boundary measure, for global surface measure, plug in "dS" instead of "ds", which is used for subdomain                                                                   boundary measures
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

##########  Calculating center of mass #########

new_V = VectorFunctionSpace(mesh,"P",1)
new_R = VectorFunctionSpace(mesh,"R",0)
new_position = Function(new_V)
new_position.assign(Expression(["x[0]","x[1]","x[2]"], element=new_V.ufl_element()))
new_c = TestFunction(new_R)
new_bulk = assemble(Constant(1.0)*dx(domain=mesh))
new_centroid = assemble(dot(new_c,new_position)*dx)
new_f = new_centroid/new_bulk
new_f_np = new_f.get_local()

center = new_f_np

print(f"center of mass = {center}",flush=True)

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

V = VectorFunctionSpace(mesh, 'P', 1)

u = Function(V) 
v = TestFunction(V)

########################  Defining function to impose DBC at center   ##################

class Gamma(SubDomain):
   def inside(self, x, on_boundary):
     return near(center[0], 1e-16) and near(center[1], 1e-16) and near(center[2], 1e-16)  # I have tried with different values for the tolerance level

####################################################################
r = 0.02
nu1 = 0.01   
mu1 = 5.0

d = mesh.geometry().dim()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
g_bact = as_tensor([[1.0+r,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])

Fe = variable(F*inv(g_bact))
C = Fe.T*Fe
J  = det(Fe)

M = cofac(F)

# Invariants of deformation tensors
Ic = tr(C)     

psi = (mu1/2)*(Ic - 3) + ((nu1)*(mu1)/(1-2*nu1))*(J-1)**2 - mu1*(ln(J))

TT = det(g_bact)*diff(psi,Fe)*(inv(g_bact).T)  # First Piola-Kirchoff stress

# Compute first variation of Pi (directional derivative about u in the direction of v)
Pi = inner(TT,grad(v))*dx_custom(5) 

r_center = Constant((0.0,0.0,0.0))

bcs = DirichletBC(V, r_center, Gamma(), "pointwise")   
# Solve variational problem
solve(Pi == 0, u , bcs,
     solver_parameters={"newton_solver": {"relative_tolerance": 5e-8,   
     "absolute_tolerance": 5e-9,"maximum_iterations": 50}})
     
File('displacement.pvd')<<u

How do I get the Dirichlet BC be imposed at a node that is closest to the center of mass? Any suggestions would be extremely useful.

Thank you.

2

See Setting Delta-function-like initial condition - #2 by dokken

3

Hi @dokken ! Thanks for the suggestion. In the link, the displacement is made 0. But I am not sure how to pass it as a Dirichlet boundary condition. Following the suggestion in the link, I tried to make the displacement 0 at the center of mass after solving the problem. But upon viewing it in Paraview, I see that the displacement is not 0 at the center of mass. In fact, the solution is never zero. It is close to 0 far away from the center of mass. Here is the modified MWE (mesh is a simple unit square mesh created with Gmsh):

msh = meshio.read("square.msh") # Simple square mesh created in Gmsh

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

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_custom = Measure("dS", domain=mesh, subdomain_data=mf) # Defining global boundary measure, for global surface measure, plug in "dS" instead of "ds", which is used for subdomain                                                                   boundary measures
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

##########  Calculating center of mass #########

new_V = VectorFunctionSpace(mesh,"P",1)
new_R = VectorFunctionSpace(mesh,"R",0)
new_position = Function(new_V)
new_position.assign(Expression(["x[0]","x[1]","x[2]"], element=new_V.ufl_element()))
new_c = TestFunction(new_R)
new_bulk = assemble(Constant(1.0)*dx(domain=mesh))
new_centroid = assemble(dot(new_c,new_position)*dx)
new_f = new_centroid/new_bulk
new_f_np = new_f.get_local()


center = Point(0.5,0.5,0.0)

# Find which cell 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 = 1e3
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]

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

V = VectorFunctionSpace(mesh, 'P', 1)

# Find dof corresponding to vertex
vtd_p = vertex_to_dof_map(V)
dof_number = vtd_p[closest_vertex]

u = Function(V) 
v = TestFunction(V)

####################################################################
r = 0.02
nu1 = 0.01   
mu1 = 5.0

d = mesh.geometry().dim()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
g_bact = as_tensor([[1.0+r,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])

Fe = variable(F*inv(g_bact))
C = Fe.T*Fe
J  = det(Fe)

M = cofac(F)

# Invariants of deformation tensors
Ic = tr(C)     

psi = (mu1/2)*(Ic - 3) + ((nu1)*(mu1)/(1-2*nu1))*(J-1)**2 - mu1*(ln(J))

TT = det(g_bact)*diff(psi,Fe)*(inv(g_bact).T)  # First Piola-Kirchoff stress

# Compute first variation of Pi (directional derivative about u in the direction of v)
Pi = inner(TT,grad(v))*dx_custom(5) 

all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)

bc = DirichletBC(V.sub(1), Constant(0.0), all_domains, 1)  # making Y-component of displacement 0
bcs = [bcs2]
# Solve variational problem
solve(Pi == 0, u , bcs,
     solver_parameters={"newton_solver": {"relative_tolerance": 5e-8,   
     "absolute_tolerance": 5e-9,"maximum_iterations": 50}})
     
u.vector()[dof_number] = 0.0  # making displacement at center of mass 0
     
u.rename('u','u')
     
File('displacement.pvd')<<u

So here are my questions:

  1. How do I pass u.vector()[dof_number]=0 as a Dirichlet Boundary condition?
  2. I am writing the center of mass manually as center = Point(0.5,0.5,0.0). I have already calculated the center of mass as new_f_np. It’s a numpy array. How do I convert new_f_np to a Point? I need that because in my original problem, the center of mass is changing over time. So I will not be able to manually calculate it and put it in a Point() format.