OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
thskewc.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| thskewc ../starter/source/output/th/thskewc.F
25!||--- called by ------------------------------------------------------
26!|| lectur ../starter/source/starter/lectur.F
27!||====================================================================
28 SUBROUTINE thskewc(
29 1 RTHBUF ,ITHGRP ,ITHBUF,X ,IXC ,IXTG ,SKEW,NTHGRP)
30
31
32C-----------------------------------------------
33C I m p l i c i t T y p e s
34C-----------------------------------------------
35#include "implicit_f.inc"
36C-----------------------------------------------
37C C o m m o n B l o c k s
38C-----------------------------------------------
39#include "param_c.inc"
40C-----------------------------------------------
41C D u m m y A r g u m e n t s
42C-----------------------------------------------
43 INTEGER NTHGRP,ITHGRP(NITHGR,*),ITHBUF(*),IXC(NIXC,*),IXTG(NIXTG,*)
44 my_real
45 . rthbuf(*), x(3,*), skew(lskew,*)
46C-----------------------------------------------
47C L o c a l V a r i a b l e s
48C-----------------------------------------------
49 INTEGER NNE,IAD,IAD2,IADR,ISK,NN,N1,N2,N3,N4,IGS,N,ITYP,K
50C REAL
51 my_real
52 . x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4,
53 . e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z,
54 . x31, y31, z31, x32, y32, z32, x21, y21, z21,
55 . x42, y42, z42, s1, s2, vx, vy, vz, v, vr, vs,
56 . suma,area
57
58C Fill table RTHBUF
59 iadr=0
60 DO n=1,nthgrp
61 ityp=ithgrp(2,n)
62 nne =ithgrp(4,n)
63 iad =ithgrp(5,n)
64 iad2=iad+3*nne
65 IF(ityp==3)THEN
66 DO k=1,nne
67 nn=ithbuf(iad)
68 isk=1+ithbuf(iad2)
69c
70 IF(isk > 1) THEN
71C Corotational Frame E1 E2 E3
72 n1=ixc(2,nn)
73 n2=ixc(3,nn)
74 n3=ixc(4,nn)
75 n4=ixc(5,nn)
76
77 x1=x(1,n1)
78 x2=x(1,n2)
79 x3=x(1,n3)
80 x4=x(1,n4)
81
82 y1=x(2,n1)
83 y2=x(2,n2)
84 y3=x(2,n3)
85 y4=x(2,n4)
86
87 z1=x(3,n1)
88 z2=x(3,n2)
89 z3=x(3,n3)
90 z4=x(3,n4)
91
92
93 x21=x2-x1
94 y21=y2-y1
95 z21=z2-z1
96 x31=x3-x1
97 y31=y3-y1
98 z31=z3-z1
99 x42=x4-x2
100 y42=y4-y2
101 z42=z4-z2
102
103 e3x=y31*z42-z31*y42
104 e3y=z31*x42-x31*z42
105 e3z=x31*y42-y31*x42
106 suma=sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
107
108 e1x = x2+x3-x1-x4
109 e1y = y2+y3-y1-y4
110 e1z = z2+z3-z1-z4
111c
112 e2x = x3+x4-x1-x2
113 e2y = y3+y4-y1-y2
114 e2z = z3+z4-z1-z2
115c
116 e3x = e1y*e2z-e1z*e2y
117 e3y = e1z*e2x-e1x*e2z
118 e3z = e1x*e2y-e1y*e2x
119
120 suma = e3x*e3x+e3y*e3y+e3z*e3z
121 suma = one/max(sqrt(suma),em20)
122 e3x = e3x*suma
123 e3y = e3y*suma
124 e3z = e3z*suma
125c
126 s1 = e1x*e1x+e1y*e1y+e1z*e1z
127 s2 = e2x*e2x+e2y*e2y+e2z*e2z
128 suma = sqrt(s1/s2)
129 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
130 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
131 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
132c
133 suma = e1x*e1x+e1y*e1y+e1z*e1z
134 suma = one/max(sqrt(suma),em20)
135 e1x = e1x*suma
136 e1y = e1y*suma
137 e1z = e1z*suma
138c
139 e2x = e3y * e1z - e3z * e1y
140 e2y = e3z * e1x - e3x * e1z
141 e2z = e3x * e1y - e3y * e1x
142
143C Project First axe of the skew
144 vx = skew(1,isk)
145 vy = skew(2,isk)
146 vz = skew(3,isk)
147
148 v =vx*e3x+vy*e3y+vz*e3z
149 vx=vx-v*e3x
150 vy=vy-v*e3y
151 vz=vz-v*e3z
152 v =sqrt(vx*vx+vy*vy+vz*vz)
153
154 vx=vx/max(v,em20)
155 vy=vy/max(v,em20)
156 vz=vz/max(v,em20)
157
158C Cos and Sin calculation
159 vr=vx*e1x+vy*e1y+vz*e1z
160 vs=vx*e2x+vy*e2y+vz*e2z
161c Save data in RTHBUF
162 ithbuf(iad2)=iadr+1
163 rthbuf(iadr+1)=vr
164 rthbuf(iadr+2)=vs
165
166 iadr=iadr+2
167 ENDIF
168 iad=iad+1
169 iad2=iad2+1
170 ENDDO
171 ELSEIF(ityp==7)THEN
172 DO k=1,nne
173 nn=ithbuf(iad)
174 isk=ithbuf(iad2)
175 IF(isk /= 0) THEN
176C Corotational Frame E1 E2 E3
177 n1=ixc(2,nn)
178 n2=ixc(3,nn)
179 n3=ixc(4,nn)
180
181 x1=x(1,n1)
182 x2=x(1,n2)
183 x3=x(1,n3)
184
185 y1=x(2,n1)
186 y2=x(2,n2)
187 y3=x(2,n3)
188
189 z1=x(3,n1)
190 z2=x(3,n2)
191 z3=x(3,n3)
192
193 x21=x2-x1
194 y21=y2-y1
195 z21=z2-z1
196 x31=x3-x1
197 y31=y3-y1
198 z31=z3-z1
199 x32=x3-x2
200 y32=y3-y2
201 z32=z3-z2
202c
203 e1x= x21
204 e1y= y21
205 e1z= z21
206 suma = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
207 e1x=e1x/suma
208 e1y=e1y/suma
209 e1z=e1z/suma
210c
211 e3x=y31*z32-z31*y32
212 e3y=z31*x32-x31*z32
213 e3z=x31*y32-y31*x32
214 suma = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
215 e3x=e3x/suma
216 e3y=e3y/suma
217 e3z=e3z/suma
218 area = half * suma
219c
220 e2x=e3y*e1z-e3z*e1y
221 e2y=e3z*e1x-e3x*e1z
222 e2z=e3x*e1y-e3y*e1x
223 suma = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
224 e2x=e2x/suma
225 e2y=e2y/suma
226 e2z=e2z/suma
227
228C Project First axe of the skew
229 vx = skew(1,isk)
230 vy = skew(2,isk)
231 vz = skew(3,isk)
232
233 v =vx*e3x+vy*e3y+vz*e3z
234 vx=vx-v*e3x
235 vy=vy-v*e3y
236 vz=vz-v*e3z
237 v =sqrt(vx*vx+vy*vy+vz*vz)
238
239 vx=vx/max(v,em20)
240 vy=vy/max(v,em20)
241 vz=vz/max(v,em20)
242C Cos and Sin calculation
243 vr=vx*e1x+vy*e1y+vz*e1z
244 vs=vx*e2x+vy*e2y+vz*e2z
245
246c Save data in RTHBUF
247 ithbuf(iad2)=iadr+1
248 rthbuf(iadr+1)=vr
249 rthbuf(iadr+2)=vs
250 iadr=iadr+2
251 ENDIF
252 iad=iad+1
253 iad2=iad2+1
254 ENDDO
255 ENDIF
256 ENDDO
257 RETURN
258 END
259
my_real
#define my_real
Definition cppsort.cpp:32
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
max
#define max(a, b)
Definition macros.h:21
thskewc
subroutine thskewc(rthbuf, ithgrp, ithbuf, x, ixc, ixtg, skew, nthgrp)
Definition thskewc.F:30