def shape_function_hex8(element_node, r, s, t):

  # r,s,t are each in the range [-1,1] with 0,0,0 in the middle of the
  # element. They are in the same directions as the R,S,T used in
  # Mesh Tools -> Refine -> Custom.

  if element_node == 1:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * -1) * (1 + t * -1)
  if element_node == 2:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * -1) * (1 + t * -1)
  if element_node == 3:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * 1) * (1 + t * -1)
  if element_node == 4:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * 1) * (1 + t * -1)
  if element_node == 5:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * -1) * (1 + t * 1)
  if element_node == 6:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * -1) * (1 + t * 1)
  if element_node == 7:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * 1) * (1 + t * 1)
  if element_node == 8:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * 1) * (1 + t * 1)


def interpolate(element_id, variable, r, s, t):
  # Only for hex8 elements.

  sum = 0
  nodes = solution.nodes(element_id)
  for element_node in range(1, len(nodes) + 1):
    shape_function = shape_function_hex8(element_node, r, s, t)
    node_id = nodes[element_node - 1]
    value = solution.node_value(variable, node_id)
    sum += shape_function * value

  return sum


assert mw.version() >= 15, "Wrong Mecway version"

# Enter parameters here
value = interpolate(element_id = 1,
                    variable = "stressxx",
                    r = 0,
                    s = -1,
                    t = 0)

mw.message(str(value))