SOLID186 — 20-node quadratic hex solid#

Kinematics. 20-node serendipity hexahedron; three translational DOFs per node (\(u_x, u_y, u_z\)); 60 DOFs per element. Also dispatches degenerate wedge (15-node) and pyramid (13-node) shapes automatically when a CDB-stored element has collapsed corners.

Stiffness. \(K_e = \int_{\Omega_e} B^\top C B\,\mathrm{d}V\) integrated on the reference cube with 3×3×3 Gauss-Legendre (full) or 2×2×2 (reduced).

Mass. Consistent \(M_e = \int \rho\, N^\top N\,\mathrm{d}V\) with 3×3×3 Gauss.

MAPDL compatibility. Reproduces SOLID186 at KEYOPT(2) = 1 (full 3×3×3 integration) and KEYOPT(2) = 0 (uniform-reduced 2×2×2 without explicit hourglass control).

Shape functions (serendipity 20-node hex)#

The 20-node brick is the quadratic serendipity extension of the trilinear SOLID185 [ZTZ2013] [Bathe2014]. Eight corner nodes sit at \((\xi_i, \eta_i, \zeta_i) \in \{\pm 1\}^3\); twelve mid-edge nodes sit where exactly one natural coordinate is zero. There are no mid-face or body-centre nodes — that’s the hallmark of the serendipity family [Cook2001].

  • Corner node \(i\) (all three \(|\xi_i^c| = 1\)):

    \[N_i(\xi, \eta, \zeta) = \tfrac{1}{8}\,(1 + \xi_i \xi)\,(1 + \eta_i \eta)\, (1 + \zeta_i \zeta)\,\bigl(\xi_i \xi + \eta_i \eta + \zeta_i \zeta - 2\bigr).\]

  • Mid-edge node on a \(\xi\)-directed edge (\(\xi_i^m = 0\)):

    \[N_i(\xi, \eta, \zeta) = \tfrac{1}{4}\,(1 - \xi^2)\,(1 + \eta_i \eta)\, (1 + \zeta_i \zeta),\]

    with the symmetric forms for \(\eta\)- and \(\zeta\)-directed mid-edges.

Jacobian, \(B\) matrix, and the isotropic constitutive \(C\) are the same six-component Voigt objects as for SOLID185 — only the shape-function and derivative entries change.

3×3×3 Gauss quadrature#

Full integration uses 27 points at natural coordinates \(\xi_g \in \{-\sqrt{3/5},\, 0,\, +\sqrt{3/5}\}\) (and similarly for \(\eta\), \(\zeta\)) with weights \(w_g\) from the tensor product of the one-dimensional 3-point rule (\(5/9,\ 8/9,\ 5/9\)). The rule is exact for polynomials up to degree 5 in each direction, which integrates \(B^\top C B\) on an undistorted hex without approximation.

Reduced ("reduced") integration drops to the 2×2×2 rule (same points as SOLID185). On a single element this produces six extra zero-energy modes (hourglass modes); in a patch of elements the modes usually don’t connect and the global \(K\) stays stable. MAPDL’s KEYOPT(2) = 0 option applies explicit hourglass-control stabilisation that femorph-solver does not yet implement, so use "full" in any production run.

Numpy walk-through#

The following listing reuses the SOLID185 pattern with the 20-node shape functions. On a unit cube with generic steel it produces the expected 6 rigid-body + 54 elastic eigenvalues:

import numpy as np

# Natural coordinates of the 20 nodes, in MAPDL / VTK order:
# corners (bottom 4, top 4), bottom-face midsides, top-face
# midsides, vertical midsides.
xi_n = np.array([
    (-1,-1,-1), (+1,-1,-1), (+1,+1,-1), (-1,+1,-1),
    (-1,-1,+1), (+1,-1,+1), (+1,+1,+1), (-1,+1,+1),
    ( 0,-1,-1), (+1, 0,-1), ( 0,+1,-1), (-1, 0,-1),
    ( 0,-1,+1), (+1, 0,+1), ( 0,+1,+1), (-1, 0,+1),
    (-1,-1, 0), (+1,-1, 0), (+1,+1, 0), (-1,+1, 0),
], dtype=float)

# 3x3x3 Gauss-Legendre.
p1d = np.array([-np.sqrt(3/5), 0.0, +np.sqrt(3/5)])
w1d = np.array([5/9, 8/9, 5/9])
gp  = np.array([(a, b, c) for a in p1d for b in p1d for c in p1d])
w   = np.array([a * b * c for a in w1d for b in w1d for c in w1d])

def shape_and_derivs(xi, eta, zeta):
    N  = np.empty(20)
    dN = np.empty((20, 3))
    for i, (xi_i, eta_i, zeta_i) in enumerate(xi_n):
        one_xi   = 1.0 + xi_i * xi
        one_eta  = 1.0 + eta_i * eta
        one_zeta = 1.0 + zeta_i * zeta
        if abs(xi_i) + abs(eta_i) + abs(zeta_i) == 3:              # corner
            combo = xi_i * xi + eta_i * eta + zeta_i * zeta - 2.0
            N[i] = 0.125 * one_xi * one_eta * one_zeta * combo
            dN[i, 0] = 0.125 * one_eta * one_zeta * (xi_i * combo + one_xi * xi_i)
            dN[i, 1] = 0.125 * one_xi  * one_zeta * (eta_i * combo + one_eta * eta_i)
            dN[i, 2] = 0.125 * one_xi  * one_eta  * (zeta_i * combo + one_zeta * zeta_i)
        elif xi_i == 0:                                            # mid-edge on ξ
            q = 1.0 - xi * xi
            N[i] = 0.25 * q * one_eta * one_zeta
            dN[i, 0] = 0.25 * (-2.0 * xi) * one_eta * one_zeta
            dN[i, 1] = 0.25 * q * eta_i * one_zeta
            dN[i, 2] = 0.25 * q * one_eta * zeta_i
        elif eta_i == 0:                                           # mid-edge on η
            q = 1.0 - eta * eta
            N[i] = 0.25 * q * one_xi * one_zeta
            dN[i, 0] = 0.25 * q * xi_i * one_zeta
            dN[i, 1] = 0.25 * (-2.0 * eta) * one_xi * one_zeta
            dN[i, 2] = 0.25 * q * one_xi * zeta_i
        else:                                                      # mid-edge on ζ
            q = 1.0 - zeta * zeta
            N[i] = 0.25 * q * one_xi * one_eta
            dN[i, 0] = 0.25 * q * xi_i * one_eta
            dN[i, 1] = 0.25 * q * one_xi * eta_i
            dN[i, 2] = 0.25 * (-2.0 * zeta) * one_xi * one_eta
    return N, dN

# Assemble K_e the same way as the SOLID185 listing: loop over Gauss
# points, compute J = dN^T X, detJ, dN/dx = J^{-1} dN/dξ, build the
# 6x60 B, accumulate K += B^T C B detJ w.

The element passes the patch test to machine precision: a linear displacement field ramps exactly through every node of any distorted hex.

Degenerate shapes — wedge and pyramid#

MAPDL stores pyramid and wedge elements as 20-slot hex connectivity with repeated node IDs at collapsed corners:

  • Wedge (15 nodes) — corner collapse K = L and O = P on adjacent vertical edges.

  • Pyramid (13 nodes) — corner collapse M = N = O = P on all four top-face nodes (the apex).

  • Degenerate tet would collapse further, but in practice MAPDL writes those as SOLID187.

mapdl_archive parses the degenerate connectivity and emits VTK_QUADRATIC_WEDGE (15 pts) or VTK_QUADRATIC_PYRAMID (13 pts) cells. femorph-solver’s SOLID186W and SOLID186P wrappers restore the full 20-slot layout, delegate \(K_e\) / \(M_e\) to the plain SOLID186 kernel, and fold the 60 × 60 result back down to 45 × 45 (wedge) or 39 × 39 (pyramid) with a constant indicator matrix \(T\):

\[K_\text{deg} = T^\top K_\text{hex}\, T.\]

This is exactly what the global assembler would compute if the element record had 20 node IDs with repeats — bit-for-bit identical to MAPDL’s degenerate-hex formulation.

Validated on a 5 166-element rotor-blade CDB (3 069 regular hex + 497 pyramid + 19 wedge + 1 581 tet): 12-mode modal parity to 0.08 % relative vs MAPDL (the residual is the plain-Gauss vs B-bar offset, not the degenerate dispatch).

Formulation flag#

Name

MAPDL KEYOPT(2)

Description

"full"

1

3×3×3 Gauss (27 points). Fully ranked single-element \(K_e\).

"reduced"

0

2×2×2 Gauss (8 points). 6 hourglass modes per element. No stabilisation yet — caller’s responsibility at mesh level.

 
m.materials[1]["_SOLID186_FORMULATION"] = "reduced"

Validation#

Rigid-body modes. 54 non-zero elastic modes + 6 rigid-body zeros on an undistorted hex ("full"); 48 + 12 (6 rigid + 6 hourglass) for "reduced".

Degenerate-shape correctness. SOLID186W on a 15-node wedge and SOLID186P on a 13-node pyramid each have exactly 6 rigid-body zero eigenvalues on the reference element.

MAPDL parity. Plain modal on a quarter-arc blade sector (fixed root, 12 lowest modes) agrees with MAPDL to 0.08 % relative.

API reference#

class femorph_solver.elements.solid186.SOLID186[source]#

Bases: ElementBase

static eel_batch(coords: ndarray, u_e: ndarray) → ndarray | None[source]#

Per-element elastic strain at every element node.

coords: (n_elem, 20, 3); u_e: (n_elem, 60). Returns (n_elem, 20, 6) — engineering-shear Voigt strain.

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

  • coords ((n_nodes, 3) float64)

  • material (mapping with MAPDL property names as keys (EX, PRXY, DENS, ...))

  • real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#

Return the element mass matrix in global DOF ordering.

class femorph_solver.elements.solid186_degenerate.SOLID186W[source]#

Bases: ElementBase

15-node degenerate SOLID186 wedge (K=L, O=P collapse).

Dedicated quadratic serendipity kernel with shape functions in area coords (ξ₁, ξ₂, ζ), stiffness at 9-pt Gauss (3-pt triangle × 3-pt line), consistent mass at 21-pt Gauss (7-pt triangle × 3-pt line).

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

  • coords ((n_nodes, 3) float64)

  • material (mapping with MAPDL property names as keys (EX, PRXY, DENS, ...))

  • real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#

Return the element mass matrix in global DOF ordering.

class femorph_solver.elements.solid186_degenerate.SOLID186P[source]#

Bases: ElementBase

13-node degenerate SOLID186 pyramid.

Default: dedicated Bedrosian apex-singular serendipity kernel integrated with a 2×2×2 Duffy collapsed-hex Gauss rule. This is bit-exact (~1 × 10⁻¹¹ rel Frobenius) against MAPDL’s ptr_em.emat pyramid block on every measured element. The 2×2×2 rule (not 3×3×3 or higher) is what MAPDL actually uses: the Bedrosian shape functions integrated with a higher-order rule diverge from MAPDL’s output by ~13% because MAPDL deliberately under-integrates in the SOLID186 URI style for pyramids just like it does for hexes.

The hex-fold wrapper (T^T · K_20-node · T) is retained as material["_SOLID186P_FORMULATION"] = "hex_fold" for back-compat and regression testing.

static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

  • coords ((n_nodes, 3) float64)

  • material (mapping with MAPDL property names as keys (EX, PRXY, DENS, ...))

  • real (1-D array of real constants)

static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#

Return the element mass matrix in global DOF ordering.

References#

[ZTZ2013]

Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th ed. Butterworth-Heinemann. https://doi.org/10.1016/C2009-0-24909-9

[Bathe2014]

Bathe, K.-J. (2014). Finite Element Procedures, 2nd ed. Prentice-Hall / Klaus-Jürgen Bathe. Ch. 5 (isoparametric elements). https://www.klausjurgenbathe.com/fepbook/

[Cook2001]

Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J. (2001). Concepts and Applications of Finite Element Analysis, 4th ed. Wiley. ISBN 978-0-471-35605-9. https://www.wiley.com/en-us/Concepts+and+Applications+of+Finite+Element+Analysis%2C+4th+Edition-p-9780471356059