From 90a0b41da9dfac3695e67c4066e58186510ff9ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Fri, 2 Oct 2026 07:01:43 -0300 Subject: [PATCH 1/6] TenNodeTetrahedron: correct -doInitDisp, responses, parameters; add verification - -doInitDisp is enabled by a non-zero value (the page said 0); show the dash consistently; it may follow $matTag or the optional body forces. - Replace the undefined |bool| substitution with |integer|. - Document stresses/strains layout (4 GPs x 6) and per-point `material $gp` parameters. - Note the fixed volume-scaling, response-buffer and parser bugs, and add a downloadable verification script (body-force and uniform-strain patch tests, per-point parameter update) with expected output. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../model/elements/TenNodeTetrahedron.rst | 91 ++++++++++-- .../example_tet10_patch.tcl | 133 ++++++++++++++++++ 2 files changed, 215 insertions(+), 9 deletions(-) create mode 100644 source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_patch.tcl diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index aced5071..bdc4d0d8 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -4,24 +4,27 @@ TenNodeTetrahedron Element ^^^^^^^^^^^^^^^^^^^^^^^^^^ -This command is used to construct an ten-node tetrahedron element object, which uses the standard isoparametric formulation. +This command is used to construct a ten-node tetrahedron element object, which uses the standard isoparametric formulation. .. admonition:: Command - **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> ** + **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value>** -.. csv-table:: +.. csv-table:: :header: "Argument", "Type", "Description" :widths: 10, 10, 40 $eleTag, |integer|, unique element object tag $node1 .. $node10, 10 |integer|, nodes of tet (ordered as shown in fig below) $matTag, |integer|, tag of nDMaterial - $b1 $b2 $b3, |listFloat|, optional: body forces in global x y z directions - <-doInitDisp $value>, |bool|, optional: consider initial displacements if $value is 0 + $b1 $b2 $b3, |listFloat|, optional: body forces in global x y z directions (default: 0.0 0.0 0.0) + -doInitDisp $value, |integer|, "optional flag; if $value is non-zero, the nodal displacements present when the element is added to the domain are stored and subtracted from all later trial displacements before strains are computed (default: 0, i.e. disabled, so ordinary total displacements are used)" -This element is based on second-order interpolation of nodal quantities, this means that the strain and stress field inside the element are linearly interpolated. Four Gauss-points inside the element are used for integration. +This element is based on second-order interpolation of nodal quantities, so the strain and stress field inside the element vary linearly over the element. Four Gauss points inside the element are used for integration. +.. note:: + + ``-doInitDisp`` may follow ``$matTag`` directly or any of the optional body forces ``$b1 $b2 $b3`` (earlier versions required all three body forces to be given first; otherwise the command failed with "WARNING: invalid double data"). .. figure:: figures/TenNodeTetrahedron/TenNodeTetrahedron.png :align: center @@ -29,13 +32,81 @@ This element is based on second-order interpolation of nodal quantities, this me TenNodeTetrahedron Element Node Numbering +Theory +"""""" + +The element uses the full quadratic (10-node) interpolation of the tetrahedron: the 4 corner nodes plus the 6 edge mid-side nodes shown above, giving a displacement field that is quadratic in the natural (barycentric) coordinates :math:`(\xi,\eta,\zeta,\lambda)`, with :math:`\lambda = 1-\xi-\eta-\zeta`, and a strain/stress field that varies linearly over the element. + +Integration uses the classical symmetric 4-point rule for tetrahedra, exact for polynomials up to degree 2 (sufficient for the quadratic displacement field). The four sample points are the permutations of the barycentric coordinates + +.. math:: + + \alpha = \frac{5+3\sqrt{5}}{20} \approx 0.58541020,\qquad + \beta = \frac{5-\sqrt{5}}{20} \approx 0.13819660 + +i.e. point :math:`k` has one barycentric coordinate equal to :math:`\alpha` and the other three equal to :math:`\beta`, each with the same natural weight :math:`w = 1/24` (so that :math:`4w = 1/6`, the volume of the reference tetrahedron in natural coordinates). For an element with straight edges, the Jacobian determinant is constant over the element and equal to :math:`6V`, where :math:`V` is the physical element volume, so that + +.. math:: + + \int_V f\, dV = \sum_{k=1}^{4} w\, f(\mathbf{x}_k)\, \det(J) = \sum_{k=1}^{4} \frac{1}{24} f(\mathbf{x}_k)\, (6V) + +reproduces the exact element volume when :math:`f \equiv 1` (see, e.g., Keast, "Moderate degree tetrahedral quadrature formulas," *Computer Methods in Applied Mechanics and Engineering*, 55(3), 1986, for this family of quadrature rules). + +Responses and parameters +"""""""""""""""""""""""" + +Valid queries to a ``TenNodeTetrahedron`` element (e.g. via :ref:`ElementRecorder ` or the ``eleResponse`` command) are: + +.. csv-table:: + :header: "Response", "Size", "Description" + :widths: 15, 10, 40 + + forces (or force), 30, "nodal resisting forces, 3 components per node, nodes 1 to 10 in order" + stresses, 24, ":math:`(\sigma_{11},\sigma_{22},\sigma_{33},\sigma_{12},\sigma_{23},\sigma_{13})` at each of the 4 Gauss points, concatenated gp 1..4" + strains, 24, ":math:`(\epsilon_{11},\epsilon_{22},\epsilon_{33},\epsilon_{12},\epsilon_{23},\epsilon_{13})` at each of the 4 Gauss points, concatenated gp 1..4" + "material $gp arg1 ...", varies, "forwarded to ``setResponse`` of the NDMaterial copy at Gauss point $gp (1 to 4); ``integrPoint $gp ...`` is accepted as a synonym" + +Each of the 4 Gauss points holds its own independent copy of the assigned ``nDMaterial`` (made with ``getCopy`` when the element is added to the domain). ``setParameter``/``updateParameter`` can target one copy or all of them: + +.. csv-table:: + :header: "Parameter", "Description" + :widths: 25, 50 + + "material $gp arg1 ...", "forwarded to ``setParameter`` of the material copy at Gauss point $gp only (1 to 4)" + "arg1 ... (no ``material`` keyword)", "forwarded to ``setParameter`` of all 4 material copies, so the same change is applied at every Gauss point" + .. note:: - The valid queries to a `TenNodeTetrahedron` element when creating an ElementRecorder object are 'forces', 'stresses,' ('strains' version > 2.2.0) and 'material $matNum matArg1 matArg2 ...' Where $matNum refers to the material object at the integration point corresponding to the node numbers in the isoparametric domain. + Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, and ``-doInitDisp`` could not be used without body forces. All three are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#XXXX `_); the example below reproduces and checks the fix. This element can only be defined after a :ref:`model` with **-ndm 3 -ndf 3** -.. admonition:: Example +Verification +"""""""""""" + +:download:`example_tet10_patch.tcl ` builds a single, deliberately distorted ``TenNodeTetrahedron`` with corners at :math:`(0,0,0)`, :math:`(2,0,0)`, :math:`(0,3,0)`, :math:`(0,0,1.5)` (volume :math:`V = 2 \cdot 3 \cdot 1.5/6 = 1.5`) and runs three checks: + +#. **Body-force patch test.** With a uniform body force :math:`b_3=-2.0` and all nodes restrained, the sum of the vertical reactions over all 10 nodes must equal :math:`-b_3 V = 3.0` exactly (any underintegration of the volume shows up directly here). The script prints the summed reaction and compares it to the expected value. +#. **Uniform-strain patch test.** Prescribing :math:`\epsilon_{xx}=10^{-3}` on all 10 nodes (linear displacement BC) must reproduce the exact uniform stress :math:`(\sigma_{11},\sigma_{22},\sigma_{33})=(1.2,\,0.4,\,0.4)` (with :math:`E=1000`, :math:`\nu=0.25`) identically at all 4 Gauss points; this also exercises the ``stresses`` element response (24 values) without heap corruption. +#. **Per-Gauss-point parameter targeting.** For each Gauss point in turn, ``parameter ... element 1 material $gp E`` followed by ``updateParameter`` doubles the stiffness at that point only; the stress at that Gauss point must double while the other three stay unchanged. + +Running it with the OpenSees interpreter gives, on the fixed element: + +.. code-block:: text + + BODYFORCE sum Rz = 3.0000000000 expected = 3.0000000000 + STRESSES (expect per GP: 1.2 0.4 0.4 0 0 0) + PATCH max abs error = 3.802e-11 + SETPARAMETER per-Gauss-point E override (E: 1000.0 -> 2000.0, scale = 2.0) + gp=1 sigma_xx before=1.200000 after=2.400000 (expected 2.400000) other-points-unchanged=yes + gp=2 sigma_xx before=1.200000 after=2.400000 (expected 2.400000) other-points-unchanged=yes + gp=3 sigma_xx before=1.200000 after=2.400000 (expected 2.400000) other-points-unchanged=yes + gp=4 sigma_xx before=1.200000 after=2.400000 (expected 2.400000) other-points-unchanged=yes + OVERALL: PASS + +On the unfixed (vanilla upstream) element the same model instead gives a body-force reaction sum of ``0.5`` (should be ``3.0``) and a work-of-reactions under the imposed :math:`\epsilon_{xx}=10^{-3}` of ``3.0e-4`` (should be ``1.8e-3``), both exactly a factor of 6 too small, confirming the volume-scaling bug above. + +.. admonition:: Example The following example constructs a TenNodeTetrahedron element with tag **1** between nodes **1, 2, 3, 4, 5, 6, 7, 8, 9, 10** with an nDMaterial of tag **1** and body forces given by varaiables **b1, b2, b3**. @@ -51,4 +122,6 @@ This element is based on second-order interpolation of nodal quantities, this me element('TenNodeTetrahedron',1, 1,2,3,4,5,6,7,8,9,10, 1, b1, b2, b3) -Code Developed by: `José Antonio Abell `_ and José Luis Larenas (UANDES). For bugs and features, start a new issue on the `OpenSees github repo `_ and tag me (@jaabell). +This command is registered in both the Tcl and the (Python-backing) generic OpenSees interpreter. + +Code Developed by: `José Antonio Abell `_ and José Luis Larenas (UANDES). For bugs and features, start a new issue on the `OpenSees github repo `_ and tag me (@jaabell). diff --git a/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_patch.tcl b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_patch.tcl new file mode 100644 index 00000000..041765ee --- /dev/null +++ b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_patch.tcl @@ -0,0 +1,133 @@ +# TenNodeTetrahedron (quadratic 10-node tet, 4-point integration) checks. +# +# Fixes exercised (branch fix/tet10-integration vs upstream/master): +# 1. Jacobian/weight scaling of the 4-point rule (commit 2f3362a28) +# 2. setParameter/updateParameter reaching all 4 material points (a9430a76d) +# 3. heap overflow in eleResponse stresses|strains (110bb4dcc) +# +# Corners: (0,0,0) (2,0,0) (0,3,0) (0,0,1.5) -> V = 2*3*1.5/6 = 1.5 +# Midside ordering (from shape functions): 5:1-2 6:2-3 7:3-1 8:1-4 9:3-4 10:2-4 +proc build {b3} { + wipe + model BasicBuilder -ndm 3 -ndf 3 + set X {{0 0 0} {2 0 0} {0 3 0} {0 0 1.5}} + set mids {{1 2} {2 3} {3 1} {1 4} {3 4} {2 4}} + for {set i 0} {$i < 4} {incr i} { eval node [expr $i+1] [lindex $X $i] } + set n 5 + foreach e $mids { + set a [lindex $X [expr [lindex $e 0]-1]]; set b [lindex $X [expr [lindex $e 1]-1]] + node $n [expr ([lindex $a 0]+[lindex $b 0])/2.] [expr ([lindex $a 1]+[lindex $b 1])/2.] [expr ([lindex $a 2]+[lindex $b 2])/2.] + incr n + } + nDMaterial ElasticIsotropic 1 1000.0 0.25 + element TenNodeTetrahedron 1 1 2 3 4 5 6 7 8 9 10 1 0.0 0.0 $b3 +} + +proc applyUniformStrainBC {} { + timeSeries Linear 2; pattern Plain 2 2 { + for {set i 1} {$i <= 10} {incr i} { + set x [lindex [nodeCoord $i] 0] + sp $i 1 [expr 1e-3*$x]; sp $i 2 0.0; sp $i 3 0.0 + } + } + constraints Penalty 1e14 1e14; numberer Plain; system FullGeneral + test NormUnbalance 1e-10 5; algorithm Linear + integrator LoadControl 1.0; analysis Static + analyze 1 +} + +set failed 0 + +# ------------------------------------------------------------------ +# 1) Body-force patch test (fix #1: Jacobian/weight scaling). +# Net reaction over all 10 nodes must equal -b3*V exactly: the +# 4-pt rule must reconstruct the true element volume. +# ------------------------------------------------------------------ +build -2.0 +timeSeries Constant 1; pattern Plain 1 1 { + for {set i 1} {$i <= 10} {incr i} { sp $i 1 0.0; sp $i 2 0.0; sp $i 3 0.0 } +} +constraints Penalty 1e14 1e14; numberer Plain; system FullGeneral +test NormUnbalance 1e-10 5; algorithm Linear +integrator LoadControl 1.0; analysis Static +analyze 1 +reactions +set Rz 0.0 +for {set i 1} {$i <= 10} {incr i} { set Rz [expr $Rz + [nodeReaction $i 3]] } +set expRz 3.0 +puts [format "BODYFORCE sum Rz = %.10f expected = %.10f" $Rz $expRz] +if {abs($Rz-$expRz) > 1e-6} { puts " -> FAIL"; set failed 1 } + +# ------------------------------------------------------------------ +# 2) Uniform-strain patch test (fix #1 again, and fix #3: eleResponse +# "stresses" used to overflow a size-6 static Vector with 24 +# doubles -- this call used to corrupt the heap / crash). +# eps_xx = 1e-3 imposed on all nodes -> exact, identical stresses +# at all 4 Gauss points. +# ------------------------------------------------------------------ +build 0.0 +applyUniformStrainBC +set s [eleResponse 1 stresses] +set exp {1.2 0.4 0.4 0 0 0} +puts "STRESSES (expect per GP: $exp)" +puts $s +set maxerr 0.0 +for {set g 0} {$g < 4} {incr g} { + for {set k 0} {$k < 6} {incr k} { + set e [expr abs([lindex $s [expr 6*$g+$k]] - [lindex $exp $k])] + if {$e > $maxerr} { set maxerr $e } + } +} +puts [format "PATCH max abs error = %.3e" $maxerr] +if {$maxerr > 1e-6} { puts " -> FAIL"; set failed 1 } + +# ------------------------------------------------------------------ +# 3) setParameter/updateParameter reaching every material point +# (fix #2: both the per-point branch and the "all points" loop +# in TenNodeTetrahedron::setParameter stopped at point 1 upstream). +# +# Each Gauss point holds an *independent copy* of the nDMaterial +# (TenNodeTetrahedron::setDomain -> theMaterial.getCopy(...)), so +# raising E only at Gauss point $gp must rescale the stress at +# that point only (strain is fixed by the displacement BC, so for +# an isotropic elastic material: sigma = E-dependent stiffness +# tensor * fixed strain -> doubling E doubles sigma at that point). +# ------------------------------------------------------------------ +set E0 1000.0 +set E1 2000.0 +set scale [expr $E1/$E0] +puts "SETPARAMETER per-Gauss-point E override (E: $E0 -> $E1, scale = $scale)" +for {set gp 1} {$gp <= 4} {incr gp} { + build 0.0 + applyUniformStrainBC + set before [eleResponse 1 stresses] + + set ptag [expr 100+$gp] + parameter $ptag element 1 material $gp E + updateParameter $ptag $E1 + + set after [eleResponse 1 stresses] + + set ok 1 + for {set g 0} {$g < 4} {incr g} { + set expScale [expr {$g == $gp-1 ? $scale : 1.0}] + for {set k 0} {$k < 6} {incr k} { + set b [lindex $before [expr 6*$g+$k]] + set a [lindex $after [expr 6*$g+$k]] + set e [expr abs($a - $expScale*$b)] + if {$e > 1e-6} { set ok 0 } + } + } + set targetBefore [lrange $before [expr 6*($gp-1)] [expr 6*($gp-1)+5]] + set targetAfter [lrange $after [expr 6*($gp-1)] [expr 6*($gp-1)+5]] + puts [format " gp=%d sigma_xx before=%.6f after=%.6f (expected %.6f) other-points-unchanged=%s" \ + $gp [lindex $targetBefore 0] [lindex $targetAfter 0] [expr $scale*[lindex $targetBefore 0]] \ + [expr {$ok ? "yes" : "NO"}]] + if {!$ok} { puts " -> FAIL (gp=$gp)"; set failed 1 } +} + +if {$failed} { + puts "OVERALL: FAIL" +} else { + puts "OVERALL: PASS" +} From 3b90ca4c0b2931137d99c61dc9e25c281261bb64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Sat, 10 Oct 2026 11:30:45 -0300 Subject: [PATCH 2/6] TenNodeTetrahedron: link the code fix PR (#1833) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- source/user/manual/model/elements/TenNodeTetrahedron.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index bdc4d0d8..c7d13a16 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -77,7 +77,7 @@ Each of the 4 Gauss points holds its own independent copy of the assigned ``nDMa .. note:: - Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, and ``-doInitDisp`` could not be used without body forces. All three are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#XXXX `_); the example below reproduces and checks the fix. + Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, and ``-doInitDisp`` could not be used without body forces. All three are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. This element can only be defined after a :ref:`model` with **-ndm 3 -ndf 3** From d7556706690c525df8d4078932460bba817a885f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Mon, 12 Oct 2026 22:16:16 -0300 Subject: [PATCH 3/6] TenNodeTetrahedron: document -numbering gmsh|abaqus Add the node-numbering table (gmsh default vs Abaqus C3D10 for mid-edge nodes 9/10), note that the figure shows the Abaqus positions, and add a verified example comparing both orders on a curved element. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../model/elements/TenNodeTetrahedron.rst | 62 ++++++++++++++++- .../example_tet10_numbering.tcl | 68 +++++++++++++++++++ 2 files changed, 127 insertions(+), 3 deletions(-) create mode 100644 source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index c7d13a16..65790eb7 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -8,17 +8,18 @@ This command is used to construct a ten-node tetrahedron element object, which u .. admonition:: Command - **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value>** + **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value> <-numbering $order>** .. csv-table:: :header: "Argument", "Type", "Description" :widths: 10, 10, 40 $eleTag, |integer|, unique element object tag - $node1 .. $node10, 10 |integer|, nodes of tet (ordered as shown in fig below) + $node1 .. $node10, 10 |integer|, "nodes of tet: corners 1-4, then mid-edge nodes 5-10 in the order selected with ``-numbering`` (see :ref:`TenNodeTetrahedronNumbering` below)" $matTag, |integer|, tag of nDMaterial $b1 $b2 $b3, |listFloat|, optional: body forces in global x y z directions (default: 0.0 0.0 0.0) -doInitDisp $value, |integer|, "optional flag; if $value is non-zero, the nodal displacements present when the element is added to the domain are stored and subtracted from all later trial displacements before strains are computed (default: 0, i.e. disabled, so ordinary total displacements are used)" + -numbering $order, |string|, "optional: ``gmsh`` (default) or ``abaqus``, the convention used for mid-edge nodes 9 and 10 (see below)" This element is based on second-order interpolation of nodal quantities, so the strain and stress field inside the element vary linearly over the element. Four Gauss points inside the element are used for integration. @@ -30,7 +31,36 @@ This element is based on second-order interpolation of nodal quantities, so the :align: center :figclass: align-center - TenNodeTetrahedron Element Node Numbering + TenNodeTetrahedron Element Node Numbering (as drawn: ``-numbering abaqus``; with the default ``gmsh`` numbering nodes 9 and 10 are swapped) + +.. _TenNodeTetrahedronNumbering: + +Node numbering +"""""""""""""" + +Corner nodes 1-4 and mid-edge nodes 5-8 are the same in both conventions. Only nodes 9 and 10 differ: + +.. csv-table:: + :header: "Node", "``gmsh`` (default)", "``abaqus`` (Abaqus C3D10, VTK quadratic tetra)" + :widths: 10, 20, 30 + + 5, edge 1-2, edge 1-2 + 6, edge 2-3, edge 2-3 + 7, edge 3-1, edge 3-1 + 8, edge 1-4, edge 1-4 + 9, edge 3-4, edge 2-4 + 10, edge 2-4, edge 3-4 + +The default matches gmsh's second-order tetrahedron, so meshes exported from gmsh can be used as they are. +For meshes in Abaqus (``C3D10``) order, add ``-numbering abaqus``. The element keeps the nodes in the order +given, so the ``forces`` response, the stiffness and mass matrices and the element connectivity seen by +recorders all follow the chosen convention. Giving Abaqus-ordered connectivity without the flag puts nodes 9 and 10 +on the wrong edges and produces wrong results without any error message. + +.. note:: + + The figure above has always shown the Abaqus position of nodes 9 and 10, while the element has always used + the gmsh one by default. Models written from the figure need ``-numbering abaqus`` (or nodes 9 and 10 swapped). Theory """""" @@ -104,6 +134,22 @@ Running it with the OpenSees interpreter gives, on the fixed element: gp=4 sigma_xx before=1.200000 after=2.400000 (expected 2.400000) other-points-unchanged=yes OVERALL: PASS +:download:`example_tet10_numbering.tcl ` checks the two +numberings against each other. It uses a curved element (mid-edge nodes moved off the edges) under an imposed +quadratic displacement field and a body force. It is defined once in gmsh order and once in Abaqus order +(nodes 9 and 10 swapped, ``-numbering abaqus``). The nodal forces must agree after swapping nodes 9 and 10, and the +Gauss-point stresses must be identical. As a control, the Abaqus connectivity without the flag must disagree. It +then repeats the body-force and uniform-strain patch tests in Abaqus order: + +.. code-block:: text + + max |force| = 1.4612e+00 + ABAQUS forces vs gmsh (9<->10 swapped): max diff = 2.220e-16 + ABAQUS stresses vs gmsh: max diff = 1.110e-16 + control: Abaqus conn WITHOUT flag: max diff = 4.009e+00 (must be large) + ABAQUS BODYFORCE sum Rz = 3.0000000000 expected 3.0 + ABAQUS PATCH max abs error = 3.802e-11 + On the unfixed (vanilla upstream) element the same model instead gives a body-force reaction sum of ``0.5`` (should be ``3.0``) and a work-of-reactions under the imposed :math:`\epsilon_{xx}=10^{-3}` of ``3.0e-4`` (should be ``1.8e-3``), both exactly a factor of 6 too small, confirming the volume-scaling bug above. .. admonition:: Example @@ -122,6 +168,16 @@ On the unfixed (vanilla upstream) element the same model instead gives a body-fo element('TenNodeTetrahedron',1, 1,2,3,4,5,6,7,8,9,10, 1, b1, b2, b3) + 3. **Abaqus (C3D10) node order** + + .. code-block:: tcl + + element TenNodeTetrahedron 1 1 2 3 4 5 6 7 8 9 10 1 -numbering abaqus + + .. code-block:: python + + element('TenNodeTetrahedron', 1, 1,2,3,4,5,6,7,8,9,10, 1, '-numbering', 'abaqus') + This command is registered in both the Tcl and the (Python-backing) generic OpenSees interpreter. Code Developed by: `José Antonio Abell `_ and José Luis Larenas (UANDES). For bugs and features, start a new issue on the `OpenSees github repo `_ and tag me (@jaabell). diff --git a/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl new file mode 100644 index 00000000..ad5add20 --- /dev/null +++ b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl @@ -0,0 +1,68 @@ +# TenNodeTetrahedron: gmsh vs Abaqus (-numbering abaqus) node order. +# Node tags: 1-4 corners, 5:1-2 6:2-3 7:3-1 8:1-4, 9 on edge 3-4, 10 on edge 2-4. +# gmsh connectivity: 1..10. Abaqus C3D10 connectivity: 1 2 3 4 5 6 7 8 10 9. +# Midside nodes are moved off the edges (curved element) so a wrong 9/10 assignment shows. +proc build {conn opts b3} { + wipe; model BasicBuilder -ndm 3 -ndf 3 + set X {{0 0 0} {2 0 0} {0 3 0} {0 0 1.5}} + set mids {{1 2} {2 3} {3 1} {1 4} {3 4} {2 4}} + set bump {{0.05 -0.1 0.02} {0.1 0.05 -0.03} {-0.08 0.04 0.06} {0.03 0.07 -0.05} {-0.06 0.09 0.04} {0.07 -0.04 0.08}} + for {set i 0} {$i < 4} {incr i} { eval node [expr $i+1] [lindex $X $i] } + set n 5 + foreach e $mids d $bump { + set a [lindex $X [expr [lindex $e 0]-1]]; set b [lindex $X [expr [lindex $e 1]-1]] + set c {} + foreach k {0 1 2} { lappend c [expr ([lindex $a $k]+[lindex $b $k])/2. + [lindex $d $k]] } + eval node $n $c; incr n + } + nDMaterial ElasticIsotropic 1 1000.0 0.25 + eval element TenNodeTetrahedron 1 $conn 1 0.0 0.0 $b3 $opts +} +# Same quadratic displacement field and body force imposed on each element; results are compared between orderings. +proc run {conn opts} { + build $conn $opts -2.0 + timeSeries Linear 1; pattern Plain 1 1 { + for {set i 1} {$i <= 10} {incr i} { + lassign [nodeCoord $i] x y z + sp $i 1 [expr 1e-3*($x + 0.3*$y*$z)]; sp $i 2 [expr 1e-3*(0.5*$x*$x - 0.2*$z)]; sp $i 3 [expr 1e-3*($y*$y + 0.1*$x)] + } + } + constraints Penalty 1e14 1e14; numberer Plain; system FullGeneral; test NormUnbalance 1e-10 5 + algorithm Linear; integrator LoadControl 1.0; analysis Static + analyze 1 + return [list [eleResponse 1 forces] [eleResponse 1 stresses]] +} +proc swap910 {f} { return [concat [lrange $f 0 23] [lrange $f 27 29] [lrange $f 24 26]] } +proc maxdiff {a b} { set m 0.0; foreach x $a y $b { set d [expr abs($x-$y)]; if {$d > $m} {set m $d} }; return $m } +proc maxabs {a} { set m 0.0; foreach x $a { if {abs($x) > $m} {set m [expr abs($x)]} }; return $m } + +lassign [run {1 2 3 4 5 6 7 8 9 10} ""] fg sg +lassign [run {1 2 3 4 5 6 7 8 10 9} "-numbering abaqus"] fa sa +lassign [run {1 2 3 4 5 6 7 8 10 9} ""] fw sw +puts [format "max |force| = %.4e" [maxabs $fg]] +puts [format "ABAQUS forces vs gmsh (9<->10 swapped): max diff = %.3e" [maxdiff $fa [swap910 $fg]]] +puts [format "ABAQUS stresses vs gmsh: max diff = %.3e" [maxdiff $sa $sg]] +puts [format "control: Abaqus conn WITHOUT flag: max diff = %.3e (must be large)" [maxdiff $fw [swap910 $fg]]] + +# patch test + body force in Abaqus order with straight edges (same check as smoke/tet10_patch.tcl) +wipe; model BasicBuilder -ndm 3 -ndf 3 +set X {{0 0 0} {2 0 0} {0 3 0} {0 0 1.5}} +set n 1; foreach p $X { eval node $n $p; incr n } +foreach e {{1 2} {2 3} {3 1} {1 4} {2 4} {3 4}} { + set a [lindex $X [expr [lindex $e 0]-1]]; set b [lindex $X [expr [lindex $e 1]-1]] + node $n [expr ([lindex $a 0]+[lindex $b 0])/2.] [expr ([lindex $a 1]+[lindex $b 1])/2.] [expr ([lindex $a 2]+[lindex $b 2])/2.]; incr n +} +nDMaterial ElasticIsotropic 1 1000.0 0.25 +element TenNodeTetrahedron 1 1 2 3 4 5 6 7 8 9 10 1 0.0 0.0 -2.0 -numbering abaqus +timeSeries Linear 1; pattern Plain 1 1 { + for {set i 1} {$i <= 10} {incr i} { sp $i 1 [expr 1e-3*[lindex [nodeCoord $i] 0]]; sp $i 2 0.0; sp $i 3 0.0 } +} +constraints Penalty 1e14 1e14; numberer Plain; system FullGeneral; test NormUnbalance 1e-10 5 +algorithm Linear; integrator LoadControl 1.0; analysis Static +analyze 1 +reactions +set Rz 0.0; for {set i 1} {$i <= 10} {incr i} { set Rz [expr $Rz + [nodeReaction $i 3]] } +puts [format "ABAQUS BODYFORCE sum Rz = %.10f expected 3.0" $Rz] +set s [eleResponse 1 stresses]; set m 0.0 +for {set g 0} {$g < 4} {incr g} { foreach k {0 1 2 3 4 5} e {1.2 0.4 0.4 0 0 0} { set d [expr abs([lindex $s [expr 6*$g+$k]]-$e)]; if {$d>$m} {set m $d} } } +puts [format "ABAQUS PATCH max abs error = %.3e" $m] From e9cbc2a4fe87aee1ec3441da590e04713bc9266b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Mon, 12 Oct 2026 22:23:08 -0300 Subject: [PATCH 4/6] TenNodeTetrahedron: document TenNodeTetrahedronStd and the parallel fixes The node-order option is now a second command, TenNodeTetrahedronStd (standard edge order), instead of a -numbering flag. Update the command syntax, numbering table, verified example and the fixes note (parallel support, debug output). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../model/elements/TenNodeTetrahedron.rst | 62 ++++++++++--------- .../example_tet10_numbering.tcl | 35 ++++++----- 2 files changed, 52 insertions(+), 45 deletions(-) diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index 65790eb7..31649a9f 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -4,22 +4,23 @@ TenNodeTetrahedron Element ^^^^^^^^^^^^^^^^^^^^^^^^^^ -This command is used to construct a ten-node tetrahedron element object, which uses the standard isoparametric formulation. +This command is used to construct a ten-node tetrahedron element object, which uses the standard isoparametric formulation. The two commands ``TenNodeTetrahedron`` and ``TenNodeTetrahedronStd`` create the same element and differ only in the order of mid-edge nodes 9 and 10 (see :ref:`TenNodeTetrahedronNumbering`). .. admonition:: Command - **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value> <-numbering $order>** + **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value>** + + **element TenNodeTetrahedronStd $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value>** .. csv-table:: :header: "Argument", "Type", "Description" :widths: 10, 10, 40 $eleTag, |integer|, unique element object tag - $node1 .. $node10, 10 |integer|, "nodes of tet: corners 1-4, then mid-edge nodes 5-10 in the order selected with ``-numbering`` (see :ref:`TenNodeTetrahedronNumbering` below)" + $node1 .. $node10, 10 |integer|, "nodes of tet: corners 1-4, then mid-edge nodes 5-10 (see :ref:`TenNodeTetrahedronNumbering` below)" $matTag, |integer|, tag of nDMaterial $b1 $b2 $b3, |listFloat|, optional: body forces in global x y z directions (default: 0.0 0.0 0.0) -doInitDisp $value, |integer|, "optional flag; if $value is non-zero, the nodal displacements present when the element is added to the domain are stored and subtracted from all later trial displacements before strains are computed (default: 0, i.e. disabled, so ordinary total displacements are used)" - -numbering $order, |string|, "optional: ``gmsh`` (default) or ``abaqus``, the convention used for mid-edge nodes 9 and 10 (see below)" This element is based on second-order interpolation of nodal quantities, so the strain and stress field inside the element vary linearly over the element. Four Gauss points inside the element are used for integration. @@ -31,7 +32,7 @@ This element is based on second-order interpolation of nodal quantities, so the :align: center :figclass: align-center - TenNodeTetrahedron Element Node Numbering (as drawn: ``-numbering abaqus``; with the default ``gmsh`` numbering nodes 9 and 10 are swapped) + TenNodeTetrahedron Element Node Numbering (as drawn: ``TenNodeTetrahedronStd``; for ``TenNodeTetrahedron`` nodes 9 and 10 are swapped) .. _TenNodeTetrahedronNumbering: @@ -41,7 +42,7 @@ Node numbering Corner nodes 1-4 and mid-edge nodes 5-8 are the same in both conventions. Only nodes 9 and 10 differ: .. csv-table:: - :header: "Node", "``gmsh`` (default)", "``abaqus`` (Abaqus C3D10, VTK quadratic tetra)" + :header: "Node", "``TenNodeTetrahedron`` (gmsh order)", "``TenNodeTetrahedronStd`` (standard order, as VTK)" :widths: 10, 20, 30 5, edge 1-2, edge 1-2 @@ -51,16 +52,19 @@ Corner nodes 1-4 and mid-edge nodes 5-8 are the same in both conventions. Only n 9, edge 3-4, edge 2-4 10, edge 2-4, edge 3-4 -The default matches gmsh's second-order tetrahedron, so meshes exported from gmsh can be used as they are. -For meshes in Abaqus (``C3D10``) order, add ``-numbering abaqus``. The element keeps the nodes in the order -given, so the ``forces`` response, the stiffness and mass matrices and the element connectivity seen by -recorders all follow the chosen convention. Giving Abaqus-ordered connectivity without the flag puts nodes 9 and 10 -on the wrong edges and produces wrong results without any error message. +``TenNodeTetrahedron`` follows gmsh's second-order tetrahedron, so meshes exported from gmsh can be used as they +are. ``TenNodeTetrahedronStd`` takes the standard order used by VTK's quadratic tetrahedron and most other +finite-element programs and mesh formats. The element keeps the nodes in the order given, so the ``forces`` response, +the stiffness and mass matrices and the connectivity seen by recorders all follow the convention of the command used. +Giving standard-order connectivity to ``TenNodeTetrahedron`` (or gmsh-order connectivity to +``TenNodeTetrahedronStd``) puts nodes 9 and 10 on the wrong edges and produces wrong results without any error +message. .. note:: - The figure above has always shown the Abaqus position of nodes 9 and 10, while the element has always used - the gmsh one by default. Models written from the figure need ``-numbering abaqus`` (or nodes 9 and 10 swapped). + The figure above has always shown the standard position of nodes 9 and 10, while ``TenNodeTetrahedron`` has + always used the gmsh one. Models written from the figure should use ``TenNodeTetrahedronStd`` (or swap nodes 9 + and 10). Theory """""" @@ -107,7 +111,7 @@ Each of the 4 Gauss points holds its own independent copy of the assigned ``nDMa .. note:: - Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, and ``-doInitDisp`` could not be used without body forces. All three are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. + Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, ``-doInitDisp`` could not be used without body forces, the element could not be used in OpenSeesSP/OpenSeesMP (it was missing from the object broker and sent only 4 of its 10 nodes), and analyses using the initial stiffness printed every shape-function value to the screen. All of these are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. This element can only be defined after a :ref:`model` with **-ndm 3 -ndf 3** @@ -135,22 +139,24 @@ Running it with the OpenSees interpreter gives, on the fixed element: OVERALL: PASS :download:`example_tet10_numbering.tcl ` checks the two -numberings against each other. It uses a curved element (mid-edge nodes moved off the edges) under an imposed -quadratic displacement field and a body force. It is defined once in gmsh order and once in Abaqus order -(nodes 9 and 10 swapped, ``-numbering abaqus``). The nodal forces must agree after swapping nodes 9 and 10, and the -Gauss-point stresses must be identical. As a control, the Abaqus connectivity without the flag must disagree. It -then repeats the body-force and uniform-strain patch tests in Abaqus order: +commands against each other. It uses a curved element (mid-edge nodes moved off the edges) under an imposed +quadratic displacement field and a body force. The element is defined once with ``TenNodeTetrahedron`` and once with +``TenNodeTetrahedronStd`` (nodes 9 and 10 swapped in the command). The nodal forces must agree after swapping nodes 9 +and 10, and the Gauss-point stresses must be identical. As a control, the standard-order connectivity given to +``TenNodeTetrahedron`` must disagree. The script then repeats the body-force and uniform-strain patch tests with +``TenNodeTetrahedronStd``: .. code-block:: text max |force| = 1.4612e+00 - ABAQUS forces vs gmsh (9<->10 swapped): max diff = 2.220e-16 - ABAQUS stresses vs gmsh: max diff = 1.110e-16 - control: Abaqus conn WITHOUT flag: max diff = 4.009e+00 (must be large) - ABAQUS BODYFORCE sum Rz = 3.0000000000 expected 3.0 - ABAQUS PATCH max abs error = 3.802e-11 + Std forces vs gmsh (9<->10 swapped): max diff = 2.220e-16 + Std stresses vs gmsh: max diff = 1.110e-16 + control: Std conn in TenNodeTetrahedron: max diff = 4.009e+00 (must be large) + Std BODYFORCE sum Rz = 3.0000000000 expected 3.0 + Std PATCH max abs error = 3.802e-11 -On the unfixed (vanilla upstream) element the same model instead gives a body-force reaction sum of ``0.5`` (should be ``3.0``) and a work-of-reactions under the imposed :math:`\epsilon_{xx}=10^{-3}` of ``3.0e-4`` (should be ``1.8e-3``), both exactly a factor of 6 too small, confirming the volume-scaling bug above. +The element also runs in OpenSeesSP/OpenSeesMP. A cube meshed with 48 elements of either command gives the same +displacements in serial and in OpenSeesSP on 2, 3 and 4 processes. .. admonition:: Example @@ -168,15 +174,15 @@ On the unfixed (vanilla upstream) element the same model instead gives a body-fo element('TenNodeTetrahedron',1, 1,2,3,4,5,6,7,8,9,10, 1, b1, b2, b3) - 3. **Abaqus (C3D10) node order** + 3. **Standard node order** (``TenNodeTetrahedronStd``) .. code-block:: tcl - element TenNodeTetrahedron 1 1 2 3 4 5 6 7 8 9 10 1 -numbering abaqus + element TenNodeTetrahedronStd 1 1 2 3 4 5 6 7 8 9 10 1 .. code-block:: python - element('TenNodeTetrahedron', 1, 1,2,3,4,5,6,7,8,9,10, 1, '-numbering', 'abaqus') + element('TenNodeTetrahedronStd', 1, 1,2,3,4,5,6,7,8,9,10, 1) This command is registered in both the Tcl and the (Python-backing) generic OpenSees interpreter. diff --git a/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl index ad5add20..cba842d2 100644 --- a/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl +++ b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl @@ -1,8 +1,8 @@ -# TenNodeTetrahedron: gmsh vs Abaqus (-numbering abaqus) node order. +# TenNodeTetrahedron (gmsh order) vs TenNodeTetrahedronStd (standard order). # Node tags: 1-4 corners, 5:1-2 6:2-3 7:3-1 8:1-4, 9 on edge 3-4, 10 on edge 2-4. -# gmsh connectivity: 1..10. Abaqus C3D10 connectivity: 1 2 3 4 5 6 7 8 10 9. +# TenNodeTetrahedron connectivity: 1..10. TenNodeTetrahedronStd connectivity: 1 2 3 4 5 6 7 8 10 9. # Midside nodes are moved off the edges (curved element) so a wrong 9/10 assignment shows. -proc build {conn opts b3} { +proc build {conn ele b3} { wipe; model BasicBuilder -ndm 3 -ndf 3 set X {{0 0 0} {2 0 0} {0 3 0} {0 0 1.5}} set mids {{1 2} {2 3} {3 1} {1 4} {3 4} {2 4}} @@ -16,11 +16,12 @@ proc build {conn opts b3} { eval node $n $c; incr n } nDMaterial ElasticIsotropic 1 1000.0 0.25 - eval element TenNodeTetrahedron 1 $conn 1 0.0 0.0 $b3 $opts + eval element $ele 1 $conn 1 0.0 0.0 $b3 } -# Same quadratic displacement field and body force imposed on each element; results are compared between orderings. -proc run {conn opts} { - build $conn $opts -2.0 +# quadratic displacement field imposed at all nodes (exactly representable by a curved Tet10 only in the +# sense of the isoparametric map; used here as a common input, results are compared between orderings) +proc run {conn ele} { + build $conn $ele -2.0 timeSeries Linear 1; pattern Plain 1 1 { for {set i 1} {$i <= 10} {incr i} { lassign [nodeCoord $i] x y z @@ -36,15 +37,15 @@ proc swap910 {f} { return [concat [lrange $f 0 23] [lrange $f 27 29] [lrange $f proc maxdiff {a b} { set m 0.0; foreach x $a y $b { set d [expr abs($x-$y)]; if {$d > $m} {set m $d} }; return $m } proc maxabs {a} { set m 0.0; foreach x $a { if {abs($x) > $m} {set m [expr abs($x)]} }; return $m } -lassign [run {1 2 3 4 5 6 7 8 9 10} ""] fg sg -lassign [run {1 2 3 4 5 6 7 8 10 9} "-numbering abaqus"] fa sa -lassign [run {1 2 3 4 5 6 7 8 10 9} ""] fw sw +lassign [run {1 2 3 4 5 6 7 8 9 10} TenNodeTetrahedron] fg sg +lassign [run {1 2 3 4 5 6 7 8 10 9} TenNodeTetrahedronStd] fa sa +lassign [run {1 2 3 4 5 6 7 8 10 9} TenNodeTetrahedron] fw sw puts [format "max |force| = %.4e" [maxabs $fg]] -puts [format "ABAQUS forces vs gmsh (9<->10 swapped): max diff = %.3e" [maxdiff $fa [swap910 $fg]]] -puts [format "ABAQUS stresses vs gmsh: max diff = %.3e" [maxdiff $sa $sg]] -puts [format "control: Abaqus conn WITHOUT flag: max diff = %.3e (must be large)" [maxdiff $fw [swap910 $fg]]] +puts [format "Std forces vs gmsh (9<->10 swapped): max diff = %.3e" [maxdiff $fa [swap910 $fg]]] +puts [format "Std stresses vs gmsh: max diff = %.3e" [maxdiff $sa $sg]] +puts [format "control: Std conn in TenNodeTetrahedron: max diff = %.3e (must be large)" [maxdiff $fw [swap910 $fg]]] -# patch test + body force in Abaqus order with straight edges (same check as smoke/tet10_patch.tcl) +# patch test + body force with TenNodeTetrahedronStd with straight edges (same check as smoke/tet10_patch.tcl) wipe; model BasicBuilder -ndm 3 -ndf 3 set X {{0 0 0} {2 0 0} {0 3 0} {0 0 1.5}} set n 1; foreach p $X { eval node $n $p; incr n } @@ -53,7 +54,7 @@ foreach e {{1 2} {2 3} {3 1} {1 4} {2 4} {3 4}} { node $n [expr ([lindex $a 0]+[lindex $b 0])/2.] [expr ([lindex $a 1]+[lindex $b 1])/2.] [expr ([lindex $a 2]+[lindex $b 2])/2.]; incr n } nDMaterial ElasticIsotropic 1 1000.0 0.25 -element TenNodeTetrahedron 1 1 2 3 4 5 6 7 8 9 10 1 0.0 0.0 -2.0 -numbering abaqus +element TenNodeTetrahedronStd 1 1 2 3 4 5 6 7 8 9 10 1 0.0 0.0 -2.0 timeSeries Linear 1; pattern Plain 1 1 { for {set i 1} {$i <= 10} {incr i} { sp $i 1 [expr 1e-3*[lindex [nodeCoord $i] 0]]; sp $i 2 0.0; sp $i 3 0.0 } } @@ -62,7 +63,7 @@ algorithm Linear; integrator LoadControl 1.0; analysis Static analyze 1 reactions set Rz 0.0; for {set i 1} {$i <= 10} {incr i} { set Rz [expr $Rz + [nodeReaction $i 3]] } -puts [format "ABAQUS BODYFORCE sum Rz = %.10f expected 3.0" $Rz] +puts [format "Std BODYFORCE sum Rz = %.10f expected 3.0" $Rz] set s [eleResponse 1 stresses]; set m 0.0 for {set g 0} {$g < 4} {incr g} { foreach k {0 1 2 3 4 5} e {1.2 0.4 0.4 0 0 0} { set d [expr abs([lindex $s [expr 6*$g+$k]]-$e)]; if {$d>$m} {set m $d} } } -puts [format "ABAQUS PATCH max abs error = %.3e" $m] +puts [format "Std PATCH max abs error = %.3e" $m] From dc127237980f1c8185bf190411611466ea2d4c4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Wed, 14 Oct 2026 18:14:05 -0300 Subject: [PATCH 5/6] TenNodeTetrahedron: note the ground-motion inertia-load fix Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- source/user/manual/model/elements/TenNodeTetrahedron.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index 31649a9f..8ff06e90 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -111,7 +111,7 @@ Each of the 4 Gauss points holds its own independent copy of the assigned ``nDMa .. note:: - Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, ``-doInitDisp`` could not be used without body forces, the element could not be used in OpenSeesSP/OpenSeesMP (it was missing from the object broker and sent only 4 of its 10 nodes), and analyses using the initial stiffness printed every shape-function value to the screen. All of these are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. + Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, ``-doInitDisp`` could not be used without body forces, the element could not be used in OpenSeesSP/OpenSeesMP (it was missing from the object broker and sent only 4 of its 10 nodes), analyses using the initial stiffness printed every shape-function value to the screen, and ground-motion inertia loads (``UniformExcitation``) were formed for 4 nodes only, giving wrong base reactions. All of these are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. This element can only be defined after a :ref:`model` with **-ndm 3 -ndf 3** From d8b764ff0da17c34a239b6a0296d3cbef0a0fa17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Thu, 15 Oct 2026 09:52:49 -0300 Subject: [PATCH 6/6] TenNodeTetrahedron: document the curved-element and robustness fixes Add a verified curved-element patch test, update the numbering example's expected output (its element is curved, so the numbers changed with the fix), and extend the fixes note. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../model/elements/TenNodeTetrahedron.rst | 24 +++++++++++++--- .../example_tet10_curved_patch.tcl | 28 +++++++++++++++++++ 2 files changed, 48 insertions(+), 4 deletions(-) create mode 100644 source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_curved_patch.tcl diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index 8ff06e90..20dd96c5 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -111,7 +111,7 @@ Each of the 4 Gauss points holds its own independent copy of the assigned ``nDMa .. note:: - Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, ``-doInitDisp`` could not be used without body forces, the element could not be used in OpenSeesSP/OpenSeesMP (it was missing from the object broker and sent only 4 of its 10 nodes), analyses using the initial stiffness printed every shape-function value to the screen, and ground-motion inertia loads (``UniformExcitation``) were formed for 4 nodes only, giving wrong base reactions. All of these are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. + Earlier versions of this element integrated the stiffness, mass and consistent body-force vectors over a volume **6 times smaller** than the true element volume (an erroneous extra division by 6 of the Jacobian determinant, on top of the `1/24` quadrature weights that already account for it), the ``stresses``/``strains`` element responses wrote 24 values (4 Gauss points x 6 components) into a statically-sized length-6 buffer, corrupting the heap on any call to ``eleResponse ... stresses`` or an equivalent recorder, ``-doInitDisp`` could not be used without body forces, the element could not be used in OpenSeesSP/OpenSeesMP (it was missing from the object broker and sent only 4 of its 10 nodes), analyses using the initial stiffness printed every shape-function value to the screen, ground-motion inertia loads (``UniformExcitation``) were formed for 4 nodes only, giving wrong base reactions, and elements with curved edges (mid-edge nodes off the edge midpoints, as in second-order meshes of curved boundaries) computed wrong strains because the shape-function derivatives used the corner coordinates instead of the full Jacobian. Smaller fixes in the same PR: an element with left-handed corner ordering (negative Jacobian) is now reported instead of silently giving negative stiffness and mass; ``eleResponse ... material`` without a Gauss-point number no longer crashes; the ``setDispInit`` parameter now works; a deactivated element no longer contributes initial stiffness; and ``print -JSON`` lists all 10 nodes. All of these are fixed in OpenSees PR "TenNodeTetrahedron fixes" (`#1833 `_); the example below reproduces and checks the fix. This element can only be defined after a :ref:`model` with **-ndm 3 -ndf 3** @@ -138,6 +138,22 @@ Running it with the OpenSees interpreter gives, on the fixed element: gp=4 sigma_xx before=1.200000 after=2.400000 (expected 2.400000) other-points-unchanged=yes OVERALL: PASS +:download:`example_tet10_curved_patch.tcl ` repeats the +uniform-strain patch test on one element with straight edges and with its mid-edge nodes moved off the edges (by about +0.1 and 0.3), for both commands. An isoparametric element reproduces a linear displacement field exactly even when +curved, so every case must give the exact stress: + +.. code-block:: text + + TenNodeTetrahedron bump=0.0 patch max stress error = 3.802e-11 |sum nodal forces| = 2.220e-16 + TenNodeTetrahedron bump=1.0 patch max stress error = 3.601e-11 |sum nodal forces| = 4.441e-16 + TenNodeTetrahedron bump=3.0 patch max stress error = 3.498e-11 |sum nodal forces| = 2.220e-16 + TenNodeTetrahedronStd bump=0.0 patch max stress error = 3.802e-11 |sum nodal forces| = 2.346e-16 + TenNodeTetrahedronStd bump=1.0 patch max stress error = 3.508e-11 |sum nodal forces| = 2.220e-16 + TenNodeTetrahedronStd bump=3.0 patch max stress error = 3.618e-11 |sum nodal forces| = 2.220e-16 + +Before the fix, the curved cases gave stress errors of 0.14-0.50 against an exact value of 1.2. + :download:`example_tet10_numbering.tcl ` checks the two commands against each other. It uses a curved element (mid-edge nodes moved off the edges) under an imposed quadratic displacement field and a body force. The element is defined once with ``TenNodeTetrahedron`` and once with @@ -148,10 +164,10 @@ and 10, and the Gauss-point stresses must be identical. As a control, the standa .. code-block:: text - max |force| = 1.4612e+00 + max |force| = 1.4035e+00 Std forces vs gmsh (9<->10 swapped): max diff = 2.220e-16 - Std stresses vs gmsh: max diff = 1.110e-16 - control: Std conn in TenNodeTetrahedron: max diff = 4.009e+00 (must be large) + Std stresses vs gmsh: max diff = 3.331e-16 + control: Std conn in TenNodeTetrahedron: max diff = 6.967e+00 (must be large) Std BODYFORCE sum Rz = 3.0000000000 expected 3.0 Std PATCH max abs error = 3.802e-11 diff --git a/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_curved_patch.tcl b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_curved_patch.tcl new file mode 100644 index 00000000..ff4afe92 --- /dev/null +++ b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_curved_patch.tcl @@ -0,0 +1,28 @@ +# Uniform-strain patch test on ONE Tet10 with straight vs curved edges (mid-edge nodes moved off the midpoints). +# Imposed u = (1e-3 x, 0, 0) on all 10 nodes. An isoparametric element reproduces any linear field exactly, +# straight or curved, so every Gauss point must give sigma = (1.2, 0.4, 0.4, 0, 0, 0) for E=1000, nu=0.25. +proc run {ele bumpScale} { + wipe; model BasicBuilder -ndm 3 -ndf 3 + set X {{0 0 0} {2 0 0} {0 3 0} {0 0 1.5}} + if {$ele == "TenNodeTetrahedron"} { set mids {{1 2} {2 3} {3 1} {1 4} {3 4} {2 4}} } else { set mids {{1 2} {2 3} {3 1} {1 4} {2 4} {3 4}} } + set bump {{0.05 -0.1 0.02} {0.1 0.05 -0.03} {-0.08 0.04 0.06} {0.03 0.07 -0.05} {-0.06 0.09 0.04} {0.07 -0.04 0.08}} + for {set i 0} {$i < 4} {incr i} { eval node [expr $i+1] [lindex $X $i] } + set n 5 + foreach e $mids d $bump { + set a [lindex $X [expr [lindex $e 0]-1]]; set b [lindex $X [expr [lindex $e 1]-1]] + set c {}; foreach k {0 1 2} { lappend c [expr ([lindex $a $k]+[lindex $b $k])/2. + $bumpScale*[lindex $d $k]] } + eval node $n $c; incr n + } + nDMaterial ElasticIsotropic 1 1000.0 0.25 + element $ele 1 1 2 3 4 5 6 7 8 9 10 1 + timeSeries Linear 1; pattern Plain 1 1 { for {set i 1} {$i <= 10} {incr i} { sp $i 1 [expr 1e-3*[lindex [nodeCoord $i] 0]]; sp $i 2 0.0; sp $i 3 0.0 } } + constraints Penalty 1e14 1e14; numberer Plain; system FullGeneral; test NormDispIncr 1e-14 5; algorithm Linear + integrator LoadControl 1.0; analysis Static; analyze 1 + set s [eleResponse 1 stresses]; set m 0.0 + for {set g 0} {$g < 4} {incr g} { foreach k {0 1 2 3 4 5} e {1.2 0.4 0.4 0 0 0} { set d [expr abs([lindex $s [expr 6*$g+$k]]-$e)]; if {$d>$m} {set m $d} } } + # element force equilibrium: sum of nodal forces must vanish + set f [eleResponse 1 forces]; set r 0.0 + foreach k {0 1 2} { set t 0.0; for {set i 0} {$i < 10} {incr i} { set t [expr $t+[lindex $f [expr 3*$i+$k]]] }; if {abs($t)>$r} {set r [expr abs($t)]} } + return [format "%-22s bump=%.1f patch max stress error = %.3e |sum nodal forces| = %.3e" $ele $bumpScale $m $r] +} +foreach ele {TenNodeTetrahedron TenNodeTetrahedronStd} { foreach s {0.0 1.0 3.0} { puts [run $ele $s] } }