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