48
49
50
58
59
60
61#include "implicit_f.inc"
62
63
64
65#include "com04_c.inc"
66#include "sphcom.inc"
67#include "units_c.inc"
68#include "scr17_c.inc"
69#include "r2r_c.inc"
70
71
72
73 INTEGER ISPHIO(NISPHIO,*), IPART(LIPART1,*),
74 . NOD2SP(*),IPARTSP(*),
75 . ITAB(*),MFI,LWASPIO,ITABM1(*)
76 INTEGER, INTENT(IN) :: NRTRANS
78 . vsphio(*),x(3,*)
79 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
80 TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB
81 TYPE(SUBMODEL_DATA) LSUBMODEL(*)
82 my_real,
DIMENSION(NTRANSF,NRTRANS),
INTENT(IN) :: rtrans
83
84
85
86 INTEGER I, N, ID, J, IDS, IDPRT, IPRT, INOD, ICEL,
87 . IPROP, ISU, NSEG,
88 . MFITMP, IVAD, LVAD, ITYPE, IDSURF, IFTEMP,
89 . IFVITS,IFDENS, IFPRES, IFENER,SKIP,
90 . IAD,IN1,IN2,IN3,IN4, IUN,NBOX,NBOY,NBOZ,NBAND,SUB_ID
91 CHARACTER MESS*40
92 CHARACTER(LEN=NCHARTITLE) :: TITR
93 CHARACTER(LEN=NCHARKEY) :: KEY
95 . bid,rhoin,pin,ein,dist,
96 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
97 . x12,y12,z12,x13,y13,z13,nn,nx,ny,nz,
98 . dbucs,xbmin,ybmin,zbmin,xbmax,ybmax,zbmax,
99 . pinfini,carl,xx(9),fcut,
100 . rhoin_unit,pin_unit,ein_unit,xl(3)
101 LOGICAL IS_AVAILABLE,FOUND
102
103 DATA mess/'SPH INLET/OUTLET DEFINITION '/
104 DATA iun/1/
105
106
107
108 INTEGER USR2SYS
110 . get_u_geo
111 EXTERNAL get_u_geo
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153 lwaspio = 0
154 nseg_io = 0
155 is_available = .false.
156
157 mfitmp = mfi
158 ivad = 1
159
160
162
163
164 i = 0
165 DO n = 1,nsphio
166 skip = 0
167
168 IF (nsubdom > 0) THEN
170 ENDIF
171
172 IF (skip == 0) THEN
173
174
175 titr = ''
178 . option_titr = titr,
179 . submodel_id = sub_id)
180
181 i = i+1
183
184
185 CALL hm_get_intv(
'Itype' ,itype ,is_available,lsubmodel)
186 CALL hm_get_intv(
'pid' ,idprt ,is_available,lsubmodel)
187 CALL hm_get_intv(
'SURF_ID' ,idsurf ,is_available,lsubmodel)
188 CALL hm_get_floatv(
'DIST' ,dist ,is_available,lsubmodel, unitab)
189
190 CALL hm_get_intv(
'node_ID1',in1 ,is_available,lsubmodel)
191 CALL hm_get_intv(
'node_ID2',in2 ,is_available,lsubmodel)
192 CALL hm_get_intv(
'node_ID3',in3 ,is_available,lsubmodel)
193 CALL hm_get_floatv(
'Fcut' ,fcut ,is_available,lsubmodel, unitab)
194 isphio(1,i) =itype
195
196 ids = 0
197 found = .false.
198 DO j = 1,npart
199 IF (ipart(4,j) == idprt) THEN
200 ids = j
201 found = .true.
202 EXIT
203 ENDIF
204 ENDDO
205 IF (.NOT.found) THEN
207 . msgtype=msgerror,
208 . anmode=aninfo,
210 . c1=titr,
211 . i2=idprt)
212 ENDIF
213 isphio(2,i) = ids
214 isphio(4,i) = ivad
215
216 IF (idsurf > 0) THEN
217
218 ids=0
219 found = .false.
220 DO j=1,nsurf
221 IF (igrsurf(j)%ID == idsurf) THEN
222 ids=j
223 found = .true.
224 EXIT
225 ENDIF
226 ENDDO
227 IF (.NOT.found) THEN
229 . msgtype=msgerror,
230 . anmode=aninfo,
232 . c1=titr,
233 . i2=idsurf)
234 ENDIF
235 isphio(3,i) =ids
236 ELSEIF ((in1 == 0).AND.(in2 == 0).AND.(in3 == 0)) THEN
237
238 isphio(12,i) =1
239 CALL hm_get_floatv(
'XM' ,xx(1) ,is_available,lsubmodel, unitab)
240 CALL hm_get_floatv(
'YM' ,xx(2) ,is_available,lsubmodel, unitab)
241 CALL hm_get_floatv(
'ZM' ,xx(3) ,is_available,lsubmodel, unitab)
242 CALL hm_get_floatv(
'XM1' ,xx(4) ,is_available,lsubmodel, unitab)
243 CALL hm_get_floatv(
'YM1' ,xx(5) ,is_available,lsubmodel, unitab)
244 CALL hm_get_floatv(
'ZM1' ,xx(6) ,is_available,lsubmodel, unitab)
245 CALL hm_get_floatv(
'XM2' ,xx(7) ,is_available,lsubmodel, unitab)
246 CALL hm_get_floatv(
'YM2' ,xx(8) ,is_available,lsubmodel, unitab)
247 CALL hm_get_floatv(
'ZM2' ,xx(9) ,is_available,lsubmodel, unitab)
248 IF(sub_id /= 0) THEN
249 xl(1:3) = xx(1:3)
250 CALL subrotpoint(xl(1),xl(2),xl(3),rtrans,sub_id,lsubmodel)
251 xx(1:3) = xl(1:3)
252 xl(1:3) = xx(4:6)
253 CALL subrotpoint(xl(1),xl(2),xl(3),rtrans,sub_id,lsubmodel)
254 xx(4:6) = xl(1:3)
255 xl(1:3) = xx(7:9)
256 CALL subrotpoint(xl(1),xl(2),xl(3),rtrans,sub_id,lsubmodel)
257 xx(7:9) = xl(1:3)
258 END IF
259 DO j = 1,9
260 vsphio(ivad+3+j) = xx(j)
261 ENDDO
262 ELSE
263
264 isphio(12,i) = 2
265 isphio(13,i) =
usr2sys(in1,itabm1,mess,
id)
266 isphio(14,i) =
usr2sys(in2,itabm1,mess,
id)
267 isphio(15,i) =
usr2sys(in3,itabm1,mess,
id)
268 ENDIF
269
270 IF (itype == 1) THEN
271 nseg=igrsurf(ids)%NSEG
272 lvad=4+2*nseg
273 ELSE
274 lvad=22
275 ENDIF
276
277 IF (itype == 1)THEN
278 CALL hm_get_intv(
'FUN_A1' ,ifdens ,is_available,lsubmodel)
279 CALL hm_get_floatv(
'R0k0' ,rhoin ,is_available,lsubmodel,unitab)
280 IF ((rhoin == zero).AND.(ifdens>0)) THEN
282 rhoin = one * rhoin_unit
283 ENDIF
284 CALL hm_get_intv(
'FUN_A6' ,ifener ,is_available,lsubmodel)
285 CALL hm_get_floatv(
'MAT_E0',ein ,is_available,lsubmodel,unitab)
286 IF ((ein == zero).AND.(ifener>0)) THEN
288 ein = one * ein_unit
289 ENDIF
290 CALL hm_get_intv(
'FUN_A3' ,ifvits ,is_available,lsubmodel)
291 ELSEIF (itype == 2) THEN
292 CALL hm_get_intv(
'FUN_A2' ,ifpres ,is_available,lsubmodel)
293 CALL hm_get_floatv(
'MAT_P0' ,pin ,is_available,lsubmodel,unitab)
294 IF ((pin == zero).AND.(ifpres>0)) THEN
296 pin = one * pin_unit
297 ENDIF
298 ELSEIF (itype == 3) THEN
299 CALL hm_get_intv(
'FUN_A4' ,ifpres ,is_available,lsubmodel)
300 CALL hm_get_floatv(
'MAT_PScale',pin ,is_available,lsubmodel, unitab)
301 IF ((pin == zero).AND.(ifpres>0)) THEN
303 pin = one * pin_unit
304 ENDIF
305 CALL hm_get_floatv(
'Lc' ,carl ,is_available,lsubmodel, unitab)
306 ENDIF
307
308 IF (itype == 1) THEN
309 isphio(5,i) = ifdens
310 isphio(7,i) = ifener
311 isphio(8,i) = ifvits
312 ELSEIF (itype == 2) THEN
313 isphio(6,i) = ifpres
314 ELSEIF (itype == 3) THEN
315 isphio(6,i) = ifpres
316 ENDIF
317
318 IF (itype == 1) THEN
319 IF ((ifdens /= 0).AND.(rhoin == zero)) rhoin = one
320 IF ((ifener /= 0).AND.(ein == zero)) ein = one
321 vsphio(ivad ) = rhoin
322 vsphio(ivad+2) = ein
323 vsphio(ivad+3) = dist
324 ELSEIF (itype == 2) THEN
325 vsphio(ivad+1) = pin
326 vsphio(ivad+3) = dist
327 vsphio(ivad+15) = fcut
328 ELSEIF (itype == 3) THEN
329 vsphio(ivad+1) = pin
330 vsphio(ivad+2) = carl
331 vsphio(ivad+3) = dist
332 vsphio(ivad+15) = fcut
333 ELSEIF (itype == 4) THEN
334 vsphio(ivad+15) = fcut
335 ENDIF
336 mfitmp = mfitmp+lvad
337 ivad = ivad+lvad
338
339 ENDIF
340 ENDDO
341 mfi = mfitmp
342
343
344
345 CALL udouble(isphio(nisphio,1),nisphio,nsphio,mess,0,bid)
346
347
348
349 DO n=1,nsphio
350 IF(isphio(1,n)==1)THEN
351 ivad=isphio(4,n)+5
352 isu =isphio(3,n)
353 nseg=igrsurf(isu)%NSEG
354 DO j=1,nseg
355 in1=igrsurf(isu)%NODES(j,1)
356 in2=igrsurf(isu)%NODES(j,2)
357 in3=igrsurf(isu)%NODES(j,3)
358 in4=igrsurf(isu)%NODES(j,4)
359 x1=x(1,in1)
360 y1=x(2,in1)
361 z1=x(3,in1)
362 x2=x(1,in2)
363 y2=x(2,in2)
364 z2=x(3,in2)
365 x3=x(1,in3)
366 y3=x(2,in3)
367 z3=x(3,in3)
368 x12=x1-x2
369 y12=y1-y2
370 z12=z1-z2
371 x13=x3-x2
372 y13=y3-y2
373 z13=z3-z2
374 nx=y12*z13-z12*y13
375 ny=z12*x13-x12*z13
376 nz=x12*y13-y12*x13
377 nn =sqrt(nx*nx+ny*ny+nz*nz)
378 IF(in4/=in3)THEN
379 x4=x(1,in4)
380 y4=x(2,in4)
381 z4=x(3,in4)
382 x12=x1-x4
383 y12=y1-y4
384 z12=z1-z4
385 x13=x3-x4
386 y13=y3-y4
387 z13=z3-z4
388 nx=y12*z13-z12*y13
389 ny=z12*x13-x12*z13
390 nz=x12*y13-y12*x13
391 nn =nn+sqrt(nx*nx+ny*ny+nz*nz)
392 ENDIF
393 vsphio(ivad)=half*nn
394 ivad=ivad+2
395 ENDDO
396 ENDIF
397 ENDDO
398
399
400
401 DO n=1,nsphio
402 IF(isphio(12,n)==0)THEN
403 ivad=isphio(4,n)
404 isu =isphio(3,n)
405 nseg=igrsurf(isu)%NSEG
406 nseg_io = nseg_io + nseg
407 dist =vsphio(ivad+3)
408 dbucs=dist
409 xbmin =1.e+20
410 ybmin =1.e+20
411 zbmin =1.e+20
412 xbmax=-1.e+20
413 ybmax=-1.e+20
414 zbmax=-1.e+20
415 DO j=1,nseg
416 in1=igrsurf(isu)%NODES(j,1)
417 in2=igrsurf(isu)%NODES(j,2)
418 in3=igrsurf(isu)%NODES(j,3)
419 in4=igrsurf(isu)%NODES(j,4)
420 x1=x(1,in1)
421 y1=x(2,in1)
422 z1=x(3,in1)
423 x2=x(1,in2)
424 y2=x(2,in2)
425 z2=x(3,in2)
426 x3=x(1,in3)
427 y3=x(2,in3)
428 z3=x(3,in3)
447 dbucs=
max(dbucs,abs(x1-x2))
448 dbucs=
max(dbucs,abs(y1-y2))
449 dbucs=
max(dbucs,abs(z1-z2))
450 dbucs=
max(dbucs,abs(x2-x3))
451 dbucs=
max(dbucs,abs(y2-y3))
452 dbucs=
max(dbucs,abs(z2-z3))
453 dbucs=
max(dbucs,abs(x3-x1))
454 dbucs=
max(dbucs,abs(y3-y1))
455 dbucs=
max(dbucs,abs(z3-z1))
456 in4=igrsurf(isu)%NODES(j,4)
457 IF(in4/=in3)THEN
458 x4=x(1,in4)
459 y4=x(2,in4)
460 z4=x(3,in4)
467 dbucs=
max(dbucs,abs(x1-x4))
468 dbucs=
max(dbucs,abs(y1-y4))
469 dbucs=
max(dbucs,abs(z1-z4))
470 dbucs=
max(dbucs,abs(x2-x4))
471 dbucs=
max(dbucs,abs(y2-y4))
472 dbucs=
max(dbucs,abs(z2-z4))
473 dbucs=
max(dbucs,abs(x3-x4))
474 dbucs=
max(dbucs,abs(y3-y4))
475 dbucs=
max(dbucs,abs(z3-z4))
476 ENDIF
477 ENDDO
478 xbmin=xbmin-dist
479 ybmin=ybmin-dist
480 zbmin=zbmin-dist
481 xbmax=xbmax+dist
482 ybmax=ybmax+dist
483 zbmax=zbmax+dist
484 nbox =
max(iun,int((xbmax-xbmin)/dbucs))
485 nboy =
max(iun,int((ybmax-ybmin)/dbucs))
486 nboz =
max(iun,int((zbmax-zbmin)/dbucs))
487 nband=
max(nbox,nboy,nboz)+1
488 ELSE
489 nband = 2
490 ENDIF
491 lwaspio=
max(lwaspio,15*numsph+6*(nband+1)+12*nseg)
492 ENDDO
493 lwaspio=
max(lwaspio,3*nsphio)
494
495
496
497 WRITE(iout,1000)
498 DO n=1,nsphio
499 ivad=isphio(4,n)
500 IF(isphio(1,n)==1)THEN
501 WRITE(iout,1100) isphio(nisphio,n),
502 . ipart(4,isphio(2,n)),igrsurf(isphio(3,n))%ID,
503 . isphio(5,n),vsphio(ivad),isphio(7,n),
504 . vsphio(ivad+2),isphio(8,n),vsphio(ivad+3)
505 ELSE
506 IF(isphio(1,n)==2)THEN
507 IF (isphio(12,n)==0) THEN
508 WRITE(iout,1200) isphio(nisphio,n),
509 . ipart(4,isphio(2,n)),igrsurf(isphio(3,n))%ID,
510 . isphio(6,n),vsphio(ivad+1),vsphio(ivad+3)
511 ELSE
512 WRITE(iout,1400) isphio(nisphio,n),
513 . ipart(4,isphio(2,n)),isphio(6,n),vsphio(ivad+1),vsphio(ivad+3)
514 ENDIF
515 ELSEIF(isphio(1,n)==3)THEN
516 IF (isphio(12,n)==0) THEN
517 WRITE(iout,1300) isphio(nisphio,n),
518 . ipart(4,isphio(2,n)),igrsurf(isphio(3,n))%ID,vsphio(ivad+3),
519 . isphio(6,n),vsphio(ivad+1),vsphio(ivad+2)
520 ELSEIF (isphio(12,n)==1) THEN
521 WRITE(iout,1500) isphio(nisphio,n),ipart(4,isphio(2,n)),
522 . isphio(6,n),vsphio(ivad+1),vsphio(ivad+2)
523 ENDIF
524 ELSEIF(isphio(1,n)==4)THEN
525 IF (isphio(12,n)==0) THEN
526 WRITE(iout,1600) isphio(nisphio,n),ipart(4,isphio(2,n)),
527 . igrsurf(isphio(3,n))%ID
528 ELSE
529 WRITE(iout,1700) isphio(nisphio,n),ipart(4,isphio(2,n))
530 ENDIF
531 ENDIF
532 IF (isphio(12,n)==1) THEN
533
534 ivad=isphio(4,n)
535 WRITE(iout,2100) vsphio(ivad+4),vsphio(ivad+5),vsphio(ivad+6),
536 . vsphio(ivad+7),vsphio(ivad+8),vsphio(ivad+9),vsphio(ivad+10),
537 . vsphio(ivad+11),vsphio(ivad+12)
538 ELSEIF (isphio(12,n)==2) THEN
539
540 WRITE(iout,2200) itab(isphio(13,n)),itab(isphio(14,n)),itab(isphio(15,n))
541 ENDIF
542 IF (vsphio(ivad+15)>em20) WRITE(iout,2300) vsphio(ivad+15)
543 ENDIF
544 ENDDO
545
546
547 RETURN
548
549 1000 FORMAT(//
550 .' SPH INLET/OUTLET CONDITIONS '/
551 .' --------------------------- ')
552 1100 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
553 . /5x ,'TYPE INLET ',
554 . /10x,'PART RELATED TO CONDITION ',i10
555 . /10x,'INLET SURFACE ',i10
556 . /10x,'DENSITY FUNCTION ',i10
557 . /10x,'SCALE FACTOR ON DENSITY FUNCTION ',1pg20.13
558 . /10x,'ENERGY FUNCTION ',i10
559 . /10x,'SCALE FACTOR ON ENERGY FUNCTION ',1pg20.13
560 . /10x,'NORMAL VELOCITY FUNCTION ',i10
561 . /10x,'WITHIN DISTANCE FOR PARTICLES SETTING ',1pg20.13)
562 1200 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
563 . /5x ,'TYPE GENERAL OUTLET ',
564 . /10x,'PART RELATED TO CONDITION ',i10
565 . /10x,'OUTLET SURFACE ',i10
566 . /10x,'PRESSURE FUNCTION ',i10
567 . /10x,'SCALE FACTOR ON PRESSURE FUNCTION ',1pg20.13
568 . /10x,'WITHOUT DISTANCE FOR PARTICLES SETTING ',1pg20.13)
569 1300 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
570 . /5x ,'TYPE SILENT BOUNDARY ',
571 . /10x,'PART RELATED TO CONDITION ',i10
572 . /10x,'OUTLET SURFACE ',i10
573 . /10x,'WITHOUT DISTANCE FOR PARTICLES SETTING ',1pg20.13
574 . /10x,'PRESSURE FUNCTION ',i10
575 . /10x,'SCALE FACTOR ON PRESSURE FUNCTION ',1pg20.13
576 . /10x,'CHARACTERISTIC LENGTH ',1pg20.13)
577 1400 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
578 . /5x ,'TYPE GENERAL OUTLET ',
579 . /10x,'PART RELATED TO CONDITION ',i10
580 . /10x,'PRESSURE FUNCTION ',i10
581 . /10x,'SCALE FACTOR ON PRESSURE FUNCTION ',1pg20.13
582 . /10x,'WITHOUT DISTANCE FOR PARTICLES SETTING ',1pg20.13)
583 1500 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
584 . /5x ,'TYPE SILENT BOUNDARY ',
585 . /10x,'PART RELATED TO CONDITION ',i10
586 . /10x,'PRESSURE FUNCTION ',i10
587 . /10x,'SCALE FACTOR ON PRESSURE FUNCTION ',1pg20.13
588 . /10x,'WITHOUT DISTANCE FOR PARTICLES SETTING ',1pg20.13
589 . /10x,'CHARACTERISTIC LENGTH ',1pg20.13)
590 1600 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
591 . /5x ,'TYPE SPH CONTROL SECTION ',
592 . /10x,'PART RELATED TO CONTROL SECTION ',i10
593 . /10x,'SURFACE ',i10)
594 1700 FORMAT(/5x ,'SPH INLET/OUTLET CONDITION ID ',i10
595 . /5x ,'TYPE SPH CONTROL SECTION ',
596 . /10x,'PART RELATED TO CONTROL SECTION ',i10)
597
598 2100 FORMAT(10x,'surface defined by coordinates ',
599 . /10X,' --> coordinates of node1 ',1PG20.13,1PG20.13,1PG20.13,
600 . /10X,' --> coordinates of node2 ',1PG20.13,1PG20.13,1PG20.13,
601 . /10X,' --> coordinates of node3 ',1PG20.13,1PG20.13,1PG20.13)
602
603 2200 FORMAT(10X,'surface defined by nodes ',I10,I10,I10)
604
605 2300 FORMAT(10X,'4-pole
butterworth corner frequency
',1PG20.13)
606
607 RETURN
subroutine butterworth(dt, freq, x2, x1, x, fx2, fx1, fx)
subroutine hm_get_floatv(name, rval, is_available, lsubmodel, unitab)
subroutine hm_get_floatv_dim(name, dim_fac, is_available, lsubmodel, unitab)
subroutine hm_get_intv(name, ival, is_available, lsubmodel)
subroutine hm_option_start(entity_type)
integer, parameter nchartitle
integer, parameter ncharkey
integer, dimension(:), allocatable tagsphio
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
integer function usr2sys(iu, itabm1, mess, id)
subroutine udouble(list, ilist, nlist, mess, ir, rlist)
subroutine subrotpoint(x, y, z, rtrans, sub_id, lsubmodel)