Dear All,
I am handling a bar in GetFEM++ framework. This post is questioning about
unexpected non-zero values in a stiffness matrix K. I believe that thinking
about it helps deeper understanding of the framework.
The problem considered is:
"""
Bar:
0----1
(point: 0, 1; convex: 0 (0-1), segment or 1d-simplex))
(2 dof (u, v) on a single point)
Coordinate system:
^ y, v, Fy
|
----> x, u, Fx
Expected for the convex 0 (0-1):
K U = F
[ 1, 0, -1, 0] [u0] [Fx0]
[ 0, 0, 0, 0] [v0] [Fy0]
[-1, 0, 1, 0] [u1] [Fx1]
[ 0, 0, 0, 0] [v1] [Fy1]
(K: stiffness matrix, U: displacement vector, F: force vector)
(A bar should have zero stiffness in the transverse direction (i.e. y-
or v-direction here).)
"""
I believe that we can handle a bar using FEM_PK(1, 1) ((1, 1): first 1 for
1d-simplex (i.e. segment), second 1 for 2-nodes on a segment).
However, obtained K has unexpected non-zero values:
"""
Obtained for the convex 0 (0-1) :
K U = F
[ 1, 0, -1, 0] [u0] [Fx0]
[ 0, 0.5, 0, -0.5] [v0] [Fy0]
[ -1, 0, 1, 0] [u1] [Fx1]
[ 0, -0.5, 0, 0.5] [v1] [Fy1]
(K22, K24, K42, K44 are expected to be 0.)
"""
Can we know why K has 0.5 values like this? This behavior is unexpected for
me (or possibly for people with structural background), but I guess that
this behavior is logically expectable and explainable. I searched
getfem-users-archives and looked into FEM_PK, GT_PK and some
implementations but could not find the reason. The reason would be about
the same interpolation of u and v if |K22|, |K24|, |K42|, |K44| were 1, but
actually they are 0.5.
Could you please give your comments? Any ideas are welcome. It really helps
understanding the framework.
Attached is the script for python3. The script explicitly obtains stiffness
matrix K using gf.asm_linear_elasticity(...). The other part is similar to
Roman's one:
[Getfem-users] Solving truss structure in GetFEM++ (via Python)
https://lists.nongnu.org/archive/html/getfem-users/2011-03/msg00009.html
*I like GetFEM++ way: separative design (GeoTrans, Mesh, Fem, Im, ...);
high and low layers (Model, asm_...); point and convex; and so on.and so
on.
Best regards,
Alex Spring
"""
A single bar problem
0----1
^ Fy
|
^ y, v
|
----> x, u
bar: stiffness E*A = 1
length L = 1
Load: Fy = 1
Constraints(displacements: u,v): u_0 = 0
v_0 = 0
"""
"""
It yields:
KU = F
[ 1.0 0.0 -1.0 0.0 ][ 0.0 ] [ 0.0 ]
[ 0.0 0.5 0.0 -0.5 ][ 0.0 ] [ 0.0 ]
[ -1.0 0.0 1.0 0.0 ][ 0.0 ] [ 0.0 ]
[ 0.0 -0.5 0.0 0.5 ][ 2.0 ] [ 1.0 ]
The bar has 1.0 stiffness in x direction.
It is because stiffness E*A = 1.
On the other hand, the bar has 0.5 stiffness in y direction.
Can we know why?
"""
import getfem as gf
import numpy as np
m = gf.Mesh('empty', 2)
# set points
pts = np.array([[0, 1],
[0, 0]])
pids = m.add_point(pts)
# set convexes
gt = gf.GeoTrans('GT_PK(1, 1)')
cvids = m.add_convex(gt, pts)
# set regions
cvfids = m.faces_from_pid(pids) # cvfids: [[0, 0], [0, 1]]
DIRICHLET_BOUNDARY = 1
m.set_region(DIRICHLET_BOUNDARY, np.array([[0], [1]]))
NEUMANN_BOUNDARY = 2
m.set_region(NEUMANN_BOUNDARY, np.array([[0], [0]]))
# set fem im
mfu = gf.MeshFem(m, 2) # 2: vector on xy-plane
mfu.set_fem(gf.Fem('FEM_PK(1, 1)')) # 2-dof on one segment
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS1D(1)'))
# handle DIRICHLET_BOUNDARY
mfd = gf.MeshFem(m, 1) #1: single value
mfd.set_fem(gf.Fem('FEM_PK(1, 0)')) # 1-dof on one segment
Hs = mfd.eval('[[1, 0], [0, 1]]')
Rs = mfd.eval('[0, 0]')
H, R = gf.asm_dirichlet(DIRICHLET_BOUNDARY, mim, mfu, mfd, Hs, Rs)
N, U0 = H.dirichlet_nullspace(R)
# handle NEUMANN_BOUNDARY
Fs = np.repeat([[0.], [1.]], mfd.nbdof(), axis=1) # Fx=0, Fy=400
F = gf.asm_boundary_source(NEUMANN_BOUNDARY, mim, mfu, mfd, Fs)
# get K: stiffness matrix
E = 1.
Nu = 0. # no interaction of deformation in x- and y-direction
A = 1.
Lambda = A * E * Nu / ((1. + Nu) * (1. - 2. * Nu))
Mu = A * E / (2. * (1. + Nu))
K = gf.asm_linear_elasticity(mim, mfu, mfd,
np.repeat([Lambda], mfd.nbdof()),
np.repeat([Mu], mfd.nbdof()), -1)
# solve
Nt = gf.Spmat('copy', N)
Nt.transpose()
KK = Nt * K * N
FF = Nt * F
P = gf.Precond('ildlt', KK)
UU = gf.linsolve_cg(KK, FF, P)
U = N * UU + U0
# print KU = F
print('KU = F')
for i in range(len(U)):
print('[', end='')
for kj in K.full()[i]:
print(str(kj).center(6), end='')
print(']', end='')
print('[', str(U[i]).center(6), '] [', str(F[i]).center(6), ']')
