30
31
32
33#include "implicit_f.inc"
34
35
36
37#include "mvsiz_p.inc"
38
39
40
41 INTEGER LFT,LLT,NEL
42
44 . dd(3,3,*),hh(2,*),sig(nel,6),dd1(3,3,*),gt(3,3,*),ht(*)
45
46
47
48 INTEGER I,J,K
49
51 . ss(3,mvsiz),st(3,mvsiz),p,sigy2,g(mvsiz),
52 . ss2,ds2,st2,
norm(mvsiz),normi,g2,tt,tv
54 . ew(3,mvsiz),tol,lamda(mvsiz),ev(3,3),kk(3,3)
55
56 tol=em5
57
58 DO i=lft,llt
59 g(i) = hh(2,i)
60 g2 = two*g(i)
61 tv=hh(1,i)
62 tt=tv+g2
63 dd(1,1,i)=tt
64 dd(2,2,i)=tt
65 dd(3,3,i)=tt
66 dd(1,2,i)=tv
67 dd(1,3,i)=tv
68 dd(2,3,i)=tv
69 dd(2,1,i)=tv
70 dd(3,1,i)=tv
71 dd(3,2,i)=tv
72 DO j=1,3
73 gt(j,j,i)=g2
74 ENDDO
75 ENDDO
76
77 DO i=lft,llt
78 IF (ht(i)>zero) THEN
79 p =-third*(sig(i,1)+sig(i,2)+sig(i,3))
80 ss(1,i)=sig(i,1)+p
81 ss(2,i)=sig(i,2)+p
82 ss(3,i)=sig(i,3)+p
83 st(1,i)=sig(i,4)
84 st(2,i)=sig(i,5)
85 st(3,i)=sig(i,6)
86 ENDIF
87 ENDDO
88
89 DO i=lft,llt
90 IF (ht(i)>zero) THEN
91 ss2=(ss(1,i)**2+ss(2,i)**2+ss(3,i)**2)*half
92 st2=st(1,i)**2+st(2,i)**2+st(3,i)**2
93 sigy2=
max(em20,three*(ss2+st2))
94 ds2=three*g(i)+ht(i)
95 norm(i)=three*g(i)/sqrt(ds2*sigy2)
96 ENDIF
97 ENDDO
98
99 DO i=lft,llt
100 IF (ht(i)>zero) THEN
101 DO j=1,3
102 ss(j,i)=ss(j,i)*
norm(i)
103 st(j,i)=st(j,i)*
norm(i)
104 ENDDO
105 ENDIF
106 ENDDO
107
108 DO i=lft,llt
109 IF (ht(i)>zero) THEN
110 DO j=1,3
111 DO k=j,3
112 dd(j,k,i)=dd(j,k,i)-ss(j,i)*ss(k,i)
113 gt(j,k,i)=gt(j,k,i)*half-st(j,i)*st(k,i)
114 dd(k,j,i)=dd(j,k,i)
115 gt(k,j,i)=gt(j,k,i)
116 ENDDO
117 ENDDO
118 DO j=1,3
119 DO k=1,3
120 dd1(j,k,i)=-ss(j,i)*st(k,i)
121 ENDDO
122 ENDDO
123 ENDIF
124 ENDDO
125
126 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB