OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
gjnt_diff.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!|| gjnt_diff ../engine/source/tools/lagmul/gjnt_diff.f
25!||--- called by ------------------------------------------------------
26!|| lag_gjnt ../engine/source/tools/lagmul/lag_gjnt.F
27!||====================================================================
28 SUBROUTINE gjnt_diff(SK ,L1 ,L2 ,L3 ,ALPHA ,
29 2 IADLL ,LLL ,JLL ,SLL ,XLL ,
30 3 X ,N0 ,N1 ,N2 ,N3 ,
31 4 NC )
32C-----------------------------------------------
33C I m p l i c i t T y p e s
34C-----------------------------------------------
35#include "implicit_f.inc"
36C-----------------------------------------------
37C D u m m y A r g u m e n t s
38C-----------------------------------------------
39 INTEGER N0,N1,N2,N3,NC, LLL(*),JLL(*),SLL(*),IADLL(*)
40 my_real
41 . X(3,*),XLL(*),SK(9),L1(3),L2(3),L3(3),ALPHA
42C-----------------------------------------------
43C L o c a l V a r i a b l e s
44C-----------------------------------------------
45 INTEGER I,II,J,JJ,IK,IC,IAD,INOD(4)
46 my_real
47 . UX(3),UY(3),UZ(3),VX(3),VY(3),VZ(3),WX(3),WY(3),WZ(3),
48 . rx(3),ry(3),rz(3),
49 . x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3,xs,ys,zs,
50 . x12,x22,x32,x42,y12,y22,y32,y42,z12,z22,z32,z42,
51 . xx,yy,zz,xxx,yyy,zzz,xy,yz,zx,xy2,yz2,zx2,
52 . lx,ly,lz,norm,unt,deut,b1,b2,b3,c1,c2,c3,det
53C-----------------------------------------------
54C NC : nombre de condition cinematique
55C IC : numero de la condition cinematique (1,NC)
56C IK :
57C I : numero global du noeud (1,NUMNOD)
58C J : direction 1,2,3,4,5,6
59C------
60C IADLL(NC) : IAD = IADLL(IC)
61C IK = IAD,IAD+1,IAD+2,...
62C LLL(LAG_NKF) : I = LLL(IK)
63C JLL(LAG_NKF) : J = JLL(IK)
64C======================================================================|
65 unt = third
66 deut= two_third
67C--- Node 1 skew
68 ux(1) = sk(1)*l1(1)+sk(4)*l1(2)+sk(7)*l1(3)
69 uy(1) = sk(2)*l1(1)+sk(5)*l1(2)+sk(8)*l1(3)
70 uz(1) = sk(3)*l1(1)+sk(6)*l1(2)+sk(9)*l1(3)
71 norm = one/sqrt(ux(1)*ux(1)+uy(1)*uy(1)+uz(1)*uz(1))
72 ux(1) = ux(1)*norm
73 uy(1) = uy(1)*norm
74 uz(1) = uz(1)*norm
75 IF (abs(ux(1))>zep99) THEN
76 wx(1) =-uz(1)
77 wy(1) = zero
78 wz(1) = ux(1)
79 ELSE
80 wx(1) = zero
81 wy(1) =-uz(1)
82 wz(1) = uy(1)
83 ENDIF
84 norm = one/sqrt(wx(1)*wx(1)+wy(1)*wy(1)+wz(1)*wz(1))
85 wx(1) = wx(1)*norm
86 wy(1) = wy(1)*norm
87 wz(1) = wz(1)*norm
88 vx(1) = wy(1)*uz(1)-wz(1)*uy(1)
89 vy(1) = wz(1)*ux(1)-wx(1)*uz(1)
90 vz(1) = wx(1)*uy(1)-wy(1)*ux(1)
91C--- Node 2 skew
92 ux(2) = sk(1)*l2(1)+sk(4)*l2(2)+sk(7)*l2(3)
93 uy(2) = sk(2)*l2(1)+sk(5)*l2(2)+sk(8)*l2(3)
94 uz(2) = sk(3)*l2(1)+sk(6)*l2(2)+sk(9)*l2(3)
95 norm = one/sqrt(ux(2)*ux(2)+uy(2)*uy(2)+uz(2)*uz(2))
96 ux(2) = ux(2)*norm
97 uy(2) = uy(2)*norm
98 uz(2) = uz(2)*norm
99 IF (abs(ux(2))>zep99) THEN
100 wx(2) =-uz(2)
101 wy(2) = zero
102 wz(2) = ux(2)
103 ELSE
104 wx(2) = zero
105 wy(2) =-uz(2)
106 wz(2) = uy(2)
107 ENDIF
108 norm = one/sqrt(wx(1)*wx(1)+wy(1)*wy(1)+wz(1)*wz(1))
109 wx(1) = wx(1)*norm
110 wy(1) = wy(1)*norm
111 wz(1) = wz(1)*norm
112 vx(2) = wy(2)*uz(2)-wz(2)*uy(2)
113 vy(2) = wz(2)*ux(2)-wx(2)*uz(2)
114 vz(2) = wx(2)*uy(2)-wy(2)*ux(2)
115C--- Node 3 skew
116 ux(3) = sk(1)*l3(1)+sk(4)*l3(2)+sk(7)*l3(3)
117 uy(3) = sk(2)*l3(1)+sk(5)*l3(2)+sk(8)*l3(3)
118 uz(3) = sk(3)*l3(1)+sk(6)*l3(2)+sk(9)*l3(3)
119 norm = one/sqrt(ux(3)*ux(3)+uy(3)*uy(3)+uz(3)*uz(3))
120 ux(3) = ux(3)*norm
121 uy(3) = uy(3)*norm
122 uz(3) = uz(3)*norm
123 IF (abs(ux(3))>zep99) THEN
124 wx(3) =-uz(3)
125 wy(3) = zero
126 wz(3) = ux(3)
127 ELSE
128 wx(3) = zero
129 wy(3) =-uz(3)
130 wz(3) = uy(3)
131 ENDIF
132 norm = one/sqrt(wx(1)*wx(1)+wy(1)*wy(1)+wz(1)*wz(1))
133 wx(1) = wx(1)*norm
134 wy(1) = wy(1)*norm
135 wz(1) = wz(1)*norm
136 vx(3) = wy(3)*uz(3)-wz(3)*uy(3)
137 vy(3) = wz(3)*ux(3)-wx(3)*uz(3)
138 vz(3) = wx(3)*uy(3)-wy(3)*ux(3)
139C---------------------------
140 x1=x(1,n1)
141 y1=x(2,n1)
142 z1=x(3,n1)
143 x2=x(1,n2)
144 y2=x(2,n2)
145 z2=x(3,n2)
146 x3=x(1,n3)
147 y3=x(2,n3)
148 z3=x(3,n3)
149 x0=(x1+x2+x3)*third
150 y0=(y1+y2+y3)*third
151 z0=(z1+z2+z3)*third
152C---------------------------
153 inod(1) = n1
154 inod(2) = n2
155 inod(3) = n3
156 inod(4) = n0
157C---------------------------
158 x1=x1-x0
159 y1=y1-y0
160 z1=z1-z0
161 x2=x2-x0
162 y2=y2-y0
163 z2=z2-z0
164 x3=x3-x0
165 y3=y3-y0
166 z3=z3-z0
167 xs=x(1,n0)-x0
168 ys=x(2,n0)-y0
169 zs=x(3,n0)-z0
170 rx(1) = x1
171 rx(2) = x2
172 rx(3) = x3
173 ry(1) = y1
174 ry(2) = y2
175 ry(3) = y3
176 rz(1) = z1
177 rz(2) = z2
178 rz(3) = z3
179C---------------------------
180 x12=x1*x1
181 x22=x2*x2
182 x32=x3*x3
183 y12=y1*y1
184 y22=y2*y2
185 y32=y3*y3
186 z12=z1*z1
187 z22=z2*z2
188 z32=z3*z3
189 xx=x12 + x22 + x32
190 yy=y12 + y22 + y32
191 zz=z12 + z22 + z32
192 xy=x1*y1 + x2*y2 + x3*y3
193 yz=y1*z1 + y2*z2 + y3*z3
194 zx=z1*x1 + z2*x2 + z3*x3
195 zzz=xx+yy
196 xxx=yy+zz
197 yyy=zz+xx
198 xy2=xy*xy
199 yz2=yz*yz
200 zx2=zx*zx
201 det=xxx*yyy*zzz - xxx*yz2 - yyy*zx2 - zzz*xy2 - 2.*xy*yz*zx
202 det=one/det
203 b1=zzz*yyy-yz2
204 b2=xxx*zzz-zx2
205 b3=yyy*xxx-xy2
206 c3=zzz*xy+yz*zx
207 c1=xxx*yz+zx*xy
208 c2=yyy*zx+xy*yz
209C---------------------------
210C --- (vx)
211 nc = nc + 1
212 iadll(nc+1) = iadll(nc) + 10
213 iad = iadll(nc) -1
214 DO jj=1,3
215 ik = iad+jj
216 lll(ik) = inod(jj)
217 jll(ik) = 1
218 sll(ik) = 0
219 xll(ik) = fourth
220 . + det*zs*(b2*rz(jj) - c1*ry(jj))
221 . - det*ys*(c1*rz(jj) - b3*ry(jj))
222 ENDDO
223 iad = iad + 3
224 DO jj=1,3
225 ik = iad+jj
226 lll(ik) = inod(jj)
227 jll(ik) = 2
228 sll(ik) = 0
229 xll(ik) = det*zs*(c1*rx(jj) - c3*rz(jj))
230 . - det*ys*(b3*rx(jj) - c2*rz(jj))
231 ENDDO
232 iad = iad + 3
233 DO jj=1,3
234 ik = iad+jj
235 lll(ik) = inod(jj)
236 jll(ik) = 3
237 sll(ik) = 0
238 xll(ik) = det*zs*(c3*ry(jj) - b2*rx(jj))
239 . - det*ys*(c2*ry(jj) - c1*rx(jj))
240 ENDDO
241 ik = iad + 4
242 lll(ik) = inod(4)
243 jll(ik) = 1
244 sll(ik) = 0
245 xll(ik) = -1.0
246C --- (vy)
247 nc = nc + 1
248 iadll(nc+1) = iadll(nc) + 10
249 iad = iadll(nc) -1
250 DO jj=1,3
251 ik = iad+jj
252 lll(ik) = inod(jj)
253 jll(ik) = 1
254 sll(ik) = 0
255 xll(ik) = det*xs*(c1*rz(jj) - b3*ry(jj))
256 . - det*zs*(c3*rz(jj) - c2*ry(jj))
257 ENDDO
258 iad = iad + 3
259 DO jj=1,3
260 ik = iad+jj
261 lll(ik) = inod(jj)
262 jll(ik) = 2
263 sll(ik) = 0
264 xll(ik) = fourth
265 . + det*xs*(b3*rx(jj) - c2*rz(jj))
266 . - det*zs*(c2*rx(jj) - b1*rz(jj))
267 ENDDO
268 iad = iad + 3
269 DO jj=1,3
270 ik = iad+jj
271 lll(ik) = inod(jj)
272 jll(ik) = 3
273 sll(ik) = 0
274 xll(ik) = det*xs*(c2*ry(jj) - c1*rx(jj))
275 . - det*zs*(b1*ry(jj) - c3*rx(jj))
276 ENDDO
277 ik = iad + 4
278 lll(ik) = inod(4)
279 jll(ik) = 2
280 sll(ik) = 0
281 xll(ik) = -one
282C --- (vz)
283 nc = nc + 1
284 iadll(nc+1) = iadll(nc) + 10
285 iad = iadll(nc) -1
286 DO jj=1,3
287 ik = iad+jj
288 lll(ik) = inod(jj)
289 jll(ik) = 1
290 sll(ik) = 0
291 xll(ik) = det*ys*(c3*rz(jj) - c2*ry(jj))
292 . - det*xs*(b2*rz(jj) - c1*ry(jj))
293 ENDDO
294 iad = iad + 3
295 DO jj=1,3
296 ik = iad+jj
297 lll(ik) = inod(jj)
298 jll(ik) = 2
299 sll(ik) = 0
300 xll(ik) = det*ys*(c2*rx(jj) - b1*rz(jj))
301 . - det*xs*(c1*rx(jj) - c3*rz(jj))
302 ENDDO
303 iad = iad + 3
304 DO jj=1,3
305 ik = iad+jj
306 lll(ik) = inod(jj)
307 jll(ik) = 3
308 sll(ik) = 0
309 xll(ik) = fourth
310 . + det*ys*(b1*ry(jj) - c3*rx(jj))
311 . - det*xs*(c3*ry(jj) - b2*rx(jj))
312 ENDDO
313 ik = iad + 4
314 lll(ik) = inod(4)
315 jll(ik) = 3
316 sll(ik) = 0
3172 xll(ik) = -one
318C
319c rotational dof of secnd
320C --- (wx)
321 nc = nc + 1
322 iadll(nc+1) = iadll(nc) + 10
323 iad = iadll(nc) -1
324 DO jj=1,3
325 ik = iad+jj
326 lll(ik) = inod(jj)
327 jll(ik) = 1
328 sll(ik) = 0
329 xll(ik) = det*(c3*rz(jj) - c2*ry(jj))
330 ENDDO
331 iad = iad + 3
332 DO jj=1,3
333 ik = iad+jj
334 lll(ik) = inod(jj)
335 jll(ik) = 2
336 sll(ik) = 0
337 xll(ik) = det*(c2*rx(jj) - b1*rz(jj))
338 ENDDO
339 iad = iad + 3
340 DO jj=1,3
341 ik = iad+jj
342 lll(ik) = inod(jj)
343 jll(ik) = 3
344 sll(ik) = 0
345 xll(ik) = det*(b1*ry(jj) - c3*rx(jj))
346 ENDDO
347 ik = iad + 4
348 lll(ik) = inod(4)
349 jll(ik) = 4
350 sll(ik) = 0
351 xll(ik) = -1.0
352C --- (wy)
353 nc = nc + 1
354 iadll(nc+1) = iadll(nc) + 10
355 iad = iadll(nc) -1
356 DO jj=1,3
357 ik = iad+jj
358 lll(ik) = inod(jj)
359 jll(ik) = 1
360 sll(ik) = 0
361 xll(ik) = det*(b2*rz(jj) - c1*ry(jj))
362 ENDDO
363 iad = iad + 3
364 DO jj=1,3
365 ik = iad+jj
366 lll(ik) = inod(jj)
367 jll(ik) = 2
368 sll(ik) = 0
369 xll(ik) = det*(c1*rx(jj) - c3*rz(jj))
370 ENDDO
371 iad = iad + 3
372 DO jj=1,3
373 ik = iad+jj
374 lll(ik) = inod(jj)
375 jll(ik) = 3
376 sll(ik) = 0
377 xll(ik) = det*(c3*ry(jj) - b2*rx(jj))
378 ENDDO
379 ik = iad + 4
380 lll(ik) = inod(4)
381 jll(ik) = 5
382 sll(ik) = 0
383 xll(ik) = -1.0
384C --- (wz)
385 nc = nc + 1
386 iadll(nc+1) = iadll(nc) + 10
387 iad = iadll(nc) -1
388 DO jj=1,3
389 ik = iad+jj
390 lll(ik) = inod(jj)
391 jll(ik) = 1
392 sll(ik) = 0
393 xll(ik) = det*(c1*rz(jj) - b3*ry(jj))
394 ENDDO
395 iad = iad + 3
396 DO jj=1,3
397 ik = iad+jj
398 lll(ik) = inod(jj)
399 jll(ik) = 2
400 sll(ik) = 0
401 xll(ik) = det*(b3*rx(jj) - c2*rz(jj))
402 ENDDO
403 iad = iad + 3
404 DO jj=1,3
405 ik = iad+jj
406 lll(ik) = inod(jj)
407 jll(ik) = 3
408 sll(ik) = 0
409 xll(ik) = det*(c2*ry(jj) - c1*rx(jj))
410 ENDDO
411 ik = iad + 4
412 lll(ik) = inod(4)
413 jll(ik) = 6
414 sll(ik) = 0
415 xll(ik) = -one
416C
417C --- Constraints
418C X local
419C
420 inod(4) = n0
421C
422 nc = nc + 1
423 iadll(nc+1) = iadll(nc) + 12
424 ik = iadll(nc)
425 lll(ik) = n1
426 jll(ik) = 4
427 sll(ik) = 0
428 xll(ik) = alpha*ux(1)
429 ik = ik+1
430 lll(ik) = n1
431 jll(ik) = 5
432 sll(ik) = 0
433 xll(ik) = alpha*uy(1)
434 ik = ik+1
435 lll(ik) = n1
436 jll(ik) = 6
437 sll(ik) = 0
438 xll(ik) = alpha*uz(1)
439 DO i=2,3
440 ik = ik+1
441 lll(ik) = inod(i)
442 jll(ik) = 4
443 sll(ik) = 0
444 xll(ik) = ux(i)
445 ik = ik+1
446 lll(ik) = inod(i)
447 jll(ik) = 5
448 sll(ik) = 0
449 xll(ik) = uy(i)
450 ik = ik+1
451 lll(ik) = inod(i)
452 jll(ik) = 6
453 sll(ik) = 0
454 xll(ik) = uz(i)
455 ENDDO
456 ik = ik+1
457 lll(ik) = n0
458 jll(ik) = 4
459 sll(ik) = 0
460 xll(ik) =-alpha*ux(1)-ux(2)-ux(3)
461 ik = ik+1
462 lll(ik) = n0
463 jll(ik) = 5
464 sll(ik) = 0
465 xll(ik) =-alpha*uy(1)-uy(2)-uy(3)
466 ik = ik+1
467 lll(ik) = n0
468 jll(ik) = 6
469 sll(ik) = 0
470 xll(ik) =-alpha*uz(1)-uz(2)-uz(3)
471C
472C Y local
473 DO i=1,3
474 nc = nc + 1
475 iadll(nc+1) = iadll(nc) + 6
476 ik = iadll(nc)
477 lll(ik) = inod(i)
478 jll(ik) = 4
479 sll(ik) = 0
480 xll(ik) = vx(i)
481 ik = ik+1
482 lll(ik) = inod(i)
483 jll(ik) = 5
484 sll(ik) = 0
485 xll(ik) = vy(i)
486 ik = ik+1
487 lll(ik) = inod(i)
488 jll(ik) = 6
489 sll(ik) = 0
490 xll(ik) = vz(i)
491 ik = ik+1
492 lll(ik) = n0
493 jll(ik) = 4
494 sll(ik) = 0
495 xll(ik) =-vx(i)
496 ik = ik+1
497 lll(ik) = n0
498 jll(ik) = 5
499 sll(ik) = 0
500 xll(ik) =-vy(i)
501 ik = ik+1
502 lll(ik) = n0
503 jll(ik) = 6
504 sll(ik) = 0
505 xll(ik) =-vz(i)
506 ENDDO
507C
508C Z local
509 DO i=1,3
510 nc = nc + 1
511 iadll(nc+1) = iadll(nc) + 6
512 ik = iadll(nc)
513 lll(ik) = inod(i)
514 jll(ik) = 4
515 sll(ik) = 0
516 xll(ik) = wx(i)
517 ik = ik+1
518 lll(ik) = inod(i)
519 jll(ik) = 5
520 sll(ik) = 0
521 xll(ik) = wy(i)
522 ik = ik+1
523 lll(ik) = inod(i)
524 jll(ik) = 6
525 sll(ik) = 0
526 xll(ik) = wz(i)
527 ik = ik+1
528 lll(ik) = n0
529 jll(ik) = 4
530 sll(ik) = 0
531 xll(ik) =-wx(i)
532 ik = ik+1
533 lll(ik) = n0
534 jll(ik) = 5
535 sll(ik) = 0
536 xll(ik) =-wy(i)
537 ik = ik+1
538 lll(ik) = n0
539 jll(ik) = 6
540 sll(ik) = 0
541 xll(ik) =-wz(i)
542 ENDDO
543C---
544 RETURN
545 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
gjnt_diff
subroutine gjnt_diff(sk, l1, l2, l3, alpha, iadll, lll, jll, sll, xll, x, n0, n1, n2, n3, nc)
Definition gjnt_diff.F:32