FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 !-------------------------------------------------------------------------------
8 contains
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
2135 end module m_heat_lib_dflux
Definition: hecmw.f90:6
This module provides subroutines for calculating distributed heat flux for various elements.
subroutine heat_dflux_242(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
subroutine heat_dflux_351(NN, XX, YY, ZZ, LTYPE, val, VECT)
subroutine heat_dflux_341(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_741(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
subroutine heat_dflux_111(NN, XX, YY, ZZ, ASECT, LTYPE, val, VECT)
subroutine heat_dflux_342(NN, XX, YY, ZZ, LTYPE, val, VECT)
subroutine heat_dflux_232(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
subroutine heat_dflux_731(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)
subroutine heat_dflux_362(NN, XX, YY, ZZ, LTYPE, val, VECT)
subroutine heat_dflux_361(NN, XX, YY, ZZ, LTYPE, val, VECT)
subroutine heat_dflux_241(NN, XX, YY, ZZ, THICK, LTYPE, val, VECT)