jfkiviaho
commented
1 month ago More importantly, it looks like something similar happens with tractions and tetrahedral elements. I recreated the analogous problem in 3D. And, again, applying the tractions to the "slanted" face of a reference tetrahedral element (the analogy of the hypotenuse for a reference triangular element) causes some unwanted scaling of the resultant nodal forces.
Here's the script I wrote.
from tacs import TACS, elements, constitutive
import numpy as np
from mpi4py import MPI
# Define two connectivities for a 12-element tetrahedral mesh of a cube. All
# that changes between the connectivities is the cell-to-vertex lists for the
# pair of tetrahedral cells that make up the right face of the cube (where the
# traction is applied). In the second connectivity, the idea (as before) is to
# permute these lists so, for each of these tetrahedral cells, the right-facing
# faces map to the larger face of the reference tetrahedral element (analogous
# to the hypotenuse of the reference triangle).
conn0 = np.array(
[
# bottom face
0, 1, 3, 8, # c0
0, 8, 3, 2, # c1
# top face
6, 7, 5, 8, # c2
6, 8, 5, 4, # c3
# front face
4, 5, 1, 8, # c4
4, 8, 1, 0, # c5
# back face
7, 6, 2, 8, # c6
7, 8, 2, 3, # c7
# left face
6, 4, 0, 8, # c8
6, 8, 0, 2, # c9
# right face
5, 7, 3, 8, # c10
5, 8, 3, 1 # c11
],
dtype=np.intc
)
conn1 = np.array(
[
# bottom face
0, 1, 3, 8, # c0
0, 8, 3, 2, # c1
# top face
6, 7, 5, 8, # c2
6, 8, 5, 4, # c3
# front face
4, 5, 1, 8, # c4
4, 8, 1, 0, # c5
# back face
7, 6, 2, 8, # c6
7, 8, 2, 3, # c7
# left face
6, 4, 0, 8, # c8
6, 8, 0, 2, # c9
# right face
8, 1, 3, 5, # c10
8, 3, 7, 5 # c11
],
dtype=np.intc
)
# The tractions are applied on elements c10 and c11. But the indices of the
# faces on which they are applied changes depending on the connectivity.
trac_face_indices0 = (0, 3) # in-plane faces
trac_face_indices1 = (2, 2) # "slanted" faces
# Define function that prints nodal forces given the connectivity and traction
# face index
comm = MPI.COMM_WORLD
def print_nodal_forces(conn, trac_face_indices):
props = constitutive.MaterialProperties(rho=1.0, E=1.0, nu=0.3, ys=1e2)
stiff = constitutive.SolidConstitutive(props, t=1.0, tNum=0)
model = elements.LinearElasticity3D(stiff)
basis = elements.LinearTetrahedralBasis()
elem = elements.Element3D(model, basis)
vars_per_node = model.getVarsPerNode()
creator = TACS.Creator(comm, vars_per_node)
num_nodes = 9
ptr = np.arange(0, 4*(12+1), 4, dtype=np.intc)
comp_ids = np.zeros(12, dtype=np.intc)
creator.setGlobalConnectivity(num_nodes, ptr, conn, comp_ids)
node_coords = np.array(
[
0.0, 0.0, 0.0, # v0
1.0, 0.0, 0.0, # v1
0.0, 1.0, 0.0, # v2
1.0, 1.0, 0.0, # v3
0.0, 0.0, 1.0, # v4
1.0, 0.0, 1.0, # v5
0.0, 1.0, 1.0, # v6
1.0, 1.0, 1.0, # v7
0.5, 0.5, 0.5 # v8
]
)
creator.setNodes(node_coords)
bcnodes = np.array([0, 2, 4, 6], dtype=np.intc)
bcdofs = np.array([0, 1, 2, 0, 1, 2, 0, 1, 0, 1], dtype=np.intc)
bcptr = np.array([0, 3, 6, 8, 10], dtype=np.intc)
creator.setBoundaryConditions(bcnodes, bcptr, bcdofs)
element_list = [elem]
creator.setElements(element_list)
assembler = creator.createTACS()
trac_vec = np.array([0.1, 0.0, 0.0], dtype=np.double)
aux_elems = TACS.AuxElements()
traction0 = elem.createElementTraction(trac_face_indices[0], trac_vec)
traction1 = elem.createElementTraction(trac_face_indices[1], trac_vec)
aux_elems.addElement(10, traction0)
aux_elems.addElement(11, traction1)
assembler.setAuxElements(aux_elems)
res = assembler.createVec()
ans = assembler.createVec()
mat = assembler.createSchurMat()
alpha = 1.0
beta = 0.0
gamma = 0.0
assembler.zeroVariables()
assembler.assembleJacobian(alpha, beta, gamma, res, mat)
res.scale(-1.0)
print('forces:')
print(res.getArray().reshape(-1,vars_per_node))
print('sum:')
print(np.sum(res.getArray()))
print()
node_vec = assembler.createNodeVec()
assembler.getNodes(node_vec)
print('nodes:')
print(node_vec.getArray().reshape(-1,3))
print()
# Print the nodal forces for the two difference cases (connectivities)
print('Case 1:')
print('-------')
print_nodal_forces(conn0, trac_face_indices0)
print('Case 2:')
print('-------')
print_nodal_forces(conn1, trac_face_indices1)And here are the results.
Case 1:
-------
forces:
[[-0. -0. -0. ]
[ 0.01666667 -0. -0. ]
[ 0.03333333 -0. -0. ]
[-0. -0. -0. ]
[-0. -0. -0. ]
[-0. -0. -0. ]
[ 0.01666667 -0. -0. ]
[ 0.03333333 -0. -0. ]
[-0. -0. -0. ]]
sum:
0.09999999999999999
nodes:
[[0. 0. 0. ]
[1. 0. 0. ]
[1. 1. 0. ]
[0.5 0.5 0.5]
[0. 1. 0. ]
[0. 1. 1. ]
[1. 1. 1. ]
[1. 0. 1. ]
[0. 0. 1. ]]
Case 2:
-------
forces:
[[-0.00000000e+00 -0.00000000e+00 -0.00000000e+00]
[ 2.88675135e-02 -0.00000000e+00 -0.00000000e+00]
[ 5.77350269e-02 -0.00000000e+00 -0.00000000e+00]
[ 1.44222201e-17 -0.00000000e+00 -0.00000000e+00]
[-0.00000000e+00 -0.00000000e+00 -0.00000000e+00]
[-0.00000000e+00 -0.00000000e+00 -0.00000000e+00]
[ 2.88675135e-02 -0.00000000e+00 -0.00000000e+00]
[ 5.77350269e-02 -0.00000000e+00 -0.00000000e+00]
[-0.00000000e+00 -0.00000000e+00 -0.00000000e+00]]
sum:
0.17320508075688776
nodes:
[[0. 0. 0. ]
[1. 0. 0. ]
[1. 1. 0. ]
[0.5 0.5 0.5]
[0. 1. 0. ]
[0. 1. 1. ]
[1. 1. 1. ]
[1. 0. 1. ]
[0. 0. 1. ]]Again, the net force is higher than it ought to be.
I was trying solve the simple 2D problem (see diagram) of a bar in pure extension when I encountered a problem with how the tractions are calculated.
The problem can be reproduced with the following script.
The script assumes that the right-hand side has a height of
1and the distributed load has value0.1. Thus, the nodal forces on the right-hand side should sum up to a net force of1 * 0.1 = 0.1. These are the results for the two different (but valid) connectivities defined in the script.In the second case, the net force is larger than
0.1. The forces appear to have been scaled bysqrt(2), the length of the hypotenuse of the reference triangle. The upshot is that, unless the mesh is prepared so that the tractions are only applied to the legs of the reference triangle element, TACS cannot correctly solve the simple problem 2D problem I illustrated above.