34
35
36
37
38
39
40
41
42
43
44
45 USE elbufdef_mod
46
47
48
49#include "implicit_f.inc"
50
51
52
53 INTEGER,INTENT(IN) :: NUMEL,NEL,NUMNOD,NPARG,NGROUP,N2D,NFT,ITYP
54 INTEGER,INTENT(IN) :: IX(NIX,NUMEL),IPARG(NPARG,NGROUP),NIX,NG
55 my_real,
INTENT(IN) :: x(3,numnod),v(3,numnod)
56 my_real,
INTENT(OUT) :: evar(nel)
57 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP), INTENT(IN), TARGET :: ELBUF_TAB
58
59
60
61 INTEGER I, , IE
62 my_real :: div , xi(8), yi(8), zi(8)
63 INTEGER :: NC(8)
64 my_real :: n(3,6), vel(6,3),
area, nx, a1, a2, dotprod(6,nel), vol(nel)
65
66
67
68
69 mlw = iparg(01,ng)
70
71 DO i=1,nel
72 dotprod(1,i)=zero
73 dotprod(2,i)=zero
74 dotprod(3,i)=zero
75 dotprod(4,i)=zero
76 dotprod(5,i)=zero
77 dotprod(6,i)=zero
78 ENDDO
79
80
81
82
83 IF(n2d == 0 .AND. ityp == 1)THEN
84 DO i=1,nel
85 ie=i+nft
86
87 nc(1)=ix(2,ie)
88 nc(2)=ix(3,ie)
89 nc(3)=ix(4,ie)
90 nc(4)=ix(5,ie)
91 nc(5)=ix(6,ie)
92 nc(6)=ix(7,ie)
93 nc(7)=ix(8,ie)
94 nc(8)=ix(9,ie)
95
96 xi(1)=x(1,nc(1))
97 yi(1)=x(2,nc(1))
98 zi(1)=x(3,nc(1))
99
100 xi(2)=x(1,nc(2))
101 yi(2)=x(2,nc(2))
102 zi(2)=x(3,nc(2))
103
104 xi(3)=x(1,nc(3))
105 yi(3)=x(2,nc(3))
106 zi(3)=x(3,nc(3))
107
108 xi(4)=x(1,nc(4))
109 yi(4)=x(2,nc(4))
110 zi(4)=x(3,nc(4))
111
112 xi(5)=x(1,nc(5))
113 yi(5)=x(2,nc(5))
114 zi(5)=x(3,nc(5))
115
116 xi(6)=x(1,nc(6))
117 yi(6)=x(2,nc(6))
118 zi(6)=x(3,nc(6))
119
120 xi(7)=x(1,nc(7))
121 yi(7)=x(2,nc(7))
122 zi(7)=x(3,nc(7))
123
124 xi(8)=x(1,nc(8))
125 yi(8)=x(2,nc(8))
126 zi(8)=x(3,nc(8))
127
128 n(1,1)=(yi(3)-yi(1))*(zi(2)-zi(4)) - (zi(3)-zi(1))*(yi(2)-yi(4))
129 n(2,1)=(zi(3)-zi(1))*(xi(2)-xi(4)) - (xi(3)-xi(1))*(zi(2)-zi(4))
130 n(3,1)=(xi(3)-xi(1))*(yi(2)-yi(4)) - (yi(3)-yi(1))*(xi(2)-xi(4))
131
132 n(1,2)=(yi(7)-yi(4))*(zi(3)-zi(8)) - (zi(7)-zi(4))*(yi(3)-yi(8))
133 n(2,2)=(zi(7)-zi(4))*(xi(3)-xi(8)) - (xi(7)-xi(4))*(zi(3)-zi(8))
134 n(3,2)=(xi(7)-xi(4))*(yi(3)-yi(8)) - (yi(7)-yi(4))*(xi(3)-xi(8))
135
136 n(1,3)=(yi(6)-yi(8))*(zi(7)-zi(5)) - (zi(6)-zi(8))*(yi(7)-yi(5))
137 n(2,3)=(zi(6)-zi(8))*(xi(7)-xi(5)) - (xi(6)-xi(8))*(zi(7)-zi(5))
138 n(3,3)=(xi(6)-xi(8))*(yi(7)-yi(5)) - (yi(6)-yi(8))*(xi(7)-xi(5))
139
140 n(1,4)=(yi(2)-yi(5))*(zi(6)-zi(1)) - (zi(2)-zi(5))*(yi(6)-yi(1))
141 n(2,4)=(zi(2)-zi(5))*(xi(6)-xi(1)) - (xi(2)-xi(5))*(zi(6)-zi(1))
142 n(3,4)=(xi(2)-xi(5))*(yi(6)-yi(1)) - (yi(2)-yi(5))*(xi(6)-xi(1))
143
144 n(1,5)=(yi(7)-yi(2))*(zi(6)-zi(3)) - (zi(7)-zi(2))*(yi(6)-yi(3))
145 n(2,5)=(zi(7)-zi(2))*(xi(6)-xi(3)) - (xi(7)-xi(2))*(zi(6)-zi(3))
146 n(3,5)=(xi(7)-xi(2))*(yi(6)-yi(3)) - (yi(7)-yi(2))*(xi(6)-xi(3))
147
148 n(1,6)=(yi(8)-yi(1))*(zi(4)-zi(5)) - (zi(8)-zi(1))*(yi(4)-yi(5))
149 n(2,6)=(zi(8)-zi(1))*(xi(4)-xi(5)) - (xi(8)-xi(1))*(zi(4)-zi(5))
150 n(3,6)=(xi(8)-xi(1))*(yi(4)-yi(5)) - (yi(8)-yi(1))*(xi(4)-xi(5))
151
152 n(1:3,1) = half * n(1:3,1)
153 n(1:3,2) = half * n(1:3,2)
154 n(1:3,3) = half * n(1:3,3)
155 n(1:3,4) = half * n(1:3,4)
156 n(1:3,5) = half * n(1:3,5)
157 n(1:3,6) = half * n(1:3,6)
158
159 vel(1,1) = fourth*( v(1,nc(1)) + v(1,nc(2)) + v(1,nc(3)) + v(1,nc(4)) )
160 vel(1,2) = fourth*( v(2,nc(1)) + v(2,nc(2)) + v(2,nc(3)) + v(2,nc(4)) )
161 vel(1,3) = fourth*( v(3,nc(1)) + v(3,nc(2)) + v(3,nc(3)) + v(3,nc(4)) )
162
163 vel(2,1) = fourth*( v(1,nc(3)) + v(1,nc(4)) + v(1,nc(7)) + v(1,nc(8)) )
164 vel(2,2) = fourth*( v(2,nc(3)) + v(2,nc(4)) + v(2,nc(7)) + v(2,nc(8)) )
165 vel(2,3) = fourth*( v(3,nc(3)) + v(3,nc(4)) + v(3,nc(7)) + v(3,nc(8)) )
166
167 vel(3,1) = fourth*( v(1,nc(5)) + v(1,nc(6)) + v(1,nc(7)) + v(1,nc(8)) )
168 vel(3,2) = fourth*( v(2,nc(5)) + v(2,nc(6)) + v(2,nc(7)) + v(2,nc(8)) )
169 vel(3,3) = fourth*( v(3,nc(5)) + v(3,nc(6)) + v(3,nc(7)) + v(3,nc(8)) )
170
171 vel(4,1) = fourth*( v(1,nc(1)) + v(1,nc(2)) + v(1,nc(5)) + v
172 vel(4,2) = fourth*( v(2,nc(1)) + v(2,nc(2)) + v(2,nc(5)) + v(2,nc(6)) )
173 vel(4,3) = fourth*( v(3,nc(1)) + v(3,nc(2)) + v(3,nc(5)) + v(3,nc(6)) )
174
175 vel(5,1) = fourth*( v(1,nc(2)) + v(1,nc(3)) + v(1,nc(6)) + v(1,nc(7)) )
176 vel(5,2) = fourth*( v(2,nc(2)) + v(2,nc(3)) + v(2,nc(6)) + v(2,nc(7)) )
177 vel(5,3) = fourth*( v(3,nc(2)) + v(3,nc(3)) + v(3,nc(6)) + v(3,nc(7)) )
178
179 vel(6,1) = fourth*( v(1,nc(1)) + v(1,nc(4)) + v(1,nc(5)) + v(1,nc(8)) )
180 vel(6,2) = fourth*( v(2,nc(1)) + v(2,nc(4)) + v(2,nc(5)) + v(2,nc(8)) )
181 vel(6,3) = fourth*( v(3,nc(1)) + v(3,nc(4)) + v(3,nc(5)) + v(3,nc(8)) )
182
183 dotprod(1,i) = vel(1,1)*n(1,1) + vel(1,2)*n(2,1) + vel(1,3)*n(3,1)
184 dotprod(2,i) = vel(2,1)*n(1,2) + vel(2,2)*n(2,2) + vel(2,3)*n(3,2)
185 dotprod(3,i) = vel(3,1)*n(1,3) + vel(3,2)*n(2,3) + vel(3,3)*n(3,3)
186 dotprod(4,i) = vel(4,1)*n(1,4) + vel(4,2)*n(2,4) + vel(4,3)*n(3,4)
187 dotprod(5,i) = vel(5,1)*n(1,5) + vel(5,2)*n(2,5) + vel(5,3)*n(3,5)
188 dotprod(6,i) = vel(6,1)*n(1,6) + vel(6,2)*n(2,6) + vel(6,3)*n(3,6)
189
190 vol(i) = elbuf_tab(ng)%GBUF%VOL(i)
191 ENDDO
192 ELSEIF(n2d > 0 .AND. ityp == 2)THEN
193 DO i=1,nel
194 ie=i+nft
195
196 nc(1)=ix(2,ie)
197 nc(2)=ix(3,ie)
198 nc(3)=ix(4,ie)
199 nc(4)=ix(5,ie)
200
201
202 xi(1)=zero
203 yi(1)=x(2,nc(1))
204 zi(1)=x(3,nc(1))
205
206 xi(2)=zero
207 yi(2)=x(2,nc(2))
208 zi(2)=x(3,nc(2))
209
210 xi(3)=zero
211 yi(3)=x(2,nc(3))
212 zi(3)=x(3,nc(3))
213
214 xi(4)=zero
215 yi(4)=x(2,nc(4))
216 zi(4)=x(3,nc(4))
217
218
219 n(1,1) = zero
220 n(2,1) = (zi(2)-zi(1))
221 n(3,1) =-(yi(2)-yi(1))
222
223 n(1,2) = zero
224 n(2,2) = (zi(3)-zi(2))
225 n(3,2) =-(yi(3)-yi(2))
226
227 n(1,3) = zero
228 n(2,3) = (zi(4)-zi(3))
229 n(3,3) =-(yi(4)-yi(3))
230
231 n(1,4) = zero
232 n(2,4) = (zi(1)-zi(4))
233 n(3,4) =-(yi(1)-yi(4))
234
235 n(1:3,5:6)=zero
236
237 vel(1,2)=half*(v(2,nc(1)) + v(2,nc(2)))
238 vel(1,3)=half*(v(3,nc(1)) + v(3,nc(2)))
239
240 vel(2,2)=half*(v(2,nc(2)) + v(2,nc(3)))
241 vel(2,3)=half*(v(3,nc(2)) + v(3,nc(3)))
242
243 vel(3,2)=half*(v(2,nc(3)) + v(2,nc(4)))
244 vel(3,3)=half*(v(3,nc(3)) + v(3,nc(4)))
245
246 vel(4,2)=half*(v(2,nc(4)) + v(2,nc(1)))
247 vel(4,3)=half*(v(3,nc(4)) + v(3,nc(1)))
248
249 dotprod(1,i) = vel(1,2)*n(2,1) + vel(1,3)*n(3,1)
250 dotprod(2,i) = vel(2,2)*n(2,2) + vel(2,3)*n(3,2)
251 dotprod(3,i) = vel(3,2)*n(2,3) + vel(3,3)*n(3,3)
252 dotprod(4,i) = vel(4,2)*n(2,4) + vel(4,3)*n(3,4)
253 dotprod(5,i) = zero
254 dotprod(6,i) = zero
255
256 IF(mlw /= 151)THEN
257 a1 = yi(2)*(zi(3)-zi(4))+yi(3)*(zi(4)-zi(2))+yi(4)*(zi(2)-zi(3))
258 a2 = yi(2)*(zi(4)-zi(1))+yi(4)*(zi(1)-zi(2))+yi(1)*(zi(2)-zi(4))
260 ELSE
261 area = elbuf_tab(ng)%GBUF%AREA(i)
262 ENDIF
264 ENDDO
265 ELSEIF(n2d > 0 .AND. ityp == 7)THEN
266 DO i=1,nel
267 ie=i+nft
268
269 nc(1)=ix(2,ie)
270 nc(2)=ix(3,ie)
271 nc(3)=ix(4,ie)
272
273 xi(1)=zero
274 yi(1)=x(2,nc(1))
275 zi(1)=x(3,nc(1))
276
277 xi(2)=zero
278 yi(2)=x(2,nc(2))
279 zi(2)=x(3,nc(2))
280
281 xi(3)=zero
282 yi(3)=x(2,nc(3))
283 zi(3)=x(3,nc(3))
284
285
286 n(1,1) = zero
287 n(2,1) = (zi(2)-zi(1))
288 n(3,1) =-(yi(2)-yi(1))
289
290 n(1,2) = zero
291 n(2,2) = (zi(3)-zi(2))
292 n(3,2) =-(yi(3)-yi(2))
293
294 n(1,3) = zero
295 n(2,3) = (zi(1)-zi(3))
296 n(3,3) =-(yi(1)-yi(3))
297
298 n(1:3,4:6)=zero
299
300 vel(1,2)=half*(v(2,nc(1)) + v(2,nc(2)))
301 vel(1,3)=half*(v(3,nc(1)) + v(3,nc(2)))
302
303 vel(2,2)=half*(v(2,nc(2)) + v(2,nc(3)))
304 vel(2,3)=half*(v(3,nc(2)) + v(3,nc(3)))
305
306 vel(3,2)=half*(v(2,nc(3)) + v(2,nc(1)))
307 vel(3,3)=half*(v(3,nc(3)) + v(3,nc(1)))
308
309 dotprod(1,i) = vel(1,2)*n(2,1) + vel(1,3)*n(
310 dotprod(2,i) = vel(2,2)*n(2,2) + vel(2,3)*n(3,2)
311 dotprod(3,i) = vel(3,2)*n(2,3) + vel(3,3)*n(3,3)
312 dotprod(4,i) = zero
313 dotprod(5,i) = zero
314 dotprod(6,i) = zero
315
316 IF(mlw /= 151)THEN
317 nx = half * ((yi(2) - yi(1)) * (zi(3) - zi(1)) - (zi(2) - zi(1)) * (yi(3) - yi(1)))
319 ELSE
320 area = elbuf_tab(ng)%GBUF%AREA(i)
321 ENDIF
323 ENDDO
324 ELSE
325 vol(1:nel) = one
326 ENDIF
327
328
329
330
331
332 DO i=1,nel
333 ie=i+nft
334
335
336
337 div = zero
338 div = div + dotprod(1,i)
339 div = div + dotprod(2,i)
340 div = div + dotprod(3,i)
341 div = div + dotprod(4,i)
342 div = div + dotprod(5,i)
343 div = div + dotprod(6,i)
344 div = div / vol(i)
345 evar(i) = div
346
347 enddo
348
subroutine area(d1, x, x2, y, y2, eint, stif0)