49 . MAS ,IXS ,PM ,X ,MSS ,
50 . DETONATORS,GEO ,VEUL ,ALE_CONNECTIVITY,IPARG ,
51 . DTELEM ,SIGI ,NEL ,SKEW ,IGEO ,
52 . STIFN ,PARTSAV ,V ,IPARTS ,IPART ,
53 . SIGSP ,NSIGI ,MSNF ,MSSF ,IPM ,
54 . IUSER ,NSIGS ,VOLNOD ,BVOLNOD ,VNS ,
55 . BNS ,WMA ,PTSOL ,BUFMAT ,MCP ,
56 . MCPS ,TEMP ,NPF ,TF ,MSSA ,
57 . STRSGLOB ,STRAGLOB ,ORTHOGLOB ,FAIL_INI ,ILOADP,
58 . FACLOAD ,RNOISE ,PERTURB ,GLOB_THERM,MAT_PARAM)
72 use element_mod ,
only : nixs
76#include "implicit_f.inc"
89#include "vect01_c.inc"
93 INTEGER NEL,NSIGI, IUSER, NSIGS
94 INTEGER IXS(NIXS,*),IPARG(*),IPARTS(*),IPART(LIPART1,*),
95 . IPM(NPROPMI,*),PTSOL(*), NPF(*),IGEO(NPROPGI,*),
96 . STRSGLOB(*),STRAGLOB(*),ORTHOGLOB(*),FAIL_INI(*),PERTURB(NPERTURB)
98 . MAS(*), PM(NPROPM,*), X(*), GEO(NPROPG,*),
99 . VEUL(LVEUL,*), DTELEM(*),SIGI(NSIGS,*),SKEW(LSKEW,*),STIFN(*),
100 . PARTSAV(20,*), V(*), MSS(8,*),
101 . SIGSP(NSIGI, *),MSNF(*), MSSF(8,*), WMA(*),
102 . VOLNOD(*), BVOLNOD(*), VNS(8,*), BNS(8,*), BUFMAT(*),
103 . mcp(*),mcps(8,*),temp(*), tf(*), mssa(*),rnoise(nperturb,*)
104 TYPE(elbuf_struct_),
TARGET :: ELBUF_STR
105 INTEGER,
INTENT(IN) :: ILOADP(SIZLOADP,*)
106 my_real,
INTENT(IN) :: FACLOAD(LFACLOAD,*)
108 TYPE(t_ale_connectivity),
INTENT(INOUT) :: ALE_CONNECTIVITY
109 type (glob_therm_) ,
intent(in) :: glob_therm
110 TYPE (MATPARAM_STRUCT_) ,
DIMENSION(NUMMAT) ,
INTENT(INOUT) :: MAT_PARAM
114 INTEGER NF1, IBID, I, NLAY,IGTYP,NLYMAX,IS,NUVAR,IREP,NCC,JHBE,
115 . idef, ip, ipang, ipthk, ippos, ipmat,ig,im,mtn0,ipid1,
116 . nptr,npts,nptt,l_pla,l_sigb
117 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
118 . ix5(mvsiz), ix6(mvsiz), ix7(mvsiz), ix8(mvsiz)
119 INTEGER MAT(MVSIZ), PID(MVSIZ), NGL(MVSIZ) , MAT0(MVSIZ)
120 CHARACTER(LEN=NCHARTITLE)::TITR1
122 . BID, FV, STI, ZI,WI
124 . v8loc(51,mvsiz),volu(mvsiz),dtx(mvsiz),vzl(mvsiz),vzq(mvsiz),
125 . x1(mvsiz),x2(mvsiz),x3
126 . x7(mvsiz),x8(mvsiz),y1(mvsiz),y2(mvsiz),y3(mvsiz),y4(mvsiz),
127 . y5(mvsiz),y6(mvsiz),y7(mvsiz),y8(mvsiz),z1(mvsiz),z2(mvsiz),
128 . z3(mvsiz),z4(mvsiz),z5(mvsiz),z6
129 . rx(mvsiz) ,ry(mvsiz) ,rz(mvsiz) ,sx(mvsiz) ,
130 . sy(mvsiz) ,sz(mvsiz) ,tx(mvsiz) ,ty(mvsiz) ,
131 . tz(mvsiz) ,e1x(mvsiz),e1y(mvsiz),e1z(mvsiz),e2x(mvsiz),
132 . e2y(mvsiz),e2z(mvsiz),e3x(mvsiz),e3y(mvsiz),e3z(mvsiz),
133 . f1x(mvsiz) ,f1y(mvsiz) ,f1z
134 . f2x(mvsiz) ,f2y(mvsiz) ,f2z(mvsiz) ,gama(6,mvsiz),
135 . rhocp(mvsiz) ,temp0(mvsiz),angle(mvsiz),dtx0(mvsiz),
136 . deltax(mvsiz), aire(mvsiz),llsh(mvsiz)
137 my_real,
DIMENSION(8,MVSIZ) :: bid8mvsiz
138 my_real,
DIMENSION(MVSIZ) :: bidmvsiz
139 my_real :: tempel(nel)
142 TYPE(g_bufel_) ,
POINTER :: GBUF
143 TYPE(BUF_LAY_) ,
POINTER :: BUFLY
144 TYPE(L_BUFEL_) ,
POINTER :: LBUF
145 TYPE(BUF_MAT_) ,
POINTER :: MBUF
148 . w_gauss(9,9),a_gauss(9,9)
157 3 0.555555555555556,0.888888888888889,0.555555555555556,
160 4 0.347854845137454,0.652145154862546,0.652145154862546,
161 4 0.347854845137454,0. ,0. ,
163 5 0.236926885056189,0.478628670499366,0.568888888888889,
164 5 0.478628670499366,0.236926885056189,0. ,
166 6 0.171324492379170,0.360761573048139,0.467913934572691,
167 6 0.467913934572691,0.360761573048139,0.171324492379170,
169 7 0.129484966168870,0.279705391489277,0.381830050505119,
170 7 0.417959183673469,0.381830050505119,0.279705391489277,
171 7 0.129484966168870,0. ,0. ,
172 8 0.101228536290376,0.222381034453374,0.313706645877887,
173 8 0.362683783378362,0.362683783378362,0.313706645877887,
174 8 0.222381034453374,0.101228536290376,0. ,
175 9 0.081274388361574,0.180648160694857,0.260610696402935,
176 9 0.312347077040003,0.330239355001260,0.312347077040003,
177 9 0.260610696402935,0.180648160694857,0.081274388361574/
182 2 -.577350269189626,0.577350269189626,0. ,
185 3 -.774596669241483,0. ,0.774596669241483,
188 4 -.861136311594053,-.339981043584856,0.339981043584856,
189 4 0.861136311594053,0. ,0. ,
191 5 -.906179845938664,-.538469310105683,0. ,
192 5 0.538469310105683,0.906179845938664,0. ,
194 6 -.932469514203152,-.661209386466265,-.238619186083197,
195 6 0.238619186083197,0.661209386466265,0.932469514203152,
197 7 -.949107912342759,-.741531185599394,-.405845151377397,
198 7 0. ,0.405845151377397,0.741531185599394,
199 7 0.949107912342759,0. ,0. ,
200 8 -.960289856497536,-.796666477413627,-.525532409916329,
201 8 -.183434642495650,0.183434642495650,0.525532409916329,
202 8 0.796666477413627,0.960289856497536,0. ,
203 9 -.968160239507626,-.836031107326636,-.613371432700590,
204 9 -.324253423403809,0. ,0.324253423403809,
205 9 0.613371432700590,0.836031107326636,0.968160239507626/
209 gbuf => elbuf_str%GBUF
210 lbuf => elbuf_str%BUFLY(1)%LBUF(1,1,1)
211 mbuf => elbuf_str%BUFLY(1)%MAT(1,1,1)
212 bufly => elbuf_str%BUFLY(1)
213 nptr = elbuf_str%NPTR
214 npts = elbuf_str%NPTS
215 nptt = elbuf_str%NPTT
216 nlay = elbuf_str%NLAY
223 IF (jcvt==1.AND.isorth/=0) jcvt=2
224 IF (igtyp /= 22) isorth = 0
229 rhocp(i) = pm(69,ixs(1,nft+i))
230 temp0(i) = pm(79,ixs(1,nft+i))
232 CALL sccoor3(x ,ixs(1,nf1),geo ,mat ,pid ,ngl ,
233 . ix1 ,ix2 ,ix3 ,ix4 ,ix5 ,ix6 ,ix7 ,ix8 ,
234 . x1 ,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,
235 . y1 ,y2 ,y3 ,y4 ,y5 ,y6 ,y7 ,y8 ,
236 . z1 ,z2 ,z3 ,z4 ,z5 ,z6 ,z7 ,z8 ,
237 . rx ,ry ,rz ,sx ,sy ,sz ,tx ,ty ,tz ,
238 . e1x ,e1y ,e1z ,e2x ,e2y ,e2z ,e3x ,e3y ,e3z ,
239 . f1x ,f1y ,f1z ,f2x ,f2y ,f2z ,temp0, temp,glob_therm%NINTEMP)
240 IF (igtyp == 21 .OR. igtyp == 22)
THEN
242 angle(i) = geo(1,pid(i))
244 CALL scmorth3(pid ,geo ,igeo ,skew ,irep ,gbuf%GAMA ,
245 . rx ,ry ,rz ,sx ,sy ,sz ,tx ,ty ,tz ,
246 . e1x ,e1y ,e1z ,e2x ,e2y ,e2z ,e3x ,e3y ,e3z ,
247 . ngl ,angle,nsigi,sigsp,nsigs,sigi ,ixs(1,nf1) ,1 ,
248 . orthoglob,ptsol,nel)
249 IF (igtyp == 22)
THEN
263 CALL scderi3(nel,gbuf%VOL,jeul,veul(1,nf1),geo ,
264 . vzl ,vzq ,ngl ,pid ,
265 . x1 ,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,
266 . y1 ,y2 ,y3 ,y4 ,y5 ,y6 ,y7 ,y8 ,
267 . z1 ,z2 ,z3 ,z4 ,z5 ,z6 ,z7 ,z8 , volu)
269 CALL sdlen3(x1 ,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,
270 . y1 ,y2 ,y3 ,y4 ,y5 ,y6 ,y7 ,y8 ,
271 . z1 ,z2 ,z3 ,z4 ,z5 ,z6 ,z7 ,z8 ,
275 . x1, x2, x3, x4, x5, x6, x7, x8,
276 . y1, y2, y3, y4, y5, y6, y7, y8,
277 . z1, z2, z3, z4, z5, z6, z7, z8)
279 IF (gbuf%IDT_TSH(i)>0)
280 . deltax(i)=
max(llsh(i),deltax(i))
284 IF (jthe == 0 .and. glob_therm%NINTEMP > 0)
THEN
286 tempel(i) = one_over_8 *(temp(ixs(2,i)) + temp(ixs(3,i))
287 . + temp(ixs(4,i)) + temp(ixs(5,i))
288 . + temp(ixs(6,i)) + temp(ixs(7,i))
289 . + temp(ixs(8,i)) + temp(ixs(9,i)))
292 tempel(1:nel) = temp0(1:nel)
296 CALL matini(pm ,ixs ,nixs ,x ,
297 . geo ,ale_connectivity ,detonators ,iparg ,
298 . sigi ,nel ,skew ,igeo ,
300 . mat ,ipm ,nsigs ,numsol ,ptsol ,
301 . ip ,ngl ,npf ,tf ,bufmat ,
302 . gbuf ,lbuf ,mbuf ,elbuf_str ,iloadp ,
303 . facload, deltax ,tempel ,mat_param )
305 IF (igtyp == 22)
CALL sczero3(gbuf%RHO,gbuf%SIG,gbuf%EINT,nel)
308 IF(jthe /=0)
CALL atheri(mat,pm,gbuf%TEMP)
313 lbuf => elbuf_str%BUFLY(is)%LBUF(1,1,1)
314 mbuf => elbuf_str%BUFLY(is)%MAT(1,1,1)
315 l_pla = elbuf_str%BUFLY(is)%L_PLA
316 l_sigb= elbuf_str%BUFLY(is)%L_SIGB
318 IF (igtyp == 22)
THEN
319 zi = geo(ippos+is,ig)
320 wi = geo(ipthk+is,ig)
325 angle(i) = geo(ipang+is,pid(i))
328 zi = a_gauss(is,nlay)
329 wi = w_gauss(is,nlay)
332 lbuf%VOL0DP(i) = half*wi*(gbuf%VOL(i)+vzl(i)*zi)
333 lbuf%VOL(i) = lbuf%VOL0DP(i)
336 .
CALL scmorth3(pid ,geo ,igeo ,skew ,irep ,lbuf%GAMA ,
337 . rx ,ry ,rz ,sx ,sy ,sz ,tx ,ty ,tz ,
338 . e1x ,e1y ,e1z ,e2x ,e2y ,e2z ,e3x ,e3y ,e3z ,
339 . ngl ,angle,nsigi,sigsp,nsigs,sigi ,ixs(1,nf1) ,is ,
340 . orthoglob,ptsol,nel)
342 IF (jthe == 0 .and. glob_therm%NINTEMP > 0)
THEN
344 tempel(i) = one_over_8 *(temp(ixs(2,i)) + temp(ixs(3,i))
345 . + temp(ixs(4,i)) + temp(ixs(5,i))
346 . + temp(ixs(6,i)) + temp(ixs(7,i))
347 . + temp(ixs(8,i)) + temp(ixs(9,i)))
350 tempel(1:nel) = temp0(1:nel)
353 CALL matini(pm ,ixs ,nixs ,x ,
354 . geo ,ale_connectivity ,detonators ,iparg ,
355 . sigi ,nel ,skew ,igeo ,
357 . mat ,ipm ,nsigs ,numsol ,ptsol,
358 . is ,ngl ,npf ,tf ,bufmat,
359 . gbuf ,lbuf ,mbuf ,elbuf_str,iloadp,
360 . facload, deltax ,tempel ,mat_param )
362 nuvar = elbuf_str%BUFLY(is)%NVAR_MAT
366 IF(mtn == 14 .OR. mtn == 12)
THEN
368 ELSEIF(mtn == 24)
THEN
370 ELSEIF(istrain == 1)
THEN
377 ELSEIF(mtn == 3.OR.mtn == 6.OR.mtn == 10.OR.
378 . mtn == 21.OR.mtn == 22.OR.mtn == 23.
386 . pm, lbuf%VOL,sigsp,sigi,lbuf%EINT,lbuf%RHO,mbuf%VAR ,
387 . lbuf%STRA,ixs ,nixs,nsigi, is, nuvar,nel,iuser,idef,
388 . nsigs,strsglob,straglob,jhbe,igtyp,x,gbuf%GAMA,
389 . mat ,lbuf%PLA,l_pla,ptsol,lbuf%SIGB,l_sigb,ipm ,
390 . bufmat,lbuf%VOL0DP)
394 CALL dtmain(geo ,pm ,ipm ,pid ,mat ,fv ,
395 . lbuf%EINT ,lbuf%TEMP ,lbuf%DELTAX ,lbuf%RK ,lbuf%RE ,bufmat, deltax, aire,
396 . volu, dtx , igeo ,igtyp)
399 . lbuf%RHO,lbuf%VOL,lbuf%OFF,lbuf%SIG,lbuf%EINT,dtx,
400 . gbuf%RHO,gbuf%VOL,gbuf%OFF,gbuf%SIG,gbuf%EINT,dtx0,
413 bid8mvsiz(1:8,1:mvsiz) = zero
414 bidmvsiz(1:mvsiz) = zero
416 . gbuf%RHO ,mas ,partsav ,x ,v ,
417 . iparts(nf1),mss(1,nf1),volu ,
418 . msnf ,mssf(1,nf1),bid ,
419 . bid ,bid8mvsiz ,wma ,rhocp ,mcp ,
420 . mcps(1,nf1),mssa ,bidmvsiz ,bidmvsiz ,gbuf%FILL,
421 . ix1, ix2, ix3, ix4, ix5, ix6, ix7, ix8)
424 CALL failini(elbuf_str,nptr,npts,nptt,nlay,
425 . ipm,sigsp,nsigi,fail_ini ,
426 . sigi,nsigs,ixs,nixs,ptsol,rnoise,perturb,bufmat)
432 CALL sbulk3(volu ,ix1 ,ncc,mat,pm ,
433 2 volnod,bvolnod,vns(1,nf1),bns(1,nf1),bid,
439 CALL dtmain(geo ,pm ,ipm ,pid ,mat ,fv ,
440 . lbuf%EINT ,lbuf%TEMP ,lbuf%DELTAX ,lbuf%RK ,lbuf%RE ,bufmat, deltax, aire,
450 IF (ixs(10,i+nft) /= 0)
THEN
451 IF (igtyp < 20 .OR. igtyp > 22)
THEN
452 ipid1=ixs(nixs-1,i+nft)
453 CALL fretitl2(titr1,igeo(npropgi-ltitr+1,ipid1),ltitr)
456 . anmode=aninfo_blind_1,
464 sti = fourth * gbuf%FILL(i) * gbuf%RHO(i) * volu(i) /
465 .
max(em20,dtx(i)*dtx(i))
466 stifn(ixs(2,i+nft))=stifn(ixs(2,i+nft))+sti
467 stifn(ixs(3,i+nft))=stifn(ixs(3,i+nft))+sti
468 stifn(ixs(4,i+nft))=stifn(ixs(4,i+nft))+sti
469 stifn(ixs(5,i+nft))=stifn(ixs(5,i+nft))+sti
470 stifn(ixs(6,i+nft))=stifn(ixs(6,i+nft))
471 stifn(ixs(7,i+nft))=stifn(ixs(7,i+nft))+sti
472 stifn(ixs(8,i+nft))=stifn(ixs(8,i+nft))+sti
473 stifn(ixs(9,i+nft))=stifn(ixs(9,i+nft))+sti
571 . X1, X2, X3, X4, X5, X6, X7, X8,
572 . Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,
573 . Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8)
577#include "implicit_f.inc"
581#include "mvsiz_p.inc"
585#include "scr17_c.inc"
590 my_real,
DIMENSION(MVSIZ),
INTENT(IN) ::
591 . X1, X2, X3, X4, X5, X6, X7, X8,
592 . Y1, Y2, , Y4, Y5, Y6, Y7, Y8,
593 . Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8
594 my_real,
DIMENSION(MVSIZ),
INTENT(OUT) :: LLSH
600 . RX(),RY(MVSIZ),RZ(MVSIZ),SX(MVSIZ),SY(MVSIZ),SZ(MVSIZ),
601 . VQ(MVSIZ,3,3), LXYZ0(3),DETA1(MVSIZ),XX,YY,ZZ,
602 . XL2(MVSIZ),(MVSIZ),XL4(MVSIZ),YL2(MVSIZ),
603 . YL3(MVSIZ),YL4(MVSIZ),ZL1(MVSIZ),AREA(MVSIZ),
604 . xn(mvsiz,4) , yn(mvsiz,4) , zn(mvsiz,4)
606 . al1,al2,ll(mvsiz),corel(2,4)
608 . x13,x24,y13,y24,l13,l24,c1,c2,thkly,posly,
609 . fac,visce,rx1,ry1,sx1,sy1,s1,fac1,fac2,faci,fac11,facdt
612 xn(i,1) = half*(x1(i)+x5(i))
613 yn(i,1) = half*(y1(i)+y5(i))
614 zn(i,1) = half*(z1(i)+z5(i))
615 xn(i,2) = half*(x2(i)+x6(i))
616 yn(i,2) = half*(y2(i)+y6(i))
617 zn(i,2) = half*(z2(i)+z6(i))
618 xn(i,3) = half*(x3(i)+x7(i))
619 yn(i,3) = half*(y3(i)+y7(i))
620 zn(i,3) = half*(z3(i)+z7(i))
621 xn(i,4) = half*(x4(i)+x8(i))
622 yn(i,4) = half*(y4(i)+y8(i))
623 zn(i,4) = half*(z4(i)+z8(i))
627 rx(i)=xn(i,2)+xn(i,3)-xn(i,1)-xn(i,4)
628 ry(i)=yn(i,2)+yn(i,3)-yn(i,1)-yn(i,4)
629 rz(i)=zn(i,2)+zn(i,3)-zn(i,1)-zn(i,4)
630 sx(i)=xn(i,3)+xn(i,4)-xn(i,1)-xn(i,2)
631 sy(i)=yn(i,3)+yn(i,4)-yn(i,1)-yn(i,2)
632 sz(i)=zn(i,3)+zn(i,4)-zn(i,1)-zn(i,2)
635 CALL clsys3(rx, ry, rz, sx, sy, sz,
636 . vq, deta1,nel ,mvsiz)
639 lxyz0(1)=fourth*(xn(i,1)+xn(i,2)+xn(i,3)+xn(i,4))
640 lxyz0(2)=fourth*(yn(i,1)+yn(i,2)+yn(i,3)+yn(i,4))
641 lxyz0(3)=fourth*(zn(i,1)+zn(i,2)+zn(i,3)+zn(i,4))
645 xl2(i)=vq(i,1,1)*xx+vq(i,2,1)*yy+vq(i,3,1)*zz
646 yl2(i)=vq(i,1,2)*xx+vq(i,2,2)*yy+vq(i,3,2)*zz
650 zl1(i)=vq(i,1,3)*xx+vq(i,2,3)*yy+vq(i,3,3)*zz
655 xl3(i)=vq(i,1,1)*xx+vq(i,2,1)*yy+vq(i,3,1)*zz
656 yl3(i)=vq(i,1,2)*xx+vq(i,2,2)*yy+vq(i,3,2)*zz
661 xl4(i)=vq(i,1,1)*xx+vq(i,2,1)*yy+vq(i,3,1)*zz
662 yl4(i)=vq(i,1,2)*xx+vq(i,2,2)*yy+vq(i,3,2)*zz
663 area(i)=fourth*deta1(i)
668 IF (idt1sol>0) facdt =four_over_3
671 lxyz0(1)=fourth*(xl2(i)+xl3(i)+xl4(i))
672 lxyz0(2)=fourth*(yl2(i)+yl3(i)+yl4(i))
674 corel(1,2)=xl2(i)-lxyz0(1)
675 corel(1,3)=xl3(i)-lxyz0(1)
676 corel(1,4)=xl4(i)-lxyz0(1)
678 corel(2,2)=yl2(i)-lxyz0(2)
679 corel(2,3)=yl3(i)-lxyz0(2)
680 corel(2,4)=yl4(i)-lxyz0(2)
681 x13=(corel(1,1)-corel(1,3))*half
682 x24=(corel(1,2)-corel(1,4))*half
683 y13=(corel(2,1)-corel(2,3))*half
684 y24=(corel(2,2)-corel(2,4))*half
689 c1 =corel(1,2)*corel(2,4)-corel(2,2)*corel(1,4)
690 c2 =corel(1,1)*corel(2,3)-corel(2,1)*corel(1,3)
691 al2 =
max(abs(c1),abs(c2))/area(i)
696 c1=sqrt(rx1*rx1+ry1*ry1)
697 c2=sqrt(sx1*sx1+sy1*sy1)
698 s1=fourth*(
max(c1,c2)/
min(c1,c2)-one)
699 fac1=
min(half,s1)+one
701 fac2=3.413*
max(zero,fac2-0.7071)
702 fac2=0.78+0.22*fac2*fac2*fac2
704 s1 = sqrt(faci*(facdt+al2)*al1)