Applying the nodal fixation

[ywan459]

Hi, I am trying to setup a simple linear elastic model for my further algorithm development. The following snippet does the read a bar model, and I am trying to fix the bottom side nodes, and pull the top size nodes in Z direction. The results looks quite wrong that the deformation shows a lot of randomness. I guess my way of applying nodal fixation could be wrong, or did I miss something important? The mesh "elastic_bar.med" is attached here.

Here is the code:

from skfem import *
from skfem.models.elasticity import *
import numpy as np

# read in an elastic bar model in med file format, meshio supports this file format
m = Mesh3D.load("elastic_bar.med")

basis = Basis(m, ElementVector(ElementHex1()))

# get coo_data object in order to access element-wise stiffness matrices
data = linear_elasticity(*lame_parameters(1e9, 0.3)).coo_data(basis)
# For future application that manipulates the element stiffness matrix
data = data.fromlocal(data.tolocal()) 
# assemble general stiffness matrix 
K = data.tocsr()

# vertices index on top of elastic bar
group_btm = [17,51,85,119,153,187,221,255,34,68,102,136,170,204,238]
# vertices index on bottom of elastic bar
group_top = [1,35,69,103,137,171,205,239,18,52, 86,120,154,188,222]

u = basis.zeros()

# get the dofs for the top nodes
pulled = basis.get_dofs(group_top)
# dofs for the bottom nodes
fixed = basis.get_dofs(group_btm)

# pull the top nodes in Z direction by 0.001
u[pulled.nodal['u^3']] = 0.001

# solve the system, 
# D for eliminated dofs, in our case, it meant fixed nodes, the bottom nodes group_btm
u = solve(*condense(K, x=u, D=fixed.flatten()))

# update mesh nodes with displacements.
sf = 1.0
m = m.translated(sf * u[basis.nodal_dofs])

# section calculating the system energy
elem_mats = data.tolocal()
obj = 0

for elem_i, elem_i_mat in enumerate(elem_mats):
     u_idx = basis.element_dofs[:,elem_i]
     u_i = u[u_idx]
     obj = obj + u_i.dot(elem_i_mat.dot(u_i))
# compare the element-wise energy calculation with the overall system calculation
# they should be identical
print(u.dot(K.dot(u)))
print(obj)

if __name__ == "__main__":
    from os.path import splitext
    from sys import argv
    
    #m.save(splitext(argv[0])[0] + '.med')
    m.save(splitext(argv[0])[0] + '.vtk')

This is the original shape, and bottom nodes shall be fixed, and top nodes were pulled.

The deformed states showed a random deformation.

elastic_bar.zip

[kinnala]

The method basis.get_dofs requires a set of facets instead of a set of nodes because it is designed to work with general finite elements that have all nodal, facet and edge DOFs.

Here you can use something like

basis.nodal_dofs[2, group_top].flatten()

to get the DOF indices of z-directional nodal displacement and

basis.nodal_dofs[:, group_top].flatten()

to get all DOF indices matching group_top.

I think *.med might not support facet sets but if you decide to use some more general format like *.msh then Mesh3D.load can read also the names of the facet sets.

[ywan459]

The method basis.get_dofs requires a set of facets instead of a set of nodes because it is designed to work with general finite elements that have all nodal, facet and edge DOFs.

Here you can use something like

basis.nodal_dofs[2, group_top].flatten()

to get the DOF indices of z-directional nodal displacement and

basis.nodal_dofs[:, group_top].flatten()

to get all DOF indices matching group_top.

I think *.med might not support facet sets but if you decide to use some more general format like *.msh then Mesh3D.load can read also the names of the facet sets.

Thanks. I have changed my code as shown below. I found there is no deformation at all. The internal energy is 0.0, and result is exactly same as original state.

from skfem import *
from skfem.io import from_meshio
from skfem.models.elasticity import *
import numpy as np

import meshio

# read in an elastic bar model in med file format, meshio supports this file format
m = Mesh3D.load("elastic_bar.med")

#mio = meshio.read("elastic_bar.med")
#m = from_meshio(mio)

basis = Basis(m, ElementVector(ElementHex1()))

# get coo_data object in order to access element-wise stiffness matrices
data = linear_elasticity(*lame_parameters(1e9, 0.3)).coo_data(basis)
# For future application that manipulates the element stiffness matrix
data = data.fromlocal(data.tolocal()) 
# assemble general stiffness matrix 
K = data.tocsr()

# vertices index on top of elastic bar
group_btm = np.array([17,51,85,119,153,187,221,255,34,68,102,136,170,204,238])-1  
# vertices index on bottom of elastic bar
group_top = np.array([1,35,69,103,137,171,205,239,18,52, 86,120,154,188,222])-1 

u = basis.zeros()

# fix the bottom nodes dofs
fixed = basis.dofs.nodal_dofs[:, group_btm]
# only get the dofs in Z direction
pulled = basis.dofs.nodal_dofs[2, group_top]
# top nodes move 0.01 in Z direction
u[pulled.flatten()] = 0.01

# solve the system, 
# D for eliminated dofs, in our case, it meant fixed nodes, the bottom nodes group_btm
u = solve(*condense(K, x=u, D=fixed.flatten()))

# update mesh nodes with displacements.
sf = 1.0
m = m.translated(sf * u[basis.nodal_dofs])

# section calculating the system energy
elem_mats = data.tolocal()
obj = 0

for elem_i, elem_i_mat in enumerate(elem_mats):
     u_idx = basis.element_dofs[:,elem_i]
     u_i = u[u_idx]
     obj = obj + u_i.dot(elem_i_mat.dot(u_i))
# compare the element-wise energy calculation with the overall system calculation
# they should be identical
print(u.dot(K.dot(u)))
print(obj)

if __name__ == "__main__":
    from os.path import splitext
    from sys import argv
    
    #m.save(splitext(argv[0])[0] + '.med')
    m.save(splitext(argv[0])[0] + '.vtk')

[kinnala]

In your code D=fixed.flatten() sets the Dirichlet nodes as only the fixed nodes so there is no load being applied. You must include all Dirichlet nodes, e.g., D=np.concatenate((fixed.flatten(), pulled.flatten())).

[ywan459]

In your code D=fixed.flatten() sets the Dirichlet nodes as only the fixed nodes so there is no load being applied. You must include all Dirichlet nodes, e.g., D=np.concatenate((fixed.flatten(), pulled.flatten())).

Yes, that's the real root cause, I did not fully understand the D parameter. Now it works. Thank you very much! I will keep working on the algorithm.