OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
s4deri3.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!|| s4deri3 ../engine/source/elements/solid/solide4/s4deri3.F
25!||--- called by ------------------------------------------------------
26!|| s4forc3 ../engine/source/elements/solid/solide4/s4forc3.F
27!||--- calls -----------------------------------------------------
28!|| schkjabt3 ../engine/source/elements/solid/solide4/schkjabt3.F
29!||====================================================================
30 SUBROUTINE s4deri3(
31 1 OFF, DET, NGL, DELTAX,
32 2 MXT, X1, X2, X3,
33 3 X4, Y1, Y2, Y3,
34 4 Y4, Z1, Z2, Z3,
35 5 Z4, PX1, PX2, PX3,
36 6 PX4, PY1, PY2, PY3,
37 7 PY4, PZ1, PZ2, PZ3,
38 8 PZ4, RX, RY, RZ,
39 9 SX, SY, SZ, TX,
40 A TY, TZ, SAV, OFFG,
41 B NEL, PM, VOLDP, ISMSTR,
42 C IFORMDT, JLAG)
43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47#include "comlock.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "mvsiz_p.inc"
52C-----------------------------------------------
53C C o m m o n B l o c k s
54C-----------------------------------------------
55#include "param_c.inc"
56#include "scr17_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER, INTENT(IN) :: JLAG
61 INTEGER, INTENT(IN) :: ISMSTR
62 INTEGER, INTENT(IN) :: IFORMDT
63 INTEGER NEL, MXT(MVSIZ)
64 DOUBLE PRECISION
65 . X1(*), X2(*), X3(*), X4(*),
66 . Y1(*), Y2(*), Y3(*), Y4(*),
67 . Z1(*), Z2(*), Z3(*), Z4(*), SAV(NEL,9),VOLDP(NEL)
68
69 my_real
70 . OFF(*), DET(*), DELTAX(*), PM(NPROPM,*),
71 . PX1(*), PX2(*), PX3(*), PX4(*),
72 . py1(*), py2(*), py3(*), py4(*),
73 . pz1(*), pz2(*), pz3(*), pz4(*),offg(*),
74 . rx(*), ry(*), rz(*), sx(*), sy(*), sz(*),tx(*), ty(*), tz(*)
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER NGL(*), I,J
79 INTEGER NNEGA,INDEX(MVSIZ)
80 DOUBLE PRECISION
81 . x41, y41, z41, x42, y42, z42, x43, y43, z43,b1dp,c1dp,d1dp
82 my_real
83 . a1, a2, a3, a4, d, areamx2,
84 . b1(mvsiz), b2(mvsiz), b3(mvsiz), b4(mvsiz),
85 . c1(mvsiz), c2(mvsiz), c3(mvsiz), c4(mvsiz),
86 . d1(mvsiz), d2(mvsiz), d3(mvsiz), d4(mvsiz),
87 . pxx, pyy, pzz, pxy, pyz, pxz, gfac, aa, bb, p, ld
88C-----------------------------------------------
89C
90 DO i=1,nel
91 x43 = x4(i) - x3(i)
92 y43 = y4(i) - y3(i)
93 z43 = z4(i) - z3(i)
94 x41 = x4(i) - x1(i)
95 y41 = y4(i) - y1(i)
96 z41 = z4(i) - z1(i)
97 x42 = x4(i) - x2(i)
98 y42 = y4(i) - y2(i)
99 z42 = z4(i) - z2(i)
100C
101 rx(i) = -x41
102 ry(i) = -y41
103 rz(i) = -z41
104 sx(i) = -x42
105 sy(i) = -y42
106 sz(i) = -z42
107C
108 tx(i) = -x43
109 ty(i) = -y43
110 tz(i) = -z43
111C
112 b1dp = y43*z42 - y42*z43
113 b1(i) = b1dp
114 b2(i) = y41*z43 - y43*z41
115 b3(i) = y42*z41 - y41*z42
116 b4(i) = -(b1(i) + b2(i) + b3(i))
117C
118 c1dp = z43*x42 - z42*x43
119 c1(i) = c1dp
120 c2(i) = z41*x43 - z43*x41
121 c3(i) = z42*x41 - z41*x42
122 c4(i) = -(c1(i) + c2(i) + c3(i))
123C
124 d1dp = x43*y42 - x42*y43
125 d1(i) = d1dp
126 d2(i) = x41*y43 - x43*y41
127 d3(i) = x42*y41 - x41*y42
128 d4(i) = -(d1(i) + d2(i) + d3(i))
129C
130 voldp(i) = (x41*b1dp + y41*c1dp + z41*d1dp)*one_over_6
131 det(i) = voldp(i)
132C
133 ENDDO
134C
135C
136 CALL schkjabt3(
137 1 off, det, ngl, offg,
138 2 nnega, index, nel, ismstr,
139 3 jlag)
140 IF (nnega>0) THEN
141 IF (ismstr==10.OR.ismstr==12) THEN
142#include "vectorize.inc"
143 DO j=1,nnega
144 i = index(j)
145 x1(i)=sav(i,1)
146 y1(i)=sav(i,4)
147 z1(i)=sav(i,7)
148 x2(i)=sav(i,2)
149 y2(i)=sav(i,5)
150 z2(i)=sav(i,8)
151 x3(i)=sav(i,3)
152 y3(i)=sav(i,6)
153 z3(i)=sav(i,9)
154 x4(i)=zero
155 y4(i)=zero
156 z4(i)=zero
157 ENDDO
158 ELSE
159#include "vectorize.inc"
160 DO j=1,nnega
161 i = index(j)
162 x1(i)=sav(i,1)
163 y1(i)=sav(i,2)
164 z1(i)=sav(i,3)
165 x2(i)=sav(i,4)
166 y2(i)=sav(i,5)
167 z2(i)=sav(i,6)
168 x3(i)=sav(i,7)
169 y3(i)=sav(i,8)
170 z3(i)=sav(i,9)
171 x4(i)=zero
172 y4(i)=zero
173 z4(i)=zero
174 ENDDO
175 END IF
176#include "vectorize.inc"
177 DO j=1,nnega
178 i = index(j)
179 x43 = x4(i) - x3(i)
180 y43 = y4(i) - y3(i)
181 z43 = z4(i) - z3(i)
182 x41 = x4(i) - x1(i)
183 y41 = y4(i) - y1(i)
184 z41 = z4(i) - z1(i)
185 x42 = x4(i) - x2(i)
186 y42 = y4(i) - y2(i)
187 z42 = z4(i) - z2(i)
188C
189 rx(i) = -x41
190 ry(i) = -y41
191 rz(i) = -z41
192 sx(i) = -x42
193 sy(i) = -y42
194 sz(i) = -z42
195 tx(i) = -x43
196 ty(i) = -y43
197 tz(i) = -z43
198C
199C
200 b1dp = y43*z42 - y42*z43
201 b1(i) = b1dp
202 b2(i) = y41*z43 - y43*z41
203 b3(i) = y42*z41 - y41*z42
204 b4(i) = -(b1(i) + b2(i) + b3(i))
205C
206 c1dp = z43*x42 - z42*x43
207 c1(i) = c1dp
208 c2(i) = z41*x43 - z43*x41
209 c3(i) = z42*x41 - z41*x42
210 c4(i) = -(c1(i) + c2(i) + c3(i))
211C
212 d1dp = x43*y42 - x42*y43
213 d1(i) = d1dp
214 d2(i) = x41*y43 - x43*y41
215 d3(i) = x42*y41 - x41*y42
216 d4(i) = -(d1(i) + d2(i) + d3(i))
217C
218 voldp(i) = (x41*b1dp + y41*c1dp + z41*d1dp)*one_over_6
219 det(i) = voldp(i)
220 offg(i) = two
221 ENDDO
222 END IF
223C
224 DO i=1,nel
225 d = one/det(i)/six
226 px1(i)=-b1(i)*d
227 py1(i)=-c1(i)*d
228 pz1(i)=-d1(i)*d
229 px2(i)=-b2(i)*d
230 py2(i)=-c2(i)*d
231 pz2(i)=-d2(i)*d
232 px3(i)=-b3(i)*d
233 py3(i)=-c3(i)*d
234 pz3(i)=-d3(i)*d
235 px4(i)=-b4(i)*d
236 py4(i)=-c4(i)*d
237 pz4(i)=-d4(i)*d
238 END DO
239
240 IF(idt1sol==0)THEN
241
242 DO i=1,nel
243 d = max(px1(i)*px1(i)+py1(i)*py1(i)+pz1(i)*pz1(i),
244 . px2(i)*px2(i)+py2(i)*py2(i)+pz2(i)*pz2(i),
245 . px3(i)*px3(i)+py3(i)*py3(i)+pz3(i)*pz3(i),
246 . px4(i)*px4(i)+py4(i)*py4(i)+pz4(i)*pz4(i))
247 deltax(i) = one / sqrt(d)
248 END DO
249
250 ELSEIF(iformdt==0)THEN
251 DO i=1,nel
252 d = px1(i)*px1(i)+py1(i)*py1(i)+pz1(i)*pz1(i)
253 . + px2(i)*px2(i)+py2(i)*py2(i)+pz2(i)*pz2(i)
254 . + px3(i)*px3(i)+py3(i)*py3(i)+pz3(i)*pz3(i)
255 . + px4(i)*px4(i)+py4(i)*py4(i)+pz4(i)*pz4(i)
256 deltax(i) = one / sqrt(d)
257 END DO
258
259 ELSEIF(iformdt==1)THEN
260
261 gfac=pm(105,mxt(1))
262 ld =two*sqrt(max(one-gfac,zero))+one
263 DO i=1,nel
264 pxx=px1(i)*px1(i)+px2(i)*px2(i)+px3(i)*px3(i)+px4(i)*px4(i)
265 pyy=py1(i)*py1(i)+py2(i)*py2(i)+py3(i)*py3(i)+py4(i)*py4(i)
266 pzz=pz1(i)*pz1(i)+pz2(i)*pz2(i)+pz3(i)*pz3(i)+pz4(i)*pz4(i)
267 pxy=px1(i)*py1(i)+px2(i)*py2(i)+px3(i)*py3(i)+px4(i)*py4(i)
268 pxz=px1(i)*pz1(i)+px2(i)*pz2(i)+px3(i)*pz3(i)+px4(i)*pz4(i)
269 pyz=py1(i)*pz1(i)+py2(i)*pz2(i)+py3(i)*pz3(i)+py4(i)*pz4(i)
270C
271 aa = -(pxx+pyy+pzz)
272 bb = (pxx*pyy+pxx*pzz+pyy*pzz-pxy**2-pxz**2-pyz**2)
273 p = bb-third*aa*aa
274 d = two*sqrt(third*max(-p,zero))-third*aa
275C
276 d=ld*d
277C
278 deltax(i) = one / sqrt(d)
279 END DO
280
281 ELSEIF(iformdt==2)THEN
282
283 gfac=pm(105,mxt(1))
284 DO i=1,nel
285 pxx=px1(i)*px1(i)+px2(i)*px2(i)+px3(i)*px3(i)+px4(i)*px4(i)
286 pyy=py1(i)*py1(i)+py2(i)*py2(i)+py3(i)*py3(i)+py4(i)*py4(i)
287 pzz=pz1(i)*pz1(i)+pz2(i)*pz2(i)+pz3(i)*pz3(i)+pz4(i)*pz4(i)
288 pxy=px1(i)*py1(i)+px2(i)*py2(i)+px3(i)*py3(i)+px4(i)*py4(i)
289 pxz=px1(i)*pz1(i)+px2(i)*pz2(i)+px3(i)*pz3(i)+px4(i)*pz4(i)
290 pyz=py1(i)*pz1(i)+py2(i)*pz2(i)+py3(i)*pz3(i)+py4(i)*pz4(i)
291C
292 aa = -(pxx+pyy+pzz)
293 bb = gfac*(pxx*pyy+pxx*pzz+pyy*pzz-pxy**2-pxz**2-pyz**2)
294 p = bb-third*aa*aa
295 d = two*sqrt(third*max(-p,zero))-third*aa
296C
297 deltax(i) = one / sqrt(d)
298 END DO
299
300 END IF
301C
302 RETURN
303C
304 1000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB',i10/)
305 2000 FORMAT(/' ZERO OR NEGATIVE VOLUME : DELETE 3D-ELEMENT NB',i10/)
306 END
s4deri3
subroutine s4deri3(off, det, ngl, deltax, mxt, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, rx, ry, rz, sx, sy, sz, tx, ty, tz, sav, offg, nel, pm, voldp, ismstr, iformdt, jlag)
Definition s4deri3.F:43
max
#define max(a, b)
Definition macros.h:21
schkjabt3
subroutine schkjabt3(off, det, ngl, offg, nnega, index, nel, ismstr, jlag)
Definition schkjabt3.F:40