189
190
191
192#include "implicit_f.inc"
193
194
195
196#include "mvsiz_p.inc"
197
198
199
200 INTEGER JFT, JLT,ISROT,NEL
201 INTEGER IXTG(NIXTG,*)
203 . x(3,*), offg(*), dr(3,*),
204 . e1x(*), e1y(*), e1z(*),
205 . e2x(*), e2y(*), e2z(*),e3x(*), e3y(*), e3z(*),
206 . xl2(*),xl3(*),yl2(*),yl3(*),
area(*),
207 . v21x(*),v31x(*),v21y(*),v31y(*),rz13(*),rz23(*),
208 . x2_t(*),x3_t(*),y2_t(*),y3_t(*)
209 DOUBLE PRECISION
210 . SMSTR(*)
211
212
213
214 INTEGER I, J, NC1, NC2, NC3,II(6),NN(3)
216 . x0g2(mvsiz),x0g3(mvsiz),y0g2(mvsiz),y0g3(mvsiz),off_l,
217 . z0g2(mvsiz),z0g3(mvsiz),axyz(mvsiz,3,3),
218 . e01x(mvsiz), e01y(mvsiz), e01z(mvsiz),
219 . e02x(mvsiz), e02y(mvsiz), e02z(mvsiz),e03x(mvsiz),
220 . e03y(mvsiz), e03z(mvsiz),x0l2(mvsiz), x0l3(mvsiz),
221 . y0l2(mvsiz),y0l3(mvsiz),sum(mvsiz),
norm,xl,yl,vr1_12,vr1_21,
222 . rlz1,rlz2,rlz3,areai,x0g32,y0g32,z0g32,dirz(mvsiz,2)
223
224 DO i=1,6
225 ii(i) = nel*(i-1)
226 ENDDO
227
228 DO i=jft,jlt
229 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
230 axyz(i,1:3,1:3)= zero
231
232 IF (isrot > 0 ) THEN
233 nn(1) = ixtg(2,i)
234 nn(2) = ixtg(3,i)
235 nn(3) = ixtg(4,i)
236 axyz(i,1,1) = dr(1,nn(1))
237 axyz(i,2,1) = dr(2,nn(1))
238 axyz(i,3,1) = dr(3,nn(1))
239 axyz(i,1,2) = dr(1,nn(2))
240 axyz(i,2,2) = dr(2,nn(2))
241 axyz(i,3,2) = dr(3,nn(2))
242 axyz(i,1,3) = dr(1,nn(3))
243 axyz(i,2,3) = dr(2,nn(3))
244 axyz(i,3,3) = dr(3,nn(3))
245 END IF
246
247 x0g2(i) = smstr(ii(1)+i)
248 y0g2(i) = smstr(ii(2)+i)
249 z0g2(i) = smstr(ii(3)+i)
250 x0g3(i) = smstr(ii(4)+i)
251 y0g3(i) = smstr(ii(5)+i)
252 z0g3(i) = smstr(ii(6)+i)
253 ENDDO
254
255 DO i=jft,jlt
256 e01x(i)= x0g2(i)
257 e01y(i)= y0g2(i)
258 e01z(i)= z0g2(i)
259 sum(i) = sqrt(e01x(i)*e01x(i)+e01y(i)*e01y(i)+e01z(i)*e01z(i))
260 e01x(i)=e01x(i)/sum(i)
261 e01y(i)=e01y(i)/sum(i)
262 e01z(i)=e01z(i)/sum(i)
263 ENDDO
264
265 DO i=jft,jlt
266 x0g32=x0g3(i)-x0g2(i)
267 y0g32=y0g3(i)-y0g2(i)
268 z0g32=z0g3(i)-z0g2(i)
269 e03x(i)=y0g3(i)*z0g32-z0g3(i)*y0g32
270 e03y(i)=z0g3(i)*x0g32-x0g3(i)*z0g32
271 e03z(i)=x0g3(i)*y0g32-y0g3(i)*x0g32
272 sum(i) = sqrt(e03x(i)*e03x(i)+e03y(i)*e03y(i)+e03z(i)*e03z(i))
273 e03x(i)=e03x(i)/sum(i)
274 e03y(i)=e03y(i)/sum(i)
275 e03z(i)=e03z(i)/sum(i)
276 area(i) = half * sum(i)
277 ENDDO
278
279 DO i=jft,jlt
280 e02x(i)=e03y(i)*e01z(i)-e03z(i)*e01y(i)
281 e02y(i)=e03z(i)*e01x(i)-e03x(i)*e01z(i)
282 e02z(i)=e03x(i)*e01y(i)-e03y(i)*e01x(i)
283 sum(i) = sqrt(e02x(i)*e02x(i)+e02y(i)*e02y(i)+e02z(i)*e02z(i))
284 e02x(i)=e02x(i)/sum(i)
285 e02y(i)=e02y(i)/sum(i)
286 e02z(i)=e02z(i)/sum(i)
287 ENDDO
288
289
290
291 DO i=jft,jlt
292 vr1_12=e01x(i)*e2x(i)+e01y(i)*e2y(i)+e01z(i)*e2z(i
293 vr1_21=e02x(i)*e1x(i)+e02y(i)*e1y(i)+e02z(i)*e1z(i)
294 dirz(i,2) = half*(vr1_12-vr1_21)
295 norm = one-dirz(i,2)*dirz(i,2)
296 dirz(i,1) = sqrt(
max(zero,
norm))
297 ENDDO
298 DO i=jft,jlt
299 x0l2(i)=e01x(i)*x0g2(i)+e01y(i)*y0g2(i)+e01z(i)*z0g2(i)
300 y0l2(i)=e02x(i)*x0g2(i)+e02y(i)*y0g2(i)+e02z(i)*z0g2(i)
301 x0l3(i)=e01x(i)*x0g3(i)+e01y(i)*y0g3(i)+e01z(i)*z0g3(i)
302 y0l3(i)=e02x(i)*x0g3(i)+e02y(i)*y0g3(i)+e02z(i)*z0g3(i)
303 ENDDO
304
305
306
307 DO i=jft,jlt
308 xl= x0l2(i)*dirz(i,1)-y0l2(i)*dirz(i,2)
309 yl= x0l2(i)*dirz(i,2)+y0l2(i)*dirz(i,1)
310 x0l2(i)=xl
311 y0l2(i)=yl
312 xl= x0l3(i)*dirz(i,1)-y0l3(i)*dirz(i,2)
313 yl= x0l3(i)*dirz(i,2)+y0l3(i)*dirz(i,1)
314 x0l3(i)=xl
315 y0l3(i)=yl
316 ENDDO
317
318 DO i=jft,jlt
319 v21x(i)=xl2(i)-x0l2(i)
320 v31x(i)=xl3(i)-x0l3(i)
321 v21y(i)=yl2(i)-y0l2(i)
322 v31y(i)=yl3(i)-y0l3(i)
323 ENDDO
324 DO i=jft,jlt
325 x2_t(i) = x0l2(i)
326 x3_t(i) = x0l3(i)
327 y2_t(i) = y0l2(i)
328 y3_t(i) = y0l3(i)
329 ENDDO
330 IF (isrot>0) THEN
331
332 DO i=jft,jlt
334 rlz1 =e3x(i)*axyz(i,1,1)+e3y(i)*axyz(i,2,1)+e3z(i)*axyz(i,3,1)
335 rlz2 =e3x(i)*axyz(i,1,2)+e3y(i)*axyz(i,2,2)+e3z(i)*axyz(i,3,2)
336 rlz3 =e3x(i)*axyz(i,1,3)+e3y(i)*axyz(i,2,3)+e3z(i)*axyz(i,3,3)
337 rz13(i)=(rlz1-rlz3)*areai
338 rz23(i)=(rlz2-rlz3)*areai
339 ENDDO
340 END IF
341
342 off_l = zero
343 DO i=jft,jlt
344 off_l =
min(off_l,offg(i))
345 ENDDO
346 IF (off_l < zero) THEN
347 DO i=jft,jlt
348 IF (offg(i) < zero) THEN
349 v21x(i) = zero
350 v31x(i) = zero
351 v21y(i) = zero
352 v31y(i) = zero
353 rz13(i) = zero
354 rz23(i) = zero
355 ENDIF
356 ENDDO
357 ENDIF
358
359 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
subroutine area(d1, x, x2, y, y2, eint, stif0)