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