OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i7curv.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!|| i7norm ../engine/source/interfaces/int07/i7curv.F
25!||--- called by ------------------------------------------------------
26!|| i7mainf ../engine/source/interfaces/int07/i7mainf.F
27!||====================================================================
28 SUBROUTINE i7norm(NRTM,IRECT,NUMNOD,X,NOD_NORMAL,
29 . NMN ,MSR )
30C-----------------------------------------------
31C I m p l i c i t T y p e s
32C-----------------------------------------------
33#include "implicit_f.inc"
34C-----------------------------------------------
35C D u m m y A r g u m e n t s
36C-----------------------------------------------
37 INTEGER NRTM,NUMNOD,IRECT(4,NRTM),NMN,MSR(*)
38C REAL
39 my_real
40 . x(3,numnod), nod_normal(3,numnod)
41C-----------------------------------------------
42C L o c a l V a r i a b l e s
43C-----------------------------------------------
44 INTEGER I ,N1,N2,N3,N4
45 my_real
46 . surfx,surfy,surfz,x13,y13,z13,x24,y24,z24
47C-----------------------------------------------
48
49C Optimizable in SPMD If Addition Flag for comm, SPMD_exchange_n
50 DO n1=1,numnod
51 nod_normal(1,n1) = zero
52 nod_normal(2,n1) = zero
53 nod_normal(3,n1) = zero
54 END DO
55
56 DO i=1,nrtm
57 n1 = irect(1,i)
58 n2 = irect(2,i)
59 n3 = irect(3,i)
60 n4 = irect(4,i)
61
62 x13 = x(1,n3) - x(1,n1)
63 y13 = x(2,n3) - x(2,n1)
64 z13 = x(3,n3) - x(3,n1)
65
66 x24 = x(1,n4) - x(1,n2)
67 y24 = x(2,n4) - x(2,n2)
68 z24 = x(3,n4) - x(3,n2)
69
70 surfx = y13*z24 - z13*y24
71 surfy = z13*x24 - x13*z24
72 surfz = x13*y24 - y13*x24
73
74 nod_normal(1,n1) = nod_normal(1,n1) + surfx
75 nod_normal(2,n1) = nod_normal(2,n1) + surfy
76 nod_normal(3,n1) = nod_normal(3,n1) + surfz
77 nod_normal(1,n2) = nod_normal(1,n2) + surfx
78 nod_normal(2,n2) = nod_normal(2,n2) + surfy
79 nod_normal(3,n2) = nod_normal(3,n2) + surfz
80 nod_normal(1,n3) = nod_normal(1,n3) + surfx
81 nod_normal(2,n3) = nod_normal(2,n3) + surfy
82 nod_normal(3,n3) = nod_normal(3,n3) + surfz
83 nod_normal(1,n4) = nod_normal(1,n4) + surfx
84 nod_normal(2,n4) = nod_normal(2,n4) + surfy
85 nod_normal(3,n4) = nod_normal(3,n4) + surfz
86 ENDDO
87
88 RETURN
89 END
90C
91!||====================================================================
92!|| i7normp ../engine/source/interfaces/int07/i7curv.F
93!||--- called by ------------------------------------------------------
94!|| i7mainf ../engine/source/interfaces/int07/i7mainf.F
95!||--- calls -----------------------------------------------------
96!|| myqsort ../common_source/tools/sort/myqsort.F
97!|| spmd_i7curvcom ../engine/source/mpi/interfaces/spmd_i7curvcom.F
98!||====================================================================
99 SUBROUTINE i7normp(NRTM ,IRECT ,NUMNOD ,X ,NOD_NORMAL,
100 . NMN ,MSR ,LENT ,MAXCC,ISDSIZ ,
101 . IRCSIZ,IAD_ELEM,FR_ELEM,ITAG )
102C-----------------------------------------------
103C I m p l i c i t T y p e s
104C-----------------------------------------------
105#include "implicit_f.inc"
106C-----------------------------------------------
107C C o m m o n B l o c k s
108C-----------------------------------------------
109#include "com01_c.inc"
110C-----------------------------------------------
111C D u m m y A r g u m e n t s
112C-----------------------------------------------
113 INTEGER NRTM,NUMNOD,NMN,MAXCC,LENT,
114 . IRECT(4,NRTM),MSR(*),
115 . IAD_ELEM(2,*),FR_ELEM(*),ISDSIZ(*),IRCSIZ(*),ITAG(*)
116C REAL
117 my_real
118 . x(3,numnod), nod_normal(3,numnod)
119C-----------------------------------------------
120C L o c a l V a r i a b l e s
121C-----------------------------------------------
122 INTEGER I ,J ,N1,N2,N3,N4, IAD, LENR, LENS, CC, ERROR,
123 . ADSKYT(0:NUMNOD+1)
124 my_real
125 . surfx,surfy,surfz,x13,y13,z13,x24,y24,z24,
126 . fskyt(3,lent), fskyt2(maxcc)
127 INTEGER :: PERM(MAXCC)
128C-----------------------------------------------
129 ADSKYT(0) = 1
130 adskyt(1) = 1
131 DO n1=1,numnod
132 adskyt(n1+1) = adskyt(n1)+itag(n1)
133 itag(n1) = adskyt(n1)
134 nod_normal(1,n1) = zero
135 nod_normal(2,n1) = zero
136 nod_normal(3,n1) = zero
137 END DO
138
139 DO i=1,nrtm
140 n1 = irect(1,i)
141 n2 = irect(2,i)
142 n3 = irect(3,i)
143 n4 = irect(4,i)
144
145 x13 = x(1,n3) - x(1,n1)
146 y13 = x(2,n3) - x(2,n1)
147 z13 = x(3,n3) - x(3,n1)
148
149 x24 = x(1,n4) - x(1,n2)
150 y24 = x(2,n4) - x(2,n2)
151 z24 = x(3,n4) - x(3,n2)
152
153 surfx = y13*z24 - z13*y24
154 surfy = z13*x24 - x13*z24
155 surfz = x13*y24 - y13*x24
156
157 iad = adskyt(n1)
158 adskyt(n1) = adskyt(n1)+1
159 fskyt(1,iad) = surfx
160 fskyt(2,iad) = surfy
161 fskyt(3,iad) = surfz
162 iad = adskyt(n2)
163 adskyt(n2) = adskyt(n2)+1
164 fskyt(1,iad) = surfx
165 fskyt(2,iad) = surfy
166 fskyt(3,iad) = surfz
167 iad = adskyt(n3)
168 adskyt(n3) = adskyt(n3)+1
169 fskyt(1,iad) = surfx
170 fskyt(2,iad) = surfy
171 fskyt(3,iad) = surfz
172 iad = adskyt(n4)
173 adskyt(n4) = adskyt(n4)+1
174 fskyt(1,iad) = surfx
175 fskyt(2,iad) = surfy
176 fskyt(3,iad) = surfz
177 END DO
178C
179 IF(nspmd>1) THEN
180 lenr = ircsiz(nspmd+1)*3+iad_elem(1,nspmd+1)-iad_elem(1,1)
181 lens = isdsiz(nspmd+1)*3+iad_elem(1,nspmd+1)-iad_elem(1,1)
182 CALL spmd_i7curvcom(iad_elem,fr_elem,adskyt,fskyt,
183 . isdsiz,ircsiz,itag ,lenr ,lens )
184 END IF
185C
186C sorting normals by packet
187C
188 DO n1 = 1, numnod
189 n2 = adskyt(n1-1)
190 n3 = adskyt(n1)-1
191 n4 = n3-n2+1
192 IF(n4>1)THEN ! cas N contribution => tri
193 DO j = 1, 3
194 DO cc = n2, n3
195 fskyt2(cc-n2+1) = fskyt(j,cc)
196 END DO
197C IF(N4>MAXCC)print*,'error cc:',n4,maxcc
198 CALL myqsort(n4,fskyt2,perm,error)
199 DO cc = n2, n3
200 nod_normal(j,n1) = nod_normal(j,n1) + fskyt2(cc-n2+1)
201 END DO
202 END DO
203 ELSEIF(n4==1)THEN ! Case 1 single contribution => direct
204 nod_normal(1,n1) = fskyt(1,n2)
205 nod_normal(2,n1) = fskyt(2,n2)
206 nod_normal(3,n1) = fskyt(3,n2)
207 END IF
208 END DO
209C
210 RETURN
211 END
212C
213!||====================================================================
214!|| i7curvsz ../engine/source/interfaces/int07/i7curv.F
215!||--- called by ------------------------------------------------------
216!|| i7mainf ../engine/source/interfaces/int07/i7mainf.F
217!||====================================================================
218 SUBROUTINE i7curvsz(NRTM,IRECT,NUMNOD,ITAG,LENT,MAXCC)
219C-----------------------------------------------
220C I m p l i c i t T y p e s
221C-----------------------------------------------
222#include "implicit_f.inc"
223C-----------------------------------------------
224C D u m m y A r g u m e n t s
225C-----------------------------------------------
226 INTEGER NRTM,NUMNOD,LENT,MAXCC,
227 . IRECT(4,NRTM),ITAG(*)
228C REAL
229C-----------------------------------------------
230C L o c a l V a r i a b l e s
231C-----------------------------------------------
232 INTEGER I ,N1, N2, N3, N4
233C-----------------------------------------------
234C
235 DO n1=1,numnod
236 itag(n1) = 0
237 END DO
238
239 DO i=1,nrtm
240 n1 = irect(1,i)
241 n2 = irect(2,i)
242 n3 = irect(3,i)
243 n4 = irect(4,i)
244 itag(n1) = itag(n1) + 1
245 itag(n2) = itag(n2) + 1
246 itag(n3) = itag(n3) + 1
247 itag(n4) = itag(n4) + 1
248 END DO
249C
250 lent = 0
251 maxcc = 0
252 DO n1=1,numnod
253 lent = lent + itag(n1)
254 maxcc = max(maxcc,itag(n1))
255 END DO
256C
257 RETURN
258 END
259
260!||====================================================================
261!|| i7cmaj ../engine/source/interfaces/int07/i7curv.f
262!||--- called by ------------------------------------------------------
263!|| i7dst3 ../engine/source/interfaces/int07/i7dst3.F
264!||====================================================================
265 SUBROUTINE i7cmaj(JLT ,CMAJ ,IRECT ,NOD_NORMAL,CAND_E,
266 2 X1 ,X2 ,X3 ,X4 ,
267 3 Y1 ,Y2 ,Y3 ,Y4 ,
268 4 Z1 ,Z2 ,Z3 ,Z4 ,
269 5 NNX1 ,NNX2 ,NNX3 ,NNX4 ,
270 6 NNY1 ,NNY2 ,NNY3 ,NNY4 ,
271 7 NNZ1 ,NNZ2 ,NNZ3 ,NNZ4 )
272c
273c + tmaj
274c / \
275c tmax/__ \
276c / \_\
277c 0-------0---> r
278c
279c
280c borne max:
281c (r+1)t1' = (r-1)t'2
282c r = - (t1' + t'2)/(t1'-t'2)
283c r = e
284c tmaj = (e+1)t1' = (e-1)t'2
285c => en 2D t = t'max => Z = (Lmax/2)t'max
286c
287c Vectors: Z, L, N
288c Z^2 = (L/2.N)^2(L/2)^2 / (N^2(L/2)^2 - (L/2.N)^2)
289c Z^2 = (L.N)^2 L^2 / 4(N^2 L^2 - (L.N)^2)
290c
291C-----------------------------------------------
292C I m p l i c i t T y p e s
293C-----------------------------------------------
294#include "implicit_f.inc"
295C-----------------------------------------------
296C G l o b a l P a r a m e t e r s
297C-----------------------------------------------
298#include "mvsiz_p.inc"
299C-----------------------------------------------
300C D u m m y A r g u m e n t s
301C-----------------------------------------------
302 INTEGER JLT ,IRECT(4,*),CAND_E(*)
303C REAL
304 my_real
305 . NNX1(MVSIZ), NNX2(MVSIZ), NNX3(MVSIZ), NNX4(MVSIZ),
306 . NNY1(MVSIZ), NNY2(MVSIZ), NNY3(MVSIZ), NNY4(MVSIZ),
307 . NNZ1(MVSIZ), NNZ2(MVSIZ), NNZ3(MVSIZ), NNZ4(MVSIZ),
308 . X1(MVSIZ), X2(MVSIZ), X3(MVSIZ), X4(MVSIZ),
309 . Y1(MVSIZ), Y2(MVSIZ), Y3(MVSIZ), Y4(MVSIZ),
310 . Z1(MVSIZ), Z2(MVSIZ), Z3(MVSIZ), Z4(MVSIZ),
311 . cmaj(mvsiz),nod_normal(3,*)
312C-----------------------------------------------
313C L o c a l V a r i a b l e s
314C-----------------------------------------------
315 INTEGER I , IX,L
316 my_real
317 . X12,Y12,Z12,X23,Y23,Z23,X34,Y34,Z34,X41,Y41,Z41,
318 . L12L12,L23L23,L34L34,L41L41,N1N1,N2N2,N3N3,N4N4,NL,
319 . N1L12N1L12,N2L12N2L12,N2L23N2L23,N3L23N3L23,
320 . N3L34N3L34,N4L34N4L34,N4L41N4L41,N1L41N1L41,AAA
321C-----------------------------------------------
322
323 DO I=1,jlt
324 l = cand_e(i)
325
326 ix=irect(1,l)
327 nnx1(i)=nod_normal(1,ix)
328 nny1(i)=nod_normal(2,ix)
329 nnz1(i)=nod_normal(3,ix)
330
331 ix=irect(2,l)
332 nnx2(i)=nod_normal(1,ix)
333 nny2(i)=nod_normal(2,ix)
334 nnz2(i)=nod_normal(3,ix)
335
336 ix=irect(3,l)
337 nnx3(i)=nod_normal(1,ix)
338 nny3(i)=nod_normal(2,ix)
339 nnz3(i)=nod_normal(3,ix)
340
341 ix=irect(4,l)
342 nnx4(i)=nod_normal(1,ix)
343 nny4(i)=nod_normal(2,ix)
344 nnz4(i)=nod_normal(3,ix)
345
346 x12 = x2(i) - x1(i)
347 y12 = y2(i) - y1(i)
348 z12 = z2(i) - z1(i)
349
350 x23 = x3(i) - x2(i)
351 y23 = y3(i) - y2(i)
352 z23 = z3(i) - z2(i)
353
354 x34 = x4(i) - x3(i)
355 y34 = y4(i) - y3(i)
356 z34 = z4(i) - z3(i)
357
358 x41 = x1(i) - x4(i)
359 y41 = y1(i) - y4(i)
360 z41 = z1(i) - z4(i)
361
362 l12l12 = x12*x12 + y12*y12 + z12*z12
363 l23l23 = x23*x23 + y23*y23 + z23*z23
364 l34l34 = x34*x34 + y34*y34 + z34*z34
365 l41l41 = x41*x41 + y41*y41 + z41*z41
366
367 n1n1 = nnx1(i)*nnx1(i)
368 + + nny1(i)*nny1(i)
369 + + nnz1(i)*nnz1(i)
370 n2n2 = nnx2(i)*nnx2(i)
371 + + nny2(i)*nny2(i)
372 + + nnz2(i)*nnz2(i)
373 n3n3 = nnx3(i)*nnx3(i)
374 + + nny3(i)*nny3(i)
375 + + nnz3(i)*nnz3(i)
376 n4n4 = nnx4(i)*nnx4(i)
377 + + nny4(i)*nny4(i)
378 + + nnz4(i)*nnz4(i)
379
380 nl = nnx1(i)*x12
381 + + nny1(i)*y12
382 + + nnz1(i)*z12
383 n1l12n1l12 = nl*nl
384 nl = nnx2(i)*x12
385 + + nny2(i)*y12
386 + + nnz2(i)*z12
387 n2l12n2l12 = nl*nl
388 nl = nnx2(i)*x23
389 + + nny2(i)*y23
390 + + nnz2(i)*z23
391 n2l23n2l23 = nl*nl
392 nl = nnx3(i)*x23
393 + + nny3(i)*y23
394 + + nnz3(i)*z23
395 n3l23n3l23 = nl*nl
396 nl = nnx3(i)*x34
397 + + nny3(i)*y34
398 + + nnz3(i)*z34
399 n3l34n3l34 = nl*nl
400 nl = nnx4(i)*x34
401 + + nny4(i)*y34
402 + + nnz4(i)*z34
403 n4l34n4l34 = nl*nl
404 nl = nnx4(i)*x41
405 + + nny4(i)*y41
406 + + nnz4(i)*z41
407 n4l41n4l41 = nl*nl
408 nl = nnx1(i)*x41
409 + + nny1(i)*y41
410 + + nnz1(i)*z41
411 n1l41n1l41 = nl*nl
412
413c Z^2 = (L.N)^2 L^2 / 4(N^2 L^2 - (L.N)^2)
414 aaa = max(n1l12n1l12*l12l12 / (n1n1*l12l12 - n1l12n1l12),
415 . n2l12n2l12*l12l12 / (n2n2*l12l12 - n2l12n2l12),
416 . n2l23n2l23*l23l23 / (n2n2*l23l23 - n2l23n2l23),
417 . n3l23n3l23*l23l23 / (n3n3*l23l23 - n3l23n3l23),
418 . n3l34n3l34*l34l34 / (n3n3*l34l34 - n3l34n3l34),
419 . n4l34n4l34*l34l34 / (n4n4*l34l34 - n4l34n4l34),
420 . n4l41n4l41*l41l41 / (n4n4*l41l41 - n4l41n4l41),
421 . n1l41n1l41*l41l41 / (n1n1*l41l41 - n1l41n1l41))
422 cmaj(i) = half*sqrt(aaa)
423
424 ENDDO
425
426 RETURN
427 END
428!||====================================================================
429!|| i7curv ../engine/source/interfaces/int07/i7curv.F
430!||--- called by ------------------------------------------------------
431!|| i20for3 ../engine/source/interfaces/int20/i20for3.F
432!|| i7for3 ../engine/source/interfaces/int07/i7for3.F
433!||--- calls -----------------------------------------------------
434!|| i7cubic ../engine/source/interfaces/int07/i7curv.F
435!||====================================================================
436 SUBROUTINE i7curv(JLT ,PENE ,N1 ,N2 ,
437 1 N3 ,GAPV ,X ,NOD_NORMAL,
438 2 IX1 ,IX2 ,IX3 ,IX4 ,
439 3 H1 ,H2 ,H3 ,H4 ,
440 4 X1 ,X2 ,X3 ,X4 ,Y1 ,
441 5 Y2 ,Y3 ,Y4 ,Z1 ,Z2 ,
442 6 Z3 ,Z4 ,XI ,YI ,ZI )
443C-----------------------------------------------
444C I m p l i c i t T y p e s
445C-----------------------------------------------
446#include "implicit_f.inc"
447C-----------------------------------------------
448C G l o b a l P a r a m e t e r s
449C-----------------------------------------------
450#include "mvsiz_p.inc"
451C-----------------------------------------------
452C D u m m y A r g u m e n t s
453C-----------------------------------------------
454 INTEGER JLT,
455 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ)
456C REAL
457 my_real
458 . X(3,*),NOD_NORMAL(3,*),
459 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
460 . PENE(MVSIZ),GAPV(MVSIZ),N1(MVSIZ),N2(MVSIZ),N3(MVSIZ),
461 . X1(MVSIZ), X2(MVSIZ), X3(MVSIZ), X4(MVSIZ),
462 . Y1(MVSIZ), Y2(MVSIZ), Y3(MVSIZ), Y4(MVSIZ),
463 . Z1(MVSIZ), Z2(MVSIZ), Z3(MVSIZ), Z4(MVSIZ),
464 . XI(MVSIZ), YI(MVSIZ), ZI(MVSIZ)
465C-----------------------------------------------
466C L o c a l V a r i a b l e s
467C-----------------------------------------------
468 INTEGER I
469 my_real
470 . xd, yd, zd,
471 . nnx1(mvsiz), nnx2(mvsiz), nnx3(mvsiz), nnx4(mvsiz),
472 . nny1(mvsiz), nny2(mvsiz), nny3(mvsiz), nny4(mvsiz),
473 . nnz1(mvsiz), nnz2(mvsiz), nnz3(mvsiz), nnz4(mvsiz),
474 . x12(mvsiz),y12(mvsiz),z12(mvsiz),
475 . x43(mvsiz),y43(mvsiz),z43(mvsiz),
476 . x23(mvsiz),y23(mvsiz),z23(mvsiz),
477 . x14(mvsiz),y14(mvsiz),z14(mvsiz),
478 . nnx12(mvsiz),nny12(mvsiz),nnz12(mvsiz),
479 . nnx43(mvsiz),nny43(mvsiz),nnz43(mvsiz),
480 . nnx23(mvsiz),nny23(mvsiz),nnz23(mvsiz),
481 . nnx14(mvsiz),nny14(mvsiz),nnz14(mvsiz),
482 . xa(mvsiz),ya(mvsiz),za(mvsiz),
483 . xb(mvsiz),yb(mvsiz),zb(mvsiz),
484 . nnxa(mvsiz),nnya(mvsiz),nnza(mvsiz),
485 . nnxb(mvsiz),nnyb(mvsiz),nnzb(mvsiz),
486 . csx(mvsiz),csy(mvsiz),csz(mvsiz),
487 . crx(mvsiz),cry(mvsiz),crz(mvsiz),
488 . r(mvsiz),s(mvsiz),tt(mvsiz),ta(mvsiz),tb(mvsiz),aaa,
489 . bbb,ccc,ddd,eee,rm,rp,sm,sp,rpmm,spmm,rppm,sppm,
490 . nx,ny,nz,dist,eps,xp,yp,zp,xx,yy,zz,xxp,yyp,zzp
491C-----------------------------------------------
492C Bicubic interpolation
493c with continuity of coordinates
494c of tangent and normal derivatives on the edges
495C-----------------------------------------------
496
497 DO i=1,jlt
498 nnx1(i)=nod_normal(1,ix1(i))
499 nny1(i)=nod_normal(2,ix1(i))
500 nnz1(i)=nod_normal(3,ix1(i))
501
502 nnx2(i)=nod_normal(1,ix2(i))
503 nny2(i)=nod_normal(2,ix2(i))
504 nnz2(i)=nod_normal(3,ix2(i))
505
506 nnx3(i)=nod_normal(1,ix3(i))
507 nny3(i)=nod_normal(2,ix3(i))
508 nnz3(i)=nod_normal(3,ix3(i))
509
510 nnx4(i)=nod_normal(1,ix4(i))
511 nny4(i)=nod_normal(2,ix4(i))
512 nnz4(i)=nod_normal(3,ix4(i))
513
514C-----------------------------------------------
515c calculation of r,s from the Hi
516C-----------------------------------------------
517 rm = (h2(i)-h1(i))/max(em20,h2(i)+h1(i))
518 rp = (h3(i)-h4(i))/max(em20,h3(i)+h4(i))
519 eps = em4
520 IF(h2(i)+h1(i)<eps) rm = rp
521 IF(h3(i)+h4(i)<eps) rp = rm
522
523 sm = (h4(i)-h1(i))/max(em20,h4(i)+h1(i))
524 sp = (h3(i)-h2(i))/max(em20,h3(i)+h2(i))
525 IF(h4(i)+h1(i)<eps) sm = sp
526 IF(h3(i)+h2(i)<eps) sp = sm
527
528 rpmm = rp - rm
529 spmm = sp - sm
530 aaa = four - rpmm*spmm
531 rppm = rp + rm
532 sppm = sp + sm
533
534 r(i)= (rppm + rppm + sppm*rpmm) / aaa
535 s(i)= (sppm + sppm + rppm*spmm) / aaa
536
537 ENDDO
538cote 12
539 CALL i7cubic(jlt,r,tt,
540 . x1,nnx1,x2,nnx2,x12,nnx12,csx,
541 . y1,nny1,y2,nny2,y12,nny12,csy,
542 . z1,nnz1,z2,nnz2,z12,nnz12,csz)
543cote 34
544 CALL i7cubic(jlt,r,tt,
545 . x4,nnx4,x3,nnx3,x43,nnx43,csx,
546 . y4,nny4,y3,nny3,y43,nny43,csy,
547 . z4,nnz4,z3,nnz3,z43,nnz43,csz)
548c ligne s
549 CALL i7cubic(jlt,s,ta,
550 . x12,nnx12,x43,nnx43,xa,nnxa,csx,
551 . y12,nny12,y43,nny43,ya,nnya,csy,
552 . z12,nnz12,z43,nnz43,za,nnza,csz)
553cote 14
554 CALL i7cubic(jlt,s,tt,
555 . x1,nnx1,x4,nnx4,x14,nnx14,crx,
556 . y1,nny1,y4,nny4,y14,nny14,cry,
557 . z1,nnz1,z4,nnz4,z14,nnz14,crz)
558cote 23
559 CALL i7cubic(jlt,s,tt,
560 . x2,nnx2,x3,nnx3,x23,nnx23,crx,
561 . y2,nny2,y3,nny3,y23,nny23,cry,
562 . z2,nnz2,z3,nnz3,z23,nnz23,crz)
563c ligne r
564 CALL i7cubic(jlt,r,tb,
565 . x14,nnx14,x23,nnx23,xb,nnxb,crx,
566 . y14,nny14,y23,nny23,yb,nnyb,cry,
567 . z14,nnz14,z23,nnz23,zb,nnzb,crz)
568
569 DO i=1,jlt
570
571C-----------------------------------------------
572c normal the cubic surface in r, s
573c cross product of the direction vectors cr and cs
574C-----------------------------------------------
575 nx = cry(i)*csz(i)-crz(i)*csy(i)
576 ny = crz(i)*csx(i)-crx(i)*csz(i)
577 nz = crx(i)*csy(i)-cry(i)*csx(i)
578 aaa = one / sqrt(nx*nx+ny*ny+nz*nz)
579 nx = nx*aaa
580 ny = ny*aaa
581 nz = nz*aaa
582
583C-----------------------------------------------
584c point calculated on the cubic surface (r,s)
585C-----------------------------------------------
586 xa(i) = half*(xa(i)+xb(i))
587 ya(i) = half*(ya(i)+yb(i))
588 za(i) = half*(za(i)+zb(i))
589
590C-----------------------------------------------
591c Vector: Calculation point -> Second node.
592C-----------------------------------------------
593 xd = xi(i) - xa(i)
594 yd = yi(i) - ya(i)
595 zd = zi(i) - za(i)
596C-----------------------------------------------
597c Distance: Cubic surface -> Second node.
598C-----------------------------------------------
599 dist = nx*xd + ny*yd + nz*zd
600
601C-----------------------------------------------
602c projection of the secondary node on the cubic surface
603C-----------------------------------------------
604 xp = xi(i) - dist*nx
605 yp = yi(i) - dist*ny
606 zp = zi(i) - dist*nz
607 xx = x43(i) - x12(i)
608 yy = y43(i) - y12(i)
609 zz = z43(i) - z12(i)
610
611 xxp = xp - x12(i)
612 yyp = yp - y12(i)
613 zzp = zp - z12(i)
614
615 bbb = xxp*xx + yyp*yy + zzp*zz
616 ccc = xx*xx + yy*yy + zz*zz
617
618C IF(BBB<ZERO.OR.BBB>CCC)THEN
619C-----------------------------------------------
620c s<-1 ou s>1
621c keep the normal and the previous penetration
622C-----------------------------------------------
623C ELSE
624 xx = x23(i) - x14(i)
625 yy = y23(i) - y14(i)
626 zz = z23(i) - z14(i)
627
628 xxp = xp - x14(i)
629 yyp = yp - y14(i)
630 zzp = zp - z14(i)
631
632 ddd = xxp*xx + yyp*yy + zzp*zz
633 eee = xx*xx + yy*yy + zz*zz
634C IF(DDD<ZERO.OR.DDD>EEE)THEN
635C-----------------------------------------------
636c r<-1 ou r>1
637c keep the normal and the previous penetration
638C-----------------------------------------------
639C ELSE
640C-----------------------------------------------
641c -1<r<+1 et -1<s<+1
642C-----------------------------------------------
643c r and can very recalculate
644c R(I) = TWO*DDD/EEE - ONE
645c S(I) = TWO*BBB/CCC - ONE
646C-----------------------------------------------
647 IF(dist > zero)THEN
648 pene(i) = gapv(i) - dist
649 n1(i) = nx
650 n2(i) = ny
651 n3(i) = nz
652 ELSE
653 pene(i) = gapv(i) + dist
654 n1(i) = -nx
655 n2(i) = -ny
656 n3(i) = -nz
657 ENDIF
658C ENDIF
659C ENDIF
660
661 ENDDO
662
663 RETURN
664 END
665!||====================================================================
666!|| i7cubic ../engine/source/interfaces/int07/i7curv.F
667!||--- called by ------------------------------------------------------
668!|| i7curv ../engine/source/interfaces/int07/i7curv.F
669!||====================================================================
670 SUBROUTINE i7cubic(JLT,R,TT,
671 . X1,NNX1,X2,NNX2,X12,NNX12,CX,
672 . Y1,NNY1,Y2,NNY2,Y12,NNY12,CY,
673 . Z1,NNZ1,Z2,NNZ2,Z12,NNZ12,CZ)
674C-----------------------------------------------
675C I m p l i c i t T y p e s
676C-----------------------------------------------
677#include "implicit_f.inc"
678C-----------------------------------------------
679C G l o b a l P a r a m e t e r s
680C-----------------------------------------------
681#include "mvsiz_p.inc"
682C-----------------------------------------------
683C D u m m y A r g u m e n t s
684C-----------------------------------------------
685 INTEGER JLT
686C REAL
687 my_real
688 . x1(mvsiz), x2(mvsiz),
689 . y1(mvsiz), y2(mvsiz),
690 . z1(mvsiz), z2(mvsiz),
691 . nnx1(mvsiz), nnx2(mvsiz),
692 . nny1(mvsiz), nny2(mvsiz),
693 . nnz1(mvsiz), nnz2(mvsiz),
694 . x12(mvsiz),y12(mvsiz),z12(mvsiz),
695 . nnx12(mvsiz),nny12(mvsiz),nnz12(mvsiz),
696 . cx(mvsiz),cy(mvsiz),cz(mvsiz),
697 . r(mvsiz),tt(mvsiz)
698C-----------------------------------------------
699C L o c a l V a r i a b l e s
700C-----------------------------------------------
701 INTEGER I
702 my_real
703 . aaa,a1,rr,dtdr,dtdr1,dtdr2,nx,ny,nz,r2,
704 . upr2,umr2,a3rr2,a2rr,
705 . erx,ery,erz,esx,esy,esz,etx,ety,etz
706c---------------------------------------------------
707c Director vector in R
708c---------------------------------------------------
709c
710c ^t
711c | ___
712c N1^ | / \ ^N2
713c \|/ \ \ 2
714c 1 0 \___-0------->r
715c -1. +1.
716c 0. R2
717c
718c t = a3 rr^3 + a2 rr^2 + a1 rr + a0
719c
720c t'r = 3 a3 rr^2 + 2 a2 rr + a1
721c
722c rr=0 t = 0 t'r = dr1
723c rr=r2 t = 0 t'r = dr2
724c
725c a0 = 0
726c 0 = a3 r2^2 + a2 r2 + a1
727c a1 = dr1
728c dr2 = 3 a3 r2^2 + 2 a2 r2 + a1
729c
730c 0 = 3 a3 r2^2 + 3 a2 r2 + 3 a1
731c dr2 = 3 a3 r2^2 + 2 a2 r2 + a1
732c a2 = - (dr2 + 2 a1)/r2
733c a3 = - (a2 r2 + a1)/r2^2
734c
735c a0 = 0
736c a1 = dr1
737c a2 = - (dr2 + 2 dr1)/r2
738c a3 = + (dr2 + dr1 )/r2^2
739c
740 DO i=1,jlt
741 erx = x2(i)-x1(i)
742 ery = y2(i)-y1(i)
743 erz = z2(i)-z1(i)
744
745 r2 = sqrt(erx*erx+ery*ery+erz*erz)
746 aaa = one / r2
747 erx = erx*aaa
748 ery = ery*aaa
749 erz = erz*aaa
750
751 etx = nnx1(i)+nnx2(i)
752 ety = nny1(i)+nny2(i)
753 etz = nnz1(i)+nnz2(i)
754
755 aaa = erx*etx + ery*ety + erz*etz
756 etx = etx - erx*aaa
757 ety = ety - ery*aaa
758 etz = etz - erz*aaa
759
760 aaa = one / sqrt(etx*etx+ety*ety+etz*etz)
761 etx = etx*aaa
762 ety = ety*aaa
763 etz = etz*aaa
764
765 esx = ety*erz-etz*ery
766 esy = etz*erx-etx*erz
767 esz = etx*ery-ety*erx
768
769c---------------------------------------------------
770c of RIV Calculation of normal (not norms)
771c---------------------------------------------------
772 dtdr1 = -(nnx1(i)*erx+nny1(i)*ery+nnz1(i)*erz)
773 . / (nnx1(i)*etx+nny1(i)*ety+nnz1(i)*etz)
774
775 dtdr2 = -(nnx2(i)*erx+nny2(i)*ery+nnz2(i)*erz)
776 . / (nnx2(i)*etx+nny2(i)*ety+nnz2(i)*etz)
777c---------------------------------------------------
778c the derivative is limited to +- 1 to restrict the size
779c of the boxes in i7tri
780c---------------------------------------------------
781 dtdr1 = dtdr1/max(one,abs(dtdr1))
782 dtdr2 = dtdr2/max(one,abs(dtdr2))
783
784 upr2 = half*(one+r(i))
785 umr2 = half*(one-r(i))
786 rr = upr2*r2
787
788c a0 = 0
789c a1 = dr1
790c a2 = - (dr2 + 2 a1)/r2
791c a3 = - (a2 r2 + a1)/r2^2
792c A1 = DTDR1
793c A2 = -(DTDR2 + TWO*DTDR1)/R2
794c A3 = (DTDR2 + DTDR1)/R2/R2
795
796c TT = ((A3*RR + A2)*RR + A1)*RR
797c DTDR = (THREE*A3*RR + TWO*A2)*RR + A1
798
799 a1 = dtdr1
800 a2rr = -(dtdr2 + dtdr1 + dtdr1)*upr2
801 a3rr2 = (dtdr2 + dtdr1)*upr2*upr2
802
803 tt(i) = ((a3rr2 + a2rr) + a1) * rr
804
805c---------------------------------------------------
806c coordonn s en r
807c---------------------------------------------------
808 x12(i) = x1(i) + rr*erx + tt(i)*etx
809 y12(i) = y1(i) + rr*ery + tt(i)*ety
810 z12(i) = z1(i) + rr*erz + tt(i)*etz
811
812 dtdr = three*a3rr2 + a2rr + a2rr + a1
813
814c---------------------------------------------------
815c Director vector in R
816c---------------------------------------------------
817 cx(i) = erx + dtdr*etx
818 cy(i) = ery + dtdr*ety
819 cz(i) = erz + dtdr*etz
820
821 nx = umr2*nnx1(i) + upr2*nnx2(i)
822 ny = umr2*nny1(i) + upr2*nny2(i)
823 nz = umr2*nnz1(i) + upr2*nnz2(i)
824
825 aaa = (nx*cx(i)+ny*cy(i)+nz*cz(i))
826 . / (cx(i)*cx(i)+cy(i)*cy(i)+cz(i)*cz(i))
827
828c---------------------------------------------------
829c normale en r
830c---------------------------------------------------
831 nnx12(i) = nx - cx(i)*aaa
832 nny12(i) = ny - cy(i)*aaa
833 nnz12(i) = nz - cz(i)*aaa
834 ENDDO
835
836 RETURN
837 END
838c
839c ^t
840c | ___
841c N1^ | / \ ^N2
842c \|/ \ \
843c 0 \___-0------->r
844c
845c 1 2
846c
847c
848c t = a3 rr^3 + a2 rr^2 + a1 rr + a0
849c
850c t'r = 3 a3 rr^2 + 2 a2 rr + a1
851c
852c rr=0 t = 0 t'r = dr1
853c rr=r2 t = 0 t'r = dr2
854c
855c a0 = 0
856c a1 = dr1
857c a2 = - (dr2 + 2 dr1)/r2
858c a3 = + (dr2 + dr1 )/r2^2
859c--------------------------------------
860c calculation tmax
861c
862c t'r = 3 a3 rr^2 + 2 a2 rr + a1 = 0
863c
864c rr = (-a2 +- sqrt(a2^2-3 a1 a3))/ 3*a3
865c rr = ((dr2 + 2 dr1)/r2
866c +- sqrt((dr2 + 2 dr1)^2/r2^2 - 3 dr1(dr2 + dr1 )/r2^2))
867c / (3(dr2 + dr1 )/r2^2)
868c rr/r2 = ((dr2 + 2 dr1)
869c +- sqrt((dr2 + 2 dr1)^2 - 3 dr1(dr2 + dr1 )))
870c / (3(dr2 + dr1))
871c rr/r2 = ((dr2 + 2 dr1) +- sqrt(dr2^2 + dr1^2 + dr1 dr2))
872c / (3(dr2 + dr1))
873c
874c R = rr/r2
875c T = t/r2
876c T = A3 R^3 + A2 R^2 + A1 R
877c A1 = dr1
878c A2 = - (dr2 + 2 dr1)
879c A3 = + (dr2 + dr1 )
880c
881c R = ((dr2 + 2 dr1) +- sqrt(dr2^2 + dr1^2 + dr1 dr2))
882c / (3(dr2 + dr1))
883c R = Rn / (3(dr2 + dr1))
884c Rn = (dr2 + 2 dr1) +- sqrt(dr2^2 + dr1^2 + dr1 dr2)
885c Rn^2 = ((dr2 + 2 dr1)^2 + (dr2^2 + dr1^2 + dr1 dr2)
886c +- 2(dr2 + 2 dr1)sqrt(dr2^2 + dr1^2 + dr1 dr2))
887c Rn^2 = (2 dr2^2 + 5 dr1^2 + 5 dr1 dr2)
888c +- 2(dr2 + 2 dr1)sqrt(dr2^2 + dr1^2 + dr1 dr2))
889c Rn^3 = [(2 dr2^2 + 5 dr1^2 + 5 dr1 dr2)
890c +- 2(dr2 + 2 dr1)sqrt(dr2^2 + dr1^2 + dr1 dr2))]
891c *[(dr2 + 2 dr1) +- sqrt(dr2^2 + dr1^2 + dr1 dr2)]
892c Rn^3 =
893c + 2 dr2^3 + 5 dr1^2 dr2 + 5 dr1 dr2^2
894c + 4 dr1 dr2^2 + 10 dr1^3 + 10 dr1^2 dr2
895c + 2 dr2^3 + 2 dr1^2 dr2 + 2 dr1 dr2^2
896c + 4 dr1 dr2^2 + 4 dr1^3 + 4 dr1^2 dr2
897c +-sqrt(dr2^2 + dr1^2 + dr1 dr2)
898c [13 dr1^2 + 4 dr2^2 + 13 dr1 dr2]
899c Rn^3 =
900c + 14 dr1^3
901c + 21 dr1^2 dr2
902c + 15 dr1 dr2^2
903c + 4 dr2^3
904c +-sqrt(dr2^2 + dr1^2 + dr1 dr2)
905c [13 dr1^2 + 4 dr2^2 + 13 dr1 dr2]
906c
907c si 0 < R < 1
908c T = (dr2 + dr1) Rn^3 / (3^3(dr2 + dr1)^3)
909c - (dr2 + 2 dr1) Rn / (3^2(dr2 + dr1)^2)
910c + dr1 Rn / (3(dr2 + dr1))
911c
912c T = Rn^3 / (27(dr2 + dr1)^2)
913c - 3(dr2 + 2 dr1) Rn^2 / (27(dr2 + dr1)^2)
914c + 9 (dr2 + dr1)dr1 Rn / (27(dr2 + dr1)^2)
915c
916c T = (Rn^3
917c - 3(dr2 + 2 dr1) Rn^2
918c + 9(dr2 + dr1)dr1 Rn ) / (27(dr2 + dr1)^2)
919c
920c
921c T (27(dr2 + dr1)^2) =
922c + 14 dr1^3
923c + 21 dr1^2 dr2
924c + 15 dr1 dr2^2
925c + 4 dr2^3
926c +-sqrt(dr2^2 + dr1^2 + dr1 dr2)
927c [13 dr1^2 + 4 dr2^2 + 13 dr1 dr2]
928c - 3(dr2 + 2 dr1)[(2 dr2^2 + 5 dr1^2 + 5 dr1 dr2)
929c +- 2(dr2 + 2 dr1)sqrt(dr2^2 + dr1^2 + dr1 dr2))]
930c + 9(dr2 + dr1)dr1 [(dr2 + 2 dr1) +- sqrt(dr2^2 + dr1^2 + dr1 dr2)]
931c
932c T (27(dr2 + dr1)^2) =
933c - 16 dr1^3
934c + 3 dr1^2 dr2
935c - 3 dr1 dr2^2
936c + 16 dr2^3
937c +-sqrt(dr2^2 + dr1^2 + dr1 dr2)
938c [13 dr1^2 + 4 dr2^2 + 13 dr1 dr2
939c - 3(dr2 + 2 dr1)2(dr2 + 2 dr1))
940c + 9(dr2 + dr1)dr1 ]
941c
942c T (27(dr2 + dr1)^2) =
943c - 16 dr1^3 + 3 dr1^2 dr2 - 3 dr1 dr2^2 + 16 dr2^3
944c -+ 2 [dr1^2 + dr2^2 + dr1dr2]^(3/2)
945c
946c T (27(dr2 + dr1)^2) =
947c - 15 (dr1^3 - dr2^3)
948c - (dr1 - dr2)^3
949c -+ 2 [dr1^2 + dr2^2 + dr1dr2]^(3/2)
950c
951c dr2=-dr1
952c T (27(-dr1+dr1)^2) =
953c - 38 dr1^3
954
955
956
957
958
my_real
#define my_real
Definition cppsort.cpp:32
i7cubic
subroutine i7cubic(jlt, r, tt, x1, nnx1, x2, nnx2, x12, nnx12, cx, y1, nny1, y2, nny2, y12, nny12, cy, z1, nnz1, z2, nnz2, z12, nnz12, cz)
Definition i7curv.F:674
i7curv
subroutine i7curv(jlt, pene, n1, n2, n3, gapv, x, nod_normal, ix1, ix2, ix3, ix4, h1, h2, h3, h4, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi)
Definition i7curv.F:443
i7cmaj
subroutine i7cmaj(jlt, cmaj, irect, nod_normal, cand_e, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, nnx1, nnx2, nnx3, nnx4, nny1, nny2, nny3, nny4, nnz1, nnz2, nnz3, nnz4)
Definition i7curv.F:272
i7curvsz
subroutine i7curvsz(nrtm, irect, numnod, itag, lent, maxcc)
Definition i7curv.F:219
i7normp
subroutine i7normp(nrtm, irect, numnod, x, nod_normal, nmn, msr, lent, maxcc, isdsiz, ircsiz, iad_elem, fr_elem, itag)
Definition i7curv.F:102
i7norm
subroutine i7norm(nrtm, irect, numnod, x, nod_normal, nmn, msr)
Definition i7curv.F:30
max
#define max(a, b)
Definition macros.h:21
myqsort
subroutine myqsort(n, a, perm, error)
Definition myqsort.F:51
spmd_i7curvcom
subroutine spmd_i7curvcom(iad_elem, fr_elem, adskyt, fskyt, isdsiz, ircsiz, itag, lenr, lens)
Definition spmd_i7curvcom.F:38