38
39
40
41
42
43
44
45
46
47
48
49
51 use element_mod , only : nixs
52
53
54
55#include "implicit_f.inc"
56
57
58
59#include "mvsiz_p.inc"
60
61
62
63#include "vect01_c.inc"
64#include "com01_c.inc"
65#include "param_c.inc"
66
67
68
69
70
71
72
73
74
75
76
77 INTEGER :: NEL
78 INTEGER :: IXS(NIXS,*),IPM(NPROPMI,*)
79 my_real :: sig(nel,6),qvis(*),x(3,*),rho(*), vol(*)
80 my_real,
INTENT(IN) :: mom(nel,3), ssp(*)
81 my_real :: n1x(*), n2x(*), n3x(*), n4x(*), n5x(*), n6x(*),
82 . n1y(*), n2y(*), n3y(*), n4y(*), n5y(*), n6y(*),
83 . n1z(*), n2z(*), n3z(*), n4z(*), n5z(*), n6z(*)
84 my_real,
DIMENSION(MVSIZ,6),
INTENT(INOUT) :: svis
85
86
87
88 INTEGER :: I, II, IMAT, IALEFVM_FLG
89 INTEGER :: NC1(MVSIZ),NC2(MVSIZ),NC3(MVSIZ),NC4(MVSIZ),NC5(MVSIZ),NC6(MVSIZ),NC7(),NC8(MVSIZ)
90
92 . s1(mvsiz), s2(mvsiz), s3(mvsiz),
93 . s4(mvsiz), s5(mvsiz), s6(mvsiz),
94 . mass
95 my_real x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
96 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
97 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz),
98 . u_n(6), surf(6),
norm(6)
99
100 INTEGER :: idbf,idbl
101 LOGICAL :: debug_outp
102
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
116
117
118
119
120
121
122 DO i=1,nel
123 s1(i) = sig(i,1) + svis(i,1) - qvis(i)
124 s2(i) = sig(i,2) + svis(i,2) - qvis(i)
125 s3(i) = sig(i,3) + svis(i,3) - qvis(i)
126 s4(i) = sig(i,4) + svis(i,4)
127 s5(i) = sig(i,5) + svis(i,5)
128 s6(i) = sig(i,6) + svis(i,6)
129 ENDDO
130
131 DO i=1,nel
132 p(i) = -third*(s1(i)+s2(i)+s3(i))
133 ENDDO
134
135
136
137
138
139 DO i=1,nel
140 ii = i + nft
141
142 nc1(i)=ixs(2,ii)
143 nc2(i)=ixs(3,ii)
144 nc3(i)=ixs(4,ii)
145 nc4(i)=ixs(5,ii)
146 nc5(i)=ixs(6,ii)
147 nc6(i)=ixs(7,ii)
148 nc7(i)=ixs(8,ii)
149 nc8(i)=ixs(9,ii)
150
151
152 x1(i)=x(1,nc1(i))
153 y1(i)=x(2,nc1(i))
154 z1(i)=x(3,nc1(i))
155
156 x2(i)=x(1,nc2(i))
157 y2(i)=x(2,nc2(i))
158 z2(i)=x(3,nc2(i))
159
160 x3(i)=x(1,nc3(i))
161 y3(i)=x(2,nc3(i))
162 z3(i)=x(3,nc3(i))
163
164 x4(i)=x(1,nc4(i))
165 y4(i)=x(2,nc4(i))
166 z4(i)=x(3,nc4(i))
167
168 x5(i)=x(1,nc5(i))
169 y5(i)=x(2,nc5(i))
170 z5(i)=x(3,nc5(i))
171
172 x6(i)=x(1,nc6(i))
173 y6(i)=x(2,nc6(i))
174 z6(i)=x(3,nc6(i))
175
176 x7(i)=x(1,nc7(i))
177 y7(i)=x(2,nc7(i))
178 z7(i)=x(3,nc7(i))
179
180 x8(i)=x(1,nc8(i))
181 y8(i)=x(2,nc8(i))
182 z8(i)=x(3,nc8(i))
183 ENDDO
184 DO i=1,nel
185
186 n1x(i)=(y3(i)-y1(i))*(z2(i)-z4(i)) - (z3(i)-z1(i))*(y2(i)-y4(i))
187 n1y(i)=(z3(i)-z1(i))*(x2(i)-x4(i)) - (x3(i)-x1(i))*(z2(i)-z4(i))
188 n1z(i)=(x3(i)-x1(i))*(y2(i)-y4(i)) - (y3(i)-y1(i))*(x2(i)-x4(i))
189
190 n2x(i)=(y7(i)-y4(i))*(z3(i)-z8(i)) - (z7(i)-z4(i))*(y3(i)-y8(i))
191 n2y(i)=(z7(i)-z4(i))*(x3(i)-x8(i)) - (x7(i)-x4(i))*(z3(i)-z8(i))
192 n2z(i)=(x7(i)-x4(i))*(y3(i)-y8(i)) - (y7(i)-y4(i))*(x3(i)-x8(i))
193
194 n3x(i)=(y6(i)-y8(i))*(z7(i)-z5(i)) - (z6(i
195 n3y(i)=(z6(i)-z8(i))*(x7(i)-x5(i)) - (x6(i)-x8(i))*(z7(i)-z5(i))
196 n3z(i)=(x6(i)-x8(i))*(y7(i)-y5(i)) - (y6(i)-y8(i))*(x7(i)-x5(i))
197
198 n4x(i)=(y2(i)-y5(i))*(z6(i)-z1(i)) - (z2(i)-z5(i))*(y6(i)-y1(i))
199 n4y(i)=(z2(i)-z5(i))*(x6(i)-x1(i)) - (x2(i)-x5(i))*(z6(i)-z1(i))
200 n4z(i)=(x2(i)-x5(i))*(y6(i)-y1(i)) - (y2(i)-y5(i))*(x6(i)-x1(i))
201
202 n5x(i)=(y7(i)-y2(i))*(z6(i)-z3(i)) - (z7(i)-z2(i))*(y6(i)-y3(i))
203 n5y(i)=(z7(i)-z2(i))*(x6(i)-x3(i)) - (x7(i)-x2(i))*(z6(i)-z3(i))
204 n5z(i)=(x7(i)-x2(i))*(y6(i)-y3(i)) - (y7(i)-y2(i))*(x6(i)-x3(i))
205
206 n6x(i)=(y8(i)-y1(i))*(z4(i)-z5(i)) - (z8(i)-z1(i))*(y4(i)-y5(i))
207 n6y(i)=(z8(i)-z1(i))*(x4(i)-x5(i)) - (x8(i)-x1(i))*(z4(i)-z5(i))
208 n6z(i)=(x8(i)-x1(i))*(y4(i)-y5(i)) - (y8(i)-y1(i))*(x4(i)-x5(i))
209 ENDDO
210
211
212
213
214
215
216
217 DO i=1, nel
218 ii = i + nft
219 mass = rho(i)*vol(i)
220 norm(1) = sqrt(n1x(i)*n1x(i) + n1y(i)*n1y(i) + n1z(i)*n1z(i))
221 norm(2) = sqrt(n2x(i)*n2x(i) + n2y(i)*n2y(i) + n2z(i)*n2z(i))
222 norm(3) = sqrt(n3x(i)*n3x(i) + n3y(i)*n3y(i) + n3z(i)*n3z(i))
223 norm(4) = sqrt(n4x(i)*n4x(i) + n4y(i)*n4y(i) + n4z(i)*n4z(i))
224 norm(5) = sqrt(n5x(i)*n5x(i) + n5y(i)*n5y(i) + n5z(i)*n5z(i))
225 norm(6) = sqrt(n6x(i)*n6x(i) + n6y(i)*n6y(i) + n6z(i)*n6z(i))
226 u_n(1) = (mom(i,1)*n1x(i) + mom(i,2)*n1y(i) + mom(i,3)*n1z(i)) / (mass*
norm(1))
227 u_n(2) = (mom(i,1)*n2x(i) + mom(i,2)*n2y(i) + mom(i,3)*n2z(i)) / (mass*
norm(2))
228 u_n(3) = (mom(i,1)*n3x(i) + mom(i,2)*n3y(i) + mom(i,3)*n3z(i)) / (mass*
norm(3))
229 u_n(4) = (mom(i,1)*n4x(i) + mom(i,2)*n4y(i) + mom(i,3)*n4z(i)) / (mass*
norm(4))
230 u_n(5) = (mom(i,1)*n5x(i) + mom(i,2)*n5y(i) + mom(i,3)*n5z(i)) / (mass*
norm(5))
231 u_n(6) = (mom(i,1)*n6x(i) + mom(i,2)*n6y(i) + mom(i,3)*n6z(i)) / (mass*
norm(6))
232 surf(1) = half*
norm(1)
233 surf(2) = half*
norm(2)
234 surf(3) = half*
norm(3)
235 surf(4) = half*
norm(4)
236 surf(5) = half*
norm(5)
237 surf(6) = half*
norm(6)
238
239
244 alefvm_buffer%F_FACE(1 ,5 ,ii) = sqrt(mom(i,1)*mom(i,1)+mom(i,2)*mom(i,2)+mom(i,3)*mom(i,3))/rho(i)/vol(i) / ssp(i)
246
248
250 ENDDO
251
252
254 debug_outp = .false.
256 do i=lft,llt
257 ii = nft + i
259 debug_outp = .true.
260 idbf = i
261 idbl = i
262 EXIT
263 endif
264 enddo
266 debug_outp=.true.
267 idbf = lft
268 idbl = llt
269 endif
270 if(debug_outp)then
271
272 print *, " |----alefvm_stress.F-----|"
273 print *, " | THREAD INFORMATION |"
274 print *, " |------------------------|"
275 print *, " NCYCLE =", ncycle
276 do i=idbf,idbl
277 ii = nft + i
278 print *, " brique=", ixs(11,nft+i)
279 write(*,fmt='(A24,1A26)') " ",
280 . "#-stress Tensor (P+VIS+Q)#"
281
282 write (*,fmt='(A,3E26.14,A)') " | ", sig(i,1),sig(i,4),sig(i,6), " |"
283 write (*,fmt='(A,3E26.14,A)') " P =| ", sig(i,4),sig(i" |"
284 write (*,fmt='(A,3E26.14,A)') " |_", sig(i,6),sig(i,5),sig(i,3), "_|"
285 write (*,fmt='(A,3E26.14,A)') " | ", svis(i,1),svis(i,4),svis(i,6)," |"
286 write (*,fmt='(A,3E26.14,A)') " VIS =| ", svis(i,4),svis(i,2),svis(i,5)," |"
287 write (*,fmt='(A,3E26.14,A)') " |_", svis(i,6),svis(i,5),svis(i,3"_|"
288 write (*,fmt='(A,3E26.14,A)') " | ", s1(i),s4(i),s6(i)," |"
289 write (*,fmt='(A,3E26.14,A)') " sigma = p+vis+q =| ", S4(I),S2(I),S5(I)," |"
290 write (*,FMT='(A,3E26.14,A)') " |_", S6(I),S5(I),S3(I),"_|"
291 write (*,FMT='(A,2E26.14)') " ",-THIRD*SUM(SIG(I,1:3)),P(I)
292 write (*,FMT='(A,1E26.14)') " q = ",QVIS(I)
293 write (*,FMT='(A,1E26.14)') " m = ",ALEFVM_Buffer%F_FACE(1 ,5 ,II)
294 write (*,FMT='(A,1E26.14)') " rho = ",ALEFVM_Buffer%F_FACE(1 ,1 ,II)
295 write (*,FMT='(A,1E26.14)') " u = ",ALEFVM_Buffer%F_FACE(1 ,5 ,II)*SSP(I)
296 write (*,FMT='(A,1E26.14)') " ssp = ",ALEFVM_Buffer%F_FACE(1 ,2 ,II)
297 write (*,FMT='(A,1E26.14)') " z = ",ALEFVM_Buffer%F_FACE(1 ,3 ,II)
298 write (*,FMT='(A,3E26.14)') " rho.u = ",MOM(I,1:3)
299 write (*,FMT='(A,3E26.14)') " v = ",VOL(I)
300 write(*,FMT='(A34,6A26)') " ",
301 . "#-------- face_1 ---------","#-------- face_2 ---------",
302 . "#-------- face_3 ---------","#-------- face_4 ---------",
303 . "#-------- face_5 ---------","#-------- face_6 --------#"
305 write (*,fmt=
'(A,8E26.14)')
" <U,n> =",
alefvm_buffer%F_FACE(2,1:6,ii)
306 print *, " "
307 enddo
308
309 endif
310 endif
311
312
313 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
type(alefvm_buffer_), target alefvm_buffer
type(alefvm_param_), target alefvm_param