diff --git a/source/user/manual/model/elements/TenNodeTetrahedron.rst b/source/user/manual/model/elements/TenNodeTetrahedron.rst index aced5071..20dd96c5 100644 --- a/source/user/manual/model/elements/TenNodeTetrahedron.rst +++ b/source/user/manual/model/elements/TenNodeTetrahedron.rst @@ -4,38 +4,177 @@ 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. 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> ** + **element TenNodeTetrahedron $eleTag $node1 $node2 $node3 $node4 $node5 $node6 $node7 $node8 $node9 $node10 $matTag <$b1 $b2 $b3> <-doInitDisp $value>** -.. csv-table:: + **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 (ordered as shown in fig 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 - <-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, so the strain and stress field inside the element vary linearly over the element. Four Gauss points inside the element are used for integration. -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. +.. 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 :figclass: align-center - TenNodeTetrahedron Element Node Numbering + TenNodeTetrahedron Element Node Numbering (as drawn: ``TenNodeTetrahedronStd``; for ``TenNodeTetrahedron`` 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", "``TenNodeTetrahedron`` (gmsh order)", "``TenNodeTetrahedronStd`` (standard order, as VTK)" + :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 + +``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 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 +"""""" + +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, ``-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** -.. 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 + +: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 +``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.4035e+00 + Std forces vs gmsh (9<->10 swapped): max diff = 2.220e-16 + 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 + +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 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 +190,16 @@ 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). + 3. **Standard node order** (``TenNodeTetrahedronStd``) + + .. code-block:: tcl + + element TenNodeTetrahedronStd 1 1 2 3 4 5 6 7 8 9 10 1 + + .. code-block:: python + + 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. + +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_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] } } 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..cba842d2 --- /dev/null +++ b/source/user/manual/model/elements/figures/TenNodeTetrahedron/example_tet10_numbering.tcl @@ -0,0 +1,69 @@ +# 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. +# 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 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}} + 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 $ele 1 $conn 1 0.0 0.0 $b3 +} +# 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 + 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} 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 "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 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 } +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 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 } +} +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 "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 "Std PATCH max abs error = %.3e" $m] 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" +}