FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
heat_LIB_DFLUX.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
8contains
9 !C************************************************************************
10 !C* DFLUX_111(NN,XX,YY,ZZ,ASECT,LTYPE,VAL,VECT)
11 !C* DFLUX_231(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
12 !C* DFLUX_232(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
13 !C* DFLUX_241(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
14 !C* DFLUX_242(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
15 !C* DFLUX_341(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
16 !C* DFLUX_342(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
17 !C* DFLUX_351(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
18 !C* DFLUX_352(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
19 !C* DFLUX_361(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
20 !C* DFLUX_362(NN,XX,YY,ZZ,LTYPE,VAL,VECT)
21 !C* DFLUX_731(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
22 !C* DFLUX_741(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT)
23 !C*--------------------------------------------------------------------*
24 subroutine heat_dflux_111(NN,XX,YY,ZZ,ASECT,LTYPE,val,VECT)
25 !C*--------------------------------------------------------------------*
26 !C***
27 !C*** SET 111 DFLUX
28 !C***
29 !C* LTYPE=0 : BODY FLUX
30 use hecmw
31 implicit real(kind=kreal)(a-h,o-z)
32 !C* I/F VARIABLES
33 integer(kind=kint) NN, LTYPE
34 real(kind=kreal) xx(nn),yy(nn),zz(nn),val,vect(nn)
35 !C*
36 !C*
37 if( ltype.EQ.0 ) then
38 asect = 1.0d0 * asect
39 else
40 asect = 0.0d0
41 endif
42
43 do i=1,nn
44 vect(i)=0.0
45 enddo
46 !C*
47 dx = xx(2) - xx(1)
48 dy = yy(2) - yy(1)
49 dz = zz(2) - zz(1)
50 al = dsqrt( dx*dx + dy*dy + dz*dz )
51 vv = asect * al
52 !
53 vect(1)= -val*vv/2.0d0
54 vect(2)= -val*vv/2.0d0
55 !
56 return
57
58 end subroutine heat_dflux_111
59 !C*--------------------------------------------------------------------*
60 subroutine heat_dflux_231(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
61 !C*--------------------------------------------------------------------*
62 !C***
63 !C*** SET 231 DFLUX
64 !C***
65 !C* LTYPE=0 : BODY FLUX
66 use hecmw
67 implicit real(kind=kreal) (a - h, o - z)
68 !C* I/F VARIABLES
69 integer(kind=kint) :: NN, LTYPE
70 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
71 real(kind=kreal) :: xg(2), wgt(2), h(4), hr(4), hs(4)
72 integer(kind=kint) :: NOD(2)
73 !*************************
74 ! GAUSS INTEGRATION POINT
75 !*************************
76 data xg/-0.5773502691896,0.5773502691896/
77 data wgt/1.0,1.0/
78 !C*
79 ivol=0
80 isuf=0
81 if( ltype.EQ.0 ) then
82 ivol=1
83 else if( ltype.GE.1 ) then
84 isuf=1
85 if(ltype.EQ.1) then
86 nod(1)=1
87 nod(2)=2
88 else if(ltype.EQ.2) then
89 nod(1)=2
90 nod(2)=3
91 else if(ltype.EQ.3) then
92 nod(1)=3
93 nod(2)=1
94 endif
95 endif
96 !** SURFACE LOAD
97 if( isuf.EQ.1 ) then
98 do i=1,nn
99 vect(i)=0.0
100 enddo
101 !** INTEGRATION OVER SURFACE
102 do lx=1,2
103 ri=xg(lx)
104 h(1)=0.5*(1.0-ri)
105 h(2)=0.5*(1.0+ri)
106 hr(1)=-0.5
107 hr(2)= 0.5
108 gx=0.0
109 gy=0.0
110 do i=1,2
111 gx=gx+hr(i)*xx(nod(i))
112 gy=gy+hr(i)*yy(nod(i))
113 enddo
114 xsum=dsqrt(gx*gx+gy*gy)*thick
115 do i=1,2
116 vect(nod(i))=vect(nod(i))-xsum*wgt(lx)*h(i)*val
117 enddo
118 enddo
119 endif
120 !** VOLUME LOAD
121 if( ivol.EQ.1 ) then
122 do i=1,nn
123 vect(i)=0.0
124 enddo
125
126 v1x=xx(2)-xx(1)
127 v1y=yy(2)-yy(1)
128 v1z=zz(2)-zz(1)
129 v2x=xx(3)-xx(1)
130 v2y=yy(3)-yy(1)
131 v2z=zz(3)-zz(1)
132 v3x= v1y*v2z-v1z*v2y
133 v3y= v1z*v2x-v1x*v2z
134 v3z= v1x*v2y-v1y*v2x
135
136 aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
137 vv=aa*thick
138 vect(1)= -val*vv/3.0
139 vect(2)= -val*vv/3.0
140 vect(3)= -val*vv/3.0
141
142 endif
143 return
144
145 end subroutine heat_dflux_231
146 !C*--------------------------------------------------------------------*
147 subroutine heat_dflux_232(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
148 !C*--------------------------------------------------------------------*
149 !C***
150 !C*** SET 232 DFLUX
151 !C***
152 !C* LTYPE=0 : BODY FLUX
153 use hecmw
154 implicit real(kind=kreal) (a - h, o - z)
155 !C* I/F VARIABLES
156 integer(kind=kint) :: NN, LTYPE
157 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
158 real(kind=kreal) :: xg(3), wgt(3), h(6), hr(3)
159 real(kind=kreal) :: hl1(6), hl2(6), hl3(6)
160 integer(kind=kint) :: NOD(3)
161 !*************************
162 ! GAUSS INTEGRATION POINT
163 !*************************
164 xg(1) = -0.7745966692
165 xg(2) = 0.0
166 xg(3) = 0.7745966692
167 wgt(1) = 0.5555555555
168 wgt(2) = 0.8888888888
169 wgt(3) = 0.5555555555
170 !C*
171 ivol=0
172 isuf=0
173 if( ltype.EQ.0 ) then
174 ivol=1
175 else if( ltype.GE.1 ) then
176 isuf=1
177 if(ltype.EQ.1) then
178 nod(1)=1
179 nod(2)=2
180 nod(3)=4
181 else if(ltype.EQ.2) then
182 nod(1)=2
183 nod(2)=3
184 nod(3)=5
185 else if(ltype.EQ.3) then
186 nod(1)=3
187 nod(2)=1
188 nod(3)=6
189 endif
190 endif
191 do i=1,nn
192 vect(i)=0.0
193 enddo
194 !** SURFACE LOAD
195 if( isuf.EQ.1 ) then
196 !** INTEGRATION OVER SURFACE
197 do lx=1,3
198 ri=xg(lx)
199 h(1)=-0.5*(1.0-ri)*ri
200 h(2)= 0.5*(1.0+ri)*ri
201 h(3)=1.0-ri*ri
202 hr(1)= ri-0.5
203 hr(2)= ri+0.5
204 hr(3)=-2.0*ri
205 gx=0.0
206 gy=0.0
207 do i=1,3
208 gx=gx+hr(i)*xx(nod(i))
209 gy=gy+hr(i)*yy(nod(i))
210 enddo
211 xsum=dsqrt(gx*gx+gy*gy)
212 wg=wgt(lx)*val*thick
213 vect(nod(1))=vect(nod(1))-h(1)*xsum*wg
214 vect(nod(2))=vect(nod(2))-h(2)*xsum*wg
215 vect(nod(3))=vect(nod(3))-h(3)*xsum*wg
216 enddo
217 endif
218 !** VOLUME LOAD
219 if( ivol.EQ.1 ) then
220
221 !* LOOP OVER ALL INTEGRATION POINTS
222 do l2=1,3
223 xl2=xg(l2)
224 x2 =(xl2+1.0)*0.5
225 do l1=1,3
226 xl1=xg(l1)
227 x1=0.5*(1.0-x2)*(xl1+1.0)
228 ! INTERPOLATION FUNCTION
229 x3=1.0-x1-x2
230 h(1)= x1*(2.0*x1-1.)
231 h(2)= x2*(2.0*x2-1.)
232 h(3)= x3*(2.0*x3-1.)
233 h(4)= 4.0*x1*x2
234 h(5)= 4.0*x2*x3
235 h(6)= 4.0*x1*x3
236 ! DERIVATIVE OF INTERPOLATION FUNCTION
237 ! FOR L1-COORDINATE
238 hl1(1)=4.0*x1-1.0
239 hl1(2)= 0.0
240 hl1(3)= 0.0
241 hl1(4)= 4.0*x2
242 hl1(5)= 0.0
243 hl1(6)= 4.0*x3
244 ! FOR L2-COORDINATE
245 hl2(1)= 0.0
246 hl2(2)= 4.0*x2-1.0
247 hl2(3)= 0.0
248 hl2(4)= 4.0*x1
249 hl2(5)= 4.0*x3
250 hl2(6)= 0.0
251 ! FOR L3-COORDINATE
252 hl3(1)= 0.0
253 hl3(2)= 0.0
254 hl3(3)= 4.0*x3-1.0
255 hl3(4)= 0.0
256 hl3(5)= 4.0*x2
257 hl3(6)= 4.0*x1
258 ! JACOBI MATRIX
259 xj11=0.0
260 xj21=0.0
261 xj12=0.0
262 xj22=0.0
263 do i=1,nn
264 xj11=xj11+(hl1(i)-hl3(i))*xx(i)
265 xj21=xj21+(hl2(i)-hl3(i))*xx(i)
266 xj12=xj12+(hl1(i)-hl3(i))*yy(i)
267 xj22=xj22+(hl2(i)-hl3(i))*yy(i)
268 enddo
269 det=xj11*xj22-xj21*xj12
270 wg=wgt(l1)*wgt(l2)*det*thick*(1.0-x2)*0.25
271 do i = 1, nn
272 vect(i) = vect(i) - h(i)*wg*val
273 enddo
274 enddo
275 enddo
276 endif
277 return
278
279 end subroutine heat_dflux_232
280 !C*--------------------------------------------------------------------*
281 subroutine heat_dflux_241(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
282 !C*--------------------------------------------------------------------*
283 !C***
284 !C*** SET 241 DFLUX
285 !C***
286 !C* LTYPE=0 : BODY FLUX
287 use hecmw
288 implicit real(kind=kreal) (a - h, o - z)
289 !C* I/F VARIABLES
290 integer(kind=kint) :: NN, LTYPE
291 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
292 real(kind=kreal) :: xg(2), wgt(2), h(4), hr(4), hs(4)
293 integer(kind=kint) :: NOD(2)
294 !*************************
295 ! GAUSS INTEGRATION POINT
296 !*************************
297 data xg/-0.5773502691896,0.5773502691896/
298 data wgt/1.0,1.0/
299 !C*
300 ivol=0
301 isuf=0
302 if( ltype.EQ.0 ) then
303 ivol=1
304 else if( ltype.GE.1 ) then
305 isuf=1
306 if(ltype.EQ.1) then
307 nod(1)=1
308 nod(2)=2
309 else if(ltype.EQ.2) then
310 nod(1)=2
311 nod(2)=3
312 else if(ltype.EQ.3) then
313 nod(1)=3
314 nod(2)=4
315 else if(ltype.EQ.4) then
316 nod(1)=4
317 nod(2)=1
318 endif
319 endif
320 do i=1,nn
321 vect(i)=0.0
322 enddo
323 !** SURFACE LOAD
324 if( isuf.EQ.1 ) then
325 !** INTEGRATION OVER SURFACE
326 do lx=1,2
327 ri=xg(lx)
328 h(1)=0.5*(1.0-ri)
329 h(2)=0.5*(1.0+ri)
330 hr(1)=-0.5
331 hr(2)= 0.5
332 gx=0.0
333 gy=0.0
334 do i=1,2
335 gx=gx+hr(i)*xx(nod(i))
336 gy=gy+hr(i)*yy(nod(i))
337 enddo
338 xsum=dsqrt(gx*gx+gy*gy)*thick
339 do i=1,2
340 vect(nod(i))=vect(nod(i))-xsum*wgt(lx)*h(i)*val
341 enddo
342 enddo
343 endif
344 !** VOLUME LOAD
345 if( ivol.EQ.1 ) then
346 do lx=1,2
347 ri=xg(lx)
348 do ly=1,2
349 si=xg(ly)
350 rp=1.0+ri
351 sp=1.0+si
352 rm=1.0-ri
353 sm=1.0-si
354 h(1)=0.25*rm*sm
355 h(2)=0.25*rp*sm
356 h(3)=0.25*rp*sp
357 h(4)=0.25*rm*sp
358 hr(1)=-.25*sm
359 hr(2)= .25*sm
360 hr(3)= .25*sp
361 hr(4)=-.25*sp
362 hs(1)=-.25*rm
363 hs(2)=-.25*rp
364 hs(3)= .25*rp
365 hs(4)= .25*rm
366 xj11=0.0
367 xj21=0.0
368 xj12=0.0
369 xj22=0.0
370 do i=1,nn
371 xj11=xj11+hr(i)*xx(i)
372 xj21=xj21+hs(i)*xx(i)
373 xj12=xj12+hr(i)*yy(i)
374 xj22=xj22+hs(i)*yy(i)
375 enddo
376 det=xj11*xj22-xj21*xj12
377 wg=wgt(lx)*wgt(ly)*det*thick
378 do i =1,nn
379 vect(i)=vect(i)-wg*h(i)*val
380 enddo
381 enddo
382 enddo
383 endif
384 return
385
386 end subroutine heat_dflux_241
387 !C*--------------------------------------------------------------------*
388 subroutine heat_dflux_242(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
389 !C*--------------------------------------------------------------------*
390 !C***
391 !C*** SET 242 DFLUX
392 !C***
393 !C* LTYPE=0 : BODY FLUX
394 use hecmw
395 implicit real(kind=kreal) (a - h, o - z)
396 !C* I/F VARIABLES
397 integer(kind=kint) :: NN, LTYPE
398 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), thick, val, vect(nn)
399 real(kind=kreal) :: xg(3), wgt(3), h(8), hr(8), hs(8)
400 integer(kind=kint) :: NOD(3)
401 !*************************
402 ! GAUSS INTEGRATION POINT
403 !*************************
404 xg(1) = -0.7745966692
405 xg(2) = 0.0
406 xg(3) = 0.7745966692
407 wgt(1) = 0.5555555555
408 wgt(2) = 0.8888888888
409 wgt(3) = 0.5555555555
410 !C*
411 ivol=0
412 isuf=0
413 if( ltype.EQ.0 ) then
414 ivol=1
415 else if( ltype.GE.1 ) then
416 isuf=1
417 if(ltype.EQ.1) then
418 nod(1)=1
419 nod(2)=2
420 nod(3)=5
421 else if(ltype.EQ.2) then
422 nod(1)=2
423 nod(2)=3
424 nod(3)=6
425 else if(ltype.EQ.3) then
426 nod(1)=3
427 nod(2)=4
428 nod(3)=7
429 else if(ltype.EQ.4) then
430 nod(1)=4
431 nod(2)=1
432 nod(3)=8
433 endif
434 endif
435 do i=1,nn
436 vect(i)=0.0
437 enddo
438 !** SURFACE LOAD
439 if( isuf.EQ.1 ) then
440 !** INTEGRATION OVER SURFACE
441 do lx=1,3
442 ri=xg(lx)
443 h(1)=-0.5*(1.0-ri)*ri
444 h(2)= 0.5*(1.0+ri)*ri
445 h(3)=1.0-ri*ri
446 hr(1)= ri-0.5
447 hr(2)= ri+0.5
448 hr(3)=-2.0*ri
449 gx=0.0
450 gy=0.0
451 do i=1,3
452 gx=gx+hr(i)*xx(nod(i))
453 gy=gy+hr(i)*yy(nod(i))
454 enddo
455 xsum=dsqrt(gx*gx+gy*gy)*thick
456 vect(nod(1))=vect(nod(1))-h(1)*wgt(lx)*val*xsum
457 vect(nod(2))=vect(nod(2))-h(2)*wgt(lx)*val*xsum
458 vect(nod(3))=vect(nod(3))-h(3)*wgt(lx)*val*xsum
459 enddo
460 endif
461 !** VOLUME LOAD
462 if( ivol.EQ.1 ) then
463 do lx=1,3
464 ri=xg(lx)
465 do ly=1,3
466 si=xg(ly)
467 rp=1.0+ri
468 sp=1.0+si
469 rm=1.0-ri
470 sm=1.0-si
471 h(1)=0.25*rm*sm*(-1.0-ri-si)
472 h(2)=0.25*rp*sm*(-1.0+ri-si)
473 h(3)=0.25*rp*sp*(-1.0+ri+si)
474 h(4)=0.25*rm*sp*(-1.0-ri+si)
475 h(5)=0.5*(1.0-ri*ri)*(1.0-si)
476 h(6)=0.5*(1.0-si*si)*(1.0+ri)
477 h(7)=0.5*(1.0-ri*ri)*(1.0+si)
478 h(8)=0.5*(1.0-si*si)*(1.0-ri)
479 hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
480 hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
481 hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
482 hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
483 hr(5)=-ri*(1.0-si)
484 hr(6)= 0.5*(1.0-si*si)
485 hr(7)=-ri*(1.0+si)
486 hr(8)=-0.5*(1.0-si*si)
487 hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
488 hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
489 hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
490 hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
491 hs(5)=-0.5*(1.0-ri*ri)
492 hs(6)=-si *(1.0+ri)
493 hs(7)= 0.5*(1.0-ri*ri)
494 hs(8)= -si*(1.0-ri)
495 xj11=0.0
496 xj21=0.0
497 xj12=0.0
498 xj22=0.0
499 do i=1,nn
500 xj11=xj11+hr(i)*xx(i)
501 xj21=xj21+hs(i)*xx(i)
502 xj12=xj12+hr(i)*yy(i)
503 xj22=xj22+hs(i)*yy(i)
504 enddo
505 det=xj11*xj22-xj21*xj12
506 wg=wgt(lx)*wgt(ly)*det*thick
507 do i =1,nn
508 vect(i)=vect(i)-wg*h(i)*val
509 enddo
510 enddo
511 enddo
512 endif
513 return
514
515 end subroutine heat_dflux_242
516 !----------------------------------------------------------------------*
517 subroutine heat_dflux_341(NN,XX,YY,ZZ,LTYPE,val,VECT)
518 !----------------------------------------------------------------------*
519 !**
520 !** SET 341 DFLUX
521 !**
522 ! BF LTYPE=0 :BODY FLUX
523 ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
524 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
525 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
526 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
527 use hecmw
528 implicit none
529 ! I/F VARIABLES
530 integer(kind=kint) :: NN, LTYPE
531 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), vect(nn)
532 ! LOCAL VARIABLES
533 real(kind=kreal) :: h(4), hr(4), hs(4), ht(4)
534 real(kind=kreal) :: xg(2), wgt(2)
535 real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
536 real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33, det, wg
537 integer(kind=kint) :: IVOL, ISUF
538 integer(kind=kint) :: NOD(4)
539 integer(kind=kint) :: L1, L2, L3, I
540 real(kind=kreal) :: val, aa
541 real(kind=kreal) :: v1x, v1y, v1z
542 real(kind=kreal) :: v2x, v2y, v2z
543 real(kind=kreal) :: v3x, v3y, v3z
544 real(kind=kreal) :: xl1, xl2, xl3
545 real(kind=kreal) :: x1, x2, x3
546 !*************************
547 ! GAUSS INTEGRATION POINT
548 !*************************
549 data xg/-0.5773502691896,0.5773502691896/
550 data wgt/1.0,1.0/
551 !
552 ivol=0
553 isuf=0
554 nod = 0
555 if( ltype.EQ.0 ) then
556 ivol=1
557 else
558 isuf=1
559 if(ltype.EQ.1) then
560 nod(1)=1
561 nod(2)=2
562 nod(3)=3
563 else if(ltype.EQ.2) then
564 nod(1)=4
565 nod(2)=2
566 nod(3)=1
567 else if(ltype.EQ.3) then
568 nod(1)=4
569 nod(2)=3
570 nod(3)=2
571 else if(ltype.EQ.4) then
572 nod(1)=4
573 nod(2)=1
574 nod(3)=3
575 endif
576 endif
577 ! CLEAR VECT
578 do i=1,nn
579 vect(i)=0.0
580 enddo
581 !** SURFACE LOAD
582 if( isuf.EQ.1 ) then
583 v1x=xx(nod(2))-xx(nod(1))
584 v1y=yy(nod(2))-yy(nod(1))
585 v1z=zz(nod(2))-zz(nod(1))
586 v2x=xx(nod(3))-xx(nod(1))
587 v2y=yy(nod(3))-yy(nod(1))
588 v2z=zz(nod(3))-zz(nod(1))
589 v3x= v1y*v2z-v1z*v2y
590 v3y= v1z*v2x-v1x*v2z
591 v3z= v1x*v2y-v1y*v2x
592 aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
593 vect(nod(1))= -val*aa/3.0
594 vect(nod(2))= -val*aa/3.0
595 vect(nod(3))= -val*aa/3.0
596 endif
597 !** VOLUME LOAD
598 if( ivol.EQ.1 ) then
599 ! LOOP FOR INTEGRATION POINTS
600 do l3=1,2
601 xl3=xg(l3)
602 x3 =(xl3+1.0)*0.5
603 do l2=1,2
604 xl2=xg(l2)
605 x2 =(1.0-x3)*(xl2+1.0)*0.5
606 do l1=1,2
607 xl1=xg(l1)
608 x1=(1.0-x2-x3)*(xl1+1.0)*0.5
609 ! INTERPOLATION FUNCTION
610 h(1)=x1
611 h(2)=x2
612 h(3)=x3
613 h(4)=1.0-x1-x2-x3
614 ! DERIVATIVE OF INTERPOLATION FUNCTION
615 ! FOR L1-COORDINATE
616 hr(1)= 1.0
617 hr(2)= 0.0
618 hr(3)= 0.0
619 hr(4)=-1.0
620 ! FOR L2-COORDINATE
621 hs(1)= 0.0
622 hs(2)= 1.0
623 hs(3)= 0.0
624 hs(4)=-1.0
625 ! FOR ZETA-COORDINATE
626 ht(1)= 0.0
627 ht(2)= 0.0
628 ht(3)= 1.0
629 ht(4)=-1.0
630 ! JACOBI MATRIX
631 xj11=hr(1)*xx(1)+hr(2)*xx(2)+hr(3)*xx(3)+hr(4)*xx(4)
632 xj21=hs(1)*xx(1)+hs(2)*xx(2)+hs(3)*xx(3)+hs(4)*xx(4)
633 xj31=ht(1)*xx(1)+ht(2)*xx(2)+ht(3)*xx(3)+ht(4)*xx(4)
634 xj12=hr(1)*yy(1)+hr(2)*yy(2)+hr(3)*yy(3)+hr(4)*yy(4)
635 xj22=hs(1)*yy(1)+hs(2)*yy(2)+hs(3)*yy(3)+hs(4)*yy(4)
636 xj32=ht(1)*yy(1)+ht(2)*yy(2)+ht(3)*yy(3)+ht(4)*yy(4)
637 xj13=hr(1)*zz(1)+hr(2)*zz(2)+hr(3)*zz(3)+hr(4)*zz(4)
638 xj23=hs(1)*zz(1)+hs(2)*zz(2)+hs(3)*zz(3)+hs(4)*zz(4)
639 xj33=ht(1)*zz(1)+ht(2)*zz(2)+ht(3)*zz(3)+ht(4)*zz(4)
640 ! DETERMINANT OF JACOBIAN
641 det=xj11*xj22*xj33 &
642 +xj12*xj23*xj31 &
643 +xj13*xj21*xj32 &
644 -xj13*xj22*xj31 &
645 -xj12*xj21*xj33 &
646 -xj11*xj23*xj32
647 !
648 do i = 1, nn
649 vect(i) = vect(i) + val*h(i)*det*(1.0-x3)*(1.0-x2-x3)*0.125
650 enddo
651
652 enddo
653 enddo
654 enddo
655 endif
656 !
657 return
658
659 end subroutine heat_dflux_341
660 !----------------------------------------------------------------------*
661 subroutine heat_dflux_342(NN,XX,YY,ZZ,LTYPE,val,VECT)
662 !----------------------------------------------------------------------*
663 !**
664 !** SET 342 DFLUX
665 !**
666 ! BF LTYPE=0 :BODY FLUX
667 ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
668 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
669 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
670 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
671 use hecmw
672 implicit none
673 ! I/F VARIABLES
674 integer(kind=kint) :: NN, LTYPE
675 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), vect(nn)
676 ! LOCAL VARIABLES
677 real(kind=kreal) :: h(10)
678 real(kind=kreal) :: hl1(10), hl2(10), hl3(10), hl4(10)
679 real(kind=kreal) :: xg(3), wgt(3)
680 real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
681 real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
682 real(kind=kreal) :: det, wg
683 integer(kind=kint) :: IVOL, ISUF
684 integer(kind=kint) :: NOD(10)
685 integer(kind=kint) :: LX, LY, LZ, I, IG1, IG2, L1, L2, L3
686 real(kind=kreal) :: val, xsum
687 real(kind=kreal) :: g1x, g1y, g1z
688 real(kind=kreal) :: g2x, g2y, g2z
689 real(kind=kreal) :: g3x, g3y, g3z
690 real(kind=kreal) :: xl1, xl2, xl3
691 real(kind=kreal) :: x1, x2, x3, x4
692 !*************************
693 ! GAUSS INTEGRATION POINT
694 !*************************
695 xg(1) = -0.7745966692
696 xg(2) = 0.0
697 xg(3) = 0.7745966692
698 wgt(1) = 0.5555555555
699 wgt(2) = 0.8888888888
700 wgt(3) = 0.5555555555
701 !
702 ivol=0
703 isuf=0
704 if( ltype.EQ.0 ) then
705 ivol=1
706 else
707 isuf=1
708 if(ltype.EQ.1) then
709 nod(1)=1
710 nod(2)=2
711 nod(3)=3
712 nod(4)=5
713 nod(5)=6
714 nod(6)=7
715 else if(ltype.EQ.2) then
716 nod(1)=4
717 nod(2)=2
718 nod(3)=1
719 nod(4)=9
720 nod(5)=5
721 nod(6)=8
722 else if(ltype.EQ.3) then
723 nod(1)=4
724 nod(2)=3
725 nod(3)=2
726 nod(4)=10
727 nod(5)=6
728 nod(6)=9
729 else if(ltype.EQ.4) then
730 nod(1)=4
731 nod(2)=1
732 nod(3)=3
733 nod(4)=8
734 nod(5)=7
735 nod(6)=10
736 endif
737 endif
738 ! CLEAR VECT
739 do i=1,10
740 vect(i)=0.0
741 enddo
742 !** SURFACE LOAD
743 if( isuf.EQ.1 ) then
744 ! INTEGRATION OVER SURFACE
745 do l2=1,3
746 xl2=xg(l2)
747 x2 =(xl2+1.0)*0.5
748 do l1=1,3
749 xl1=xg(l1)
750 x1=0.5*(1.0-x2)*(xl1+1.0)
751 ! INTERPOLATION FUNCTION
752 x3=1.0-x1-x2
753 h(1)= x1*(2.0*x1-1.)
754 h(2)= x2*(2.0*x2-1.)
755 h(3)= x3*(2.0*x3-1.)
756 h(4)= 4.0*x1*x2
757 h(5)= 4.0*x2*x3
758 h(6)= 4.0*x1*x3
759 ! DERIVATIVE OF INTERPOLATION FUNCTION
760 ! FOR L1-COORDINATE
761 hl1(1)=4.0*x1-1.0
762 hl1(2)= 0.0
763 hl1(3)= 0.0
764 hl1(4)= 4.0*x2
765 hl1(5)= 0.0
766 hl1(6)= 4.0*x3
767 ! FOR L2-COORDINATE
768 hl2(1)= 0.0
769 hl2(2)= 4.0*x2-1.0
770 hl2(3)= 0.0
771 hl2(4)= 4.0*x1
772 hl2(5)= 4.0*x3
773 hl2(6)= 0.0
774 ! FOR L3-COORDINATE
775 hl3(1)= 0.0
776 hl3(2)= 0.0
777 hl3(3)= 4.0*x3-1.0
778 hl3(4)= 0.0
779 hl3(5)= 4.0*x2
780 hl3(6)= 4.0*x1
781 ! JACOBI MATRIX
782 g1x=0.0
783 g1y=0.0
784 g1z=0.0
785 g2x=0.0
786 g2y=0.0
787 g2z=0.0
788 do i=1,6
789 g1x=g1x+(hl1(i)-hl3(i))*xx(nod(i))
790 g1y=g1y+(hl1(i)-hl3(i))*yy(nod(i))
791 g1z=g1z+(hl1(i)-hl3(i))*zz(nod(i))
792 g2x=g2x+(hl2(i)-hl3(i))*xx(nod(i))
793 g2y=g2y+(hl2(i)-hl3(i))*yy(nod(i))
794 g2z=g2z+(hl2(i)-hl3(i))*zz(nod(i))
795 enddo
796 g3x=g1y*g2z-g1z*g2y
797 g3y=g1z*g2x-g1x*g2z
798 g3z=g1x*g2y-g1y*g2x
799 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
800 g3x=g3x/xsum
801 g3y=g3y/xsum
802 g3z=g3z/xsum
803 ! JACOBI MATRIX
804 xj11=g1x
805 xj12=g1y
806 xj13=g1z
807 xj21=g2x
808 xj22=g2y
809 xj23=g2z
810 xj31=g3x
811 xj32=g3y
812 xj33=g3z
813 !DETERMINANT OF JACOBIAN
814 det=xj11*xj22*xj33 &
815 +xj12*xj23*xj31 &
816 +xj13*xj21*xj32 &
817 -xj13*xj22*xj31 &
818 -xj12*xj21*xj33 &
819 -xj11*xj23*xj32
820 wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25
821 do i=1,6
822 vect(nod(i))=vect(nod(i))-val*wg*h(i)
823 enddo
824 !
825 enddo
826 enddo
827 endif
828 !** VOLUME LOAD
829 if( ivol.EQ.1 ) then
830 ! LOOP FOR INTEGRATION POINTS
831 !DETERMINANT OF JACOBIAN
832 do l3=1,3
833 xl3=xg(l3)
834 x3 =(xl3+1.0)*0.5
835 do l2=1,3
836 xl2=xg(l2)
837 x2 =(1.0-x3)*(xl2+1.0)*0.5
838 do l1=1,3
839 xl1=xg(l1)
840 x1=(1.0-x2-x3)*(xl1+1.0)*0.5
841 !C* INTERPOLATION FUNCTION
842 x4=1.0-x1-x2-x3
843 h(1)=x1*(2.0*x1-1.0)
844 h(2)=x2*(2.0*x2-1.0)
845 h(3)=x3*(2.0*x3-1.0)
846 h(4)=x4*(2.0*x4-1.0)
847 h(5)=4.0*x1*x2
848 h(6)=4.0*x2*x3
849 h(7)=4.0*x1*x3
850 h(8)=4.0*x1*x4
851 h(9)=4.0*x2*x4
852 h(10)=4.0*x3*x4
853 !C* DERIVATIVE OF INTERPOLATION FUNCTION
854 !C* FOR L1-COORDINATE
855 hl1(1)= 4.0*x1-1.0
856 hl1(2)= 0.0
857 hl1(3)= 0.0
858 hl1(4)= 0.0
859 hl1(5)= 4.0*x2
860 hl1(6)= 0.0
861 hl1(7)= 4.0*x3
862 hl1(8)= 4.0*x4
863 hl1(9)= 0.0
864 hl1(10)= 0.0
865 !C* FOR L2-COORDINATE
866 hl2(1)= 0.0
867 hl2(2)= 4.0*x2-1.0
868 hl2(3)= 0.0
869 hl2(4)= 0.0
870 hl2(5)= 4.0*x1
871 hl2(6)= 4.0*x3
872 hl2(7)= 0.0
873 hl2(8)= 0.0
874 hl2(9)= 4.0*x4
875 hl2(10)= 0.0
876 !C* FOR L3-COORDINATE
877 hl3(1)= 0.0
878 hl3(2)= 0.0
879 hl3(3)= 4.0*x3-1.0
880 hl3(4)= 0.0
881 hl3(5)= 0.0
882 hl3(6)= 4.0*x2
883 hl3(7)= 4.0*x1
884 hl3(8)= 0.0
885 hl3(9)= 0.0
886 hl3(10)= 4.0*x4
887 !C* FOR L4-COORDINATE
888 hl4(1)= 0.0
889 hl4(2)= 0.0
890 hl4(3)= 0.0
891 hl4(4)= 4.0*x4-1.0
892 hl4(5)= 0.0
893 hl4(6)= 0.0
894 hl4(7)= 0.0
895 hl4(8)= 4.0*x1
896 hl4(9)= 4.0*x2
897 hl4(10)= 4.0*x3
898 !C* JACOBI MATRIX
899 xj11=0.0
900 xj21=0.0
901 xj31=0.0
902 xj12=0.0
903 xj22=0.0
904 xj32=0.0
905 xj13=0.0
906 xj23=0.0
907 xj33=0.0
908 do i=1,nn
909 xj11=xj11+(hl1(i)-hl4(i))*xx(i)
910 xj21=xj21+(hl2(i)-hl4(i))*xx(i)
911 xj31=xj31+(hl3(i)-hl4(i))*xx(i)
912 xj12=xj12+(hl1(i)-hl4(i))*yy(i)
913 xj22=xj22+(hl2(i)-hl4(i))*yy(i)
914 xj32=xj32+(hl3(i)-hl4(i))*yy(i)
915 xj13=xj13+(hl1(i)-hl4(i))*zz(i)
916 xj23=xj23+(hl2(i)-hl4(i))*zz(i)
917 xj33=xj33+(hl3(i)-hl4(i))*zz(i)
918 enddo
919 det=xj11*xj22*xj33 &
920 +xj12*xj23*xj31 &
921 +xj13*xj21*xj32 &
922 -xj13*xj22*xj31 &
923 -xj12*xj21*xj33 &
924 -xj11*xj23*xj32
925 det=-det
926 wg=det*wgt(l1)*wgt(l2)*wgt(l3)*(1.-x3)*(1.0-x2-x3)*0.125
927 do i = 1, nn
928 vect(i) = vect(i) - val*wg*h(i)
929 enddo
930 !
931 enddo
932 enddo
933 enddo
934 endif
935 !
936 return
937
938 end subroutine heat_dflux_342
939 !----------------------------------------------------------------------*
940 subroutine heat_dflux_351(NN,XX,YY,ZZ,LTYPE,val,VECT)
941 !----------------------------------------------------------------------*
942 !**
943 !** SET 351 DFLUX
944 !**
945 ! BF LTYPE=0 :BODY FLUX
946 ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
947 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
948 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
949 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
950 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
951 use hecmw
952 implicit real(kind=kreal)(a-h,o-z)
953 ! I/F VARIABLES
954 integer(kind=kint) :: NN, LTYPE
955 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
956 ! LOCAL VARIABLES
957 real(kind=kreal) :: h(6), hr(6), hs(6), ht(6), pl(6)
958 real(kind=kreal) :: xg(2), wgt(2), xg1(3), xg2(3), wgt1(3)
959 integer(kind=kint) :: NOD(4)
960 !*************************
961 ! GAUSS INTEGRATION POINT
962 !*************************
963 data xg / -0.5773502691896,0.5773502691896 /
964 data wgt / 1.0, 1.0 /
965 data xg1 / 0.6666666667, 0.16666666667, 0.16666666667 /
966 data xg2 / 0.1666666667, 0.66666666667, 0.16666666667 /
967 data wgt1/ 0.1666666667, 0.16666666667, 0.16666666667 /
968 !
969 ! SELCTION OF LOAD TYPE
970 ivol=0
971 isuf=0
972 if( ltype.EQ.0 ) then
973 ivol=1
974 else
975 isuf=1
976 if( ltype.EQ.1 ) then
977 nod(1) = 1
978 nod(2) = 2
979 nod(3) = 3
980 isuf=2
981 else if( ltype.EQ.2 ) then
982 nod(1) = 6
983 nod(2) = 5
984 nod(3) = 4
985 isuf=2
986 else if( ltype.EQ.3 ) then
987 nod(1) = 4
988 nod(2) = 5
989 nod(3) = 2
990 nod(4) = 1
991 else if( ltype.EQ.4 ) then
992 nod(1) = 5
993 nod(2) = 6
994 nod(3) = 3
995 nod(4) = 2
996 else if( ltype.EQ.5 ) then
997 nod(1) = 6
998 nod(2) = 4
999 nod(3) = 1
1000 nod(4) = 3
1001 endif
1002 endif
1003 do i=1,nn
1004 vect(i)=0.0
1005 enddo
1006 !
1007 !** SURFACE LOAD
1008 !
1009 if( isuf.EQ.1 ) then
1010 do ig2=1,2
1011 si=xg(ig2)
1012 do ig1=1,2
1013 ri=xg(ig1)
1014 !
1015 h(1)=0.25*(1.0-ri)*(1.0-si)
1016 h(2)=0.25*(1.0+ri)*(1.0-si)
1017 h(3)=0.25*(1.0+ri)*(1.0+si)
1018 h(4)=0.25*(1.0-ri)*(1.0+si)
1019 hr(1)=-.25*(1.0-si)
1020 hr(2)= .25*(1.0-si)
1021 hr(3)= .25*(1.0+si)
1022 hr(4)=-.25*(1.0+si)
1023 hs(1)=-.25*(1.0-ri)
1024 hs(2)=-.25*(1.0+ri)
1025 hs(3)= .25*(1.0+ri)
1026 hs(4)= .25*(1.0-ri)
1027 !
1028 g1x=0.0
1029 g1y=0.0
1030 g1z=0.0
1031 g2x=0.0
1032 g2y=0.0
1033 g2z=0.0
1034 !
1035 do i=1,4
1036 g1x=g1x+hr(i)*xx(nod(i))
1037 g1y=g1y+hr(i)*yy(nod(i))
1038 g1z=g1z+hr(i)*zz(nod(i))
1039 g2x=g2x+hs(i)*xx(nod(i))
1040 g2y=g2y+hs(i)*yy(nod(i))
1041 g2z=g2z+hs(i)*zz(nod(i))
1042 enddo
1043 !
1044 g3x=g1y*g2z-g1z*g2y
1045 g3y=g1z*g2x-g1x*g2z
1046 g3z=g1x*g2y-g1y*g2x
1047 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1048 do i=1,4
1049 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1050 enddo
1051 !
1052 enddo
1053 enddo
1054 elseif( isuf.EQ.2 ) then
1055 v1x=xx(nod(2))-xx(nod(1))
1056 v1y=yy(nod(2))-yy(nod(1))
1057 v1z=zz(nod(2))-zz(nod(1))
1058 v2x=xx(nod(3))-xx(nod(1))
1059 v2y=yy(nod(3))-yy(nod(1))
1060 v2z=zz(nod(3))-zz(nod(1))
1061 v3x= v1y*v2z-v1z*v2y
1062 v3y= v1z*v2x-v1x*v2z
1063 v3z= v1x*v2y-v1y*v2x
1064 aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
1065 vect(nod(1))= -val*aa/3.0
1066 vect(nod(2))= -val*aa/3.0
1067 vect(nod(3))= -val*aa/3.0
1068 endif
1069 !
1070 !** VOLUME LOAD
1071 !
1072 if( ivol.EQ.1 ) then
1073 do i=1,6
1074 vect(i)=0.0
1075 enddo
1076 !
1077 do lz=1,2
1078 zi=xg(lz)
1079 do l12=1,3
1080 x1=xg1(l12)
1081 x2=xg2(l12)
1082
1083 h(1)=x1*(1.0-zi)*0.5
1084 h(2)=x2*(1.0-zi)*0.5
1085 h(3)=(1.0-x1-x2)*(1.0-zi)*0.5
1086 h(4)=x1*(1.0+zi)*0.5
1087 h(5)=x2*(1.0+zi)*0.5
1088 h(6)=(1.0-x1-x2)*(1.0+zi)*0.5
1089 hr(1)= (1.0-zi)*0.5
1090 hr(2)= 0.0
1091 hr(3)=-(1.0-zi)*0.5
1092 hr(4)= (1.0+zi)*0.5
1093 hr(5)= 0.0
1094 hr(6)=-(1.0+zi)*0.5
1095 hs(1)= 0.0
1096 hs(2)= (1.0-zi)*0.5
1097 hs(3)= -(1.0-zi)*0.5
1098 hs(4)= 0.0
1099 hs(5)= (1.0+zi)*0.5
1100 hs(6)= -(1.0+zi)*0.5
1101 ht(1)= -0.5*x1
1102 ht(2)= -0.5*x2
1103 ht(3)= -0.5*(1.0-x1-x2)
1104 ht(4)= 0.5*x1
1105 ht(5)= 0.5*x2
1106 ht(6)= 0.5*(1.0-x1-x2)
1107 !
1108 xj11=hr(1)*xx(1)+hr(2)*xx(2)+hr(3)*xx(3)+hr(4)*xx(4) &
1109 +hr(5)*xx(5)+hr(6)*xx(6)
1110 xj21=hs(1)*xx(1)+hs(2)*xx(2)+hs(3)*xx(3)+hs(4)*xx(4) &
1111 +hs(5)*xx(5)+hs(6)*xx(6)
1112 xj31=ht(1)*xx(1)+ht(2)*xx(2)+ht(3)*xx(3)+ht(4)*xx(4) &
1113 +ht(5)*xx(5)+ht(6)*xx(6)
1114 xj12=hr(1)*yy(1)+hr(2)*yy(2)+hr(3)*yy(3)+hr(4)*yy(4) &
1115 +hr(5)*yy(5)+hr(6)*yy(6)
1116 xj22=hs(1)*yy(1)+hs(2)*yy(2)+hs(3)*yy(3)+hs(4)*yy(4) &
1117 +hs(5)*yy(5)+hs(6)*yy(6)
1118 xj32=ht(1)*yy(1)+ht(2)*yy(2)+ht(3)*yy(3)+ht(4)*yy(4) &
1119 +ht(5)*yy(5)+ht(6)*yy(6)
1120 xj13=hr(1)*zz(1)+hr(2)*zz(2)+hr(3)*zz(3)+hr(4)*zz(4) &
1121 +hr(5)*zz(5)+hr(6)*zz(6)
1122 xj23=hs(1)*zz(1)+hs(2)*zz(2)+hs(3)*zz(3)+hs(4)*zz(4) &
1123 +hs(5)*zz(5)+hs(6)*zz(6)
1124 xj33=ht(1)*zz(1)+ht(2)*zz(2)+ht(3)*zz(3)+ht(4)*zz(4) &
1125 +ht(5)*zz(5)+ht(6)*zz(6)
1126 !*DETERMINANT OF JACOBIAN
1127 det=xj11*xj22*xj33 &
1128 +xj12*xj23*xj31 &
1129 +xj13*xj21*xj32 &
1130 -xj13*xj22*xj31 &
1131 -xj12*xj21*xj33 &
1132 -xj11*xj23*xj32
1133 wg=wgt1(l12)*wgt(lz)*det
1134 do i=1,nn
1135 vect(i)=vect(i)-val*h(i)*wg
1136 enddo
1137 !
1138 enddo
1139 enddo
1140 endif
1141 !
1142 return
1143
1144 end subroutine heat_dflux_351
1145 !----------------------------------------------------------------------*
1146 subroutine heat_dflux_352(NN,XX,YY,ZZ,LTYPE,val,VECT)
1147 !----------------------------------------------------------------------*
1148 !**
1149 !** SET 352 DFLUX
1150 !**
1151 ! BF LTYPE=0 :BODY FLUX
1152 ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
1153 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
1154 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
1155 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
1156 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
1157 use hecmw
1158 implicit real(kind=kreal)(a-h,o-z)
1159 ! I/F VARIABLES
1160 integer(kind=kint) :: NN, LTYPE
1161 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
1162 ! LOCAL VARIABLES
1163 real(kind=kreal) :: h(15), hr(15), hs(15), ht(15), pl(15)
1164 real(kind=kreal) :: xg(3), wgt(3)
1165 integer(kind=kint) :: NOD(8)
1166 !*************************
1167 ! GAUSS INTEGRATION POINT
1168 !*************************
1169 xg(1) = -0.7745966692
1170 xg(2) = 0.0
1171 xg(3) = 0.7745966692
1172 wgt(1) = 0.5555555555
1173 wgt(2) = 0.8888888888
1174 wgt(3) = 0.5555555555
1175 !
1176 ! SELCTION OF LOAD TYPE
1177 !
1178 ivol=0
1179 isuf=0
1180 if( ltype.EQ.0 ) then
1181 ivol=1
1182 else
1183 isuf=1
1184 if( ltype.EQ.1 ) then
1185 nod(1) = 1
1186 nod(2) = 2
1187 nod(3) = 3
1188 nod(4) = 7
1189 nod(5) = 8
1190 nod(6) = 9
1191 isuf=2
1192 else if( ltype.EQ.2 ) then
1193 nod(1) = 6
1194 nod(2) = 5
1195 nod(3) = 4
1196 nod(4) = 11
1197 nod(5) = 10
1198 nod(6) = 12
1199 isuf=2
1200 else if( ltype.EQ.3 ) then
1201 nod(1) = 4
1202 nod(2) = 5
1203 nod(3) = 2
1204 nod(4) = 1
1205 nod(5) = 10
1206 nod(6) = 14
1207 nod(7) = 7
1208 nod(8) = 13
1209 else if( ltype.EQ.4 ) then
1210 nod(1) = 5
1211 nod(2) = 6
1212 nod(3) = 3
1213 nod(4) = 2
1214 nod(5) = 11
1215 nod(6) = 15
1216 nod(7) = 8
1217 nod(8) = 14
1218 else if( ltype.EQ.5 ) then
1219 nod(1) = 6
1220 nod(2) = 4
1221 nod(3) = 1
1222 nod(4) = 3
1223 nod(5) = 12
1224 nod(6) = 13
1225 nod(7) = 9
1226 nod(8) = 15
1227 endif
1228 endif
1229 do i=1,nn
1230 vect(i)=0.0
1231 enddo
1232 !
1233 !** SURFACE LOAD
1234 !
1235 if( isuf.EQ.1 ) then
1236 do ig2=1,3
1237 si=xg(ig2)
1238 do ig1=1,3
1239 ri=xg(ig1)
1240 rp=1.0+ri
1241 sp=1.0+si
1242 rm=1.0-ri
1243 sm=1.0-si
1244 h(1)=0.25*rm*sm*(-1.0-ri-si)
1245 h(2)=0.25*rp*sm*(-1.0+ri-si)
1246 h(3)=0.25*rp*sp*(-1.0+ri+si)
1247 h(4)=0.25*rm*sp*(-1.0-ri+si)
1248 h(5)=0.5*(1.0-ri*ri)*(1.0-si)
1249 h(6)=0.5*(1.0-si*si)*(1.0+ri)
1250 h(7)=0.5*(1.0-ri*ri)*(1.0+si)
1251 h(8)=0.5*(1.0-si*si)*(1.0-ri)
1252 hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1253 hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1254 hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
1255 hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
1256 hr(5)=-ri*(1.0-si)
1257 hr(6)= .5*(1.0-si*si)
1258 hr(7)=-ri*(1.0+si)
1259 hr(8)=-.5*(1.0-si*si)
1260 hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1261 hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1262 hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
1263 hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
1264 hs(5)=-.5*(1.0-ri*ri)
1265 hs(6)=-si*(1.0+ri)
1266 hs(7)= .5*(1.0-ri*ri)
1267 hs(8)=-si*(1.0-ri)
1268 g1x=0.0
1269 g1y=0.0
1270 g1z=0.0
1271 g2x=0.0
1272 g2y=0.0
1273 g2z=0.0
1274 g3x=0.0
1275 g3y=0.0
1276 g3z=0.0
1277 do i=1,8
1278 g1x=g1x+hr(i)*xx(nod(i))
1279 g1y=g1y+hr(i)*yy(nod(i))
1280 g1z=g1z+hr(i)*zz(nod(i))
1281 g2x=g2x+hs(i)*xx(nod(i))
1282 g2y=g2y+hs(i)*yy(nod(i))
1283 g2z=g2z+hs(i)*zz(nod(i))
1284 g3x=g3x+ht(i)*xx(nod(i))
1285 g3y=g3y+ht(i)*yy(nod(i))
1286 g3z=g3z+ht(i)*zz(nod(i))
1287 enddo
1288 g3x=g1y*g2z-g1z*g2y
1289 g3y=g1z*g2x-g1x*g2z
1290 g3z=g1x*g2y-g1y*g2x
1291 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1292 do i=1,8
1293 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1294 enddo
1295 enddo
1296 enddo
1297 elseif( isuf.EQ.2 ) then
1298 do ig2=1,3
1299 si=xg(ig2)
1300 do ig1=1,3
1301 ri=xg(ig1)
1302 rp=1.0+ri
1303 sp=1.0+si
1304 rm=1.0-ri
1305 sm=1.0-si
1306 h(1)=0.25*rm*sm*(-1.0-ri-si)
1307 h(2)=0.25*rp*sm*(-1.0+ri-si)
1308 h(3)=0.5*sp*si
1309 h(4)=0.5*(1.0-ri*ri)*(1.0-si)
1310 h(5)=0.5*(1.0-si*si)*(1.0+ri)
1311 h(6)=0.5*(1.0-si*si)*(1.0-ri)
1312 hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1313 hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1314 hr(3)= 0.0
1315 hr(4)=-ri*(1.0-si)
1316 hr(5)= 0.5*(1.0-si*si)
1317 hr(6)=-0.5*(1.0-si*si)
1318 hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1319 hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1320 hs(3)=0.5*(1.0+2.0*si)
1321 hs(4)=-0.5*(1.0-ri*ri)
1322 hs(5)=-si *(1.0+ri)
1323 hs(6)=-si*(1.0-ri)
1324 g1x=0.0
1325 g1y=0.0
1326 g1z=0.0
1327 g2x=0.0
1328 g2y=0.0
1329 g2z=0.0
1330 g3x=0.0
1331 g3y=0.0
1332 g3z=0.0
1333 do i=1, 6
1334 g1x=g1x+hr(i)*xx(nod(i))
1335 g1y=g1y+hr(i)*yy(nod(i))
1336 g1z=g1z+hr(i)*zz(nod(i))
1337 g2x=g2x+hs(i)*xx(nod(i))
1338 g2y=g2y+hs(i)*yy(nod(i))
1339 g2z=g2z+hs(i)*zz(nod(i))
1340 g3x=g3x+ht(i)*xx(nod(i))
1341 g3y=g3y+ht(i)*yy(nod(i))
1342 g3z=g3z+ht(i)*zz(nod(i))
1343 enddo
1344 g3x=g1y*g2z-g1z*g2y
1345 g3y=g1z*g2x-g1x*g2z
1346 g3z=g1x*g2y-g1y*g2x
1347 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1348 do i=1,6
1349 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1350 enddo
1351 enddo
1352 enddo
1353 endif
1354 !
1355 !** VOLUME LOAD
1356 !
1357 if( ivol.EQ.1 ) then
1358 do lx=1,3
1359 ri=xg(lx)
1360 do ly=1,3
1361 si=xg(ly)
1362 do lz=1,3
1363 ti=xg(lz)
1364 rp=1.0+ri
1365 sp=1.0+si
1366 tp=1.0+ti
1367 rm=1.0-ri
1368 sm=1.0-si
1369 tm=1.0-ti
1370 ! INTERPOLATION FUNCTION
1371 h(1)=-0.125*rm*sm*tm*(2.0+ri+si+ti)
1372 h(2)=-0.125*rp*sm*tm*(2.0-ri+si+ti)
1373 h(3)=-0.25*sp*tm*(1.0-si+ti)
1374 h(4)=-0.125*rm*sm*tp*(2.0+ri+si-ti)
1375 h(5)=-0.125*rp*sm*tp*(2.0-ri+si-ti)
1376 h(6)=0.25*sp*tp*(si+ti-1.0)
1377 h(7)=0.25*(1.0-ri**2)*sm*tm
1378 h(8)=0.25*rp*(1.0-si**2)*tm
1379 h(9)=0.25*rm*(1.0-si**2)*tm
1380 h(10)=0.25*(1.0-ri**2)*sm*tp
1381 h(11)=0.25*rp*(1.0-si**2)*tp
1382 h(12)=0.25*rm*(1.0-si**2)*tp
1383 h(13)=0.25*rm*sm*(1.0-ti**2)
1384 h(14)=0.25*rp*sm*(1.0-ti**2)
1385 h(15)=0.5*sp*(1.0-ti**2)
1386 ! DERIVATIVE OF INTERPOLATION FUNCTION
1387 ! FOR R-COORDINATE
1388 hr(1)=-0.125*rm*sm*tm+0.125*sm*tm*(2.0+ri+si+ti)
1389 hr(2)=+0.125*rp*sm*tm-0.125*sm*tm*(2.0-ri+si+ti)
1390 hr(3)=0.0
1391 hr(4)=-0.125*rm*sm*tp+0.125*sm*tp*(2.0+ri+si-ti)
1392 hr(5)=+0.125*rp*sm*tp-0.125*sm*tp*(2.0-ri+si-ti)
1393 hr(6)=0.0
1394 hr(7)=-0.50*ri*sm*tm
1395 hr(8)=+0.25*(1.0-si**2)*tm
1396 hr(9)=-0.25*(1.0-si**2)*tm
1397 hr(10)=-0.50*ri*sm*tp
1398 hr(11)=+0.25*(1.0-si**2)*tp
1399 hr(12)=-0.25*(1.0-si**2)*tp
1400 hr(13)=-0.25*sm*(1.0-ti**2)
1401 hr(14)=+0.25*sm*(1.0-ti**2)
1402 hr(15)=0.0
1403 ! FOR S-COORDINATE
1404 hs(1)=-0.125*rm*sm*tm+0.125*rm*tm*(2.0+ri+si+ti)
1405 hs(2)=-0.125*rp*sm*tm+0.125*rp*tm*(2.0-ri+si+ti)
1406 hs(3)=+0.25*tm*(2.0*si-ti)
1407 hs(4)=-0.125*rm*sm*tp+0.125*rm*tp*(2.0+ri+si-ti)
1408 hs(5)=-0.125*rp*sm*tp+0.125*rp*tp*(2.0-ri+si-ti)
1409 hs(6)=+0.25*tp*(2.0*si+ti)
1410 hs(7)=-0.25*(1.0-ri**2)*tm
1411 hs(8)=-0.50*rp*si*tm
1412 hs(9)=-0.50*rm*si*tm
1413 hs(10)=-0.25*(1.0-ri**2)*tp
1414 hs(11)=-0.50*rp*si*tp
1415 hs(12)=-0.50*rm*si*tp
1416 hs(13)=-0.25*rm*(1.0-ti**2)
1417 hs(14)=-0.25*rp*(1.0-ti**2)
1418 hs(15)=+0.50*(1.0-ti**2)
1419 ! FOR T-COORDINATE
1420 ht(1)=-0.125*rm*sm*tm+0.125*rm*sm*(2.0+ri+si+ti)
1421 ht(2)=-0.125*rp*sm*tm+0.125*rp*sm*(2.0-ri+si+ti)
1422 ht(3)=+0.25*sp*(2.0*ti-si)
1423 ht(4)=+0.125*rm*sm*tp-0.125*rm*sm*(2.0+ri+si-ti)
1424 ht(5)=+0.125*rp*sm*tp-0.125*rp*sm*(2.0-ri+si-ti)
1425 ht(6)=+0.25*sp*(si+2.0*ti)
1426 ht(7)=-0.25*(1.0-ri**2)*sm
1427 ht(8)=-0.25*rp*(1.0-si**2)
1428 ht(9)=-0.25*rm*(1.0-si**2)
1429 ht(10)=0.25*(1.0-ri**2)*sm
1430 ht(11)=0.25*rp*(1.0-si**2)
1431 ht(12)=0.25*rm*(1.0-si**2)
1432 ht(13)=-0.5*rm*sm*ti
1433 ht(14)=-0.5*rp*sm*ti
1434 ht(15)=-sp*ti
1435 ! JACOBI MATRIX
1436 xj11=0.0
1437 xj21=0.0
1438 xj31=0.0
1439 xj12=0.0
1440 xj22=0.0
1441 xj32=0.0
1442 xj13=0.0
1443 xj23=0.0
1444 xj33=0.0
1445 do i=1,nn
1446 xj11=xj11+hr(i)*xx(i)
1447 xj21=xj21+hs(i)*xx(i)
1448 xj31=xj31+ht(i)*xx(i)
1449 xj12=xj12+hr(i)*yy(i)
1450 xj22=xj22+hs(i)*yy(i)
1451 xj32=xj32+ht(i)*yy(i)
1452 xj13=xj13+hr(i)*zz(i)
1453 xj23=xj23+hs(i)*zz(i)
1454 xj33=xj33+ht(i)*zz(i)
1455 enddo
1456 !DETERMINANT OF JACOBIAN
1457 det=xj11*xj22*xj33 &
1458 +xj12*xj23*xj31 &
1459 +xj13*xj21*xj32 &
1460 -xj13*xj22*xj31 &
1461 -xj12*xj21*xj33 &
1462 -xj11*xj23*xj32
1463 wg=det*wgt(lx)*wgt(ly)*wgt(lz)
1464 do i=1,nn
1465 vect(i)=vect(i) - val*wg*h(i)
1466 enddo
1467 enddo
1468 enddo
1469 enddo
1470 endif
1471 return
1472
1473 end subroutine heat_dflux_352
1474 !----------------------------------------------------------------------*
1475 subroutine heat_dflux_361(NN,XX,YY,ZZ,LTYPE,val,VECT)
1476 !----------------------------------------------------------------------*
1477 !**
1478 !** SET 361 DFLUX
1479 !**
1480 ! BF LTYPE=0 :BODY FLUX
1481 ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
1482 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
1483 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
1484 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
1485 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
1486 ! S6 LTYPE=6 :FLUX IN NORMAL-DIRECTION FOR FACE-6
1487 use hecmw
1488 implicit none
1489 ! I/F VARIABLES
1490 integer(kind=kint) :: NN, LTYPE
1491 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
1492 ! LOCAL VARIABLES
1493 real(kind=kreal) :: h(8), hr(8), hs(8), ht(8), pl(8)
1494 real(kind=kreal) :: xg(2), wgt(2)
1495 real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
1496 real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33, det, wg
1497 integer(kind=kint) :: IVOL, ISUF
1498 integer(kind=kint) :: NOD(4)
1499 integer(kind=kint) :: IG1, IG2, LX, LY, LZ, I
1500 real(kind=kreal) :: vx, vy, vz
1501 real(kind=kreal) :: g1x, g1y, g1z
1502 real(kind=kreal) :: g2x, g2y, g2z
1503 real(kind=kreal) :: g3x, g3y, g3z
1504 real(kind=kreal) :: xsum
1505 !*************************
1506 ! GAUSS INTEGRATION POINT
1507 !*************************
1508 data xg/-0.5773502691896,0.5773502691896/
1509 data wgt/1.0,1.0/
1510 !
1511 ! SELCTION OF LOAD TYPE
1512 !
1513 ivol=0
1514 isuf=0
1515 if( ltype.EQ.0 ) then
1516 ivol=1
1517 else
1518 isuf=1
1519 if( ltype.EQ.1 ) then
1520 nod(1) = 1
1521 nod(2) = 2
1522 nod(3) = 3
1523 nod(4) = 4
1524 else if( ltype.EQ.2 ) then
1525 nod(1) = 8
1526 nod(2) = 7
1527 nod(3) = 6
1528 nod(4) = 5
1529 else if( ltype.EQ.3 ) then
1530 nod(1) = 5
1531 nod(2) = 6
1532 nod(3) = 2
1533 nod(4) = 1
1534 else if( ltype.EQ.4 ) then
1535 nod(1) = 6
1536 nod(2) = 7
1537 nod(3) = 3
1538 nod(4) = 2
1539 else if( ltype.EQ.5 ) then
1540 nod(1) = 7
1541 nod(2) = 8
1542 nod(3) = 4
1543 nod(4) = 3
1544 else if( ltype.EQ.6 ) then
1545 nod(1) = 8
1546 nod(2) = 5
1547 nod(3) = 1
1548 nod(4) = 4
1549 endif
1550 endif
1551 ! CLEAR VECT
1552 do i=1,nn
1553 vect(i)=0.0
1554 enddo
1555 !
1556 !** SURFACE LOAD
1557 !
1558 if( isuf.EQ.1 ) then
1559 ! INTEGRATION OVER SURFACE
1560 do ig2=1,2
1561 si=xg(ig2)
1562 do ig1=1,2
1563 ri=xg(ig1)
1564 h(1)=0.25*(1.0-ri)*(1.0-si)
1565 h(2)=0.25*(1.0+ri)*(1.0-si)
1566 h(3)=0.25*(1.0+ri)*(1.0+si)
1567 h(4)=0.25*(1.0-ri)*(1.0+si)
1568 hr(1)=-.25*(1.0-si)
1569 hr(2)= .25*(1.0-si)
1570 hr(3)= .25*(1.0+si)
1571 hr(4)=-.25*(1.0+si)
1572 hs(1)=-.25*(1.0-ri)
1573 hs(2)=-.25*(1.0+ri)
1574 hs(3)= .25*(1.0+ri)
1575 hs(4)= .25*(1.0-ri)
1576 g1x=0.0
1577 g1y=0.0
1578 g1z=0.0
1579 g2x=0.0
1580 g2y=0.0
1581 g2z=0.0
1582 do i=1,4
1583 g1x=g1x+hr(i)*xx(nod(i))
1584 g1y=g1y+hr(i)*yy(nod(i))
1585 g1z=g1z+hr(i)*zz(nod(i))
1586 g2x=g2x+hs(i)*xx(nod(i))
1587 g2y=g2y+hs(i)*yy(nod(i))
1588 g2z=g2z+hs(i)*zz(nod(i))
1589 enddo
1590 g3x=g1y*g2z-g1z*g2y
1591 g3y=g1z*g2x-g1x*g2z
1592 g3z=g1x*g2y-g1y*g2x
1593 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1594 do i=1,4
1595 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1596 enddo
1597 enddo
1598 enddo
1599 endif
1600 !
1601 !** VOLUME LOAD
1602 !
1603 if( ivol.EQ.1 ) then
1604 do i=1,nn
1605 pl(i)=0.0
1606 enddo
1607 ! LOOP FOR INTEGRATION POINTS
1608 do lx=1,2
1609 ri=xg(lx)
1610 do ly=1,2
1611 si=xg(ly)
1612 do lz=1,2
1613 ti=xg(lz)
1614 rp=1.0+ri
1615 sp=1.0+si
1616 tp=1.0+ti
1617 rm=1.0-ri
1618 sm=1.0-si
1619 tm=1.0-ti
1620 ! INTERPOLATION FUNCTION
1621 h(1)=0.125*rm*sm*tm
1622 h(2)=0.125*rp*sm*tm
1623 h(3)=0.125*rp*sp*tm
1624 h(4)=0.125*rm*sp*tm
1625 h(5)=0.125*rm*sm*tp
1626 h(6)=0.125*rp*sm*tp
1627 h(7)=0.125*rp*sp*tp
1628 h(8)=0.125*rm*sp*tp
1629 ! DERIVATIVE OF INTERPOLATION FUNCTION
1630 ! FOR R-COORDINATE
1631 hr(1)=-.125*sm*tm
1632 hr(2)= .125*sm*tm
1633 hr(3)= .125*sp*tm
1634 hr(4)=-.125*sp*tm
1635 hr(5)=-.125*sm*tp
1636 hr(6)= .125*sm*tp
1637 hr(7)= .125*sp*tp
1638 hr(8)=-.125*sp*tp
1639 ! FOR S-COORDINATE
1640 hs(1)=-.125*rm*tm
1641 hs(2)=-.125*rp*tm
1642 hs(3)= .125*rp*tm
1643 hs(4)= .125*rm*tm
1644 hs(5)=-.125*rm*tp
1645 hs(6)=-.125*rp*tp
1646 hs(7)= .125*rp*tp
1647 hs(8)= .125*rm*tp
1648 ! FOR T-COORDINATE
1649 ht(1)=-.125*rm*sm
1650 ht(2)=-.125*rp*sm
1651 ht(3)=-.125*rp*sp
1652 ht(4)=-.125*rm*sp
1653 ht(5)= .125*rm*sm
1654 ht(6)= .125*rp*sm
1655 ht(7)= .125*rp*sp
1656 ht(8)= .125*rm*sp
1657 ! JACOBI MATRIX
1658 xj11=0.0
1659 xj21=0.0
1660 xj31=0.0
1661 xj12=0.0
1662 xj22=0.0
1663 xj32=0.0
1664 xj13=0.0
1665 xj23=0.0
1666 xj33=0.0
1667 do i=1,nn
1668 xj11=xj11+hr(i)*xx(i)
1669 xj21=xj21+hs(i)*xx(i)
1670 xj31=xj31+ht(i)*xx(i)
1671 xj12=xj12+hr(i)*yy(i)
1672 xj22=xj22+hs(i)*yy(i)
1673 xj32=xj32+ht(i)*yy(i)
1674 xj13=xj13+hr(i)*zz(i)
1675 xj23=xj23+hs(i)*zz(i)
1676 xj33=xj33+ht(i)*zz(i)
1677 enddo
1678 !DETERMINANT OF JACOBIAN
1679 det=xj11*xj22*xj33 &
1680 +xj12*xj23*xj31 &
1681 +xj13*xj21*xj32 &
1682 -xj13*xj22*xj31 &
1683 -xj12*xj21*xj33 &
1684 -xj11*xj23*xj32
1685 wg=wgt(lx)*wgt(ly)*wgt(lz)*det
1686 do i=1,nn
1687 vect(i)=vect(i)-h(i)*wg*val
1688 enddo
1689 enddo
1690 enddo
1691 enddo
1692 endif
1693 !*
1694 return
1695
1696 end subroutine heat_dflux_361
1697 !----------------------------------------------------------------------*
1698 subroutine heat_dflux_362(NN,XX,YY,ZZ,LTYPE,val,VECT)
1699 !----------------------------------------------------------------------*
1700 !**
1701 !** SET 362 DFLUX
1702 !**
1703 ! BF LTYPE=0 :BODY FLUX
1704 ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1
1705 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2
1706 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3
1707 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4
1708 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5
1709 ! S6 LTYPE=6 :FLUX IN NORMAL-DIRECTION FOR FACE-6
1710 use hecmw
1711 implicit none
1712 ! I/F VARIABLES
1713 integer(kind=kint) :: NN, LTYPE
1714 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
1715 ! LOCAL VARIABLES
1716 real(kind=kreal) :: h(20), hr(20), hs(20), ht(20)
1717 real(kind=kreal) :: xg(3), wgt(3)
1718 real(kind=kreal) :: ri, si, ti, rp, sp, tp, rm, sm, tm
1719 real(kind=kreal) :: xj11, xj21, xj31, xj12, xj22, xj32, xj13, xj23, xj33
1720 real(kind=kreal) :: det, wg
1721 integer(kind=kint) :: IVOL, ISUF
1722 integer(kind=kint) :: NOD(8)
1723 integer(kind=kint) :: IG1, IG2, LX, LY, LZ, I
1724 real(kind=kreal) :: vx, vy, vz
1725 real(kind=kreal) :: g1x, g1y, g1z
1726 real(kind=kreal) :: g2x, g2y, g2z
1727 real(kind=kreal) :: g3x, g3y, g3z
1728 real(kind=kreal) :: xsum
1729 !*************************
1730 ! GAUSS INTEGRATION POINT
1731 !*************************
1732 xg(1) = -0.7745966692
1733 xg(2) = 0.0
1734 xg(3) = 0.7745966692
1735 wgt(1) = 0.5555555555
1736 wgt(2) = 0.8888888888
1737 wgt(3) = 0.5555555555
1738 !
1739 ! SELCTION OF LOAD TYPE
1740 !
1741 ivol=0
1742 isuf=0
1743 if( ltype.EQ.0 ) then
1744 ivol=1
1745 else
1746 isuf=1
1747 if( ltype.EQ.1 ) then
1748 nod(1) = 1
1749 nod(2) = 2
1750 nod(3) = 3
1751 nod(4) = 4
1752 nod(5) = 9
1753 nod(6) = 10
1754 nod(7) = 11
1755 nod(8) = 12
1756 else if( ltype.EQ.2 ) then
1757 nod(1) = 8
1758 nod(2) = 7
1759 nod(3) = 6
1760 nod(4) = 5
1761 nod(5) = 15
1762 nod(6) = 14
1763 nod(7) = 13
1764 nod(8) = 16
1765 else if( ltype.EQ.3 ) then
1766 nod(1) = 5
1767 nod(2) = 6
1768 nod(3) = 2
1769 nod(4) = 1
1770 nod(5) = 13
1771 nod(6) = 18
1772 nod(7) = 9
1773 nod(8) = 17
1774 else if( ltype.EQ.4 ) then
1775 nod(1) = 6
1776 nod(2) = 7
1777 nod(3) = 3
1778 nod(4) = 2
1779 nod(5) = 14
1780 nod(6) = 19
1781 nod(7) = 10
1782 nod(8) = 18
1783 else if( ltype.EQ.5 ) then
1784 nod(1) = 7
1785 nod(2) = 8
1786 nod(3) = 4
1787 nod(4) = 3
1788 nod(5) = 15
1789 nod(6) = 20
1790 nod(7) = 11
1791 nod(8) = 19
1792 else if( ltype.EQ.6 ) then
1793 nod(1) = 8
1794 nod(2) = 5
1795 nod(3) = 1
1796 nod(4) = 4
1797 nod(5) = 16
1798 nod(6) = 17
1799 nod(7) = 12
1800 nod(8) = 20
1801 endif
1802 endif
1803 ! CLEAR VECT
1804 do i=1,nn
1805 vect(i)=0.0
1806 enddo
1807 !
1808 !** SURFACE LOAD
1809 !
1810 if( isuf.EQ.1 ) then
1811 ! INTEGRATION OVER SURFACE
1812 do ig2=1,3
1813 si=xg(ig2)
1814 do ig1=1,3
1815 ri=xg(ig1)
1816 rp=1.0+ri
1817 sp=1.0+si
1818 rm=1.0-ri
1819 sm=1.0-si
1820 h(1)=0.25*rm*sm*(-1.0-ri-si)
1821 h(2)=0.25*rp*sm*(-1.0+ri-si)
1822 h(3)=0.25*rp*sp*(-1.0+ri+si)
1823 h(4)=0.25*rm*sp*(-1.0-ri+si)
1824 h(5)=0.5*(1.0-ri*ri)*(1.0-si)
1825 h(6)=0.5*(1.0-si*si)*(1.0+ri)
1826 h(7)=0.5*(1.0-ri*ri)*(1.0+si)
1827 h(8)=0.5*(1.0-si*si)*(1.0-ri)
1828 hr(1)=-.25*sm*(-1.0-ri-si)-0.25*rm*sm
1829 hr(2)= .25*sm*(-1.0+ri-si)+0.25*rp*sm
1830 hr(3)= .25*sp*(-1.0+ri+si)+0.25*rp*sp
1831 hr(4)=-.25*sp*(-1.0-ri+si)-0.25*rm*sp
1832 hr(5)=-ri*(1.0-si)
1833 hr(6)= .5*(1.0-si*si)
1834 hr(7)=-ri*(1.0+si)
1835 hr(8)=-.5*(1.0-si*si)
1836 hs(1)=-.25*rm*(-1.0-ri-si)-0.25*rm*sm
1837 hs(2)=-.25*rp*(-1.0+ri-si)-0.25*rp*sm
1838 hs(3)= .25*rp*(-1.0+ri+si)+0.25*rp*sp
1839 hs(4)= .25*rm*(-1.0-ri+si)+0.25*rm*sp
1840 hs(5)=-.5*(1.0-ri*ri)
1841 hs(6)=-si*(1.0+ri)
1842 hs(7)= .5*(1.0-ri*ri)
1843 hs(8)=-si*(1.0-ri)
1844 g1x=0.0
1845 g1y=0.0
1846 g1z=0.0
1847 g2x=0.0
1848 g2y=0.0
1849 g2z=0.0
1850 do i=1,8
1851 g1x=g1x+hr(i)*xx(nod(i))
1852 g1y=g1y+hr(i)*yy(nod(i))
1853 g1z=g1z+hr(i)*zz(nod(i))
1854 g2x=g2x+hs(i)*xx(nod(i))
1855 g2y=g2y+hs(i)*yy(nod(i))
1856 g2z=g2z+hs(i)*zz(nod(i))
1857 enddo
1858 g3x=g1y*g2z-g1z*g2y
1859 g3y=g1z*g2x-g1x*g2z
1860 g3z=g1x*g2y-g1y*g2x
1861 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1862 do i=1,8
1863 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1864 enddo
1865 enddo
1866 enddo
1867 endif
1868 !
1869 !** VOLUME LOAD
1870 !
1871 if( ivol.EQ.1 ) then
1872 ! LOOP FOR INTEGRATION POINTS
1873 do lx=1,3
1874 ri=xg(lx)
1875 do ly=1,3
1876 si=xg(ly)
1877 do lz=1,3
1878 ti=xg(lz)
1879 rp=1.0+ri
1880 sp=1.0+si
1881 tp=1.0+ti
1882 rm=1.0-ri
1883 sm=1.0-si
1884 tm=1.0-ti
1885 ! INTERPOLATION FUNCTION
1886 h(1)=-0.125*rm*sm*tm*(2.0+ri+si+ti)
1887 h(2)=-0.125*rp*sm*tm*(2.0-ri+si+ti)
1888 h(3)=-0.125*rp*sp*tm*(2.0-ri-si+ti)
1889 h(4)=-0.125*rm*sp*tm*(2.0+ri-si+ti)
1890 h(5)=-0.125*rm*sm*tp*(2.0+ri+si-ti)
1891 h(6)=-0.125*rp*sm*tp*(2.0-ri+si-ti)
1892 h(7)=-0.125*rp*sp*tp*(2.0-ri-si-ti)
1893 h(8)=-0.125*rm*sp*tp*(2.0+ri-si-ti)
1894 h(9)=0.25*(1.0-ri**2)*sm*tm
1895 h(10)=0.25*rp*(1.0-si**2)*tm
1896 h(11)=0.25*(1.0-ri**2)*sp*tm
1897 h(12)=0.25*rm*(1.0-si**2)*tm
1898 h(13)=0.25*(1.0-ri**2)*sm*tp
1899 h(14)=0.25*rp*(1.0-si**2)*tp
1900 h(15)=0.25*(1.0-ri**2)*sp*tp
1901 h(16)=0.25*rm*(1.0-si**2)*tp
1902 h(17)=0.25*rm*sm*(1.0-ti**2)
1903 h(18)=0.25*rp*sm*(1.0-ti**2)
1904 h(19)=0.25*rp*sp*(1.0-ti**2)
1905 h(20)=0.25*rm*sp*(1.0-ti**2)
1906 ! DERIVATIVE OF INTERPOLATION FUNCTION
1907 ! FOR R-COORDINATE
1908 hr(1)=-0.125*rm*sm*tm+0.125*sm*tm*(2.0+ri+si+ti)
1909 hr(2)=+0.125*rp*sm*tm-0.125*sm*tm*(2.0-ri+si+ti)
1910 hr(3)=+0.125*rp*sp*tm-0.125*sp*tm*(2.0-ri-si+ti)
1911 hr(4)=-0.125*rm*sp*tm+0.125*sp*tm*(2.0+ri-si+ti)
1912 hr(5)=-0.125*rm*sm*tp+0.125*sm*tp*(2.0+ri+si-ti)
1913 hr(6)=+0.125*rp*sm*tp-0.125*sm*tp*(2.0-ri+si-ti)
1914 hr(7)=+0.125*rp*sp*tp-0.125*sp*tp*(2.0-ri-si-ti)
1915 hr(8)=-0.125*rm*sp*tp+0.125*sp*tp*(2.0+ri-si-ti)
1916 hr(9 )=-0.50*ri*sm*tm
1917 hr(10)=+0.25*(1.0-si**2)*tm
1918 hr(11)=-0.50*ri*sp*tm
1919 hr(12)=-0.25*(1.0-si**2)*tm
1920 hr(13)=-0.50*ri*sm*tp
1921 hr(14)=+0.25*(1.0-si**2)*tp
1922 hr(15)=-0.50*ri*sp*tp
1923 hr(16)=-0.25*(1.0-si**2)*tp
1924 hr(17)=-0.25*sm*(1.0-ti**2)
1925 hr(18)=+0.25*sm*(1.0-ti**2)
1926 hr(19)=+0.25*sp*(1.0-ti**2)
1927 hr(20)=-0.25*sp*(1.0-ti**2)
1928 ! FOR S-COORDINATE
1929 hs(1)=-0.125*rm*sm*tm+0.125*rm*tm*(2.0+ri+si+ti)
1930 hs(2)=-0.125*rp*sm*tm+0.125*rp*tm*(2.0-ri+si+ti)
1931 hs(3)=+0.125*rp*sp*tm-0.125*rp*tm*(2.0-ri-si+ti)
1932 hs(4)=+0.125*rm*sp*tm-0.125*rm*tm*(2.0+ri-si+ti)
1933 hs(5)=-0.125*rm*sm*tp+0.125*rm*tp*(2.0+ri+si-ti)
1934 hs(6)=-0.125*rp*sm*tp+0.125*rp*tp*(2.0-ri+si-ti)
1935 hs(7)=+0.125*rp*sp*tp-0.125*rp*tp*(2.0-ri-si-ti)
1936 hs(8)=+0.125*rm*sp*tp-0.125*rm*tp*(2.0+ri-si-ti)
1937 hs(9)=-0.25*(1.0-ri**2)*tm
1938 hs(10)=-0.50*rp*si*tm
1939 hs(11)=+0.25*(1.0-ri**2)*tm
1940 hs(12)=-0.50*rm*si*tm
1941 hs(13)=-0.25*(1.0-ri**2)*tp
1942 hs(14)=-0.50*rp*si*tp
1943 hs(15)=+0.25*(1.0-ri**2)*tp
1944 hs(16)=-0.50*rm*si*tp
1945 hs(17)=-0.25*rm*(1.0-ti**2)
1946 hs(18)=-0.25*rp*(1.0-ti**2)
1947 hs(19)=+0.25*rp*(1.0-ti**2)
1948 hs(20)=+0.25*rm*(1.0-ti**2)
1949 ! FOR T-COORDINATE
1950 ht(1)=-0.125*rm*sm*tm+0.125*rm*sm*(2.0+ri+si+ti)
1951 ht(2)=-0.125*rp*sm*tm+0.125*rp*sm*(2.0-ri+si+ti)
1952 ht(3)=-0.125*rp*sp*tm+0.125*rp*sp*(2.0-ri-si+ti)
1953 ht(4)=-0.125*rm*sp*tm+0.125*rm*sp*(2.0+ri-si+ti)
1954 ht(5)=+0.125*rm*sm*tp-0.125*rm*sm*(2.0+ri+si-ti)
1955 ht(6)=+0.125*rp*sm*tp-0.125*rp*sm*(2.0-ri+si-ti)
1956 ht(7)=+0.125*rp*sp*tp-0.125*rp*sp*(2.0-ri-si-ti)
1957 ht(8)=+0.125*rm*sp*tp-0.125*rm*sp*(2.0+ri-si-ti)
1958 ht(9)=-0.25*(1.0-ri**2)*sm
1959 ht(10)=-0.25*rp*(1.0-si**2)
1960 ht(11)=-0.25*(1.0-ri**2)*sp
1961 ht(12)=-0.25*rm*(1.0-si**2)
1962 ht(13)=0.25*(1.0-ri**2)*sm
1963 ht(14)=0.25*rp*(1.0-si**2)
1964 ht(15)=0.25*(1.0-ri**2)*sp
1965 ht(16)=0.25*rm*(1.0-si**2)
1966 ht(17)=-0.5*rm*sm*ti
1967 ht(18)=-0.5*rp*sm*ti
1968 ht(19)=-0.5*rp*sp*ti
1969 ht(20)=-0.5*rm*sp*ti
1970 ! JACOBI MATRIX
1971 xj11=0.0
1972 xj21=0.0
1973 xj31=0.0
1974 xj12=0.0
1975 xj22=0.0
1976 xj32=0.0
1977 xj13=0.0
1978 xj23=0.0
1979 xj33=0.0
1980 do i=1,nn
1981 xj11=xj11+hr(i)*xx(i)
1982 xj21=xj21+hs(i)*xx(i)
1983 xj31=xj31+ht(i)*xx(i)
1984 xj12=xj12+hr(i)*yy(i)
1985 xj22=xj22+hs(i)*yy(i)
1986 xj32=xj32+ht(i)*yy(i)
1987 xj13=xj13+hr(i)*zz(i)
1988 xj23=xj23+hs(i)*zz(i)
1989 xj33=xj33+ht(i)*zz(i)
1990 enddo
1991 !DETERMINANT OF JACOBIAN
1992 det=xj11*xj22*xj33 &
1993 +xj12*xj23*xj31 &
1994 +xj13*xj21*xj32 &
1995 -xj13*xj22*xj31 &
1996 -xj12*xj21*xj33 &
1997 -xj11*xj23*xj32
1998 wg=wgt(lx)*wgt(ly)*wgt(lz)*det
1999 do i=1,nn
2000 vect(i)=vect(i)-h(i)*wg*val
2001 enddo
2002 enddo
2003 enddo
2004 enddo
2005 endif
2006 !*
2007 return
2008
2009 end subroutine heat_dflux_362
2010 !----------------------------------------------------------------------*
2011 subroutine heat_dflux_731(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
2012 !----------------------------------------------------------------------*
2013 !**
2014 !** SET S3 DFLUX
2015 !**
2016 ! BF LTYPE=0 :BODY FLUX
2017 ! S1 LTYPE=1 :SURFACE FLUX
2018 use hecmw
2019 implicit real(kind=kreal) (a - h, o - z)
2020 ! I/F VARIABLES
2021 integer(kind=kint) :: NN, LTYPE
2022 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
2023 !
2024 if( ltype.EQ.0 ) then
2025 thick = thick
2026 elseif( ltype.EQ.1 ) then
2027 thick = 1.0d0
2028 else
2029 thick = 0.0d0
2030 endif
2031
2032 do i=1,nn
2033 vect(i)=0.0
2034 enddo
2035 !
2036 v1x=xx(2)-xx(1)
2037 v1y=yy(2)-yy(1)
2038 v1z=zz(2)-zz(1)
2039 v2x=xx(3)-xx(1)
2040 v2y=yy(3)-yy(1)
2041 v2z=zz(3)-zz(1)
2042 v3x= v1y*v2z-v1z*v2y
2043 v3y= v1z*v2x-v1x*v2z
2044 v3z= v1x*v2y-v1y*v2x
2045 aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
2046 vect(1)= -val*aa*thick/3.0
2047 vect(2)= -val*aa*thick/3.0
2048 vect(3)= -val*aa*thick/3.0
2049 !
2050 return
2051
2052 end subroutine heat_dflux_731
2053 !----------------------------------------------------------------------*
2054 subroutine heat_dflux_741(NN,XX,YY,ZZ,THICK,LTYPE,val,VECT)
2055 !----------------------------------------------------------------------*
2056 !**
2057 !** SET S4 DFLUX
2058 !**
2059 ! BF LTYPE=0 : BODY FLUX
2060 ! S1 LTYPE=1 : SURFACE FLUX
2061 use hecmw
2062 implicit real(kind=kreal) (a - h, o - z)
2063 ! I/F VARIABLES
2064 integer(kind=kint) :: NN, LTYPE
2065 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
2066 ! LOCAL VARIABLES
2067 real(kind=kreal) :: h(4), hr(4), hs(4), ht(4), pl(4)
2068 real(kind=kreal) :: xg(2), wgt(2)
2069 !*************************
2070 ! GAUSS INTEGRATION POINT
2071 !*************************
2072 data xg / -0.5773502691896,0.5773502691896 /
2073 data wgt / 1.0, 1.0 /
2074 !
2075 ! SELCTION OF LOAD TYPE
2076 !
2077 if ( ltype.EQ.0 ) then
2078 thick = 1.0d0 * thick
2079 elseif( ltype.EQ.1 ) then
2080 thick = 1.0d0
2081 else
2082 thick = 0.0d0
2083 endif
2084 !
2085 do i=1,nn
2086 vect(i)=0.0
2087 enddo
2088 !
2089 do ig2=1,2
2090 si=xg(ig2)
2091 do ig1=1,2
2092 ri=xg(ig1)
2093 !
2094 h(1)=0.25*(1.0-ri)*(1.0-si)
2095 h(2)=0.25*(1.0+ri)*(1.0-si)
2096 h(3)=0.25*(1.0+ri)*(1.0+si)
2097 h(4)=0.25*(1.0-ri)*(1.0+si)
2098 hr(1)=-.25*(1.0-si)
2099 hr(2)= .25*(1.0-si)
2100 hr(3)= .25*(1.0+si)
2101 hr(4)=-.25*(1.0+si)
2102 hs(1)=-.25*(1.0-ri)
2103 hs(2)=-.25*(1.0+ri)
2104 hs(3)= .25*(1.0+ri)
2105 hs(4)= .25*(1.0-ri)
2106 !
2107 g1x=0.0
2108 g1y=0.0
2109 g1z=0.0
2110 g2x=0.0
2111 g2y=0.0
2112 g2z=0.0
2113 do i=1,nn
2114 g1x=g1x+hr(i)*xx(i)
2115 g1y=g1y+hr(i)*yy(i)
2116 g1z=g1z+hr(i)*zz(i)
2117 g2x=g2x+hs(i)*xx(i)
2118 g2y=g2y+hs(i)*yy(i)
2119 g2z=g2z+hs(i)*zz(i)
2120 enddo
2121 !
2122 g3x=g1y*g2z-g1z*g2y
2123 g3y=g1z*g2x-g1x*g2z
2124 g3z=g1x*g2y-g1y*g2x
2125 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
2126 do i=1,nn
2127 vect(i)=vect(i)-xsum*val*wgt(ig1)*wgt(ig2)*h(i)*thick
2128 enddo
2129 !
2130 enddo
2131 enddo
2132 return
2133
2134 end subroutine heat_dflux_741
2135end module m_heat_lib_dflux
Definition: hecmw.f90:6
This module provides subroutines for calculating distributed heat flux for various elements.
subroutine heat_dflux_232(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_341(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_362(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_352(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_231(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_342(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_351(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_241(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_731(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_242(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_111(nn, xx, yy, zz, asect, ltype, val, vect)
subroutine heat_dflux_741(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_361(nn, xx, yy, zz, ltype, val, vect)