hello, let’s consider a 3D problem, such as elasticity problem of a 3D beam. My question is how to apply load on an edge (line).
I can use MeshValueFuntion and measure ds to apply loads on a surface. And in 2D problem, that also works. But that seems not work on 3D edge. It would be better if you have an example, thanks in advance!
What in particular does not work when you use meshfunction for edges and corresponding measure?
ok，I use gmsh to build a 8x1x1 beam model ''temp3d.msh ‘’ and convert it into xdmf format. One side of the beam (pysicalgroup1) is fixed and another side is under a uniformed load T=(0,0,-10) on a edge (pysicalgroup6). I put my geo file and code below.
There are so many character in msh file that it exceeds the limitation. I put the geo here, gmsh $MeshFormat is 4.1 0 8
// Gmsh project created on Sat Oct 24 12:37:33 2020
//+
Point(1) = {0, 0, 0, 0.2};
//+
Point(2) = {0, 0, 1, 0.2};
//+
Point(3) = {0, 1, 1, 0.2};
//+
Point(4) = {0, 1, 0, 0.2};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {2, 3, 4, 1};
//+
Plane Surface(1) = {1};
//+
Extrude {5, 0, 0} {
Surface{1};
}
//+
Extrude {3, 0, 0} {
Surface{26};
}
//+
Physical Surface("bc") = {1};
//+
Physical Surface("mid") = {26};
//+
Physical Surface("right") = {48};
//+
Physical Volume("material1") = {1};
//+
Physical Volume("material2") = {2};
//+
Physical Curve("edge") = {29};
//+
Physical Point("p1") = {16};
//+
Physical Point("p2") = {20};
//+
Physical Surface("suf") = {47};
I use jupyter notebook and copy my code here:
import meshio
import numpy as np
msh=meshio.read("temp3d.msh")
tetra_cells=[]
triangle_cells=[]
line_cells=[]
vertex_cells=[]
for cell in msh.cells:
if cell.type == "triangle":
if len(triangle_cells)==0:
triangle_cells = cell.data
else:
triangle_cells = np.vstack([triangle_cells, cell.data])
elif cell.type == "tetra":
if len(tetra_cells)==0:
tetra_cells = cell.data
else:
tetra_cells = np.vstack([tetra_cells, cell.data])
elif cell.type == "line":
if len(line_cells)==0:
line_cells = cell.data
else:
line_cells = np.vstack([line_cells,cell.data])
elif cell.type == "vertex":
if len(vertex_cells)==0:
vertex_cells = cell.data
else:
vertex_cells = np.vstack([vertex_cells,cell.data])
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "line":
line_data = msh.cell_data_dict["gmsh:physical"][key]
elif key == "vertex":
vertex_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points,
cells={"tetra": tetra_cells},
cell_data={"volumn":[tetra_data]}
)
triangle_mesh =meshio.Mesh(points=msh.points,
cells={"triangle": triangle_cells},
cell_data={"surface":[triangle_data]}
)
line_mesh = meshio.Mesh(points=msh.points,
cells={"line":line_cells},
cell_data={"edge":[line_data]}
)
vertex_mesh = meshio.Mesh(points=msh.points,
cells={"vertex":vertex_cells},
cell_data={"point":[vertex_data]}
)
meshio.write("mesh_3d.xdmf", tetra_mesh)
meshio.write("surface_3d.xdmf", triangle_mesh)
meshio.write("edge_3d.xdmf", line_mesh)
meshio.write("point_3d.xdmf", vertex_mesh)
from dolfin import *
import matplotlib.pyplot as plt
mesh = Mesh()
with XDMFFile("mesh_3d.xdmf") as infile:
infile.read(mesh)
#plot(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("edge_3d.xdmf") as infile:
infile.read(mvc,"edge")
edge = cpp.mesh.MeshFunctionSizet(mesh, mvc)
V = VectorFunctionSpace(mesh,"P",1)
bc1 = DirichletBC(V,Constant((0,0,0)),surface,1)
bc = [bc1]
u = TrialFunction(V)
v = TestFunction(V)
from ufl import nabla_div
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)
alpha = Constant(1e-5)
def sigma(u):
return lambda_*nabla_div(u)*Identity(3) + 2*mu*epsilon(u)
# Define a
a = inner(sigma(u), epsilon(v))*dx
# Define L
ds=Measure("ds",mesh,6,subdomain_data=edge)
T=Constant((0,0,-10))
f = Constant((0,0,0))
L=dot(T,v)*ds+dot(f,v)*dx
u = Function(V)
solve(a == L, u, bc)
plt.figure()
p=plot(u*200)
plt.colorbar(p)
u_magnitude = sqrt(dot(u, u))
plt.figure()
p=plot(u_magnitude)
plt.colorbar(p)
I can run the codes and it will not raise error.But the solution seems wrong. The end displacement should be around FL^3/3EI=0.4 , which is quite different from fencis solution.The result of code is as followed.
Dear @Mh.X,
As your mesh it quite simplistic, it would have been alot easier to reproduce your issue with the built in meshes:
import numpy as np
from dolfin import *
mesh = BoxMesh(Point(0,0,0), Point(8,1,1), 80,10,10)
edge_function = MeshFunction("size_t", mesh, mesh.topology().dim()-2, 0)
class Edge(SubDomain):
def inside(self, x, on_boundary):
return np.isclose(x[0], 8) and np.isclose(x[1], 1) and on_boundary
Edge().mark(edge_function,1)
class Surface(SubDomain):
def inside(self, x, on_boundary):
return np.isclose(x[0], 0)
V = VectorFunctionSpace(mesh,"P",1)
bc1 = DirichletBC(V,Constant((0,0,0)),Surface())
ds=Measure("ds",mesh,subdomain_data=edge_function, subdomain_id=1)
Unfortunately, edge integration is currently not supported in dolfin. There has been some preliminary work by @cdaversin, but it is unfortunately in a publishable state. To approximate an edge traction, you could integrate over the surface, where your traction function has the corresponding value at the edge dofs, and zero elsewhere. At a fine mesh this should give an approximation to an edge traction
Hi, dokken. It is a pity that edge intergration not suppported. I will consider your advice. thank you!
Hi Mr. dokken. I want to know that has the current fenicsx v0.7.3 supported the edge integration. And if not, would you give me some suggestions about applying the edge uniform loading in 3D problem.