OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i16crit.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!|| i16crit ../engine/source/interfaces/int16/i16crit.F
25!||--- called by ------------------------------------------------------
26!|| i16main ../engine/source/interfaces/int16/i16main.F
27!||--- calls -----------------------------------------------------
28!|| i10box ../engine/source/interfaces/int16/i16crit.F
29!|| i16box ../engine/source/interfaces/int16/i16crit.F
30!|| i20box ../engine/source/interfaces/int16/i16crit.F
31!|| i8box ../engine/source/interfaces/int16/i16crit.F
32!||--- uses -----------------------------------------------------
33!|| element_mod ../common_source/modules/elements/element_mod.F90
34!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
35!||====================================================================
36 SUBROUTINE i16crit(
37 1 X ,NSV ,NELEM ,NSN ,EMINX ,
38 2 NME ,ITASK ,XSAV ,IXS ,IXS16 ,
39 3 IXS20 ,IXS10 ,V ,A ,XMSRG,
40 4 XSLVG )
41C-----------------------------------------------
42 USE intbufdef_mod
43 use element_mod , only : nixs
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "task_c.inc"
59 COMMON /i16tmp/SIZE
60 my_real
61 . SIZE
62C-----------------------------------------------
63C D u m m y A r g u m e n t s
64C-----------------------------------------------
65 my_real, DIMENSION(7) :: xmsrg
66 my_real, DIMENSION(7) :: xslvg
67 INTEGER NSN,ITASK,NME,
68 . NSV(*),NELEM(*),IXS(NIXS,*),IXS16(8,*),IXS20(12,*),
69 . ixs10(6,*)
70 my_real
71 . x(3,*),v(3,*),a(3,*),xsav(3,*),eminx(6,*)
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER NSNF,NMEF,NSNL,NMEL,I, J,I16,I20,LFT16,LLT16,
76 . LFT20,LLT20,LFT8,LLT8,LFT10,LLT10,I8,I10,
77 . INDEX16(MVSIZ),INDEX20(MVSIZ),INDEX8(MVSIZ),INDEX10(MVSIZ)
78 my_real
79 . xmsr(6),xslv(6),size_t ,xx,yy,zz
80C-----------------------------------------------
81C S o u r c e L i n e s
82C-----------------------------------------------
83 nsnf = 1 + itask*nsn / nthread
84 nsnl = (itask+1)*nsn / nthread
85 nmef = 1 + itask*nme / nthread
86 nmel = (itask+1)*nme / nthread
87C--------------------------------------------------------------
88C 0- criterion calculation to know if we need to sort or not
89C--------------------------------------------------------------
90 xslv(1) = -ep30
91 xslv(2) = -ep30
92 xslv(3) = -ep30
93 xslv(4) = ep30
94 xslv(5) = ep30
95 xslv(6) = ep30
96 xmsr(1) = -ep30
97 xmsr(2) = -ep30
98 xmsr(3) = -ep30
99 xmsr(4) = ep30
100 xmsr(5) = ep30
101 xmsr(6) = ep30
102C
103 size_t = zero
104C
105 DO i=nsnf,nsnl
106 j=nsv(i)
107 xx=x(1,j)+dt2*(v(1,j)+dt12*a(1,j))
108 yy=x(2,j)+dt2*(v(2,j)+dt12*a(2,j))
109 zz=x(3,j)+dt2*(v(3,j)+dt12*a(3,j))
110 xslv(1)=max(xslv(1),xx-xsav(1,j))
111 xslv(2)=max(xslv(2),yy-xsav(2,j))
112 xslv(3)=max(xslv(3),zz-xsav(3,j))
113 xslv(4)=min(xslv(4),xx-xsav(1,j))
114 xslv(5)=min(xslv(5),yy-xsav(2,j))
115 xslv(6)=min(xslv(6),zz-xsav(3,j))
116 END DO
117C------------------------------------
118C calculation of element bounds
119C------------------------------------
120 DO i=nmef,nmel
121 eminx(1,i) = ep30
122 eminx(2,i) = ep30
123 eminx(3,i) = ep30
124 eminx(4,i) = -ep30
125 eminx(5,i) = -ep30
126 eminx(6,i) = -ep30
127 ENDDO
128C
129 lft16=1
130 llt16=0
131 lft20=1
132 llt20=0
133 lft8 =1
134 llt8 =0
135 lft10=1
136 llt10=0
137 DO i=nmef,nmel
138 i8 = nelem(i)
139 i10 = i8-numels8
140 i20 = i10-numels10
141 i16 = i20-numels20
142 IF(i16>=1.AND.i16<=numels16)THEN
143 llt16=llt16+1
144 index16(llt16)=i
145 IF(llt16==mvsiz-1)THEN
146 CALL i16box(
147 1 lft16,llt16 ,nelem,eminx,nmef ,nmel ,
148 2 x ,v ,a ,ixs ,ixs16,size_t,
149 3 xmsr ,index16,xsav )
150 llt16=0
151 ENDIF
152 ELSEIF(i20>=1.AND.i20<=numels20)THEN
153 llt20=llt20+1
154 index20(llt20)=i
155 IF(llt20==mvsiz-1)THEN
156 CALL i20box(
157 1 lft20,llt20 ,nelem,eminx,nmef ,nmel ,
158 2 x ,v ,a ,ixs ,ixs20,size_t,
159 3 xmsr ,index20,xsav )
160 llt20=0
161 ENDIF
162 ELSEIF(i10>=1)THEN
163 llt10=llt10+1
164 index10(llt10)=i
165 IF(llt10==mvsiz-1)THEN
166 CALL i10box(
167 1 lft10,llt10 ,nelem,eminx,nmef ,nmel ,
168 2 x ,v ,a ,ixs ,ixs10,size_t,
169 3 xmsr ,index10,xsav )
170 llt10=0
171 ENDIF
172 ELSEIF(i8>=1)THEN
173 llt8=llt8+1
174 index8(llt8)=i
175 IF(llt8==mvsiz-1)THEN
176 CALL i8box(
177 1 lft8 ,llt8 ,nelem,eminx,nmef ,nmel ,
178 2 x ,v ,a ,ixs ,size_t,
179 3 xmsr ,index8 ,xsav )
180 llt8=0
181 ENDIF
182 ENDIF
183 END DO
184 IF(llt16>0)CALL i16box(
185 1 lft16,llt16 ,nelem,eminx,nmef ,nmel ,
186 2 x ,v ,a ,ixs ,ixs16,size_t,
187 3 xmsr ,index16,xsav )
188 IF(llt20>0)CALL i20box(
189 1 lft20,llt20 ,nelem,eminx,nmef ,nmel ,
190 2 x ,v ,a ,ixs ,ixs20,size_t,
191 3 xmsr ,index20,xsav )
192 IF(llt8>0)CALL i8box(
193 1 lft8 ,llt8 ,nelem,eminx,nmef ,nmel ,
194 2 x ,v ,a ,ixs ,size_t,
195 3 xmsr ,index8 ,xsav )
196 IF(llt10>0)CALL i10box(
197 1 lft10,llt10 ,nelem,eminx,nmef ,nmel ,
198 2 x ,v ,a ,ixs ,ixs10,size_t,
199 3 xmsr ,index10,xsav )
200C
201#include "lockon.inc"
202 xslvg(1)=max(xslvg(1),xslv(1))
203 xslvg(2)=max(xslvg(2),xslv(2))
204 xslvg(3)=max(xslvg(3),xslv(3))
205 xslvg(4)=min(xslvg(4),xslv(4))
206 xslvg(5)=min(xslvg(5),xslv(5))
207 xslvg(6)=min(xslvg(6),xslv(6))
208 xmsrg(1)=max(xmsrg(1),xmsr(1))
209 xmsrg(2)=max(xmsrg(2),xmsr(2))
210 xmsrg(3)=max(xmsrg(3),xmsr(3))
211 xmsrg(4)=min(xmsrg(4),xmsr(4))
212 xmsrg(5)=min(xmsrg(5),xmsr(5))
213 xmsrg(6)=min(xmsrg(6),xmsr(6))
214 SIZE = SIZE + size_t
215#include "lockoff.inc"
216C
217 RETURN
218 END
219!||====================================================================
220!|| i16box ../engine/source/interfaces/int16/i16crit.f
221!||--- called by ------------------------------------------------------
222!|| i16crit ../engine/source/interfaces/int16/i16crit.F
223!|| i17crit ../engine/source/interfaces/int17/i17crit.F
224!||--- uses -----------------------------------------------------
225!|| element_mod ../common_source/modules/elements/element_mod.F90
226!||====================================================================
227 SUBROUTINE i16box(LFT ,LLT ,NELEM,EMINX,NMEF,NMEL,
228 2 X ,V ,A ,IXS ,IXS16,SIZE,
229 3 XMSR,INDEX,XSAV )
230 use element_mod , only : nixs
231C-----------------------------------------------
232C I m p l i c i t T y p e s
233C-----------------------------------------------
234#include "implicit_f.inc"
235C-----------------------------------------------
236C C o m m o n B l o c k s
237C-----------------------------------------------
238#include "com04_c.inc"
239#include "com08_c.inc"
240C-----------------------------------------------
241C D u m m y A r g u m e n t s
242C-----------------------------------------------
243 INTEGER LFT ,LLT,NMEF,NMEL,
244 . IXS(NIXS,*),IXS16(8,*),NELEM(*),INDEX(*)
245C REAL
246 my_real
247 . X(3,*),V(3,*),A(3,*),EMINX(6,*),SIZE,XMSR(*),XSAV(3,*)
248C-----------------------------------------------
249C L o c a l V a r i a b l e s
250C-----------------------------------------------
251 INTEGER I,J,L,NE,IFACE,IDIR,IPERM(8,2),N16
252 my_real
253 . an,ax,bn,bx,cn,cx,dx,dn,d4,d8,x1,x2,x3,x4,
254 . x9,x10,x11,x12,xc
255 DATA iperm / 2, 3, 4, 5, 1, 2, 3, 4,
256 2 6, 7, 8, 9, 5, 6, 7, 8/
257C
258C-----------------------------------------------
259C/*
260C
261C ( 7)8==============(14)6=============( 6)7
262C //| //|
263C // | //||
264C // | // ||
265C // | // ||
266C // | // ||
267C (15)7 | (13)5 ||
268C // | // ||
269C // ( 3)4--------------(10)2------//-----( 2)3
270C // / // //
271C // / // //
272C // / // //
273C ( 8)9==============(16)8=============( 5)6 //
274C || / || //
275C || (11)3 (C) || ( 9)1
276C || / || //
277C || / || //
278C || / || //
279C || / ||//
280C ||/ ||/
281C ( 4)5==============(12)4=============( 1)2
282C
283C*/
284C---------------------------------------------------------
285C MONODIM
286C---------------------------------------------------------
287C
288C x
289C / \
290C / (2)
291C (3)
292C /
293C /
294C (1)
295C
296C N1 = -0.5 (1-r)r dN1/dr = r - 0.5
297C N2 = 0.5 (1+r)r dN2/dr = r + 0.5
298C N3 = (1-r^2) dN3/dr = - 2 r
299C
300C x = N1 x1 + N2 x2 + N3 x3
301C x = -0.5 (1-r)r x1 + 0.5 (1+r)r x2 + (1-r^2) x3
302C 2 x = (x1 + x2 - 2 x3) r^2 + (x2-x1) r + 2 x3
303C
304C 0) Point XMAX
305C
306C dx/dr = (r - 0.5) x1 + (r + 0.5) x2 - 2 r x3 = 0
307C r = 0.5 (x1 - x2) / (x1 + x2 - 2 x3)
308C
309C 2 x (x1 + x2 - 2 x3) = (x1 + x2 - 2 x3)^2 r^2
310C + (x2-x1)(x1 + x2 - 2 x3) r
311C + 2 (x1 + x2 - 2 x3)x3
312C
313C 2 x (x1 + x2 - 2 x3) = - 0.25 (x1 - x2)^2
314C + 2 (x1 + x2 - 2 x3)x3
315C
316C x = x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
317C
318C------------------------------------------------------------
319C solution 0 => x < x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
320C si x3 -> (x1 + x2)/ 2 x -> infini
321C------------------------------------------------------------
322C
323C 1) search for the xmax point between 1 and 2
324C
325C si r > 1
326C => x = x2
327C
328C si r < -1
329C => x = x1
330C
331C si -1 < r < 1
332C => -1/4 < 0.125 (x1 - x2) / (x1 + x2 - 2 x3) < 1/4
333C
334C (x - x3)/(x1 - x2) = - 0.125 (x1 - x2) / (x1 + x2 - 2 x3)
335C -1/4 < (x - x3)/ < 1/4
336C
337C si x2 > x1 x3 -1/4 (x2 - x1) < x < x3 + 1/4 (x2 - x1)
338C si x2 < x1 x3 -1/4 (x2 - x1) > x > x3 + 1/4 (x2 - x1)
339C
340C => x3 - 1/4 |x2 - x1| < x < x3 + 1/4 |x2 - x1|
341C
342C------------------------------------------------------------
343C solution 1 => x < max (x1 , x2 , x3 + 1/4 |x2 - x1|)
344C------------------------------------------------------------
345C
346C 2) search for the most unfavorable position of x3
347C
348C x = x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
349C dx/dx3 = 1 - 0.25 (x1 - x2)^2 / (x1 + x2 - 2 x3)^2 = 0
350C (2 x3 - x1 - x2 )^2 = 0.25 (x2 - x1)^2
351C if x2 > x1 and x3 > (x1 + x2)/2 or
352C si x2 < x1 et x3 < (x1 + x2)/2
353C (2 x3 - x1 - x2 ) = 0.5 (x2 - x1)
354C x3 = (x1 + x2)/2 + (x2 - x1)/4
355C if x2 > x1 and x3 < (x1 + x2)/2 or
356C si x2 < x1 et x3 > (x1 + x2)/2
357C no solution
358C
359C si x3 < (x1 + x2)/2 + (x2 - x1)/4
360C => x < max (x1 , x2 ) et
361C => x> min (x1, x2) => Verified by solution 1 and 2
362C
363C si x3 = (x1 + x2)/2 + (x2 - x1)/4
364C x = x3 + 0.125 (x1 - x2)^2 / (2 x3 - x1 - x2 )
365C x = (x1 + x2)/2 + (x2 - x1)/4
366C + 0.125 (x1 - x2)^2 / (2 ((x1 + x2)/2 + (x2 - x1)/4) - x1 - x2)
367C x = (x1 + x2)/2 + (x2 - x1)/4
368C + 0.25 (x2 - x1)^2 / (x2 - x1)
369C si x2 > x1 => x = x2
370C si x2 < x1 => x = x1
371C
372C et x < x3 + 1/4 |x2 - x1|
373C
374C si x3 = x2
375C x = x3 + 0.125 (x1 - x2)^2 / (2 x3 - x1 - x2 )
376C x = x3 + 0.125 (x1 - x2)^2 / (x2 - x1)
377C x = x3 +- 1/8 |x1 - x2| = x2 +- 1/8 |x1 - x2|
378C
379C si (x1 + x2)/2 + (x2 - x1)/4 < x3 < max(x1,x2)
380C x < max(x1,x2) + 1/8 |x1 - x2|
381C et x < x3 + 1/4 |x1 - x2|
382C
383C si x3 > max(x1,x2)
384C x < x3 + 1/8 |x1 - x2|
385C------------------------------------------------------------
386C solution 2 => x < max (x1 , x2 , x3) + 1/8 |x2 - x1|
387C------------------------------------------------------------
388C solution 1 : x < max (x1 , x2 , x3 + 1/4 |x2 - x1|)
389C------------------------------------------------------------
390C => x <min (solution 1, solution 2)
391C------------------------------------------------------------
392C
393C 3) Exact solution (0 terminal solution by solution 1)
394C
395C solution 0 :
396C x = x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
397C solution 1 :
398C x < max (x1 , x2 , x3 + 1/4 |x2 - x1|)
399C------------------------------------------------------------
400C Solution 3 => x <min (solution 1, solution 0)
401C------------------------------------------------------------
402C s = (x1+x2)/2 d = |x2-x1|
403C-------------------------------------------------------------------------
404C x3 x min x max
405C-------------------------------------------------------------------------
406C -inf < x3 < s - d/4 x3 + d^2 / 16(x3-s) max(x1,x2)
407C s - d/4 < x3 < s + d/4 min(x1,x2) max(x1,x2)
408C s + d/4 < x3 < +inf min(x1,x2) x3 + d^2 / 16(x3-s)
409C-------------------------------------------------------------------------
410C---------------------------------------------------------
411C BIDIM ( pas resolu )
412C---------------------------------------------------------
413C
414C-----------------------------------------------
415C i ri si ti Ni
416C--------------------------------------------------------------------
417C 1 -1 -1 -1 1/4(1-r)(1-t)(-r-t-1)
418C 2 -1 -1 +1 1/4(1-r)(1+t)(-r+t-1)
419C 3 +1 -1 +1 1/4(1+r)(1+t)(+r+t-1)
420C 4 +1 -1 -1 1/4(1+r)(1-t)(+r-t-1)
421C 9 -1 -1 0 1/2(1-t^2)(1-r)
422C 10 0 -1 +1 1/2(1-r^2) (1+t)
423C 11 +1 -1 0 1/2(1-t^2)(1+r)
424C 12 0 -1 -1 1/2(1-r^2) (1-t)
425C
426C x = N1 x1 + N2 x2 + N3 x3 + N4 x4
427C + N9 x9 + N10 x10 + N11 x11 + N12 x12
428C
429C 0) Point XMAX
430C
431C dx/dr = -1/4(1-t)(-2r-t) x1
432C -1/4(1+t)(-2r+t) x2
433C +1/4(1+t)(+2r+t) x3
434C +1/4(1-t)(+2r-t) x4
435C -1/2(1-t^2) x9
436C -r(1+t) x10
437C +1/2(1-t^2) x11
438C -r(1-t) x12 = 0
439C dx/dt = -1/4(1-r)(-2t-r) x1
440C +1/4(1-r)(+2t-r) x2
441C +1/4(1+r)(+2t+r) x3
442C +1/4(1+r)(-2t+r) x4
443C -t(1-r) x9
444C +1/2(1-r^2) x10
445C -t(1+r) x11
446C -1/2(1-r^2) x12 = 0
447C------------------------------------
448C calculation of element bounds
449C------------------------------------
450 DO iface=1,2
451C-----------------------------------------------------------------------
452C Face 1 2 3 4 or 5 6 7 8
453C-----------------------------------------------------------------------
454 DO idir=1,3
455C-----------------------------------------------------------------------
456C x y or z
457C-----------------------------------------------------------------------
458 DO l=lft,llt
459 i = index(l)
460 ne = nelem(i)
461 n16= ne - numels8 - numels10 - numels20
462C
463 j = ixs(iperm(1,iface),ne)
464 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
465 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
466 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
467 j = ixs(iperm(2,iface),ne)
468 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
469 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
470 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
471 j = ixs(iperm(3,iface),ne)
472 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
473 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
474 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
475 j = ixs(iperm(4,iface),ne)
476 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
477 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
478 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
479 j = ixs16(iperm(5,iface),n16)
480 x9 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
481 xmsr(idir) =max(xmsr(idir) ,x9-xsav(idir,j))
482 xmsr(idir+3)=min(xmsr(idir+3),x9-xsav(idir,j))
483 j = ixs16(iperm(6,iface),n16)
484 x10= x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
485 xmsr(idir) =max(xmsr(idir) ,x10-xsav(idir,j))
486 xmsr(idir+3)=min(xmsr(idir+3),x10-xsav(idir,j))
487 j = ixs16(iperm(7,iface),n16)
488 x11= x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
489 xmsr(idir) =max(xmsr(idir) ,x11-xsav(idir,j))
490 xmsr(idir+3)=min(xmsr(idir+3),x11-xsav(idir,j))
491 j = ixs16(iperm(8,iface),n16)
492 x12= x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
493 xmsr(idir) =max(xmsr(idir) ,x12-xsav(idir,j))
494 xmsr(idir+3)=min(xmsr(idir+3),x12-xsav(idir,j))
495C
496 xc = half*(x9+x10+x11+x12) - fourth*(x1+x2+x3+x4)
497C
498 d4 = fourth * abs(x1-x2)
499 an = min( x1 , x2 , x9-d4 )
500 ax = max( x1 , x2 , x9+d4 )
501C
502 d4 = fourth * abs(x3-x4)
503 bn = min( x3 , x4 , x11-d4 )
504 bx = max( x3 , x4 , x11+d4 )
505C
506 d4 = fourth * abs(x12-x10)
507 cn = min( x12 , x10 , xc-d4 )
508 cx = max( x12 , x10 , xc+d4 )
509C
510 d8 = one_over_8 * max( ax-bn , bx-an )
511 d4 = d8 + d8
512 dn = max(min( an , bn , cn-d4 ),min(an , bn , cn) - d8 )
513 dx = min(max( ax , bx , cx+d4 ),max( ax , bx , cx) + d8 )
514C
515 eminx(idir,i) = min( eminx(idir,i) , dn )
516 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
517C
518 SIZE = SIZE + dx - dn
519C
520 ENDDO
521 ENDDO
522 ENDDO
523C--------------------------------------------------------------
524C
525 RETURN
526 END
527!||====================================================================
528!|| i20box ../engine/source/interfaces/int16/i16crit.F
529!||--- called by ------------------------------------------------------
530!|| i16crit ../engine/source/interfaces/int16/i16crit.F
531!|| i17crit ../engine/source/interfaces/int17/i17crit.F
532!||--- uses -----------------------------------------------------
533!|| element_mod ../common_source/modules/elements/element_mod.F90
534!||====================================================================
535 SUBROUTINE i20box(LFT ,LLT ,NELEM,EMINX,NMEF,NMEL,
536 2 X ,V ,A ,IXS ,IXS20,SIZE,
537 3 XMSR,INDEX,XSAV )
538 use element_mod , only : nixs
539C-----------------------------------------------
540C I m p l i c i t T y p e s
541C-----------------------------------------------
542#include "implicit_f.inc"
543C-----------------------------------------------
544C C o m m o n B l o c k s
545C-----------------------------------------------
546#include "com04_c.inc"
547#include "com08_c.inc"
548C-----------------------------------------------
549C D u m m y A r g u m e n t s
550C-----------------------------------------------
551 INTEGER LFT ,LLT,NMEF,NMEL,
552 . IXS(NIXS,*),IXS20(12,*),NELEM(*),INDEX(*)
553C REAL
554 my_real
555 . X(3,*),V(3,*),A(3,*),EMINX(6,*),SIZE,XMSR(*),XSAV(3,*)
556C-----------------------------------------------
557C L o c a l V a r i a b l e s
558C-----------------------------------------------
559 INTEGER I,J,L,NE,IDIR,N20
560 my_real
561 . AN12,AX12,AN34,AX34,AN56,AX56,AN78,AX78,CN,CX,DX,DN,D4,D8,
562 . x1,x2,x3,x4,x5,x6,x7,x8,
563 . x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,xc
564C------------------------------------
565C calculation of element bounds
566C------------------------------------
567C-----------------------------------------------------------------------
568C Face 1 2 3 4 or 5 6 7 8
569C-----------------------------------------------------------------------
570 DO idir=1,3
571C-----------------------------------------------------------------------
572C x y or z
573C-----------------------------------------------------------------------
574 DO l=lft,llt
575 i = index(l)
576 ne = nelem(i)
577 n20= ne - numels8 - numels10
578C-----------------------------------------------------------------------
579 j = ixs(2,ne)
580 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
581 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
582 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
583 j = ixs(3,ne)
584 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
585 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
586 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
587 j = ixs(4,ne)
588 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
589 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
590 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
591 j = ixs(5,ne)
592 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
593 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
594 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
595 j = ixs(6,ne)
596 x5 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
597 xmsr(idir) =max(xmsr(idir) ,x5-xsav(idir,j))
598 xmsr(idir+3)=min(xmsr(idir+3),x5-xsav(idir,j))
599 j = ixs(7,ne)
600 x6 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
601 xmsr(idir) =max(xmsr(idir) ,x6-xsav(idir,j))
602 xmsr(idir+3)=min(xmsr(idir+3),x6-xsav(idir,j))
603 j = ixs(8,ne)
604 x7 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
605 xmsr(idir) =max(xmsr(idir) ,x7-xsav(idir,j))
606 xmsr(idir+3)=min(xmsr(idir+3),x7-xsav(idir,j))
607 j = ixs(9,ne)
608 x8 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
609 xmsr(idir) =max(xmsr(idir) ,x8-xsav(idir,j))
610 xmsr(idir+3)=min(xmsr(idir+3),x8-xsav(idir,j))
611C
612 j = ixs20(1,n20)
613 IF(j/=0)THEN
614 x9 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
615 xmsr(idir) =max(xmsr(idir) ,x9-xsav(idir,j))
616 xmsr(idir+3)=min(xmsr(idir+3),x9-xsav(idir,j))
617 ELSE
618 x9=0.5*(x(idir,ixs(2,ne))+x(idir,ixs(3,ne)))
619 ENDIF
620 j = ixs20(2,n20)
621 IF(j/=0)THEN
622 x10 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
623 xmsr(idir) =max(xmsr(idir) ,x10-xsav(idir,j))
624 xmsr(idir+3)=min(xmsr(idir+3),x10-xsav(idir,j))
625 ELSE
626 x10=0.5*(x(idir,ixs(3,ne))+x(idir,ixs(4,ne)))
627 ENDIF
628 j = ixs20(3,n20)
629 IF(j/=0)THEN
630 x11 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
631 xmsr(idir) =max(xmsr(idir) ,x11-xsav(idir,j))
632 xmsr(idir+3)=min(xmsr(idir+3),x11-xsav(idir,j))
633 ELSE
634 x11=0.5*(x(idir,ixs(4,ne))+x(idir,ixs(5,ne)))
635 ENDIF
636 j = ixs20(4,n20)
637 IF(j/=0)THEN
638 x12 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
639 xmsr(idir) =max(xmsr(idir) ,x12-xsav(idir,j))
640 xmsr(idir+3)=min(xmsr(idir+3),x12-xsav(idir,j))
641 ELSE
642 x12=0.5*(x(idir,ixs(5,ne))+x(idir,ixs(2,ne)))
643 ENDIF
644 j = ixs20(5,n20)
645 IF(j/=0)THEN
646 x13 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
647 xmsr(idir) =max(xmsr(idir) ,x13-xsav(idir,j))
648 xmsr(idir+3)=min(xmsr(idir+3),x13-xsav(idir,j))
649 ELSE
650 x13=0.5*(x(idir,ixs(2,ne))+x(idir,ixs(6,ne)))
651 ENDIF
652 j = ixs20(6,n20)
653 IF(j/=0)THEN
654 x14 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
655 xmsr(idir) =max(xmsr(idir) ,x14-xsav(idir,j))
656 xmsr(idir+3)=min(xmsr(idir+3),x14-xsav(idir,j))
657 ELSE
658 x14=0.5*(x(idir,ixs(3,ne))+x(idir,ixs(6,ne)))
659 ENDIF
660 j = ixs20(7,n20)
661 IF(j/=0)THEN
662 x15 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
663 xmsr(idir) =max(xmsr(idir) ,x15-xsav(idir,j))
664 xmsr(idir+3)=min(xmsr(idir+3),x15-xsav(idir,j))
665 ELSE
666 x15=0.5*(x(idir,ixs(4,ne))+x(idir,ixs(8,ne)))
667 ENDIF
668 j = ixs20(8,n20)
669 IF(j/=0)THEN
670 x16 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
671 xmsr(idir) =max(xmsr(idir) ,x16-xsav(idir,j))
672 xmsr(idir+3)=min(xmsr(idir+3),x16-xsav(idir,j))
673 ELSE
674 x16=0.5*(x(idir,ixs(5,ne))+x(idir,ixs(9,ne)))
675 ENDIF
676 j = ixs20(9,n20)
677 IF(j/=0)THEN
678 x17 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
679 xmsr(idir) =max(xmsr(idir) ,x17-xsav(idir,j))
680 xmsr(idir+3)=min(xmsr(idir+3),x17-xsav(idir,j))
681 ELSE
682 x17=0.5*(x(idir,ixs(6,ne))+x(idir,ixs(7,ne)))
683 ENDIF
684 j = ixs20(10,n20)
685 IF(j/=0)THEN
686 x18 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
687 xmsr(idir) =max(xmsr(idir) ,x18-xsav(idir,j))
688 xmsr(idir+3)=min(xmsr(idir+3),x18-xsav(idir,j))
689 ELSE
690 x18=0.5*(x(idir,ixs(7,ne))+x(idir,ixs(8,ne)))
691 ENDIF
692 j = ixs20(11,n20)
693 IF(j/=0)THEN
694 x19 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
695 xmsr(idir) =max(xmsr(idir) ,x19-xsav(idir,j))
696 xmsr(idir+3)=min(xmsr(idir+3),x19-xsav(idir,j))
697 ELSE
698 x19=0.5*(x(idir,ixs(8,ne))+x(idir,ixs(9,ne)))
699 ENDIF
700 j = ixs20(12,n20)
701 IF(j/=0)THEN
702 x20 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
703 xmsr(idir) =max(xmsr(idir) ,x20-xsav(idir,j))
704 xmsr(idir+3)=min(xmsr(idir+3),x20-xsav(idir,j))
705 ELSE
706 x20=0.5*(x(idir,ixs(6,ne))+x(idir,ixs(9,ne)))
707 ENDIF
708C
709C-----------------------------------------------------------------------
710C Face 1 2 3 4
711C-----------------------------------------------------------------------
712 xc = half*(x9+x10+x11+x12) - fourth*(x1+x2+x3+x4)
713C
714 d4 = fourth * abs(x1-x2)
715 an12 = min( x1 , x2 , x9-d4 )
716 ax12 = max( x1 , x2 , x9+d4 )
717C
718 d4 = fourth * abs(x3-x4)
719 an34 = min( x3 , x4 , x11-d4 )
720 ax34 = max( x3 , x4 , x11+d4 )
721C
722 d4 = fourth * abs(x12-x10)
723 cn = min( x12 , x10 , xc-d4 )
724 cx = max( x12 , x10 , xc+d4 )
725C
726 d8 = one_over_8 * max( ax12-an34 , ax34-an12 )
727 d4 = d8 + d8
728 dn = max(min(an12 , an34 , cn-d4 ),
729 . min(an12 , an34 , cn) - d8 )
730 dx = min(max(ax12 , ax34 , cx+d4 ),
731 . max(ax12 , ax34 , cx) + d8 )
732C
733 eminx(idir,i) = min( eminx(idir,i) , dn )
734 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
735C-----------------------------------------------------------------------
736C Face 5 6 7 8
737C-----------------------------------------------------------------------
738 xc = half*(x17+x18+x19+x20) - fourth*(x5+x6+x7+x8)
739C
740 d4 = fourth * abs(x5-x6)
741 an56 = min( x5 , x6 , x17-d4 )
742 ax56 = max( x5 , x6 , x17+d4 )
743C
744 d4 = fourth * abs(x7-x8)
745 an78 = min( x7 , x8 , x19-d4 )
746 ax78 = max( x7 , x8 , x19+d4 )
747C
748 d4 = fourth * abs(x20-x18)
749 cn = min( x20 , x18 , xc-d4 )
750 cx = max( x20 , x18 , xc+d4 )
751C
752 d8 = one_over_8 * max( ax56-an78 , ax78-an56 )
753 d4 = d8 + d8
754 dn = max(min(an56 , an78 , cn-d4 ),
755 . min(an56 , an78 , cn) - d8 )
756 dx = min(max(ax56 , ax78 , cx+d4 ),
757 . max(ax56 , ax78 , cx) + d8 )
758C
759 eminx(idir,i) = min( eminx(idir,i) , dn )
760 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
761C-----------------------------------------------------------------------
762C Face 1 2 6 5
763C-----------------------------------------------------------------------
764 xc = half*(x9+x14+x17+x13) - fourth*(x1+x2+x6+x5)
765C
766 d4 = fourth * abs(x13-x14)
767 cn = min( x13 , x14 , xc-d4 )
768 cx = max( x13 , x14 , xc+d4 )
769C
770 d8 = one_over_8 * max( ax12-an56 , ax56-an12 )
771 d4 = d8 + d8
772 dn = max(min(an12 , an56 , cn-d4 ),
773 . min(an12 , an56 , cn) - d8 )
774 dx = min(max(ax12 , ax56 , cx+d4 ),
775 . max(ax12 , ax56 , cx) + d8 )
776C
777 eminx(idir,i) = min( eminx(idir,i) , dn )
778 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
779C-----------------------------------------------------------------------
780C Face 3 4 8 7
781C-----------------------------------------------------------------------
782 xc = half*(x11+x15+x19+x16) - fourth*(x3+x4+x8+x7)
783C
784 d4 = fourth * abs(x16-x15)
785 cn = min( x15 , x16 , xc-d4 )
786 cx = max( x15 , x16 , xc+d4 )
787C
788 d8 = one_over_8 * max( ax34-an78 , ax78-an34 )
789 d4 = d8 + d8
790 dn = max(min(an34 , an78 , cn-d4 ),
791 . min(an34 , an78 , cn) - d8 )
792 dx = min(max(ax34 , ax78 , cx+d4 ),
793 . max(ax34 , ax78 , cx) + d8 )
794C
795 eminx(idir,i) = min( eminx(idir,i) , dn )
796 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
797C-----------------------------------------------------------------------
798C Face 4 1 5 8
799C-----------------------------------------------------------------------
800 xc = half*(x12+x13+x20+x16) - fourth*(x4+x1+x5+x8)
801C
802 d4 = fourth * abs(x4-x1)
803 an12 = min( x4 , x1 , x12-d4 )
804 ax12 = max( x4 , x1 , x12+d4 )
805C
806 d4 = fourth * abs(x8-x5)
807 an34 = min( x8 , x5 , x20-d4 )
808 ax34 = max( x8 , x5 , x20+d4 )
809C
810 d4 = fourth * abs(x16-x13)
811 cn = min( x16 , x13 , xc-d4 )
812 cx = max( x16 , x13 , xc+d4 )
813C
814 d8 = one_over_8 * max( ax12-an34 , ax34-an12 )
815 d4 = d8 + d8
816 dn = max(min(an12 , an34 , cn-d4 ),
817 . min(an12 , an34 , cn) - d8 )
818 dx = min(max(ax12 , ax34 , cx+d4 ),
819 . max(ax12 , ax34 , cx) + d8 )
820C
821 eminx(idir,i) = min( eminx(idir,i) , dn )
822 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
823C-----------------------------------------------------------------------
824C Face 3 2 6 7
825C-----------------------------------------------------------------------
826 xc = half*(x10+x14+x18+x15) - fourth*(x3+x2+x6+x7)
827C
828 d4 = fourth * abs(x3-x2)
829 an12 = min( x3 , x2 , x10-d4 )
830 ax12 = max( x3 , x2 , x10+d4 )
831C
832 d4 = fourth * abs(x7-x6)
833 an34 = min( x7 , x6 , x18-d4 )
834 ax34 = max( x7 , x6 , x18+d4 )
835C
836 d4 = fourth * abs(x15-x14)
837 cn = min( x15 , x14 , xc-d4 )
838 cx = max( x15 , x14 , xc+d4 )
839C
840 d8 = one_over_8* max( ax12-an34 , ax34-an12 )
841 d4 = d8 + d8
842 dn = max(min(an12 , an34 , cn-d4 ),
843 . min(an12 , an34 , cn) - d8 )
844 dx = min(max(ax12 , ax34 , cx+d4 ),
845 . max(ax12 , ax34 , cx) + d8 )
846C
847 eminx(idir,i) = min( eminx(idir,i) , dn )
848 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
849C-----------------------------------------------------------------------
850 SIZE = SIZE + dx - dn
851C
852 ENDDO
853 ENDDO
854C--------------------------------------------------------------
855C
856 RETURN
857 END
858!||====================================================================
859!|| i10box ../engine/source/interfaces/int16/i16crit.F
860!||--- called by ------------------------------------------------------
861!|| i16crit ../engine/source/interfaces/int16/i16crit.F
862!||--- uses -----------------------------------------------------
863!|| element_mod ../common_source/modules/elements/element_mod.F90
864!||====================================================================
865 SUBROUTINE i10box(LFT ,LLT ,NELEM,EMINX,NMEF,NMEL,
866 2 X ,V ,A ,IXS ,IXS10,SIZE,
867 3 XMSR,INDEX,XSAV )
868 use element_mod , only : nixs
869C-----------------------------------------------
870C I m p l i c i t T y p e s
871C-----------------------------------------------
872#include "implicit_f.inc"
873C-----------------------------------------------
874C C o m m o n B l o c k s
875C-----------------------------------------------
876#include "com04_c.inc"
877#include "com08_c.inc"
878C-----------------------------------------------
879C D u m m y A r g u m e n t s
880C-----------------------------------------------
881 INTEGER LFT ,LLT,NMEF,NMEL,
882 . IXS(NIXS,*),IXS10(6,*),NELEM(*),INDEX(*)
883C REAL
884 my_real
885 . X(3,*),V(3,*),A(3,*),EMINX(6,*),SIZE,XMSR(*),XSAV(3,*)
886C-----------------------------------------------
887C L o c a l V a r i a b l e s
888C-----------------------------------------------
889 INTEGER I,J,L,NE,IDIR,N10
890 my_real
891 .
892 . x1,x2,x3,x4,x5,x6,x7,x8,
893 . x9,x10,xc,xx,xn
894C------------------------------------
895C calculation of element bounds
896C------------------------------------
897C-----------------------------------------------------------------------
898C Face 1 2 3 4 or 5 6 7 8
899C-----------------------------------------------------------------------
900 DO idir=1,3
901C-----------------------------------------------------------------------
902C x y or z
903C-----------------------------------------------------------------------
904 DO l=lft,llt
905 i = index(l)
906 ne = nelem(i)
907 n10= ne - numels8
908C-----------------------------------------------------------------------
909 j = ixs(2,ne)
910 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
911 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
912 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
913 j = ixs(4,ne)
914 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
915 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
916 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
917 j = ixs(6,ne)
918 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
919 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
920 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
921 j = ixs(7,ne)
922 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
923 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
924 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
925C
926 j = ixs10(1,n10)
927 x5 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
928 xmsr(idir) =max(xmsr(idir) ,x5-xsav(idir,j))
929 xmsr(idir+3)=min(xmsr(idir+3),x5-xsav(idir,j))
930 j = ixs10(2,n10)
931 x6 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
932 xmsr(idir) =max(xmsr(idir) ,x6-xsav(idir,j))
933 xmsr(idir+3)=min(xmsr(idir+3),x6-xsav(idir,j))
934 j = ixs10(3,n10)
935 x7 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
936 xmsr(idir) =max(xmsr(idir) ,x7-xsav(idir,j))
937 xmsr(idir+3)=min(xmsr(idir+3),x7-xsav(idir,j))
938 j = ixs10(4,n10)
939 x8 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
940 xmsr(idir) =max(xmsr(idir) ,x8-xsav(idir,j))
941 xmsr(idir+3)=min(xmsr(idir+3),x8-xsav(idir,j))
942 j = ixs10(5,n10)
943 x9 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
944 xmsr(idir) =max(xmsr(idir) ,x9-xsav(idir,j))
945 xmsr(idir+3)=min(xmsr(idir+3),x9-xsav(idir,j))
946 j = ixs10(6,n10)
947 x10 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
948 xmsr(idir) =max(xmsr(idir) ,x10-xsav(idir,j))
949 xmsr(idir+3)=min(xmsr(idir+3),x10-xsav(idir,j))
950C-----------------------------------------------------------------------
951 xx=max(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,x9,x10)
952 xn=min(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,x9,x10)
953 eminx(idir,i) = min( eminx(idir,i) , xn )
954 eminx(idir+3,i) = max( eminx(idir+3,i), xx )
955C-----------------------------------------------------------------------
956C Face 1 2 3 4
957C-----------------------------------------------------------------------
958 xc = (two*(x5+x6+x7) - (x1+x2+x3))* third
959 eminx(idir,i) = min( eminx(idir,i) , xc )
960 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
961C
962 xc = (two*(x5+x8+x9) - (x1+x2+x4))*third
963 eminx(idir,i) = min( eminx(idir,i) , xc )
964 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
965C
966 xc = (two*(x6+x9+x10) - (x2+x3+x4)) * third
967 eminx(idir,i) = min( eminx(idir,i) , xc )
968 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
969C
970 xc = (two*(x7+x8+x10) - (x3+x1+x4)) * third
971 eminx(idir,i) = min( eminx(idir,i) , xc )
972 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
973C
974C-----------------------------------------------------------------------
975C
976 ENDDO
977 ENDDO
978C--------------------------------------------------------------
979C
980 RETURN
981 END
982!||====================================================================
983!|| i8box ../engine/source/interfaces/int16/i16crit.F
984!||--- called by ------------------------------------------------------
985!|| i16crit ../engine/source/interfaces/int16/i16crit.F
986!||--- uses -----------------------------------------------------
987!|| element_mod ../common_source/modules/elements/element_mod.F90
988!||====================================================================
989 SUBROUTINE i8box(LFT ,LLT ,NELEM,EMINX,NMEF,NMEL,
990 2 X ,V ,A ,IXS ,SIZE,
991 3 XMSR,INDEX,XSAV )
992 use element_mod , only : nixs
993C-----------------------------------------------
994C I m p l i c i t T y p e s
995C-----------------------------------------------
996#include "implicit_f.inc"
997C-----------------------------------------------
998C C o m m o n B l o c k s
999C-----------------------------------------------
1000#include "com08_c.inc"
1001C-----------------------------------------------
1002C D u m m y A r g u m e n t s
1003C-----------------------------------------------
1004 INTEGER LFT ,LLT,NMEF,NMEL,
1005 . IXS(NIXS,*),NELEM(*),INDEX(*)
1006C REAL
1007 my_real
1008 . X(3,*),V(3,*),A(3,*),EMINX(6,*),SIZE,XMSR(*),XSAV(3,*)
1009C-----------------------------------------------
1010C L o c a l V a r i a b l e s
1011C-----------------------------------------------
1012 INTEGER I,J,L,NE,IDIR
1013 my_real
1014 . DX,DN,
1015 . X1,X2,X3,X4,X5,X6,X7,X8
1016C------------------------------------
1017C calculation of element bounds
1018C------------------------------------
1019C-----------------------------------------------------------------------
1020C Face 1 2 3 4 or 5 6 7 8
1021C-----------------------------------------------------------------------
1022 DO idir=1,3
1023C-----------------------------------------------------------------------
1024C x y or z
1025C-----------------------------------------------------------------------
1026 DO l=lft,llt
1027 i = index(l)
1028 ne = nelem(i)
1029C-----------------------------------------------------------------------
1030 j = ixs(2,ne)
1031 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1032 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
1033 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
1034 j = ixs(3,ne)
1035 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1036 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
1037 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
1038 j = ixs(4,ne)
1039 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1040 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
1041 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
1042 j = ixs(5,ne)
1043 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1044 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
1045 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
1046 j = ixs(6,ne)
1047 x5 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1048 xmsr(idir) =max(xmsr(idir) ,x5-xsav(idir,j))
1049 xmsr(idir+3)=min(xmsr(idir+3),x5-xsav(idir,j))
1050 j = ixs(7,ne)
1051 x6 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1052 xmsr(idir) =max(xmsr(idir) ,x6-xsav(idir,j))
1053 xmsr(idir+3)=min(xmsr(idir+3),x6-xsav(idir,j))
1054 j = ixs(8,ne)
1055 x7 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1056 xmsr(idir) =max(xmsr(idir) ,x7-xsav(idir,j))
1057 xmsr(idir+3)=min(xmsr(idir+3),x7-xsav(idir,j))
1058 j = ixs(9,ne)
1059 x8 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1060 xmsr(idir) =max(xmsr(idir) ,x8-xsav(idir,j))
1061 xmsr(idir+3)=min(xmsr(idir+3),x8-xsav(idir,j))
1062C
1063
1064 dx=max(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 )
1065 dn=min(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 )
1066C
1067 eminx(idir,i) = min( eminx(idir,i) , dn )
1068 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
1069C
1070 SIZE = SIZE + dx - dn
1071C
1072 ENDDO
1073 ENDDO
1074C--------------------------------------------------------------
1075C
1076 RETURN
1077 END
my_real
#define my_real
Definition cppsort.cpp:32
i16crit
subroutine i16crit(x, nsv, nelem, nsn, eminx, nme, itask, xsav, ixs, ixs16, ixs20, ixs10, v, a, xmsrg, xslvg)
Definition i16crit.F:41
i16box
subroutine i16box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs16, size, xmsr, index, xsav)
Definition i16crit.F:230
i20box
subroutine i20box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs20, size, xmsr, index, xsav)
Definition i16crit.F:538
i10box
subroutine i10box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs10, size, xmsr, index, xsav)
Definition i16crit.F:868
i8box
subroutine i8box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, size, xmsr, index, xsav)
Definition i16crit.F:992
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21