I wish to apply 2 periodic boundary conditions on a square plate.
On the left and right, the pbc will apply in all directions.
On the top and bottom, The pbc should apply in the direction perpendicular to the force getting applied.
For e.g., If I’m applying tension upwards on the top boundary, the bottom boundary should only mimic the displacements of the top boundary in the X direction (sideways)
My code is currently as follows
import meshio
import dolfin
from dolfin import XDMFFile
import numpy as np
from dolfin import *
import matplotlib as plt
msh = meshio.read("/mnt/c/Users/ABC/sfepy/simplerve.msh")
for key in msh.cell_data_dict["gmsh:geometrical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:geometrical"][key]
else:
line_data = msh.cell_data_dict["gmsh:geometrical"][key]
for cell in msh.cells:
if cell.type == "line":
tetra_cells = cell.data
if cell.type == "triangle":
triangle_cells = cell.data
def create_mesh(mesh, cell_type, prune_z=True):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points[:,:2], cells={cell_type: cells}, cell_data={"gmsh:physical":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
meshio.write("mesh_1.xdmf", triangle_mesh)
line_mesh = create_mesh(msh, "line", prune_z=False)
meshio.write("mf.xdmf", line_mesh)
from mpi4py import MPI as mpi
from scipy.spatial.distance import cdist
import numpy as np
mesh = Mesh()
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh_1.xdmf") as infile:
infile.read(mesh)
infile.read(mvc_2d, "gmsh:physical")
cell_markers = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc_2d, "gmsh:physical")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
force = Constant((1.0, 1.0))
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds",domain=mesh, subdomain_data=boundary_markers)
tol_0 = 1E-14
tol_1 = 0.99999999999
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
# return True if on left or bottom boundary
return bool(((near(x[0],0) and on_boundary) or (near(x[1],0) and on_boundary)) and (not((near(x[0],1) and near (x[1],0) and on_boundary) or (x[0],0) and near (x[1],1)and on_boundary)))
def map(self, x, y,):
if (near(x[0],1) and near (x[1],1)): # if on top-right corner
y[0] = x[0]
y[1] = x[1] -1
elif(near(x[0],1)): # if on right boundary
y[0] = x[0] -1
y[1] = x[1]
elif(near(x[1],1)): #top boundary
y[0] = x[0]
y[1] = x[1] -1
else:
y[0] = x[0]
y[1] = x[1] -1
# Create periodic boundary condition
pbc = PeriodicBoundary()
W = VectorFunctionSpace(mesh, "P", 2, constrained_domain=pbc)
lam_values = [0.005,121] #Lame Modulus
mu_values = [0.007, 81.9] #Shear Modulus
mu = Function(W)
lam = Function(W)
for cell_no in range(len(cell_markers.array())):
subdomain_no = cell_markers.array()[cell_no]
#print(subdomain_no)
subdomain_no = int(subdomain_no-1)
#print(subdomain_no)
mu.vector()[cell_no] = mu_values[subdomain_no]
lam.vector()[cell_no] = lam_values[subdomain_no]
def clamped_boundary_bottom(x, tol):
return (x[1] < tol_0)
def free_boundary_top(x, tol):
return (x[1] > tol_1)
bc_bottom_fix = DirichletBC(W.sub(1), Constant((0.0)), clamped_boundary_bottom)
bc_top_fix = DirichletBC(W.sub(1), Constant((0.1)), free_boundary_top)
bcs = [bc_bottom_fix, bc_top_fix,]
def epsilon(u):
return 0.5*(nabla_grad(u) + (nabla_grad(u).T))
from ufl import nabla_div
def sigma_matrix(u):
return lam_values[0]*nabla_div(u)*Identity(d) + 2*mu_values[0]*epsilon(u)
def sigma_fiber(u):
return lam_values[1]*nabla_div(u)*Identity(d) + 2*mu_values[1]*epsilon(u)
u = TrialFunction(W)
d = u.geometric_dimension() # space dimension
v = TestFunction(W)
f = Constant((0,0)) #Body Force
T = Constant((0.0, 0.0)) #Traction Load
a = (inner(sigma_matrix(u), epsilon(v))*dx(1)) + (inner(sigma_fiber(u), epsilon(v))*dx(2))
L = ((dot(f, v)))*dx(1) + (dot(f, v))*dx(2) + (dot(T, v))*ds
u = Function(W)
solve(a == L, u, bcs)
The result I get is as follows
Which is not what I want (The material undergoes zero contraction)
I am aware regarding dolfinx_mpc, but I’m facing an issue with hdf5 (I have posted this on github).
Is there a way to correctly define the pbcs as I want?
My mesh - Mesh