76
77
78
79#include "implicit_f.inc"
80#include "comlock.inc"
81
82
83
84#include "mvsiz_p.inc"
85
86
87
88#include "units_c.inc"
89#include "scr17_c.inc"
90#include "impl1_c.inc"
91
92
93
94 INTEGER, INTENT(IN) :: NEL
95 INTEGER, INTENT(IN) :: MTN
96 INTEGER, INTENT(IN) :: ISMSTR
97 INTEGER, INTENT(IN) :: JHBE
98 INTEGER ICP,IDEG(*)
99
101 . off(*),offg(*),vol(*),ksi,eta,zeta,wi,
102 . pxc1(*), pxc2(*), pxc3(*), pxc4(*),
103 . pyc1(*), pyc2(*), pyc3(*), pyc4(*),
104 . pzc1(*), pzc2(*), pzc3(*), pzc4(*),
105 . px1h1(*), px1h2(*), px1h3(*),px1h4(*),
106 . px2h1(*), px2h2(*), px2h3(*),px2h4(*),
107 . px3h1(*), px3h2(*), px3h3(*),px3h4(*),
108 . px4h1(*), px4h2(*), px4h3(*),px4h4(*),
109 . hx(mvsiz,4), hy(mvsiz,4), hz(mvsiz,4),
110 . px1(*), px2(*), px3(*), px4(*),
111 . px5(*), px6(*), px7(*), px8(*),
112 . py1(*), py2(*), py3(*), py4(*),
113 . py5(*), py6(*), py7(*), py8(*),
114 . pz1(*), pz2(*), pz3(*), pz4(*),
115 . pz5(*), pz6(*), pz7(*), pz8(*),
116 . pxy1(*),pxy2(*),pxy3(*),pxy4(*),
117 . pxy5(*),pxy6(*),pxy7(*),pxy8(*),
118 . pyx1(*),pyx2(*),pyx3(*),pyx4(*),
119 . pyx5(*),pyx6(*),pyx7(*),pyx8(*),
120 . pxz1(*),pxz2(*),pxz3(*),pxz4(*),
121 . pxz5(*),pxz6(*),pxz7(*),pxz8(*),
122 . pzx1(*),pzx2(*),pzx3(*),pzx4(*),
123 . pzx5(*),pzx6(*),pzx7(*),pzx8(*),
124 . pyz1(*),pyz2(*),pyz3(*),pyz4(*),
125 . pyz5(*),pyz6(*),pyz7(*),pyz8(*),
126 . pzy1(*),pzy2(*),pzy3(*),pzy4(*),
127 . pzy5(*),pzy6(*),pzy7(*),pzy8(*),
128 . bxy1(*),bxy2(*),bxy3(*),bxy4(*),
129 . bxy5(*),bxy6(*),bxy7(*),bxy8(*),
130 . byx1(*),byx2(*),byx3(*),byx4(*),
131 . byx5(*),byx6(*),byx7(*),byx8(*),
132 . bxz1(*),bxz2(*),bxz3(*),bxz4(*),
133 . bxz5(*),bxz6(*),bxz7(*),bxz8(*),
134 . bzx1(*),bzx2(*),bzx3(*),bzx4(*),
135 . bzx5(*),bzx6(*),bzx7(*),bzx8(*),
136 . byz1(*),byz2(*),byz3(*),byz4(*),
137 . byz5(*),byz6(*),byz7(*),byz8(*),
138 . bzy1(*),bzy2(*),bzy3(*),bzy4(*),
139 . bzy5(*),bzy6(*),bzy7(*),bzy8(*),
140 . cj1(*),cj2(*),cj3(*),
141 . cj4(*),cj5(*),cj6(*),
142 . cj7(*),cj8(*),cj9(*),
143 . jac4,jac5,jac6,
144 . jac7,jac8,jac9,
145 . smax(*),deltax(*),nu(*)
146 DOUBLE PRECISION
147 . VOLDP(MVSIZ),DETDP
148
149
150
151 INTEGER NGL(*), I, J ,ICOR,ep
152
153
155 . det(mvsiz) ,dett , nu1(mvsiz),
156 . jac1,jac2,jac3,
157 . jaci1, jaci2, jaci3,
158 . jaci4, jaci5, jaci6,
159 . jaci7, jaci8, jaci9,
160 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),
161 . jac_38_29(mvsiz), jac_19_37(mvsiz), jac_27_18(mvsiz),
162 . jac_26_35(mvsiz), jac_34_16(mvsiz), jac_15_24(mvsiz),
163 . jaci12(mvsiz), jaci45(mvsiz), jaci78(mvsiz),
164 . d1x(mvsiz) , d2x(mvsiz) , d3x(mvsiz) , d4x(mvsiz) ,
165 . d1y(mvsiz) , d2y(mvsiz) , d3y(mvsiz) , d4y(mvsiz) ,
166 . d1z(mvsiz) , d2z(mvsiz) , d3z(mvsiz) , d4z(mvsiz) ,
167 . xg1(mvsiz), xg2(mvsiz), xg3(mvsiz), xg4(mvsiz),
168 . yg1(mvsiz), yg2(mvsiz), yg3(mvsiz), yg4(mvsiz),
169 . zg1(mvsiz), zg2(mvsiz), zg3(mvsiz), zg4(mvsiz),
170 . f1,f2,f3,xs,xas,ys,yas,zs,zas,cs,cas ,bxhi,byhi,bzhi
171
172
173
174
175
176
177 DO i=1,nel
178 jac1=cj1(i)+hx(i,3)*eta+(hx(i,2)+hx(i,4)*eta)*zeta
179 jac2=cj2(i)+hy(i,3)*eta+(hy(i,2)+hy(i,4)*eta)*zeta
180 jac3=cj3(i)+hz(i,3)*eta+(hz(i,2)+hz(i,4)*eta)*zeta
181
182 jac4=cj4(i)+hx(i,1)*zeta+(hx(i,3)+hx(i,4)*zeta)*ksi
183 jac5=cj5(i)+hy(i,1)*zeta+(hy(i,3)+hy(i,4)*zeta)*ksi
184 jac6=cj6(i)+hz(i,1)*zeta+(hz(i,3)+hz(i,4)*zeta)*ksi
185
186 jac7=cj7(i)+hx(i,2)*ksi+(hx(i,1)+hx(i,4)*ksi)*eta
187 jac8=cj8(i)+hy(i,2)*ksi+(hy(i,1)+hy(i,4)*ksi)*eta
188 jac9=cj9(i)+hz(i,2)*ksi+(hz(i,1)+hz(i,4)*ksi)*eta
189
190 jac_59_68(i)=jac5*jac9-jac6*jac8
191 jac_67_49(i)=jac6*jac7-jac4*jac9
192 jac_38_29(i)=(-jac2*jac9+jac3*jac8)
193 jac_19_37(i)=( jac1*jac9-jac3*jac7)
194 jac_27_18(i)=(-jac1*jac8+jac2*jac7)
195 jac_26_35(i)=( jac2*jac6-jac3*jac5)
196 jac_34_16(i)=(-jac1*jac6+jac3*jac4)
197 jac_15_24(i)=( jac1*jac5-jac2*jac4)
198 jac_48_57(i)=jac4*jac8-jac5*jac7
199
200 detdp=one_over_512*(jac1*jac_59_68(i)+jac2*jac_67_49(i)+jac3*jac_48_57(i))
201 det(i) = detdp
202 voldp(i) = wi * detdp
203 vol(i) = voldp(i)
204 ENDDO
205
206 icor = 0
207 DO i=1,nel
208 off(i)=offg(i)
209 IF(off(i)==zero)THEN
210 det(i)=one
211 IF (vol(i)<=zero) vol(i)=one
212 ELSEIF (vol(i)<=zero ) THEN
213 vol(i)= em20
214 off(i) =zero
215 icor=1
216 ENDIF
217 ENDDO
218
219 IF (icor>0.AND.impl_s>0) THEN
220 DO i=1,nel
221 IF(vol(i)<=zero.AND.off(i)/=zero)THEN
222 vol(i)= em20
223 off(i) =zero
224 IF (imp_chk>0) THEN
225#include "lockon.inc"
226 WRITE(iout ,2001) ngl(i)
227#include "lockoff.inc"
228 idel7nok = 1
229 imp_ir = imp_ir + 1
230 ELSEIF (imconv==1) THEN
231#include "lockon.inc"
232 WRITE(istdo,2000) ngl(i)
233 WRITE(iout ,2000) ngl(i)
234#include "lockoff.inc"
235 idel7nok = 1
236 ENDIF
237 ENDIF
238 ENDDO
239 END IF
240
241 f1 = eta*zeta
242 f2 = ksi*zeta
243 f3 = ksi*eta
244 DO i=1,nel
245 dett=one_over_512/det(i)
246 jaci1=dett*jac_59_68(i)
247 jaci4=dett*jac_67_49(i)
248 jaci7=dett*jac_48_57(i)
249 jaci2=dett*jac_38_29(i)
250 jaci5=dett*jac_19_37(i)
251 jaci8=dett*jac_27_18(i)
252 jaci3=dett*jac_26_35(i)
253 jaci6=dett*jac_34_16(i)
254 jaci9=dett*jac_15_24(i)
255
256 nu1(i) = nu(i)/(one - nu(i))
257
258 d1x(i)=jaci3*eta +jaci2*zeta
259 d2x(i)=jaci1*zeta +jaci3*ksi
260 d3x(i)=jaci2*ksi +jaci1*eta
261 d4x(i)=jaci1*f1+jaci2*f2+jaci3*f3
262
263 d1y(i)=jaci6*eta +jaci5*zeta
264 d2y(i)=jaci4*zeta +jaci6*ksi
265 d3y(i)=jaci5*ksi +jaci4*eta
266 d4y(i)=jaci4*f1+jaci5*f2+jaci6*f3
267
268 d1z(i)=jaci9*eta +jaci8*zeta
269 d2z(i)=jaci7*zeta +jaci9*ksi
270 d3z(i)=jaci8*ksi +jaci7*eta
271 d4z(i)=jaci7*f1+jaci8*f2+jaci9*f3
272 ENDDO
273
274 DO i=1,nel
275 xg1(i) = px1h1(i)*d1x(i)
276 xg2(i) = px1h2(i)*d2x(i)
277 xg3(i) = px1h3(i)*d3x(i)
278 xg4(i) = px1h4(i)*d4x(i)
279 yg1(i) = px1h1(i)*d1y(i)
280 yg2(i) = px1h2(i)*d2y(i)
281 yg3(i) = px1h3(i)*d3y(i)
282 yg4(i) = px1h4(i)*d4y(i)
283 zg1(i) = px1h1(i)*d1z(i)
284 zg2(i) = px1h2(i)*d2z(i)
285 zg3(i) = px1h3(i)*d3z(i)
286 zg4(i) = px1h4(i)*d4z(i)
287 ENDDO
288
289 DO i=1,nel
290 xs =d1x(i)+d2x(i)+d3x(i)
291 xas =d4x(i)+xg1(i)+xg2(i)+xg3(i)+xg4(i)-pxc1(i)
292 px1(i)=xs-xas
293 px7(i)=xs+xas
294 ys =d1y(i)+d2y(i)+d3y(i)
295 yas =d4y(i)+yg1(i)+yg2(i)+yg3(i)+yg4(i)-pyc1(i)
296 py1(i)=ys-yas
297 py7(i)=ys+yas
298 zs =d1z(i)+d2z(i)+d3z(i)
299 zas =d4z(i)+zg1(i)+zg2(i)+zg3(i)+zg4(i)-pzc1(i)
300 pz1(i)=zs-zas
301 pz7(i)=zs+zas
302
303 ys =d1y(i)+d2y(i)
304 yas =yg1(i)+yg2(i)-pyc1(i)
305 pxy1(i)=ys-yas
306 pxy7(i)=ys+yas
307 xs =d1x(i)+d2x(i)
308 xas =xg1(i)+xg2(i)-pxc1(i)
309 pyx1(i)=xs-xas
310 pyx7(i)=xs+xas
311 zs =d1z(i)+d3z(i)
312 zas =zg1(i)+zg3(i)-pzc1(i)
313 pxz1(i)=zs-zas
314 pxz7(i)=zs+zas
315 xs =d1x(i)+d3x(i)
316 xas =xg1(i)+xg3(i)-pxc1(i)
317 pzx1(i)=xs-xas
318 pzx7(i)=xs+xas
319 zs =d2z(i)+d3z(i)
320 zas =zg2(i)+zg3(i)-pzc1(i)
321 pyz1(i)=zs-zas
322 pyz7(i)=zs+zas
323 ys =d2y(i)+d3y(i)
324 yas =yg2(i)+yg3(i)-pyc1(i)
325 pzy1(i)=ys-yas
326 pzy7(i)=ys+yas
327 ENDDO
328
329 IF (icp /= 1 .AND. icp /= 11) THEN
330 DO i=1,nel
331 cas =-nu(i)*(d4x(i)+xg1(i)+xg4(i))
332 xas =-nu1(i)*xg3(i)+cas
333 bxy1(i)=-xas
334 bxy7(i)=xas
335 xas =-nu1(i)*xg2(i)+cas
336 bxz1(i)=-xas
337 bxz7(i)=xas
338 cas =-nu(i)*(d4y(i)+yg2(i)+xg4(i))
339
340
341 yas =-nu1(i)*yg3(i)+cas
342 byx1(i)=-yas
343 byx7(i)=yas
344 yas =-nu1(i)*yg1(i)+cas
345 byz1(i)=-yas
346 byz7(i)=yas
347 cas =-nu(i)*(d4z(i)+zg3(i)+zg4(i))
348 zas =-nu1(i)*zg2(i)+cas
349 bzx1(i)=-zas
350 bzx7(i)=zas
351 zas =-nu1(i)*zg1(i)+cas
352 bzy1(i)=-zas
353 bzy7(i)=zas
354 ENDDO
355 ENDIF
356
357 DO i=1,nel
358 xg1(i) = px2h1(i)*d1x(i)
359 xg2(i) = px2h2(i)*d2x(i)
360 xg3(i) = px2h3(i)*d3x(i)
361 xg4(i) = px2h4(i)*d4x(i)
362 yg1(i) = px2h1(i)*d1y(i)
363 yg2(i) = px2h2(i)*d2y(i)
364 yg3(i) = px2h3(i)*d3y(i)
365 yg4(i) = px2h4(i)*d4y(i)
366 zg1(i) = px2h1(i)*d1z(i)
367 zg2(i) = px2h2(i)*d2z(i)
368 zg3(i) = px2h3(i)*d3z(i)
369 zg4(i) = px2h4(i)*d4z(i)
370 ENDDO
371
372 DO i=1,nel
373 xs =d1x(i)-d2x(i)-d3x(i)
374 xas =-d4x(i)+xg1(i)+xg2(i)+xg3(i)+xg4(i)-pxc2(i)
375 px2(i)=xs-xas
376 px8(i)=xs+xas
377 ys =d1y(i)-d2y(i)-d3y(i)
378 yas =-d4y(i)+yg1(i)+yg2(i)+yg3(i)+yg4(i)-pyc2(i)
379 py2(i)=ys-yas
380 py8(i)=ys+yas
381 zs =d1z(i)-d2z(i)-d3z(i)
382 zas =-d4z(i)+zg1(i)+zg2(i)+zg3(i)+zg4(i)-pzc2(i)
383 pz2(i)=zs-zas
384 pz8(i)=zs+zas
385
386 ys =d1y(i)-d2y(i)
387 yas =yg1(i)+yg2(i)-pyc2(i)
388 pxy2(i)=ys-yas
389 pxy8(i)=ys+yas
390 xs =d1x(i)-d2x(i)
391 xas =xg1(i)+xg2(i)-pxc2(i)
392 pyx2(i)=xs-xas
393 pyx8(i)=xs+xas
394 zs =d1z(i)-d3z(i)
395 zas =zg1(i)+zg3(i)-pzc2(i)
396 pxz2(i)=zs-zas
397 pxz8(i)=zs+zas
398 xs =d1x(i)-d3x(i)
399 xas =xg1(i)+xg3(i)-pxc2(i)
400 pzx2(i)=xs-xas
401 pzx8(i)=xs+xas
402 zs =-d2z(i)-d3z(i)
403 zas =zg2(i)+zg3(i)-pzc2(i)
404 pyz2(i)=zs-zas
405 pyz8(i)=zs+zas
406 ys =-d2y(i)-d3y(i)
407 yas =yg2(i)+yg3(i)-pyc2(i)
408 pzy2(i)=ys-yas
409 pzy8(i)=ys+yas
410 ENDDO
411
412 IF (icp /= 1 .AND. icp /= 11) THEN
413 DO i=1,nel
414 cas =-nu(i)*(-d4x(i)+xg1(i)+xg4(i))
415 xas =-nu1(i)*xg3(i)+cas
416 bxy2(i)=-xas
417 bxy8(i)=xas
418 xas =-nu1(i)*xg2(i)+cas
419 bxz2(i)=-xas
420 bxz8(i)=xas
421 cas =-nu(i)*(-d4y(i)+yg2(i)+xg4(i))
422 yas =-nu1(i)*yg3(i)+cas
423 byx2(i)=-yas
424 byx8(i)=yas
425 yas =-nu1(i)*yg1(i)+cas
426 byz2(i)=-yas
427 byz8(i)=yas
428 cas =-nu(i)*(-d4z(i)+zg3(i)+zg4(i))
429 zas =-nu1(i)*zg2(i)+cas
430 bzx2(i)=-zas
431 bzx8(i)=zas
432 zas =-nu1(i)*zg1(i)+cas
433 bzy2(i)=-zas
434 bzy8(i)=zas
435 ENDDO
436 ENDIF
437
438
439 DO i=1,nel
440 xg1(i) = px3h1(i)*d1x(i)
441 xg2(i) = px3h2(i)*d2x(i)
442 xg3(i) = px3h3(i)*d3x(i)
443 xg4(i) = px3h4(i)*d4x(i)
444 yg1(i) = px3h1(i)*d1y(i)
445 yg2(i) = px3h2(i)*d2y(i)
446 yg3(i) = px3h3(i)*d3y(i)
447 yg4(i) = px3h4(i)*d4y(i)
448 zg1(i) = px3h1(i)*d1z(i)
449 zg2(i) = px3h2(i)*d2z(i)
450 zg3(i) = px3h3(i)*d3z(i)
451 zg4(i) = px3h4(i)*d4z(i)
452 ENDDO
453
454 DO i=1,nel
455 xs =-d1x(i)-d2x(i)+d3x(i)
456 xas =d4x(i)+xg1(i)+xg2(i)+xg3(i)+xg4(i)-pxc3(i)
457 px3(i)=xs-xas
458 px5(i)=xs+xas
459 ys =-d1y(i)-d2y(i)+d3y(i)
460 yas =d4y(i)+yg1(i)+yg2(i)+yg3(i)+yg4(i)-pyc3(i)
461 py3(i)=ys-yas
462 py5(i)=ys+yas
463 zs =-d1z(i)-d2z(i)+d3z(i)
464 zas =d4z(i)+zg1(i)+zg2(i)+zg3(i)+zg4(i)-pzc3(i)
465 pz3(i)=zs-zas
466 pz5(i)=zs+zas
467
468 ys =-d1y(i)-d2y(i)
469 yas =yg1(i)+yg2(i)-pyc3(i)
470 pxy3(i)=ys-yas
471 pxy5(i)=ys+yas
472 xs =-d1x(i)-d2x(i)
473 xas =xg1(i)+xg2(i)-pxc3(i)
474 pyx3(i)=xs-xas
475 pyx5(i)=xs+xas
476 zs =-d1z(i)+d3z(i)
477 zas =zg1(i)+zg3(i)-pzc3(i)
478 pxz3(i)=zs-zas
479 pxz5(i)=zs+zas
480 xs =-d1x(i)+d3x(i)
481 xas =xg1(i)+xg3(i)-pxc3(i)
482 pzx3(i)=xs-xas
483 pzx5(i)=xs+xas
484 zs =-d2z(i)+d3z(i)
485 zas =zg2(i)+zg3(i)-pzc3(i)
486 pyz3(i)=zs-zas
487 pyz5(i)=zs+zas
488 ys =-d2y(i)+d3y(i)
489 yas =yg2(i)+yg3(i)-pyc3(i)
490 pzy3(i)=ys-yas
491 pzy5(i)=ys+yas
492 ENDDO
493
494 IF (icp /= 1 .AND. icp /= 11) THEN
495 DO i=1,nel
496 cas =-nu(i)*(d4x(i)+xg1(i)+xg4(i))
497 xas =-nu1(i)*xg3(i)+cas
498 bxy3(i)=-xas
499 bxy5(i)=xas
500 xas =-nu1(i)*xg2(i)+cas
501 bxz3(i)=-xas
502 bxz5(i)=xas
503 cas =-nu(i)*(d4y(i)+yg2(i)+xg4(i))
504 yas =-nu1(i)*yg3(i)+cas
505 byx3(i)=-yas
506 byx5(i)=yas
507 yas =-nu1(i)*yg1(i)+cas
508 byz3(i)=-yas
509 byz5(i)=yas
510 cas =-nu(i)*(d4z(i)+zg3(i)+zg4(i))
511 zas =-nu1(i)*zg2(i)+cas
512 bzx3(i)=-zas
513 bzx5(i)=zas
514 zas =-nu1(i)*zg1(i)+cas
515 bzy3(i)=-zas
516 bzy5(i)=zas
517 ENDDO
518 ENDIF
519
520 DO i=1,nel
521 xg1(i) = px4h1(i)*d1x(i)
522 xg2(i) = px4h2(i)*d2x(i)
523 xg3(i) = px4h3(i)*d3x(i)
524 xg4(i) = px4h4(i)*d4x(i)
525 yg1(i) = px4h1(i)*d1y(i)
526 yg2(i) = px4h2(i)*d2y(i)
527 yg3(i) = px4h3(i)*d3y(i)
528 yg4(i) = px4h4(i)*d4y(i)
529 zg1(i) = px4h1(i)*d1z(i)
530 zg2(i) = px4h2(i)*d2z(i)
531 zg3(i) = px4h3(i)*d3z(i)
532 zg4(i) = px4h4(i)*d4z(i)
533 ENDDO
534
535 DO i=1,nel
536 xs =-d1x(i)+d2x(i)-d3x(i)
537 xas =-d4x(i)+xg1(i)+xg2(i)+xg3(i)+xg4(i)-pxc4(i)
538 px4(i)=xs-xas
539 px6(i)=xs+xas
540 ys =-d1y(i)+d2y(i)-d3y(i)
541 yas =-d4y(i)+yg1(i)+yg2(i)+yg3(i)+yg4(i)-pyc4(i)
542 py4(i)=ys-yas
543 py6(i)=ys+yas
544 zs =-d1z(i)+d2z(i)-d3z(i)
545 zas =-d4z(i)+zg1(i)+zg2(i)+zg3(i)+zg4(i)-pzc4(i)
546 pz4(i)=zs-zas
547 pz6(i)=zs+zas
548
549 ys =-d1y(i)+d2y(i)
550 yas =yg1(i)+yg2(i)-pyc4(i)
551 pxy4(i)=ys-yas
552 pxy6(i)=ys+yas
553 xs =-d1x(i)+d2x(i)
554 xas =xg1(i)+xg2(i)-pxc4(i)
555 pyx4(i)=xs-xas
556 pyx6(i)=xs+xas
557 zs =-d1z(i)-d3z(i)
558 zas =zg1(i)+zg3(i)-pzc4(i)
559 pxz4(i)=zs-zas
560 pxz6(i)=zs+zas
561 xs =-d1x(i)-d3x(i)
562 xas =xg1(i)+xg3(i)-pxc4(i)
563 pzx4(i)=xs-xas
564 pzx6(i)=xs+xas
565 zs =d2z(i)-d3z(i)
566 zas =zg2(i)+zg3(i)-pzc4(i)
567 pyz4(i)=zs-zas
568 pyz6(i)=zs+zas
569 ys =d2y(i)-d3y(i)
570 yas =yg2(i)+yg3(i)-pyc4(i)
571 pzy4(i)=ys-yas
572 pzy6(i)=ys+yas
573 ENDDO
574
575 IF (icp /= 1 .AND. icp /= 11) THEN
576 DO i=1,nel
577 cas =-nu(i)*(-d4x(i)+xg1(i)+xg4(i))
578 xas =-nu1(i)*xg3(i)+cas
579 bxy4(i)=-xas
580 bxy6(i)=xas
581 xas =-nu1(i)*xg2(i)+cas
582 bxz4(i)=-xas
583 bxz6(i)=xas
584 cas =-nu(i)*(-d4y(i)+yg2(i)+xg4(i))
585 yas =-nu1(i)*yg3(i)+cas
586 byx4(i)=-yas
587 byx6(i)=yas
588 yas =-nu1(i)*yg1(i)+cas
589 byz4(i)=-yas
590 byz6(i)=yas
591 cas =-nu(i)*(-d4z(i)+zg3(i)+zg4(i))
592 zas =-nu1(i)*zg2(i)+cas
593 bzx4(i)=-zas
594 bzx6(i)=zas
595 zas =-nu1(i)*zg1(i)+cas
596 bzy4(i)=-zas
597 bzy6(i)=zas
598 ENDDO
599
600 DO i=1,nel
601 cs =-nu(i)*d1x(i)
602 xs =-nu1(i)*d3x(i)
603 bxy1(i)=bxy1(i)+cs+xs
604 bxy2(i)=bxy2(i)+cs-xs
605 bxy3(i)=bxy3(i)-cs+xs
606 bxy4(i)=bxy4(i)-cs-xs
607 bxy5(i)=bxy5(i)-cs+xs
608 bxy6(i)=bxy6(i)-cs-xs
609 bxy7(i)=bxy7(i)+cs+xs
610 bxy8(i)=bxy8(i)+cs-xs
611 xs =-nu1(i)*d2x(i)
612 bxz1(i)=bxz1(i)+cs+xs
613 bxz2(i)=bxz2(i)+cs-xs
614 bxz3(i)=bxz3(i)-cs-xs
615 bxz4(i)=bxz4(i)-cs+xs
616 bxz5(i)=bxz5(i)-cs-xs
617 bxz6(i)=bxz6(i)-cs+xs
618 bxz7(i)=bxz7(i)+cs+xs
619 bxz8(i)=bxz8(i)+cs-xs
620 ENDDO
621 DO i=1,nel
622 cs =-nu(i)*d2y(i)
623 ys =-nu1(i)*d3y(i)
624 byx1(i)=byx1(i)+cs+ys
625 byx2(i)=byx2(i)-cs-ys
626 byx3(i)=byx3(i)-cs+ys
627 byx4(i)=byx4(i)+cs-ys
628 byx5(i)=byx5(i)-cs+ys
629 byx6(i)=byx6(i)+cs-ys
630 byx7(i)=byx7(i)+cs+ys
631 byx8(i)=byx8(i)-cs-ys
632 ys =-nu1(i)*d1y(i)
633 byz1(i)=byz1(i)+cs+ys
634 byz2(i)=byz2(i)-cs+ys
635 byz3(i)=byz3(i)-cs-ys
636 byz4(i)=byz4(i)+cs-ys
637 byz5(i)=byz5(i)-cs-ys
638 byz6(i)=byz6(i)+cs-ys
639 byz7(i)=byz7(i)+cs+ys
640 byz8(i)=byz8(i)-cs+ys
641 ENDDO
642 DO i=1,nel
643 cs =-nu(i)*d3z(i)
644 zs =-nu1(i)*d2z(i)
645 bzx1(i)=bzx1(i)+cs+zs
646 bzx2(i)=bzx2(i)-cs-zs
647 bzx3(i)=bzx3(i)+cs-zs
648 bzx4(i)=bzx4(i)-cs+zs
649 bzx5(i)=bzx5(i)+cs-zs
650 bzx6(i)=bzx6(i)-cs+zs
651 bzx7(i)=bzx7(i)+cs+zs
652 bzx8(i)=bzx8(i)-cs-zs
653 zs =-nu1(i)*d1z(i)
654 bzy1(i)=bzy1(i)+cs+zs
655 bzy2(i)=bzy2(i)-cs+zs
656 bzy3(i)=bzy3(i)+cs-zs
657 bzy4(i)=bzy4(i)-cs-zs
658 bzy5(i)=bzy5(i)+cs-zs
659 bzy6(i)=bzy6(i)-cs-zs
660 bzy7(i)=bzy7(i)+cs+zs
661 bzy8(i)=bzy8(i)-cs+zs
662 ENDDO
663 ENDIF
664 IF (idts6==0) THEN
665 DO i=1,nel
666 deltax(i)=
min(deltax(i),hundred28*det(i)*smax(i))
667 ENDDO
668 ELSE
669 DO i=1,nel
670 IF (ideg(i)==0) deltax(i)=
min(deltax(i),hundred28*det(i)*smax(i))
671 ENDDO
672 END IF
673
674 IF ((impl_s >0.AND.jhbe /= 14.AND.jhbe < 20).OR.
675 . ((ismstr==10.OR.ismstr==12).AND.
676 . (icp==1.OR.icp==11.OR.mtn==99.OR.mtn==1).AND.jhbe==14)) THEN
677 DO i=1,nel
678 pxy1(i)=py1(i)
679 pxy2(i)=py2(i)
680 pxy3(i)=py3(i)
681 pxy4(i)=py4(i)
682 pxy5(i)=py5(i)
683 pxy6(i)=py6(i)
684 pxy7(i)=py7(i)
685 pxy8(i)=py8(i)
686
687 pyx1(i)=px1(i)
688 pyx2(i)=px2(i)
689 pyx3(i)=px3(i)
690 pyx4(i)=px4(i)
691 pyx5(i)=px5(i)
692 pyx6(i)=px6(i)
693 pyx7(i)=px7(i)
694 pyx8(i)=px8(i)
695 ENDDO
696 DO i=1,nel
697 pyz1(i)=pz1(i)
698 pyz2(i)=pz2(i)
699 pyz3(i)=pz3(i)
700 pyz4(i)=pz4(i)
701 pyz5(i)=pz5(i)
702 pyz6(i)=pz6(i)
703 pyz7(i)=pz7(i)
704 pyz8(i)=pz8(i)
705
706 pzy1(i)=py1(i)
707 pzy2(i)=py2(i)
708 pzy3(i)=py3(i)
709 pzy4(i)=py4(i)
710 pzy5(i)=py5(i)
711 pzy6(i)=py6(i)
712 pzy7(i)=py7(i)
713 pzy8(i)=py8(i)
714 ENDDO
715 DO i=1,nel
716 pxz1(i)=pz1(i)
717 pxz2(i)=pz2(i)
718 pxz3(i)=pz3(i)
719 pxz4(i)=pz4(i)
720 pxz5(i)=pz5(i)
721 pxz6(i)=pz6(i)
722 pxz7(i)=pz7(i)
723 pxz8(i)=pz8(i)
724
725 pzx1(i)=px1(i)
726 pzx2(i)=px2(i)
727 pzx3(i)=px3(i)
728 pzx4(i)=px4(i)
729 pzx5(i)=px5(i)
730 pzx6(i)=px6(i)
731 pzx7(i)=px7(i)
732 pzx8(i)=px8(i)
733 ENDDO
734 ENDIF
735
736 IF ((ismstr==10.OR.ismstr==12).AND.mtn==1.AND.jhbe==14.AND.
737 . (icp/=1.AND.icp/=11)) THEN
738 DO i=1,nel
739 bxy1(i)=zero
740 bxz1(i)=zero
741 byx1(i)=zero
742 byz1(i)=zero
743 bzx1(i)=zero
744 bzy1(i)=zero
745 bxy2(i)=zero
746 bxz2(i)=zero
747 byx2(i)=zero
748 byz2(i)=zero
749 bzx2(i)=zero
750 bzy2(i)=zero
751 bxy3(i)=zero
752 bxz3(i)=zero
753 byx3(i)=zero
754 byz3(i)=zero
755 bzx3(i)=zero
756 bzy3(i)=zero
757 bxy4(i)=zero
758 bxz4(i)=zero
759 byx4(i)=zero
760 byz4(i)=zero
761 bzx4(i)=zero
762 bzy4(i)=zero
763 bxy5(i)=zero
764 bxz5(i)=zero
765 byx5(i)=zero
766 byz5(i)=zero
767 bzx5(i)=zero
768 bzy5(i)=zero
769 bxy6(i)=zero
770 bxz6(i)=zero
771 byx6(i)=zero
772 byz6(i)=zero
773 bzx6(i)=zero
774 bzy6(i)=zero
775 bxy7(i)=zero
776 bxz7(i)=zero
777 byx7(i)=zero
778 byz7(i)=zero
779 bzx7(i)=zero
780 bzy7(i)=zero
781 bxy8(i)=zero
782 bxz8(i)=zero
783 byx8(i)=zero
784 byz8(i)=zero
785 bzx8(i)=zero
786 bzy8(i)=zero
787 ENDDO
788 ENDIF
789
790
791
792
793 IF (icp==11) THEN
794 DO i=1,nel
795 bxhi=third*(px1(i)-pxc1(i))
796 byhi=third*(py1(i)-pyc1(i))
797 bzhi=third*(pz1(i)-pzc1(i))
798 px1(i)=pxc1(i)+two*bxhi
799 py1(i)=pyc1(i)+two*byhi
800 pz1(i)=pzc1(i)+two*bzhi
801 bxy1(i)=-bxhi
802 bxz1(i)=-bxhi
803 byx1(i)=-byhi
804 byz1(i)=-byhi
805 bzx1(i)=-bzhi
806 bzy1(i)=-bzhi
807 ENDDO
808 DO i=1,nel
809 bxhi=third*(px2(i)-pxc2(i))
810 byhi=third*(py2(i)-pyc2(i))
811 bzhi=third*(pz2(i)-pzc2(i))
812 px2(i)=pxc2(i)+two*bxhi
813 py2(i)=pyc2(i)+two*byhi
814 pz2(i)=pzc2(i)+two*bzhi
815 bxy2(i)=-bxhi
816 bxz2(i)=-bxhi
817 byx2(i)=-byhi
818 byz2(i)=-byhi
819 bzx2(i)=-bzhi
820 bzy2(i)=-bzhi
821 ENDDO
822 DO i=1,nel
823 bxhi=third*(px3(i)-pxc3(i))
824 byhi=third*(py3(i)-pyc3(i))
825 bzhi=third*(pz3(i)-pzc3(i))
826 px3(i)=pxc3(i)+two*bxhi
827 py3(i)=pyc3(i)+two*byhi
828 pz3(i)=pzc3(i)+two*bzhi
829 bxy3(i)=-bxhi
830 bxz3(i)=-bxhi
831 byx3(i)=-byhi
832 byz3(i)=-byhi
833 bzx3(i)=-bzhi
834 bzy3(i)=-bzhi
835 ENDDO
836 DO i=1,nel
837 bxhi=third*(px4(i)-pxc4(i))
838 byhi=third*(py4(i)-pyc4(i))
839 bzhi=third*(pz4(i)-pzc4(i))
840 px4(i)=pxc4(i)+two*bxhi
841 py4(i)=pyc4(i)+two*byhi
842 pz4(i)=pzc4(i)+two*bzhi
843 bxy4(i)=-bxhi
844 bxz4(i)=-bxhi
845 byx4(i)=-byhi
846 byz4(i)=-byhi
847 bzx4(i)=-bzhi
848 bzy4(i)=-bzhi
849 ENDDO
850 DO i=1,nel
851 bxhi=third*(px5(i)+pxc3(i))
852 byhi=third*(py5(i)+pyc3(i))
853 bzhi=third*(pz5(i)+pzc3(i))
854 px5(i)=-pxc3(i)+two*bxhi
855 py5(i)=-pyc3(i)+two*byhi
856 pz5(i)=-pzc3(i)+two*bzhi
857 bxy5(i)=-bxhi
858 bxz5(i)=-bxhi
859 byx5(i)=-byhi
860 byz5(i)=-byhi
861 bzx5(i)=-bzhi
862 bzy5(i)=-bzhi
863 ENDDO
864 DO i=1,nel
865 bxhi=third*(px6(i)+pxc4(i))
866 byhi=third*(py6(i)+pyc4(i))
867 bzhi=third*(pz6(i)+pzc4(i))
868 px6(i)=-pxc4(i)+two*bxhi
869 py6(i)=-pyc4(i)+two*byhi
870 pz6(i)=-pzc4(i)+two*bzhi
871 bxy6(i)=-bxhi
872 bxz6(i)=-bxhi
873 byx6(i)=-byhi
874 byz6(i)=-byhi
875 bzx6(i)=-bzhi
876 bzy6(i)=-bzhi
877 ENDDO
878 DO i=1,nel
879 bxhi=third*(px7(i)+pxc1(i))
880 byhi=third*(py7(i)+pyc1(i))
881 bzhi=third*(pz7(i)+pzc1(i))
882 px7(i)=-pxc1(i)+two*bxhi
883 py7(i)=-pyc1(i)+two*byhi
884 pz7(i)=-pzc1(i)+two*bzhi
885 bxy7(i)=-bxhi
886 bxz7(i)=-bxhi
887 byx7(i)=-byhi
888 byz7(i)=-byhi
889 bzx7(i)=-bzhi
890 bzy7(i)=-bzhi
891 ENDDO
892 DO i=1,nel
893 bxhi=third*(px8(i)+pxc2(i))
894 byhi=third*(py8(i)+pyc2(i))
895 bzhi=third*(pz8(i)+pzc2(i))
896 px8(i)=-pxc2(i)+two*bxhi
897 py8(i)=-pyc2(i)+two*byhi
898 pz8(i)=-pzc2(i)+two*bzhi
899 bxy8(i)=
900 bxz8(i)=-bxhi
901 byx8(i)=-byhi
902 byz8(i)=-byhi
903 bzx8(i)=-bzhi
904 bzy8(i)=-bzhi
905 ENDDO
906 ENDIF
907 2000 FORMAT(/' ZERO OR NEGATIVE SUB-VOLUME : DELETE 3D-ELEMENT NB',
908 . i10/)
909 RETURN
910 2001 FORMAT(/' ZERO OR NEGATIVE SOLID SUB-VOLUME : ELEMENT NB:',
911 . i10/)