40
41
42
43
44
45
46
47
48
49
50
51
53 USE elbufdef_mod
55 use element_mod , only : nixs
56
57
58
59#include "implicit_f.inc"
60
61
62
63#include "mvsiz_p.inc"
64
65
66
67#include "vect01_c.inc"
68#include "com01_c.inc"
69#include "param_c.inc"
70
71
72
73
74
75
76
77
78
79
80
81 INTEGER :: NEL
82 INTEGER :: IXS(NIXS,*),IPM(NPROPMI,*)
83 my_real :: sig(nel,6),qvis(nel),rho(nel), vol(nel), iad22(*)
84 my_real :: n1x(*), n2x(*), n3x(*), n4x(*), n5x(*), n6x(*),
85 . n1y(*), n2y(*), n3y(*), n4y(*), n5y(*), n6y
86 . n1z(*), n2z(*), n3z(*), n4z
87 . mom(nel,3), ssp(nel), isgn, v(3), w(3), vmw(3)
88 my_real,
DIMENSION(MVSIZ,6),
INTENT(INOUT) :: svis
89
90
91
92 INTEGER :: I, II, K, IMAT, IALEFVM_FLG, IB, NIN, IBM, NGM, IDLOCM, ICUT,MCELL,NUM
93 INTEGER :: ISECND, ID, IBs, ICELLs, Js, NumSECND
94
96 . s1(mvsiz) , s2(mvsiz) , s3(mvsiz),
97 . s4(mvsiz) , s5(mvsiz) , s6(mvsiz),
99 my_real :: n0x(14), n0y(14), n0z(14),mass,u_n(0:6)
100
101 INTEGER :: idbf,idbl,ICELL,NBCUT,NCELL,,IDEB,IFIN
102 LOGICAL :: debug_outp
103
104
105
107 imat = ixs(1,1+nft)
108 ialefvm_flg = ipm(251,imat)
109 IF(ialefvm_flg <= 1)RETURN
110
111
112
113
114
115 nin = 1
116
117
118
119
120
121
122
123
124
125 DO i=1,nel
126 s1(i) = sig(i,1) + svis(i,1) - qvis(i)
127 s2(i) = sig(i,2) + svis(i,2) - qvis(i)
128 s3(i) = sig(i,3) + svis(i,3) - qvis(i)
129 s4(i) = sig(i,4) + svis(i,4)
130 s5(i) = sig(i,5) + svis(i,5)
131 s6(i) = sig(i,6) + svis(i,6)
132 ENDDO
133 DO i=1,nel
134 p(i) = -third*(s1(i)+s2(i)+s3(i))
135 ENDDO
136
137 DO i=1,nel
138 ib = nint(iad22(i))
139 IF(ib<=0) cycle
148 brick_list(nin,ib)%MACH = sqrt(mom(i,1)*mom(i,1)+mom(i,2)*mom(i,2)+mom(i,3)*mom(i,3))/rho(i)/vol(i)
149 enddo
150
151
152
153
154 DO
155 ii = i + nft
156 ib = nint(iad22(i))
157 IF(ib<=0)cycle
158
162
166
170
174
178
182 ENDDO
183
184
185
186
187
188
189 DO i=1, nel
190
191 ii = i + nft
192 ib = nint(iad22(i))
193 IF(ib<=0)cycle
196 ncell = nbcut
197 icell = mcell
198 brick_list(nin,ib)%POLY(icell)%FACE(1:6)%U_N = zero
199
200 IF(nbcut>0)THEN
201 s(1:6) =
brick_list(nin,ib)%POLY(icell)%FACE(1:6)%Surf
203 ibm =
brick_list(nin,ib)%POLY(icell)%WhereIsMain(4)
206 isgn = one
207 IF(icell==9)isgn = -one
208 brick_list(nin,ib)%POLY(icell)%FACE0%U_N(1:9) = zero
209 DO k=1,nbcut
210 icut = k
215 norm = sqrt(n0x(k)**2+n0y(k)**2+n0z(k)**2)
216 n0x(k) = n0x(k) /
norm
217 n0y(k) = n0y(k) /
norm
218 n0z(k) = n0z(k) /
norm
219 mass = rho(i)*vol(i)
220 v(1:3) = mom(i,1:3) / mass
221 w(1:3) =
brick_list(nin,ib)%PCUT(icut)%VEL(1:3)
222 vmw(1) = v(1)-w(1)
223 vmw(2) = v(2)-w(2)
224 vmw(3) = v(3)-w(3)
225 u_n(0) = isgn*(vmw(1)*n0x(k) + vmw(2)*n0y(k) +
226 brick_list(nin,ib)%POLY(icell)%FACE0%U_N(k) = u_n(0)
227 ENDDO
228 ENDIF
229
230 mass = rho(i)*vol(i)
231 u_n(1) = (mom(i,1)*n1x(i) + mom(i,2)*n1y(i) + mom(i,3)*n1z(i)) / (mass)
232 u_n(2) = (mom(i,1)*n2x(i) + mom(i,2)*n2y(i) + mom(i,3)*n2z(i)) / (mass)
233 u_n(3) = (mom(i,1)*n3x(i) + mom(i,2)*n3y(i) + mom(i,3)*n3z(i)) / (mass)
234 u_n(4) = (mom(i,1)*n4x(i) + mom(i,2)*n4y(i) + mom(i,3)*n4z(i)) / (mass)
235 u_n(5) = (mom(i,1)*n5x(i) + mom(i,2)*n5y(i) + mom(i,3)*n5z(i)) / (mass)
236 u_n(6) = (mom(i,1)*n6x(i) + mom(i,2)*n6y(i) + mom(i,3)*n6z(i)) / (mass)
237 brick_list(nin,ib)%POLY(icell)%FACE(1)%U_N = u_n(1)
238 brick_list(nin,ib)%POLY(icell)%FACE(2)%U_N = u_n(2)
239 brick_list(nin,ib)%POLY(icell)%FACE(3)%U_N = u_n(3)
240 brick_list(nin,ib)%POLY(icell)%FACE(4)%U_N = u_n(4)
241 brick_list(nin,ib)%POLY(icell)%FACE(5)%U_N = u_n(5)
242 brick_list(nin,ib)%POLY(icell)%FACE(6)%U_N = u_n(6)
243
245 DO isecnd=1,numsecnd
246 ibs =
brick_list(nin,ib)%SecndList%IBV(isecnd)
247 icells =
brick_list(nin,ib)%SecndList%ICELLv(isecnd)
248 js =
brick_list(nin,ibs)%POLY(icells)%WhereIsMain(1)
250
251
252 s(1:6) =
brick_list(nin,ibs)%POLY(icells)%FACE(1:6)%Surf
253 s(0) =
brick_list(nin,ibs)%PCUT(icells)%SCUT(1)
254 ibm =
brick_list(nin,ibs)%POLY(icells)%WhereIsMain(4)
255 isgn = one
256 ideb = icells
257 ifin = icells
258 IF(icells==9)THEN
259 isgn = -one
260 ideb = 1
261 ifin = nbcutv
262 ENDIF
263 brick_list(nin,ibs)%POLY(icells)%FACE0%U_N(1:9)= zero
264 mass = rho(i)*vol(i)
265 DO k=ideb,ifin
266 icut = k
267 s0(k) =
brick_list(nin,ibs)%PCUT(icut)%Scut(1)
271 norm = sqrt(n0x(k)**2+n0y(k)**2+n0z(k)**2)
272 n0x(k) = n0x(k) /
norm
273 n0y(k) = n0y(k) /
norm
274 n0z(k) = n0z(k) /
norm
275 v(1:3) = mom(i,1:3) / mass
277 vmw(1) = v(1)-w(1)
278 vmw(2) = v(2)-w(2)
279 vmw(3) = v(3)-w(3)
280 u_n(0) = isgn*(vmw(1)*n0x(k) + vmw(2)*n0y
281 brick_list(nin,ibs)%POLY(icells)%FACE0%U_N(k) = u_n(0)
282 enddo
283
284 u_n(1) = (mom(i,1)*n1x(i) + mom
285 u_n(2) = (mom(i,1)*n2x(i) + mom(i,2)*n2y(i) + mom(i,3)*n2z(i)) / mass
286 u_n(3) = (mom(i,1)*n3x(i
287 u_n(4) = (mom(i,1)*n4x(i) + mom(i,2)*n4y(i) +
288 u_n(5) = (mom(i,1)*n5x(i) + mom(i,2)*n5y(i) + mom(i,3)*n5z(i)) / mass
289 u_n(6) = (mom(i,1)*n6x(i) + mom(i,2)*n6y(i) + mom(i,3)*n6z(i)) / mass
290 brick_list(nin,ibs)%POLY(icells)%FACE(1)%U_N = u_n(1)
291 brick_list(nin,ibs)%POLY(icells)%FACE(2)%U_N = u_n(2)
292 brick_list(nin,ibs)%POLY(icells)%FACE(3)%U_N = u_n(3)
293 brick_list(nin,ibs)%POLY(icells)%FACE(4)%U_N = u_n(4)
294 brick_list(nin,ibs)%POLY(icells)%FACE(5)%U_N = u_n(5)
295 brick_list(nin,ibs)%POLY(icells)%FACE(6)%U_N = u_n(6)
296 ENDDO
297 enddo
298
299
300
301
303 debug_outp = .false.
305 do i=lft,llt
306 ii = nft + i
308 debug_outp = .true.
309 idbf = i
310 idbl = i
311 EXIT
312 endif
313 enddo
315 debug_outp=.true.
316 idbf = lft
317 idbl = llt
318 endif
319 if(debug_outp)then
320
321 print *, " |--alefvm_stress_int22.F---|"
322 print *, " | THREAD INFORMATION |"
323 print *, " |--------------------------|"
324 print *, " NCYCLE =", ncycle
325 do i=idbf,idbl
326 ii = nft + i
327 ib = nint(iad22(i))
328 if(ib<=0)cycle
330 IF(nbcut>0)THEN
332 IF(icell==0)cycle
333
334 print *, " brique=", ixs(11,nft+i), "icell=", icell
335 print *,
" vold=",
brick_list(nin,ib)%vold_scell
336 print *,
" vnew=",
brick_list(nin,ib)%vnew_scell
337 write(*,fmt='(A24,1A26)') " ",
338 . "#-stress Tensor (P+VIS+Q)#"
339
340 write (*,fmt='(A,3E26.14,A)') " | ", sig(i,1),sig(i," |"
341 write (*,fmt='(A,3E26.14,A)') " P =| ", sig(i,4),sig(i,2),sig(i" |"
342 write (*,fmt='(A,3E26.14,A)') " |_", sig(i,6),sig(i,5),sig(i,3), "_|"
343 write (*,fmt='(A,3E26.14,A)') " | ", svis(i" |"
344 write (*,fmt='(A,3E26.14,A)') " VIS =| ", svis(i,4),svis(i,2),svis(i,5)," |"
345 write (*,fmt='(A,3E26.14,A)') " |_", svis(i"_|"
346 write (*,fmt='(A,3E26.14,A)') " | ", s1(i" |"
347 write (*,fmt='(A,3E26.14,A)') " SIGMA = P+VIS+Q =| ", s4(i),s2" |"
348 write (*,fmt='(A,3E26.14,A)') " |_", s6(i),s5"_|"
349 write (*,fmt='(A,2E26.14)') " p = ",-third
350 write (*,fmt='(A,1E26.14)') " q = ",qvis(i)
351 write(*,fmt='(A34,8A26)') " ",
352 . "#-------- face_1 ---------",
353 . "#-------- face_2 ---------","#-------- face_3 ---------"
354 . "#-------- face_4 ---------","#-------- face_5 ---------"
355 . "#-------- face_6 --------#"
356 write (*,fmt=
'(A,8E26.14)')
" <U,N> =",
brick_list(nin,ib)%POLY(icell)%FACE(1:6)%U_N
357 write (*,fmt=
'(A,9E26.14)')
" U0N0 =",
brick_list(nin,ib)%POLY
358 ELSE
359 print *, " brique=", ixs(11,nft+i)"uncut"
360 print *,
" vold=",
brick_list(nin,ib)%vold_scell
362 write(*,fmt='(A24,1A26)') " ",
363 . "#-stress Tensor (P+VIS+Q)#"
364
365 write (*,fmt='(A,3E26.14,A)') " | ", sig(i,1),sig(i,4),sig(i,6), " |"
366 write (*,fmt='(A,3E26.14,A)') " P =| ", sig(i,4),sig(i" |"
367 write (*,fmt='(A,3E26.14,A)') " |_", sig(i,6),sig(i,5),sig(i,3), "_|"
368 write (*,fmt='(A,3E26.14,A)') " | ", svis(i,1),svis(i,4),svis(i,6)," |"
369 write (*,fmt='(A,3E26.14,A)') " VIS =| ", svis(i,4)" |"
370 write (*,fmt='(A,3E26.14,A)') " |_", svis(i"_|"
371 write (*,fmt='(a,3e26.14,a)') " | ", S1(I),S4(I),S6(I)," |"
372 write (*,FMT='(a,3e26.14,a)') " SIGMA = P+VIS+Q =| ", S4(I),S2(I),S5(I)," |"
373 write (*,FMT='(a,3e26.14,a)') " |_", S6(I),S5(I),S3(I),"_|"
374 write (*,FMT='(a,2e26.14)') " p = ",-THIRD*SUM(SIG(I,1:3)),P(I)
375 write (*,FMT='(a,1e26.14)') " q = ",QVIS(I)
376 write(*,FMT='(a34,6a26)') " ",
377 . "#-------- face_1 ---------","#-------- face_2 ---------",
378 . "#-------- face_3 ---------","#-------- face_4 ---------",
379 . "#-------- face_5 ---------","#-------- face_6 --------#"
380 write (*,FMT='(a,1e26.14)') " <U,N> =", BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(1:6)%U_N
381 ENDIF!(NBCUT>0)
382
383 NUM = BRICK_LIST(NIN,IB)%SecndList%Num
384 IBM = IB
385 DO ISECND=1,NUM
386 IB = BRICK_LIST(NIN,IBM)%SecndList%IBV(ISECND)
387 ICELL = BRICK_LIST(NIN,IBM)%SecndList%ICELLv(ISECND)
388 id = BRICK_LIST(NIN,IB)%ID
389 print *, " secnd=", ixs(11,id), "icell=", ICELL
390 write(*,FMT='(a24,1a26)') " ",
391 . "#-stress Tensor (P+VIS+Q)#"
392
393 write (*,FMT='(a,3e26.14,a)') " | ", SIG(I,1),SIG(I,4),SIG(I,6), " |"
394 write (*,FMT='(a,3e26.14,a)') " P =| ", SIG(I,4),SIG(I,2),SIG(I,5), " |"
395 write (*,FMT='(a,3e26.14,a)') " |_", SIG(I,6),SIG(I,5),SIG(I,3), "_|"
396 write (*,FMT='(a,3e26.14,a)') " | ", SVIS(I,1),SVIS(I,4),SVIS(I,6)," |"
397 write (*,FMT='(a,3e26.14,a)') " VIS =| ", SVIS(I,4),SVIS(I,2),SVIS(I,5)," |"
398 write (*,FMT='(a,3e26.14,a)') " |_", SVIS(I,6),SVIS(I,5),SVIS(I,3),"_|"
399 write (*,FMT='(a,3e26.14,a)') " | ", S1(I),S4(I),S6(I)," |"
400 write (*,FMT='(a,3e26.14,a)') " SIGMA = P+VIS+Q =| ", S4(I),S2(I),S5(I)," |"
401 write (*,FMT='(a,3e26.14,a)') " |_", S6(I),S5(I),S3(I),"_|"
402 write (*,FMT='(a,2e26.14)') " p = ",-THIRD*SUM(SIG(I,1:3)),P(I)
403 write (*,FMT='(a,1e26.14)') " q = ",QVIS(I)
404 write(*,FMT='(a34,8a26)') " ",
405 . "#-------- face_1 ---------",
406 . "#-------- face_2 ---------","#-------- face_3 ---------",
407 . "#-------- face_4 ---------","#-------- face_5 ---------",
408 . "#-------- face_6 --------#"
409 write (*,FMT='(a,8e26.14)') " <U,N> =", BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE(1:6)%U_N
410 write (*,FMT='(a,9e26.14)') " U0N0 =", BRICK_LIST(NIN,IB)%POLY(ICELL)%FACE0%U_N(1:9)
411 ENDDO !next ICELL
412
413
414
415
416
417
418 print *, " "
419 enddo!next i
420!#!include "lockoff.inc"
421 endif
422 endif
423 !-----------------------------------------!
424
425 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
type(alefvm_param_), target alefvm_param
type(brick_entity), dimension(:,:), allocatable, target brick_list