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 ,J ,N1,N2,N3,N4
45 my_real
46 . surfx,surfy,surfz,x13,y13,z13,x24,y24,z24
47C-----------------------------------------------
48
49C optimisable en spmd si ajout flag pour routine de 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 tri par packet des normales
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 ! cas 1 seule 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 ,J ,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 vecteurs: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,a,b
491C-----------------------------------------------
492C Interpolation bicubique
493c avec continuit des coordonn s,
494c des d riv s tagente et normale sur les bords
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 calcul de r,s partir des 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 normale la surface cubique en r,s
573c produit vectoriel des vecteurs directeurs cr et 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 calcul sur la surface cubique (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 vecteur: point calcul -> noeud second.
592C-----------------------------------------------
593 xd = xi(i) - xa(i)
594 yd = yi(i) - ya(i)
595 zd = zi(i) - za(i)
596C-----------------------------------------------
597c distance: surface cubique -> noeud second.
598C-----------------------------------------------
599 dist = nx*xd + ny*yd + nz*zd
600
601C-----------------------------------------------
602c projection du noeud second. sur la surface cubique
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 on conserve la normale et la p n tration ancienne
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 on conserve la normale et la p n tration ancienne
638C-----------------------------------------------
639C ELSE
640C-----------------------------------------------
641c -1<r<+1 et -1<s<+1
642C-----------------------------------------------
643c r et s peuvent tres recalcul s
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,J,K
702 my_real
703 . aaa,a1,a2,a3,rr,dtdr,dtdr1,dtdr2,nr,ns,nt,nx,ny,nz,r2,
704 . upr2,umr2,a3rr2,a2rr,
705 . erx,ery,erz,esx,esy,esz,etx,ety,etz
706c---------------------------------------------------
707c vecteur directeur en 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 d riv calcul e partir des normales(non norm es)
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 la d riv est born e +- 1 pour limiter la taille
779c des boites dans 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 vecteur directeur en 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 calcul 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
i7dst3
subroutine i7dst3(ix3, ix4, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, x0, y0, z0, nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3, nx4, ny4, nz4, p1, p2, p3, p4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, last)
Definition i7dst3.F:46