Hello,
I am facing an issue where if I change the mesh a little, I get nan. Continuing from previous post https://fenicsproject.discourse.group/t/no-solution-for-simple-cantilever-beam-with-beam-element/14079
If I change the mesh I get nan, even if I don’t apply any load at the new part.
My mesh1 and mesh2 looks like this
Mesh1 works but mesh2 doesn’t.
Here is my MWE
from dolfin import *
import meshio
from ufl import Jacobian, diag
parameters['allow_extrapolation'] = True
msh = meshio.read("gmsh/t2/x-section_triangle_line.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points,
cells={"line": msh.cells_dict['line']}))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
thick = Constant(1)
width = thick/3
E = Constant(1)
nu = Constant(0.3)
G = E/2/(1+nu)
rho = Constant(1)
g = Constant(1)
S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = Constant(5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
f = Constant(-rho*g*S)
def tangent(mesh):
t = Jacobian(mesh)
return as_vector([t[0, 0], t[1, 0], t[2, 0]])/sqrt(inner(t, t))
def tgrad(u):
return dot(grad(u), t)
def generalized_strains(u):
(w, theta) = split(u)
return as_vector([dot(tgrad(w), t),
dot(tgrad(w), a1)-dot(theta, a2),
dot(tgrad(w), a2)+dot(theta, a1),
dot(tgrad(theta), t),
dot(tgrad(theta), a1),
dot(tgrad(theta), a2)])
def generalized_stresses(u):
return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))
t = tangent(mesh)
ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
def left_end(x, on_boundary):
return near(x[0], 0)
def right_end(x, on_boundary):
return near(x[0], 10, 1e-5) and near(x[1], 0, 1e-5)
Ue = VectorElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, Ue*Ue)
M = VectorFunctionSpace(mesh, "CG", 1, dim=3)
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), left_end)
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
AutoSubDomain(left_end).mark(facets, 1)
AutoSubDomain(right_end).mark(facets, 2)
dS = Measure("dS", domain=mesh, subdomain_data=facets, metadata = {"quadrature_scheme":"default","quadrature_degree": 3})
Sig = generalized_stresses(du)
Eps = generalized_strains(u_)
dx_shear = dx(scheme="default", metadata={"quadrature_scheme":"default", "quadrature_degree": 1})
k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + \
(Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
#
l_form = inner(f, 2*avg(u_[2]))*dS(2)
u = Function(W)
solve(k_form == l_form, u, bc)
print(0.5*assemble(inner(generalized_strains(u), generalized_stresses(u))*dx))
Mesh1 is a triangle which works
$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
3 3 0 0
1 0 0 0 0
2 10 0 0 0
3 5 5 0 0
1 0 0 0 10 0 0 0 2 1 -2
2 5 0 0 10 5 0 0 2 2 -3
3 0 0 0 5 5 0 0 2 3 -1
$EndEntities
$Nodes
6 3 1 3
0 1 0 1
1
0 0 0
0 2 0 1
2
10 0 0
0 3 0 1
3
5 5 0
1 1 0 0
1 2 0 0
1 3 0 0
$EndNodes
$Elements
6 6 1 6
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
1 1 1 1
4 1 2
1 2 1 1
5 2 3
1 3 1 1
6 3 1
$EndElements
Mesh2 is a triangle with a line out of plane. Which gives nan
$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
4 4 0 0
1 0 0 0 0
2 5 0 0 0
3 2.5 2.5 0 0
4 2.5 2.5 2 0
1 0 0 0 5 0 0 0 2 1 -2
2 2.5 0 0 5 2.5 0 0 2 2 -3
3 0 0 0 2.5 2.5 0 0 2 3 -1
4 2.5 2.5 0 2.5 2.5 2 0 2 3 -4
$EndEntities
$Nodes
8 4 1 4
0 1 0 1
1
0 0 0
0 2 0 1
2
5 0 0
0 3 0 1
3
2.5 2.5 0
0 4 0 1
4
2.5 2.5 2
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
$EndNodes
$Elements
8 8 1 8
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
0 4 15 1
4 4
1 1 1 1
5 1 2
1 2 1 1
6 2 3
1 3 1 1
7 3 1
1 4 1 1
8 3 4
$EndElements