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!||--- uses -----------------------------------------------------
28!||====================================================================
29 SUBROUTINE thskewc(
30 1 RTHBUF ,ITHGRP ,ITHBUF,X ,IXC ,IXTG ,SKEW,NTHGRP)
31 use element_mod , only : nixc,nixtg
32
33C-----------------------------------------------
34C I m p l i c i t T y p e s
35C-----------------------------------------------
36#include "implicit_f.inc"
37C-----------------------------------------------
38C C o m m o n B l o c k s
39C-----------------------------------------------
40#include "param_c.inc"
41C-----------------------------------------------
42C D u m m y A r g u m e n t s
43C-----------------------------------------------
44 INTEGER NTHGRP,ITHGRP(NITHGR,*),ITHBUF(*),IXC(NIXC,*),IXTG(NIXTG,*)
45 my_real
46 . rthbuf(*), x(3,*), skew(lskew,*)
47C-----------------------------------------------
48C L o c a l V a r i a b l e s
49C-----------------------------------------------
50 INTEGER NNE, IAD, IAD2, IADR, ISK, NN, N1, N2, N3, N4, N, ITYP, K
51C REAL
52 my_real
53 . x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4,
54 . e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z,
55 . x31, y31, z31, x32, y32, z32, x21, y21, z21,
56 . x42, y42, z42, s1, s2, vx, vy, vz, v, vr, vs,
57 . suma,area
58
59C Fill table RTHBUF
60 iadr=0
61 DO n=1,nthgrp
62 ityp=ithgrp(2,n)
63 nne =ithgrp(4,n)
64 iad =ithgrp(5,n)
65 iad2=iad+3*nne
66 IF(ityp==3)THEN
67 DO k=1,nne
68 nn=ithbuf(iad)
69 isk=1+ithbuf(iad2)
70c
71 IF(isk > 1) THEN
72C Corotational Frame E1 E2 E3
73 n1=ixc(2,nn)
74 n2=ixc(3,nn)
75 n3=ixc(4,nn)
76 n4=ixc(5,nn)
77
78 x1=x(1,n1)
79 x2=x(1,n2)
80 x3=x(1,n3)
81 x4=x(1,n4)
82
83 y1=x(2,n1)
84 y2=x(2,n2)
85 y3=x(2,n3)
86 y4=x(2,n4)
87
88 z1=x(3,n1)
89 z2=x(3,n2)
90 z3=x(3,n3)
91 z4=x(3,n4)
92
93
94 x21=x2-x1
95 y21=y2-y1
96 z21=z2-z1
97 x31=x3-x1
98 y31=y3-y1
99 z31=z3-z1
100 x42=x4-x2
101 y42=y4-y2
102 z42=z4-z2
103
104 e3x=y31*z42-z31*y42
105 e3y=z31*x42-x31*z42
106 e3z=x31*y42-y31*x42
107 suma=sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
108
109 e1x = x2+x3-x1-x4
110 e1y = y2+y3-y1-y4
111 e1z = z2+z3-z1-z4
112c
113 e2x = x3+x4-x1-x2
114 e2y = y3+y4-y1-y2
115 e2z = z3+z4-z1-z2
116c
117 e3x = e1y*e2z-e1z*e2y
118 e3y = e1z*e2x-e1x*e2z
119 e3z = e1x*e2y-e1y*e2x
120
121 suma = e3x*e3x+e3y*e3y+e3z*e3z
122 suma = one/max(sqrt(suma),em20)
123 e3x = e3x*suma
124 e3y = e3y*suma
125 e3z = e3z*suma
126c
127 s1 = e1x*e1x+e1y*e1y+e1z*e1z
128 s2 = e2x*e2x+e2y*e2y+e2z*e2z
129 suma = sqrt(s1/s2)
130 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
131 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
132 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
133c
134 suma = e1x*e1x+e1y*e1y+e1z*e1z
135 suma = one/max(sqrt(suma),em20)
136 e1x = e1x*suma
137 e1y = e1y*suma
138 e1z = e1z*suma
139c
140 e2x = e3y * e1z - e3z * e1y
141 e2y = e3z * e1x - e3x * e1z
142 e2z = e3x * e1y - e3y * e1x
143
144C Project First axe of the skew
145 vx = skew(1,isk)
146 vy = skew(2,isk)
147 vz = skew(3,isk)
148
149 v =vx*e3x+vy*e3y+vz*e3z
150 vx=vx-v*e3x
151 vy=vy-v*e3y
152 vz=vz-v*e3z
153 v =sqrt(vx*vx+vy*vy+vz*vz)
154
155 vx=vx/max(v,em20)
156 vy=vy/max(v,em20)
157 vz=vz/max(v,em20)
158
159C Cos and Sin calculation
160 vr=vx*e1x+vy*e1y+vz*e1z
161 vs=vx*e2x+vy*e2y+vz*e2z
162c Save data in RTHBUF
163 ithbuf(iad2)=iadr+1
164 rthbuf(iadr+1)=vr
165 rthbuf(iadr+2)=vs
166
167 iadr=iadr+2
168 ENDIF
169 iad=iad+1
170 iad2=iad2+1
171 ENDDO
172 ELSEIF(ityp==7)THEN
173 DO k=1,nne
174 nn=ithbuf(iad)
175 isk=ithbuf(iad2)
176 IF(isk /= 0) THEN
177C Corotational Frame E1 E2 E3
178 n1=ixc(2,nn)
179 n2=ixc(3,nn)
180 n3=ixc(4,nn)
181
182 x1=x(1,n1)
183 x2=x(1,n2)
184 x3=x(1,n3)
185
186 y1=x(2,n1)
187 y2=x(2,n2)
188 y3=x(2,n3)
189
190 z1=x(3,n1)
191 z2=x(3,n2)
192 z3=x(3,n3)
193
194 x21=x2-x1
195 y21=y2-y1
196 z21=z2-z1
197 x31=x3-x1
198 y31=y3-y1
199 z31=z3-z1
200 x32=x3-x2
201 y32=y3-y2
202 z32=z3-z2
203c
204 e1x= x21
205 e1y= y21
206 e1z= z21
207 suma = sqrt(e1x*e1x+e1y*e1y+e1z*e1z)
208 e1x=e1x/suma
209 e1y=e1y/suma
210 e1z=e1z/suma
211c
212 e3x=y31*z32-z31*y32
213 e3y=z31*x32-x31*z32
214 e3z=x31*y32-y31*x32
215 suma = sqrt(e3x*e3x+e3y*e3y+e3z*e3z)
216 e3x=e3x/suma
217 e3y=e3y/suma
218 e3z=e3z/suma
219 area = half * suma
220c
221 e2x=e3y*e1z-e3z*e1y
222 e2y=e3z*e1x-e3x*e1z
223 e2z=e3x*e1y-e3y*e1x
224 suma = sqrt(e2x*e2x+e2y*e2y+e2z*e2z)
225 e2x=e2x/suma
226 e2y=e2y/suma
227 e2z=e2z/suma
228
229C Project First axe of the skew
230 vx = skew(1,isk)
231 vy = skew(2,isk)
232 vz = skew(3,isk)
233
234 v =vx*e3x+vy*e3y+vz*e3z
235 vx=vx-v*e3x
236 vy=vy-v*e3y
237 vz=vz-v*e3z
238 v =sqrt(vx*vx+vy*vy+vz*vz)
239
240 vx=vx/max(v,em20)
241 vy=vy/max(v,em20)
242 vz=vz/max(v,em20)
243C Cos and Sin calculation
244 vr=vx*e1x+vy*e1y+vz*e1z
245 vs=vx*e2x+vy*e2y+vz*e2z
246
247c Save data in RTHBUF
248 ithbuf(iad2)=iadr+1
249 rthbuf(iadr+1)=vr
250 rthbuf(iadr+2)=vs
251 iadr=iadr+2
252 ENDIF
253 iad=iad+1
254 iad2=iad2+1
255 ENDDO
256 ENDIF
257 ENDDO
258 RETURN
259 END
260
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
starter
program starter
Definition starter.F:39
thskewc
subroutine thskewc(rthbuf, ithgrp, ithbuf, x, ixc, ixtg, skew, nthgrp)
Definition thskewc.F:31