Skip to content

Commit 42ded19

Browse files
claudioperezclaudioperez
authored
Merge pull request #1 from Facundo-Pfeffer/master
Constrain Documentation
2 parents 8744f88 + d8d2583 commit 42ded19

3 files changed

Lines changed: 321 additions & 19 deletions

File tree

‎source/user/manual/model/constraints/constrain.rst‎

Lines changed: 215 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,37 +1,234 @@
1-
1+
.. _constrain:
22

33
constrain
44
^^^^^^^^^
55

6+
The ``constrain`` command imposes a multi-point constraint between two nodes.
7+
When no rotation is specified, it behaves like :ref:`equalDOF`, constraining the degrees-of-freedom at the constrained node to equal those at the retained node.
8+
When a rotation vector is provided, the constraint enforces a rotated relationship between the nodes through a rotation matrix.
69

710
.. tabs::
811

9-
.. tab:: Python (RT)
10-
12+
.. tab:: Python
13+
1114
.. py:method:: Model.constrain(rNode, cNode, rotate=None)
12-
13-
:param rNode: tag of the retained :ref:`node`
14-
:type rNode: |integer|
15-
:param cNode: tag of the constrained node
16-
:type cNode: |integer|
17-
:param rotate: optional rotation vector
18-
:type rotate: |list| of |float|
1915
16+
Impose a multi-point constraint between a retained node and a constrained node.
17+
18+
:param int|tuple rNode: integer tag identifying the retained node, or tuple (rNode, cNode) for both nodes
19+
:param int cNode: integer tag identifying the constrained node (required if rNode is not a tuple)
20+
:param list rotate: optional rotation vector [v1, v2, v3] representing axis-angle in radians (see :ref:`Rotation Vector <constrain-rotation-vector>` below). Only available for 3D problems (ndm=3, ndf=6).
21+
:return: None
2022

2123
.. tab:: Tcl
2224

23-
.. function:: constrain rNode cNode <-rotate vector>
25+
.. function:: constrain $rNodeTag $cNodeTag <-rotate {v1 v2 v3}>
2426

25-
:param rNode: tag of the retained node
26-
:type rNode: |integer|
27-
:param cNode: tag of the constrained node
28-
:type cNode: |integer|
27+
.. csv-table::
28+
:header: "Argument", "Type", "Description"
29+
:widths: 10, 10, 40
2930

31+
$rNodeTag, |integer|, integer tag identifying the retained node (*rNode*)
32+
$cNodeTag, |integer|, integer tag identifying the constrained node (*cNode*)
33+
-rotate {v1 v2 v3}, |optional|, rotation vector representing axis-angle in radians (see :ref:`Rotation Vector <constrain-rotation-vector>` below). Only available for 3D problems (ndm=3, ndf=6).
3034

31-
Theory
35+
36+
Theory
3237
------
3338

34-
``code``. *italic* **bold**.
39+
The ``constrain`` command creates a multi-point constraint that enforces the relationship :math:`\boldsymbol{u}_c = \boldsymbol{C}_{cr} \boldsymbol{u}_r`, where :math:`\boldsymbol{u}_c` is the response vector at the constrained node, :math:`\boldsymbol{u}_r` is the response vector at the retained node, and :math:`\boldsymbol{C}_{cr}` is the constraint matrix.
40+
41+
.. _constrain-without-rotation:
42+
43+
Without Rotation
44+
~~~~~~~~~~~~~~~~
45+
46+
When no rotation is specified, the constraint matrix is the identity matrix, resulting in behavior equivalent to :ref:`equalDOF`: :math:`\boldsymbol{C}_{cr} = \boldsymbol{I}`.
47+
This form of the constraint works for both 2D and 3D problems.
48+
49+
With Rotation
50+
~~~~~~~~~~~~~
51+
52+
When a rotation vector :math:`\boldsymbol{v} = [v_1, v_2, v_3]` is provided, the rotation matrix :math:`\boldsymbol{R}` is computed using `Rodrigues' rotation formula <https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#:~:text=In%20terms%20of%20the%20matrix%20exponential%2C>`_ as the matrix exponential :math:`\boldsymbol{R} = \exp(\theta \boldsymbol{K})`, where :math:`\theta` is the magnitude of the rotation vector and :math:`\boldsymbol{K}` is the skew-symmetric matrix corresponding to the normalized rotation axis.
53+
54+
.. _constrain-rotation-vector:
55+
56+
Rotation Vector
57+
~~~~~~~~~~~~~~~
58+
59+
The rotation vector uses the **axis-angle representation**: the normalized direction defines the rotation axis, and the magnitude is the rotation angle in radians. The rotation follows the **right-hand rule**: a positive angle rotates counterclockwise when viewed from the tip of the axis vector.
60+
61+
In 3D models, the components :math:`[v_1, v_2, v_3]` correspond to rotations about the global X, Y, and Z axes respectively:
62+
- :math:`[\theta, 0, 0]` rotates about the X-axis
63+
- :math:`[0, \theta, 0]` rotates about the Y-axis
64+
- :math:`[0, 0, \theta]` rotates about the Z-axis
65+
66+
The coordinate system and degrees-of-freedom conventions are documented in :ref:`Conventions <conventions>`.
67+
68+
For 3D problems (ndm=3, ndf=6), the constraint matrix is:
69+
70+
.. math::
71+
72+
\boldsymbol{C}_{cr} = \begin{bmatrix}
73+
\boldsymbol{R} & \boldsymbol{0} \\
74+
\boldsymbol{0} & \boldsymbol{R}
75+
\end{bmatrix}
3576
36-
Example
77+
where the upper-left block applies the rotation matrix :math:`\boldsymbol{R}` to translational degrees-of-freedom (1-3), and the lower-right block applies the same rotation matrix to rotational degrees-of-freedom (4-6).
78+
79+
.. note::
80+
81+
**3D limitation:** The rotation option is only implemented for 3D problems (ndm=3, ndf=6) because Rodrigues' formula requires a 3D rotation vector, and the constraint matrix applies the rotation to all six degrees-of-freedom.
82+
83+
**2D alternative:** For 2D problems requiring rotated constraints (e.g., skewed supports), use the penalty method with :ref:`zeroLength <zeroLength>` elements, as demonstrated in the :ref:`Penalty Method for Skewed Supports in 2D <zeroLength-penalty-2d>` example.
84+
85+
**Rotational DOF approximation:** The same rotation matrix :math:`\boldsymbol{R}` is applied to both translational and rotational DOFs. While this is a proper coordinate transformation for translational DOFs, it is an approximation for rotational DOFs (treating rotations as vectors) that is valid for small rotations but may have limitations for large finite rotations.
86+
87+
Examples
3788
--------
89+
90+
Basic Constraint (No Rotation)
91+
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
92+
93+
This example demonstrates the basic usage of ``constrain`` without rotation, corresponding to the :ref:`Without Rotation <constrain-without-rotation>` theory section above. When no rotation is specified, the constraint matrix is the identity matrix :math:`\boldsymbol{C}_{cr} = \boldsymbol{I}`, making this equivalent to :ref:`equalDOF`.
94+
95+
The following example constrains node **5** to have the same degrees-of-freedom as node **100**. This works for both 2D and 3D problems:
96+
97+
.. tabs::
98+
.. tab:: Tcl
99+
100+
.. code-block:: none
101+
102+
constrain 100 5
103+
104+
.. tab:: Python
105+
106+
.. code-block:: python
107+
108+
model.constrain(100, 5)
109+
110+
Constraint with Rotation
111+
~~~~~~~~~~~~~~~~~~~~~~~~~
112+
113+
The following example constrains node **5** to node **100** with a rotation about the Y-axis. This requires a 3D model (ndm=3, ndf=6):
114+
115+
.. tabs::
116+
.. tab:: Tcl
117+
118+
.. code-block:: none
119+
120+
constrain 100 5 -rotate {0 -0.6435 0}
121+
122+
.. tab:: Python
123+
124+
.. code-block:: python
125+
126+
angle = 0.5 # Rotation angle in radians (about 28.6 degrees)
127+
# Negative sign rotates clockwise about Y-axis (right-hand rule)
128+
model.constrain(100, 5, rotate=[0, -angle, 0])
129+
130+
Frame with Skewed Support
131+
~~~~~~~~~~~~~~~~~~~~~~~~~~
132+
133+
This example demonstrates using ``constrain`` to model a frame with a skewed support in 3D, where the support is not aligned with the global coordinate system.
134+
135+
.. code-block:: python
136+
:linenos:
137+
138+
import xara
139+
from xara.units.english import inch, kip, ksi
140+
from shps.rotor import exp
141+
142+
# Geometry
143+
H = 144.0 * inch
144+
W = 144.0 * inch
145+
146+
# Section properties
147+
width = 12.0 * inch
148+
height = 12.0 * inch
149+
Iz = (width * height**3) / 12.0
150+
Iy = (height * width**3) / 12.0
151+
A = width * height
152+
J = (width * height**3) / 3.0
153+
154+
# Material
155+
E = 29000 * ksi
156+
G = 11500 * ksi
157+
nu = 0.30
158+
159+
# Load
160+
P = -10.0 * kip
161+
angle = 0.5 # Rotation angle in radians (about 28.6 degrees)
162+
# Negative sign rotates clockwise about Y-axis (right-hand rule),
163+
# rotating the support coordinate system clockwise from vertical
164+
165+
# Compute rotation matrix for reaction transformation
166+
R = exp([0.0, -angle, 0.0])
167+
168+
# Model (3D required for rotation)
169+
model = xara.Model('basic', ndm=3, ndf=6)
170+
171+
# Nodes
172+
model.node(1, (0.0, 0.0, 0.0)) # Fixed support
173+
model.node(2, (0.0, 0.0, H)) # Top-left
174+
model.node(3, (W/2, 0.0, H)) # Load point
175+
model.node(4, (W, 0.0, H)) # Top-right
176+
model.node(5, (W, 0.0, 0.0)) # Skewed support
177+
178+
# Fixed support at node 1
179+
model.fix(1, (1, 1, 1, 1, 1, 1))
180+
181+
# Ground node for skewed support
182+
ground_node = 100
183+
model.node(ground_node, (W, 0.0, 0.0))
184+
model.fix(ground_node, (0, 1, 1, 0, 0, 0)) # Fixed in Y and Z, free in X
185+
186+
# Apply rotated constraint: node 5 is constrained to ground_node with rotation
187+
model.constrain((ground_node, 5), rotate=[0, -angle, 0])
188+
189+
# Frame section and elements
190+
model.material('ElasticIsotropic', 1, E, nu)
191+
model.section("ElasticFrame", 1, E=E, G=G, A=A, Iy=Iy, Iz=Iz, J=J)
192+
model.geomTransf("Linear", 1, (0, 1, 0)) # Columns
193+
model.geomTransf("Linear", 2, (0, 1, 0)) # Beams
194+
195+
model.element('PrismFrame', 1, (1, 2), section=1, transform=1)
196+
model.element('PrismFrame', 2, (2, 3), section=1, transform=2)
197+
model.element('PrismFrame', 3, (3, 4), section=1, transform=2)
198+
model.element('PrismFrame', 4, (4, 5), section=1, transform=1)
199+
200+
# Load
201+
model.timeSeries('Constant', 1)
202+
model.pattern('Plain', 1, 1)
203+
model.load(3, (0.0, 0.0, P, 0.0, 0.0, 0.0))
204+
205+
# Analysis
206+
model.system('BandGeneral')
207+
model.numberer('RCM')
208+
model.constraints('Transformation')
209+
model.integrator('LoadControl', 1.0)
210+
model.algorithm('Newton')
211+
model.test('Energy', 1e-8, 10)
212+
model.analysis('Static')
213+
model.analyze(1)
214+
215+
# Results
216+
model.reactions()
217+
218+
# Reaction at fixed support (node 1)
219+
Fz_jt1 = model.nodeReaction(1, 3) / kip
220+
print(f"Reaction at node 1 (Z-direction): {Fz_jt1:.3f} kip")
221+
222+
# Reaction at skewed support (node 5) - transform to rotated coordinate system
223+
R4 = model.nodeReaction(5)[:3] # Get translational reactions (X, Y, Z)
224+
r4 = R.T @ R4 / kip # Transform to rotated coordinate system
225+
F3_jt4 = r4[2] # Z-component in rotated coordinate system
226+
print(f"Reaction at skewed support (node 5, normal direction): {F3_jt4:.3f} kip")
227+
228+
References
229+
----------
230+
231+
* Cook, R.D., Malkus, D.S., Plesha, M. E., and Witt, R. J., "Concepts and Applications of Finite Element Analysis," 4th edition, John Wiley and Sons publishers, 2002.
232+
233+
* `Rodrigues' rotation formula <https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula>`_ - Wikipedia article on Rodrigues' rotation formula.
234+

‎source/user/manual/model/constraints/index.rst‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,14 +4,15 @@ Constraints
44
***********
55

66
*Multi-point constraints* are constraints that allow the user to define the relationship between the response of a set of the degrees-of-freedom at one node (the constrained node) in relation to the response of the degrees-of-freedom at another node (the retained node).
7-
In structural analysis MP_Constraints are used to enforce rigid-diapgragm constraints or equal constraints (response of degrees-of-freedom at two nodes move the same)
7+
In structural analysis MP_Constraints are used to enforce rigid-diaphragm constraints, equal constraints (response of degrees-of-freedom at two nodes move the same), or rotated constraints (response of degrees-of-freedom at two nodes related through a rotation matrix).
88

99

1010
.. toctree::
1111
:maxdepth: 1
1212

1313
constrain
1414
equalDOF
15+
constrain
1516
diaphragm
1617
rigidLink
1718

‎source/user/manual/model/elements/other/zeroLength.rst‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,4 +63,108 @@ Element **2** has as its end nodes **4** and **5**, has only one material **1**
6363
model.element("zeroLength",2,4,5,"-mat",1,"-dir",1,"-orient",1,1,0,-1,1,0)
6464
model.element("zeroLength",3,5,6,"-mat",1,"-dir",1,"-doRayleigh",1)
6565
66+
.. _zeroLength-penalty-2d:
67+
68+
Penalty Method for Skewed Supports in 2D
69+
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
70+
71+
ZeroLength elements can be used with the penalty method to model rotated constraints in 2D problems, such as skewed supports. This approach is particularly useful when the :ref:`constrain <constrain>` command's rotation option is not available (it is only implemented for 3D problems).
72+
73+
The penalty method uses a very stiff spring (penalty stiffness) oriented in the desired constraint direction to approximate a rigid constraint. The zeroLength element connects a ground node (with fixed boundary conditions) to the constrained node, with the element's local coordinate system oriented to match the desired constraint direction.
74+
75+
The following example demonstrates modeling a 2D frame with a skewed support using the penalty method with a zeroLength element:
76+
77+
.. code-block:: python
78+
:linenos:
79+
80+
import xara
81+
from xara.units.english import inch, kip, ksi
82+
from math import atan2
83+
from shps.rotor import exp
84+
85+
# Geometry
86+
H = 144.0 * inch
87+
W = 144.0 * inch
88+
89+
# Section properties
90+
width = 12.0 * inch
91+
height = 12.0 * inch
92+
Iz = (width * height**3) / 12.0
93+
A = width * height
94+
95+
# Material
96+
E = 29000 * ksi
97+
nu = 0.30
98+
99+
# Load
100+
P = -10.0 * kip
101+
angle = atan2(0.6, 0.8) # Rotation angle for skewed support
102+
103+
# Model (2D)
104+
model = xara.Model('basic', ndm=2, ndf=3)
105+
106+
# Nodes
107+
model.node(1, (0.0, 0.0)) # Fixed support
108+
model.node(2, (0.0, H)) # Top-left
109+
model.node(3, (W/2, H)) # Load point
110+
model.node(4, (W, H)) # Top-right
111+
model.node(5, (W, 0.0)) # Skewed support
112+
113+
# Fixed support at node 1
114+
model.fix(1, (1, 1, 1))
115+
116+
# Skewed support at node 5 using penalty method
117+
cos_theta = 0.8
118+
sin_theta = 0.6
119+
120+
R = exp([0.0, 0.0, angle])
121+
ground_node = 100
122+
model.node(ground_node, (W, 0.0))
123+
model.fix(ground_node, (1, 1, 1))
124+
125+
# Penalty stiffness (very large to approximate rigid constraint)
126+
k_penalty = 1.0e10 * kip / inch
127+
model.uniaxialMaterial('Elastic', 1001, k_penalty)
128+
129+
# Normal direction for constraint (perpendicular to support)
130+
nx, ny = R[:2, 0]
131+
# Tangential direction (along support)
132+
tx, ty = cos_theta, sin_theta
133+
134+
# Create zeroLength element with oriented local axes
135+
# The element constrains motion in the normal direction (dir=2, local y-axis)
136+
model.element('zeroLength', 1001, (ground_node, 5),
137+
mat=1001,
138+
dir=2, # Constrain in local y-direction (normal to support)
139+
orient=(nx, ny, 0.0, tx, ty, 0.0))
140+
141+
# Frame section and elements
142+
model.material('ElasticIsotropic', 1, E, nu)
143+
model.section("ElasticFrame", 1, E=E, G=E/(2*(1+nu)), A=A, Iy=0, Iz=Iz, J=0)
144+
model.geomTransf("Linear", 1)
145+
146+
model.element('PrismFrame', 1, (1, 2), section=1, transform=1)
147+
model.element('PrismFrame', 2, (2, 3), section=1, transform=1)
148+
model.element('PrismFrame', 3, (3, 4), section=1, transform=1)
149+
model.element('PrismFrame', 4, (4, 5), section=1, transform=1)
150+
151+
# Load
152+
model.timeSeries('Constant', 1)
153+
model.pattern('Plain', 1, 1)
154+
model.load(3, (0.0, P, 0.0))
155+
156+
# Analysis
157+
model.system('BandGeneral')
158+
model.numberer('RCM')
159+
model.constraints('Transformation')
160+
model.integrator('LoadControl', 1.0)
161+
model.algorithm('Newton')
162+
model.test('Energy', 1e-10, 10)
163+
model.analysis('Static')
164+
model.analyze(1)
165+
166+
.. note::
167+
168+
The penalty stiffness should be chosen large enough to approximate a rigid constraint, but not so large that it causes numerical conditioning problems. A typical value is 1.0e10 times the characteristic stiffness of the structure. The orientation vectors define the local coordinate system of the zeroLength element, allowing the constraint to be applied in the desired direction (normal to the skewed support).
169+
66170
Code Developed by: |glf|

0 commit comments

Comments
 (0)