Single 2D Truss mesh

mesh
1

Hello everybody,

I am trying to solve a single truss element problem defined by two points A and B. In each point I have two DoF i.e. u1,u2,u3,u4.
I tried to use the tutorial for 2d beams see and an old thread `but I could not solve my problem defining the correct mesh and boundary condition.

So far I have this code:

import numpy as np
from dolfin import *

# Material variables
E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# Define Single Truss mesh from point A to point B with DoF 2 in each point
#
#

mesh = Mesh(mesh)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

# Define Dirichlet boundary i.e. left point A is fixed
#
#

def eps(v):
    return sym(grad(v))

def sigma(u):
    d = u.geometric_dimension()
    I = Identity(d)
    return lmbda * tr(eps(u)) * I + 2 * mu * eps(u)


# Solution variables
dx = Measure("dx")
du = TrialFunction(V)
U = Function(V)
u = U.vector()
u_ = TestFunction(V)
d = U.geometric_dimension()
I = Identity(d)
f = Constant((0, 0))

# Assembly system
a = inner(sigma(du), eps(u_)) * dx
L = inner(f, u_) * dx
K, f_vec = assemble_system(a, L, bc)

# External force on point B
f_vec[3] = 0
f_vec[4] = -1

# Solve system
solve(K, u, f_vec)
2

To create your mesh, you can use external software such as pygmsh which is a python wrapper for gmsh, or use the gmsh GUI directly.
There are many examples on how to export meshes from gmsh to dolfin, see: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio .
For setting boundary conditions at nodes, see for instance: Pointsource on specific points python

3

This does not answer my question. The problem I have is the definition of a truss element. For example I just need a mesh which consists of one Element and both nodes have DoF 2 i.e. displacemnt in direction x and y. Do I really need an external software?

4

Should the element be quadrialteral?
Does this do the trick?

from dolfin import *
mesh = UnitSquareMesh.create(MPI.comm_world, 1, 1,
                             CellType.Type.quadrilateral)
VE = VectorElement("DQ", quadrilateral, 1)
V = FunctionSpace(mesh, VE)
5

Unfortunatly no. I am searching for a method to define something like: Unbenannt

The blue line is one Element defined by two points. The left side is fixed by Dirichlet boundary. On the right point a force is acting in the y-direction. Therefore I would expect a displacement vector [u1_x, u1_y, u2_x, u2_y]. Where u1 is the displacement of the left point and u2 of the right point.

6

That is just an interval mesh with a 2D vector space:

from dolfin import *
L = 2
mesh = IntervalMesh(1, 0, L)
VE = VectorElement("CG", interval, 1, dim=2)
V = FunctionSpace(mesh, VE)

you will then have two dofs at 0 and two dofs at L.

1 Like
7

Exactly what I searched for. Just a further question. How would I define the mesh if I have something like this:

Unbenannt

8

Hi,
you can use the mesh Editor

mesh = Mesh()
ed = MeshEditor()
ed.open(mesh, "interval", 1, 2)
ed.init_vertices(3)
ed.add_vertex(0, [0., 0.])
ed.add_vertex(1, [1.,1])
ed.add_vertex(2, [0, 1])
ed.init_cells(2)
ed.add_cell(0, [0, 1])
ed.add_cell(1, [1, 2])
ed.close()
1 Like
9

I get an error calculating epsilon with

def eps(v):
    return 1/2 * (nabla_grad(v) + nabla_grad(v).T)

and

 def eps(v):
        return sym(nabla_grad(v))

The error message is: Can’t add expressions with different shapes. I can resolve this issue by just returning nabla_grad(v). I guess this has something to do with taking the transpose of a vector and adding them.

10

Without a minimal working example I can’t give you much advise.

11

If I understood correctly you just want to solve a simple truss problem. The axial strain measure must therefore be computed from the tangential part of the displacement field in the local frame of a truss member \epsilon = d(\underline{u}\cdot\underline{t})/ds which is scalar. You cannot consider the full symmetrized gradient of the displacement (this is for 2D/3D solids).
You can adapt this demo on 3D beam elements which involves all 6 internal forces of a beam (axial force, 2 bending moments, 2 shear forces and one torsion moment). You can then adapt it ignoring rotations and all forces/strains except the axial force/strain.

12

MY working example looks like:

import numpy as np
from dolfin import *

# Material variables
L = 2; W = 0.0000001
E = Constant(3e5)

# Define Single Truss mesh from point A to point B with DoF 2 in each point

mesh = IntervalMesh(1, 0, L)
VE = VectorElement("CG", interval, 1, dim=2)
V = FunctionSpace(mesh, VE)
mesh = Mesh(mesh)

# Define Dirichlet boundary
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)

def eps(v):
    return nabla_grad(v)

def sig(u):
    return E * eps(u)


# Solution variables
dx = Measure("dx")
du = TrialFunction(V)
U = Function(V)
u = U.vector()
v = TestFunction(V)

f = Constant((0, 0))
T = Constant((0, -10))
a = inner(sig(du), eps(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Assembly system
K, f_vec = assemble_system(a, L, bc)
print(np.array(f_vec))

# Solve system
solve(K, u, f_vec)

u = np.array(u)
print(u)
13

Oh that is nice. I will try to work with this too. Thank you :slight_smile:

14

Dear @loop ,

Did you manage to solve your problem?

I would like to see the code with this two-bar truss!

15

Hi! I would like to ask you did you solve also the problem with the two beams? I am trying to find similar way to do it in 3D.

16

To be honest, I switched to Ferrite.jl

17

Updates on this topic in FEniCSx

1 Like