Error in traction calculation for simplicial elements

[jfkiviaho]

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.

from tacs import TACS, elements, constitutive, functions

import numpy as np
from mpi4py import MPI
'''
Mesh used in the demonstration:


    v2                 v3
       + • • • • • • +   
       • •           •
       •   •         •
       •     •   f1  •
       •  f0   •     •
       •         •   •
       •           • •
       + • • • • • • +   
    v0                 v1              


Reference triangle element:

      v2
       +
       • •
       •   • e1
    e2 •     •
       •       •
       + • • • • +
      v0   e0    v1
'''
# Define two connectivities, one in which the right-facing edge of element f1
# maps to a leg of the reference triangle (e0 or e2) and one in which it maps
# to the hypotenuse of the reference triangle (e1). Note: the only difference
# in the connectivities is that the face-vertex list for element f1 is
# cyclically permuted. The relative orientation of the elements is preserved.
# So orientation is not at issue.
conn0 = np.array(
    [
        0, 1, 2,  # f0
        1, 3, 2   # f1
    ],
    dtype=np.intc
)

conn1 = np.array(
    [
        0, 1, 2,  # f0
        2, 1, 3   # f1
    ],
    dtype=np.intc
)


# The traction is applied on element f2. But the index of the face on which it
# is applied changes depending on the connectivity.
trac_face_index0 = 0
trac_face_index1 = 1


# Define function that prints nodal forces given the connectivity and traction
# face index
comm = MPI.COMM_WORLD
def print_nodal_forces(conn, trac_face_index):
    #props = constitutive.MaterialProperties(rho=1.0, E=1.0, nu=0.3, ys=1.0e2)
    props = constitutive.MaterialProperties(rho=2570.0, E=1.0, nu=0.3, ys=350e6)
    stiff = constitutive.PlaneStressConstitutive(props)
    model = elements.LinearElasticity2D(stiff)
    basis = elements.LinearTriangleBasis()
    elem = elements.Element2D(model, basis)

    vars_per_node = model.getVarsPerNode()
    creator = TACS.Creator(comm, vars_per_node)

    num_nodes = 4
    ptr = np.array([0, 3, 6], dtype=np.intc)
    comp_ids = np.zeros(2, 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
        ]
    )
    creator.setNodes(node_coords)

    bcnodes = np.array([0, 2], dtype=np.intc)
    bcptr = np.array([0, 2, 3], dtype=np.intc)
    bcdofs = np.array([0, 1, 0], 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], dtype=np.double)
    aux_elems = TACS.AuxElements()
    traction = elem.createElementTraction(trac_face_index, trac_vec)
    aux_elems.addElement(1, traction)
    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()

    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_nodal_forces(conn0, trac_face_index0)
print_nodal_forces(conn1, trac_face_index1)

The script assumes that the right-hand side has a height of 1 and the distributed load has value 0.1. Thus, the nodal forces on the right-hand side should sum up to a net force of 1 * 0.1 = 0.1. These are the results for the two different (but valid) connectivities defined in the script.

Case 1
------
forces:
[[-0.   -0.  ]
 [ 0.05 -0.  ]
 [-0.   -0.  ]
 [ 0.05 -0.  ]]

nodes:
[[0. 0. 0.]
 [1. 0. 0.]
 [0. 1. 0.]
 [1. 1. 0.]]


Case 2
------
forces:
[[-0.         -0.        ]
 [ 0.07071068 -0.        ]
 [-0.         -0.        ]
 [ 0.07071068 -0.        ]]

nodes:
[[0. 0. 0.]
 [1. 0. 0.]
 [0. 1. 0.]
 [1. 1. 0.]]

In the second case, the net force is larger than 0.1. The forces appear to have been scaled by sqrt(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.

[jfkiviaho]

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.