Hi
I was trying to convert the 3D beam found here FEniCSx.
I am facing 2 issues.
- I was able to convert the code. But if i change the loading condition to concentrated load FEniCSx gives no solution for complicated mesh (For simple mesh works fine).
- the
.xdmffile written can’t be opened by paraview forMeshtagshowever openingmeshworks fine. And if i try writing thewrite_functionmy Jupyter notebook crashes and restarts.
For the First problem, If i use simple mesh everything works (concentrated load or self-weight). If i use a little complicated mesh self-weight problem works fine but concentrated one gives no solution.
I installed FEniCSx with conda and dolfinx version is 0.7.3
My MWE
# static 1D Cantilever Beam with a rectangular cross section
import matplotlib.pyplot as plt
import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
locate_dofs_topological,Constant,Expression, locate_dofs_geometrical, assemble_scalar)
from dolfinx.io import VTKFile, gmshio, XDMFFile, VTXWriter, VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities_boundary,create_interval ,locate_entities, meshtags
from dolfinx import mesh
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, dS, ds, Jacobian, as_vector, sqrt, cross, dot, avg,
diag,VectorElement, split, Measure)
import meshio
import gmsh
## Complicated mesh
msh = meshio.read("../fenics_nlopt/examples/Truss/gmsh/t3/test.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points = msh.points, cells={"line": msh.cells_dict['line']}))
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
## Simple mesh
#length = 10.0
#N = 40 # number of elements
#points = np.zeros((N+1, 3))
#points[:, 0] = np.linspace(0, length, N+1)
#cells = np.array([[i, i+1] for i in range(N)])
#meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
def tangent(domain):
t = Jacobian(domain)
return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))
t = tangent(domain)
ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
thick = Constant(domain,0.3)
width = thick/3
E = Constant(domain,70e3)
nu = Constant(domain,0.3)
G = E/2/(1+nu)
rho = Constant(domain,2.7e-3)
g = Constant(domain,9.81)
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(domain,5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
Ue = VectorElement("CG", domain.ufl_cell(), 1, dim=3)
W = FunctionSpace(domain, Ue*Ue)
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)
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))
Sig = generalized_stresses(du)
Eps = generalized_strains(u_)
def left(x):
return(np.isclose(x[0], 0, 1))
def load(x):
return np.isclose(x[0], 10, 1) #and np.isclose(x[1], 0, 12)
boundaries = [(1, left),
(2, load)]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
W0, submap0 = W.sub(0).collapse()
W1, submap1 = W.sub(1).collapse()
left_dofs_x = locate_dofs_geometrical((W.sub(0),W0),left )
left_dofs_r = locate_dofs_geometrical((W.sub(1),W1),left)
#right_dofs = locate_dofs_geometrical(W, right)
zero = Function(W0)
with zero.vector.localForm() as zero_local:
zero_local.set(0.0)
bcs = [
dirichletbc(zero, left_dofs_x, W.sub(0)),
dirichletbc(zero, left_dofs_r, W.sub(1))]
ds = Measure("ds", domain=domain, subdomain_data=facet_tag, metadata = {"quadrature_degree": 4})
dS = Measure("dS", domain=domain, subdomain_data=facet_tag,metadata = {"quadrature_degree": 4})
#dx = Measure("dx",domain=domain, subdomain_data=subdomains, metadata = {"quadrature_degree": 4})
dx_shear = dx(scheme="default",metadata={ "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(Constant(domain,-1.), 2*avg(u_[1]))*dS(2)#Works for simple mesh not for complicated one
#l_form = inner(Constant(domain,-1.), u_[1])*dx# Works for both
u = Function(W)
problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()
assemble_scalar((form(inner(uh,uh)*dx)))
with XDMFFile(domain.comm, "results/tag.xdmf", "w" ,encoding=XDMFFile.Encoding.HDF5) as xdmf:
xdmf.write_meshtags(facet_tag, domain.geometry)
with XDMFFile(domain.comm, "results/fun.xdmf", "w" ) as xdmf:
xdmf.write_function(uh)
the test.msh is below
$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
6 11 0 0
1 0 0 0 0
2 10 0 0 0
3 10 10 0 0
4 0 10 0 0
6 20 0 0 0
10 20 10 0 0
1 0 0 0 10 0 0 0 2 1 -2
2 10 0 0 10 10 0 0 2 2 -3
3 0 10 0 10 10 0 0 2 3 -4
4 0 0 0 0 10 0 0 2 4 -1
5 0 0 0 10 10 0 0 2 2 -4
6 0 0 0 10 10 0 0 2 3 -1
7 10 0 0 20 0 0 0 2 2 -6
8 20 0 0 20 10 0 0 2 6 -10
9 10 10 0 20 10 0 0 2 10 -3
11 10 0 0 20 10 0 0 2 6 -3
12 10 0 0 20 10 0 0 2 10 -2
$EndEntities
$Nodes
17 6 1 6
0 1 0 1
1
0 0 0
0 2 0 1
2
10 0 0
0 3 0 1
3
10 10 0
0 4 0 1
4
0 10 0
0 6 0 1
5
20 0 0
0 10 0 1
6
20 10 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 5 0 0
1 6 0 0
1 7 0 0
1 8 0 0
1 9 0 0
1 11 0 0
1 12 0 0
$EndNodes
$Elements
17 17 1 17
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
0 4 15 1
4 4
0 6 15 1
5 5
0 10 15 1
6 6
1 1 1 1
7 1 2
1 2 1 1
8 2 3
1 3 1 1
9 3 4
1 4 1 1
10 4 1
1 5 1 1
11 2 4
1 6 1 1
12 3 1
1 7 1 1
13 2 5
1 8 1 1
14 5 6
1 9 1 1
15 6 3
1 11 1 1
16 5 3
1 12 1 1
17 6 2
$EndElements
