31 implicit real(kind=kreal)(a-h,o-z)
33 integer(kind=kint) NN, LTYPE
34 real(kind=kreal) xx(nn),yy(nn),zz(nn),val,vect(nn)
50 al = dsqrt( dx*dx + dy*dy + dz*dz )
53 vect(1)= -val*vv/2.0d0
54 vect(2)= -val*vv/2.0d0
67 implicit real(kind=kreal) (a - h, o - z)
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)
76 data xg/-0.5773502691896,0.5773502691896/
83 else if( ltype.GE.1 )
then
88 else if(ltype.EQ.2)
then
91 else if(ltype.EQ.3)
then
111 gx=gx+hr(i)*xx(nod(i))
112 gy=gy+hr(i)*yy(nod(i))
114 xsum=dsqrt(gx*gx+gy*gy)*thick
116 vect(nod(i))=vect(nod(i))-xsum*wgt(lx)*h(i)*val
136 aa=0.5*dsqrt( v3x*v3x + v3y*v3y + v3z*v3z )
154 implicit real(kind=kreal) (a - h, o - z)
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)
164 xg(1) = -0.7745966692
167 wgt(1) = 0.5555555555
168 wgt(2) = 0.8888888888
169 wgt(3) = 0.5555555555
173 if( ltype.EQ.0 )
then
175 else if( ltype.GE.1 )
then
181 else if(ltype.EQ.2)
then
185 else if(ltype.EQ.3)
then
199 h(1)=-0.5*(1.0-ri)*ri
200 h(2)= 0.5*(1.0+ri)*ri
208 gx=gx+hr(i)*xx(nod(i))
209 gy=gy+hr(i)*yy(nod(i))
211 xsum=dsqrt(gx*gx+gy*gy)
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
227 x1=0.5*(1.0-x2)*(xl1+1.0)
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)
269 det=xj11*xj22-xj21*xj12
270 wg=wgt(l1)*wgt(l2)*det*thick*(1.0-x2)*0.25
272 vect(i) = vect(i) - h(i)*wg*val
288 implicit real(kind=kreal) (a - h, o - z)
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)
297 data xg/-0.5773502691896,0.5773502691896/
302 if( ltype.EQ.0 )
then
304 else if( ltype.GE.1 )
then
309 else if(ltype.EQ.2)
then
312 else if(ltype.EQ.3)
then
315 else if(ltype.EQ.4)
then
335 gx=gx+hr(i)*xx(nod(i))
336 gy=gy+hr(i)*yy(nod(i))
338 xsum=dsqrt(gx*gx+gy*gy)*thick
340 vect(nod(i))=vect(nod(i))-xsum*wgt(lx)*h(i)*val
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)
376 det=xj11*xj22-xj21*xj12
377 wg=wgt(lx)*wgt(ly)*det*thick
379 vect(i)=vect(i)-wg*h(i)*val
395 implicit real(kind=kreal) (a - h, o - z)
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)
404 xg(1) = -0.7745966692
407 wgt(1) = 0.5555555555
408 wgt(2) = 0.8888888888
409 wgt(3) = 0.5555555555
413 if( ltype.EQ.0 )
then
415 else if( ltype.GE.1 )
then
421 else if(ltype.EQ.2)
then
425 else if(ltype.EQ.3)
then
429 else if(ltype.EQ.4)
then
443 h(1)=-0.5*(1.0-ri)*ri
444 h(2)= 0.5*(1.0+ri)*ri
452 gx=gx+hr(i)*xx(nod(i))
453 gy=gy+hr(i)*yy(nod(i))
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
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
484 hr(6)= 0.5*(1.0-si*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)
493 hs(7)= 0.5*(1.0-ri*ri)
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)
505 det=xj11*xj22-xj21*xj12
506 wg=wgt(lx)*wgt(ly)*det*thick
508 vect(i)=vect(i)-wg*h(i)*val
530 integer(kind=kint) :: NN, LTYPE
531 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), vect(nn)
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
549 data xg/-0.5773502691896,0.5773502691896/
555 if( ltype.EQ.0 )
then
563 else if(ltype.EQ.2)
then
567 else if(ltype.EQ.3)
then
571 else if(ltype.EQ.4)
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))
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
605 x2 =(1.0-x3)*(xl2+1.0)*0.5
608 x1=(1.0-x2-x3)*(xl1+1.0)*0.5
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)
649 vect(i) = vect(i) + val*h(i)*det*(1.0-x3)*(1.0-x2-x3)*0.125
674 integer(kind=kint) :: NN, LTYPE
675 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), vect(nn)
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
695 xg(1) = -0.7745966692
698 wgt(1) = 0.5555555555
699 wgt(2) = 0.8888888888
700 wgt(3) = 0.5555555555
704 if( ltype.EQ.0 )
then
715 else if(ltype.EQ.2)
then
722 else if(ltype.EQ.3)
then
729 else if(ltype.EQ.4)
then
750 x1=0.5*(1.0-x2)*(xl1+1.0)
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))
799 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
820 wg=wgt(l1)*wgt(l2)*det*(1.0-x2)*0.25
822 vect(nod(i))=vect(nod(i))-val*wg*h(i)
837 x2 =(1.0-x3)*(xl2+1.0)*0.5
840 x1=(1.0-x2-x3)*(xl1+1.0)*0.5
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)
926 wg=det*wgt(l1)*wgt(l2)*wgt(l3)*(1.-x3)*(1.0-x2-x3)*0.125
928 vect(i) = vect(i) - val*wg*h(i)
952 implicit real(kind=kreal)(a-h,o-z)
954 integer(kind=kint) :: NN, LTYPE
955 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
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)
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 /
972 if( ltype.EQ.0 )
then
976 if( ltype.EQ.1 )
then
981 else if( ltype.EQ.2 )
then
986 else if( ltype.EQ.3 )
then
991 else if( ltype.EQ.4 )
then
996 else if( ltype.EQ.5 )
then
1009 if( isuf.EQ.1 )
then
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)
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))
1047 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1049 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
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
1072 if( ivol.EQ.1 )
then
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
1097 hs(3)= -(1.0-zi)*0.5
1100 hs(6)= -(1.0+zi)*0.5
1103 ht(3)= -0.5*(1.0-x1-x2)
1106 ht(6)= 0.5*(1.0-x1-x2)
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)
1127 det=xj11*xj22*xj33 &
1133 wg=wgt1(l12)*wgt(lz)*det
1135 vect(i)=vect(i)-val*h(i)*wg
1158 implicit real(kind=kreal)(a-h,o-z)
1160 integer(kind=kint) :: NN, LTYPE
1161 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
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)
1169 xg(1) = -0.7745966692
1171 xg(3) = 0.7745966692
1172 wgt(1) = 0.5555555555
1173 wgt(2) = 0.8888888888
1174 wgt(3) = 0.5555555555
1180 if( ltype.EQ.0 )
then
1184 if( ltype.EQ.1 )
then
1192 else if( ltype.EQ.2 )
then
1200 else if( ltype.EQ.3 )
then
1209 else if( ltype.EQ.4 )
then
1218 else if( ltype.EQ.5 )
then
1235 if( isuf.EQ.1 )
then
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
1257 hr(6)= .5*(1.0-si*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)
1266 hs(7)= .5*(1.0-ri*ri)
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))
1291 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1293 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1297 elseif( isuf.EQ.2 )
then
1306 h(1)=0.25*rm*sm*(-1.0-ri-si)
1307 h(2)=0.25*rp*sm*(-1.0+ri-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
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)
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))
1347 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1349 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1357 if( ivol.EQ.1 )
then
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)
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)
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)
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)
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)
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
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)
1457 det=xj11*xj22*xj33 &
1463 wg=det*wgt(lx)*wgt(ly)*wgt(lz)
1465 vect(i)=vect(i) - val*wg*h(i)
1490 integer(kind=kint) :: NN, LTYPE
1491 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
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
1508 data xg/-0.5773502691896,0.5773502691896/
1515 if( ltype.EQ.0 )
then
1519 if( ltype.EQ.1 )
then
1524 else if( ltype.EQ.2 )
then
1529 else if( ltype.EQ.3 )
then
1534 else if( ltype.EQ.4 )
then
1539 else if( ltype.EQ.5 )
then
1544 else if( ltype.EQ.6 )
then
1558 if( isuf.EQ.1 )
then
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)
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))
1593 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1595 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1603 if( ivol.EQ.1 )
then
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)
1679 det=xj11*xj22*xj33 &
1685 wg=wgt(lx)*wgt(ly)*wgt(lz)*det
1687 vect(i)=vect(i)-h(i)*wg*val
1713 integer(kind=kint) :: NN, LTYPE
1714 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
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
1732 xg(1) = -0.7745966692
1734 xg(3) = 0.7745966692
1735 wgt(1) = 0.5555555555
1736 wgt(2) = 0.8888888888
1737 wgt(3) = 0.5555555555
1743 if( ltype.EQ.0 )
then
1747 if( ltype.EQ.1 )
then
1756 else if( ltype.EQ.2 )
then
1765 else if( ltype.EQ.3 )
then
1774 else if( ltype.EQ.4 )
then
1783 else if( ltype.EQ.5 )
then
1792 else if( ltype.EQ.6 )
then
1810 if( isuf.EQ.1 )
then
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
1833 hr(6)= .5*(1.0-si*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)
1842 hs(7)= .5*(1.0-ri*ri)
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))
1861 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
1863 vect(nod(i))=vect(nod(i))-xsum*val*wgt(ig1)*wgt(ig2)*h(i)
1871 if( ivol.EQ.1 )
then
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)
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)
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)
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
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)
1992 det=xj11*xj22*xj33 &
1998 wg=wgt(lx)*wgt(ly)*wgt(lz)*det
2000 vect(i)=vect(i)-h(i)*wg*val
2019 implicit real(kind=kreal) (a - h, o - z)
2021 integer(kind=kint) :: NN, LTYPE
2022 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
2024 if( ltype.EQ.0 )
then
2026 elseif( ltype.EQ.1 )
then
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
2062 implicit real(kind=kreal) (a - h, o - z)
2064 integer(kind=kint) :: NN, LTYPE
2065 real(kind=kreal) :: xx(nn), yy(nn), zz(nn), val, vect(nn)
2067 real(kind=kreal) :: h(4), hr(4), hs(4), ht(4), pl(4)
2068 real(kind=kreal) :: xg(2), wgt(2)
2072 data xg / -0.5773502691896,0.5773502691896 /
2073 data wgt / 1.0, 1.0 /
2077 if ( ltype.EQ.0 )
then
2078 thick = 1.0d0 * thick
2079 elseif( ltype.EQ.1 )
then
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)
2125 xsum=dsqrt(g3x**2+g3y**2+g3z**2)
2127 vect(i)=vect(i)-xsum*val*wgt(ig1)*wgt(ig2)*h(i)*thick
This module provides subroutines for calculating distributed heat flux for various elements.
subroutine heat_dflux_232(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_341(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_362(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_352(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_231(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_342(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_351(nn, xx, yy, zz, ltype, val, vect)
subroutine heat_dflux_241(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_731(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_242(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_111(nn, xx, yy, zz, asect, ltype, val, vect)
subroutine heat_dflux_741(nn, xx, yy, zz, thick, ltype, val, vect)
subroutine heat_dflux_361(nn, xx, yy, zz, ltype, val, vect)