174 use element_mod , only : nixc
175
176
177
178#include "implicit_f.inc"
179
180
181
182#include "mvsiz_p.inc"
183
184
185
186 INTEGER JFT, JLT,NEL
187 INTEGER IXC(NIXC,*)
189 . x(3,*), offg(*), dr(3,*),
190 . e1x(*), e1y(*), e1z(*),
191 . e2x(*), e2y(*), e2z(*),e3x(*), e3y(*), e3z(*),
192 . smstr(*),
area(*),px1(*),px2(*),py1(*),py2(*),
193 . v13x(*),v24x(*),v13y(*),v24y(*)
194
195
196
197 INTEGER , II(9),NN(4)
199 . x0g2(mvsiz),x0g3(mvsiz),x0g4(mvsiz),y0g2(mvsiz),
200 . y0g3(mvsiz),y0g4(mvsiz),z0g2(mvsiz),z0g3(mvsiz),z0g4(mvsiz),
201 . off_l,vg13(3),vg24(3),
202 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),yl3(mvsiz),
203 . yl4(mvsiz),uxyz(mvsiz,3,4),a_i(mvsiz)
204
205 DO i=1,9
206 ii(i) = nel*(i-1)
207 ENDDO
208
209 DO i=jft,jlt
210 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
211 uxyz(i,1:3,1:4)= zero
212 nn(1)=ixc(2,i)
213 nn(2)=ixc(3,i)
214 nn(3)=ixc(4,i)
215 nn(4)=ixc(5,i)
216 x0g2(i) = x(1,nn(2))-x(1,nn(1))
217 y0g2(i) = x(2,nn(2))-x(2,nn(1))
218 z0g2(i) = x(3,nn(2))-x(3,nn(1))
219 x0g3(i) = x(1,nn(3))-x(1,nn(1))
220 y0g3(i) = x(2,nn(3))-x(2,nn(1))
221 z0g3(i) = x(3,nn(3))-x(3,nn(1))
222 x0g4(i) = x(1,nn(4))-x(1,nn(1))
223 y0g4(i) = x(2,nn(4))-x(2,nn(1))
224 z0g4(i) = x(3,nn(4))-x(3,nn(1))
225 IF(abs(offg(i))==two)THEN
226 uxyz(i,1,2) = x0g2(i)-smstr(ii(1)+i)
227 uxyz(i,2,2) = y0g2(i)-smstr(ii(2)+i)
228 uxyz(i,3,2) = z0g2(i)-smstr(ii(3)+i)
229 uxyz(i,1,3) = x0g3(i)-smstr(ii(4)+i)
230 uxyz(i,2,3) = y0g3(i)-smstr(ii(5)+i)
231 uxyz(i,3,3) = z0g3(i)-smstr(ii(6)+i)
232 uxyz(i,1,4) = x0g4(i)-smstr(ii(7)+i)
233 uxyz(i,2,4) = y0g4(i)-smstr(ii(8)+i)
234 uxyz(i,3,4) = z0g4(i)-smstr(ii(9)+i)
235
236 x0g2(i) = smstr(ii(1)+i)
237 y0g2(i) = smstr(ii(2)+i)
238 z0g2(i) = smstr(ii(3)+i)
239 x0g3(i) = smstr(ii(4)+i)
240 y0g3(i) = smstr(ii(5)+i)
241 z0g3(i) = smstr(ii(6)+i)
242 x0g4(i) = smstr(ii(7)+i)
243 y0g4(i) = smstr(ii(8)+i)
244 z0g4(i) = smstr(ii(9)+i)
245
246 ELSE
247 smstr(ii(1)+i)= x0g2(i)
248 smstr(ii(2)+i)= y0g2(i)
249 smstr(ii(3)+i)= z0g2(i)
250 smstr(ii(4)+i)= x0g3(i)
251 smstr(ii(5)+i)= y0g3(i)
252 smstr(ii(6)+i)= z0g3(i)
253 smstr(ii(7)+i)= x0g4(i)
254 smstr(ii(8)+i)= y0g4(i)
255 smstr(ii(9)+i)= z0g4(i)
256 ENDIF
257 ENDDO
258
259 DO i=jft,jlt
260 xl2(i)=e1x(i)*x0g2(i)+e1y(i)*y0g2(i)+e1z(i)*z0g2(i)
261 xl3(i)=e1x(i)*x0g3(i)+e1y(i)*y0g3(i)+e1z(i)*z0g3(i)
262 xl4(i)=e1x(i)*x0g4(i)+e1y(i)*y0g4(i)+e1z(i)*z0g4(i)
263 yl2(i)=e2x(i)*x0g2(i)+e2y(i)*y0g2(i)+e2z(i)*z0g2(i)
264 yl3(i)=e2x(i)*x0g3(i)+e2y(i)*y0g3(i)+e2z(i)*z0g3(i)
265 yl4(i)=e2x(i)*x0g4(i)+e2y(i)*y0g4(i)+e2z(i)*z0g4(i)
266 ENDDO
267 DO i=jft,jlt
268 px1(i)= half*(yl2(i)-yl4(i))
269 py1(i)= half*(xl4(i)-xl2(i))
270 px2(i)= half* yl3(i)
271 py2(i)=-half* xl3(i)
272 area(i)=
max(two*(py2(i)*px1(i)-py1(i)*px2(i)),em20)
273 a_i(i) = one /
area(i)
274 ENDDO
275
276 DO i=jft,jlt
277 vg13(1)=uxyz(i,1,1)-uxyz(i,1,3)
278 vg24(1)=uxyz(i,1,2)-uxyz(i,1,4)
279 vg13(2)=uxyz(i,2,1)-uxyz(i,2,3)
280 vg24(2)=uxyz(i,2,2)-uxyz(i,2,4)
281 vg13(3)=uxyz(i,3,1)-uxyz(i,3,3)
282 vg24(3)=uxyz(i,3,2)-uxyz(i,3,4)
283
284 v13x(i)=e1x(i)*vg13(1)+e1y(i)*vg13(2)+e1z(i)*vg13(3)
285 v24x(i)=e1x(i)*vg24(1)+e1y(i)*vg24(2)+e1z(i)*vg24(3)
286 v13y(i)=e2x(i)*vg13(1)+e2y(i)*vg13(2)+e2z(i)*vg13(3)
287 v24y(i)=e2x(i)*vg24(1)+e2y(i)*vg24(2)+e2z(i)*vg24(3)
288
289 px1(i)=a_i(i)*px1(i)
290 py1(i)=a_i(i)*py1(i)
291 px2(i)=a_i(i)*px2(i)
292 py2(i)=a_i(i)*py2(i)
293 ENDDO
294
295 off_l = zero
296 DO i=jft,jlt
297 off_l =
min(off_l,offg(i))
298 ENDDO
299 IF (off_l < zero) THEN
300 DO i=jft,jlt
301 IF (offg(i) < zero) THEN
302 v13x(i)=zero
303 v24x(i)=zero
304 v13y(i)=zero
305 v24y(i)=zero
306 ENDIF
307 ENDDO
308 ENDIF
309
310 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)