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