32
33
34
36
37
38
39#include "implicit_f.inc"
40
41
42
43#include "mvsiz_p.inc"
44
45
46
47#include "vect01_c.inc"
48#include "tabsiz_c.inc"
49
50
51
52
53
54
55
56
57
58
59
60 INTEGER,INTENT(IN) :: IXS(NIXS,SIXS/NIXS)
62 my_real,
INTENT(INOUT) :: grad(6,*)
63 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
64
65
66
67 INTEGER I, II, IE, IV1, IV2, IV3, IV4, IV5, IV6, IAD2
69 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
70 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz), z1(mvsiz),
71 . z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz), xc(mvsiz), yc(mvsiz),
72 . zc(mvsiz), n1x(mvsiz), n2x(mvsiz), n3x(mvsiz), n4x(mvsiz), n5x(mvsiz), n6x(mvsiz), n1y(mvsiz),
73 . n2y(mvsiz), n3y(mvsiz), n4y(mvsiz), n5y(mvsiz), n6y(mvsiz), n1z(mvsiz), n2z(mvsiz), n3z(mvsiz),
74 . n4z(mvsiz), n5z(mvsiz), n6z(mvsiz), dd1(mvsiz), dd2(mvsiz), dd3(mvsiz), dd4(mvsiz), dd5(mvsiz),
75 . dd6(mvsiz), d1x(mvsiz), d2x(mvsiz), d3x(mvsiz), d4x(mvsiz), d5x(mvsiz), d6x(mvsiz), d1y(mvsiz),
76 . d2y(mvsiz), d3y(mvsiz), d4y(mvsiz), d5y(mvsiz), d6y(mvsiz), d1z(mvsiz), d2z(mvsiz), d3z(mvsiz),
77 . d4z(mvsiz), d5z(mvsiz), d6z(mvsiz)
78
79
80
81 DO i=lft,llt
82 ii=i+nft
83 x1(i)=x(1,ixs(2,ii))
84 y1(i)=x(2,ixs(2,ii))
85 z1(i)=x(3,ixs(2,ii))
86
87 x2(i)=x(1,ixs(3,ii))
88 y2(i)=x(2,ixs(3,ii))
89 z2(i)=x(3,ixs(3,ii))
90
91 x3(i)=x(1,ixs(4,ii))
92 y3(i)=x(2,ixs(4,ii))
93 z3(i)=x(3,ixs(4,ii))
94
95 x4(i)=x(1,ixs(5,ii))
96 y4(i)=x(2,ixs(5,ii))
97 z4(i)=x(3,ixs(5,ii))
98
99 x5(i)=x(1,ixs(6,ii))
100 y5(i)=x(2,ixs(6,ii))
101 z5(i)=x(3,ixs(6,ii))
102
103 x6(i)=x(1,ixs(7,ii))
104 y6(i)=x(2,ixs(7,ii))
105 z6(i)=x(3,ixs(7,ii))
106
107 x7(i)=x(1,ixs(8,ii))
108 y7(i)=x(2,ixs(8,ii))
109 z7(i)=x(3,ixs(8,ii))
110
111 x8(i)=x(1,ixs(9,ii))
112 y8(i)=x(2,ixs(9,ii))
113 z8(i)=x(3,ixs(9,ii))
114 ENDDO
115
116
117
118 DO i=lft,llt
119 n1x(i) = (y3(i)-y1(i))*(z2(i)-z4(i)) - (z3(i)-z1(i))*(y2(i)-y4(i))
120 n1y(i) = (z3(i)-z1(i))*(x2(i)-x4(i)) - (x3(i)-x1(i))*(z2(i)-z4(i))
121 n1z(i) = (x3(i)-x1(i))*(y2(i)-y4(i)) - (y3(i)-y1(i))*(x2(i)-x4(i))
122
123 n2x(i) = (y7(i)-y4(i))*(z3(i)-z8(i)) - (z7(i)-z4(i))*(y3(i)-y8(i))
124 n2y(i) = (z7(i)-z4(i))*(x3(i)-x8(i)) - (x7(i)-x4(i))*(z3(i)-z8(i))
125 n2z(i) = (x7(i)-x4(i))*(y3(i)-y8(i)) - (y7(i)-y4(i))*(x3(i)-x8(i))
126
127 n3x(i) = (y6(i)-y8(i))*(z7(i)-z5(i)) - (z6(i)-z8(i))*(y7(i)-y5(i))
128 n3y(i) = (z6(i)-z8(i))*(x7(i)-x5(i)) - (x6(i)-x8(i))*(z7(i)-z5(i))
129 n3z(i) = (x6(i)-x8(i))*(y7(i)-y5(i)) - (y6(i)-y8(i))*(x7(i)-x5(i))
130
131 n4x(i) = (y2(i)-y5(i))*(z6(i)-z1(i)) - (z2(i)-z5(i))*(y6(i)-y1(i))
132 n4y(i) = (z2(i)-z5(i))*(x6(i)-x1(i)) - (x2(i)-x5(i))*(z6(i)-z1(i))
133 n4z(i) = (x2(i)-x5(i))*(y6(i)-y1(i)) - (y2(i)-y5(i))*(x6(i)-x1(i))
134
135 n5x(i) = (y7(i)-y2(i))*(z6(i)-z3(i)) - (z7(i)-z2(i))*(y6(i)-y3(i))
136 n5y(i) = (z7(i)-z2(i))*(x6(i)-x3(i)) - (x7(i)-x2(i))*(z6(i)-z3(i))
137 n5z(i) = (x7(i)-x2(i))*(y6(i)-y3(i)) - (y7(i)-y2(i))*(x6(i)-x3(i))
138
139 n6x(i) = (y8(i)-y1(i))*(z4(i)-z5(i)) - (z8(i)-z1(i))*(y4(i)-y5(i))
140 n6y(i) = (z8(i)-z1(i))*(x4(i)-x5(i)) - (x8(i)-x1(i))*(z4(i)-z5(i))
141 n6z(i) = (x8(i)-x1(i))*(y4(i)-y5(i)) - (y8(i)-y1(i))*(x4(i)-x5(i))
142
143 xc(i) = (x1(i)+x2(i)+x3(i)+x4(i)+x5(i)+x6(i)+x7(i)+x8(i))
144 yc(i) = (y1(i)+y2(i)+y3(i)+y4(i)+y5(i)+y6(i)+y7(i)+y8(i))
145 zc(i) = (z1(i)+z2(i)+z3(i)+z4(i)+z5(i)+z6(i)+z7(i)+z8(i))
146 ENDDO
147
148
149
150 DO i=lft,llt
151 ie =nft+i
152 iad2 = ale_connectivity%ee_connect%iad_connect(ie)
153 iv1=ale_connectivity%ee_connect%connected(iad2 + 1 - 1)
154 IF(iv1 <= 0) iv1=ie
155 d1x(i) = - xc(i)
156 . +x(1,ixs(2,iv1))+x(1,ixs(3,iv1))+x(1,ixs(4,iv1))+x(1,ixs(5,iv1))
157 . +x(1,ixs(6,iv1))+x(1,ixs(7,iv1))+x
158 d1y(i) = - yc(i)
159 . +x(2,ixs(2,iv1))+x(2,ixs(3,iv1))+x(2,ixs(4,iv1))+x(2,ixs(5,iv1))
160 . +x(2,ixs(6,iv1))+x(2,ixs(7,iv1))+x(2,ixs(8,iv1))+x(2,ixs(9,iv1))
161 d1z(i) = - zc(i)
162 . +x(3,ixs(2,iv1))+x(3,ixs(3,iv1))+x(3,ixs(4,iv1))+x(3,ixs(5,iv1))
163 . +x(3,ixs(6,iv1))+x(3,ixs(7,iv1))+x(3,ixs(8,iv1))+x(3,ixs(9,iv1))
164 ENDDO
165 DO i=lft,llt
166 ie =nft+i
167 iad2 = ale_connectivity%ee_connect%iad_connect(ie)
168 iv2=ale_connectivity%ee_connect%connected(iad2 + 2 - 1)
169 IF(iv2 <= 0) iv2=ie
170 d2x(i) = - xc(i)
171 . +x(1,ixs(2,iv2))+x(1,ixs(3,iv2))+x(1,ixs(4,iv2))+x(1,ixs(5,iv2))
172 . +x(1,ixs(6,iv2))+x(1,ixs(7,iv2))+x(1,ixs(8,iv2))+x(1,ixs(9,iv2))
173 d2y(i) = - yc(i)
174 . +x(2,ixs(2,iv2))+x(2,ixs(3,iv2))+x(2,ixs(4,iv2))+x(2,ixs(5,iv2))
175 . +x(2,ixs(6,iv2))+x(2,ixs(7,iv2))+x(2,ixs(8,iv2))+x(2,ixs(9,iv2))
176 d2z(i) = - zc(i)
177 . +x(3,ixs(2,iv2))+x(3,ixs(3,iv2))+x(3,ixs(4,iv2))+x(3,ixs(5,iv2))
178 . +x(3,ixs(6,iv2))+x(3,ixs(7,iv2))+x(3,ixs(8,iv2))+x(3,ixs(9,iv2))
179 ENDDO
180 DO i=lft,llt
181 ie =nft+i
182 iad2 = ale_connectivity%ee_connect%iad_connect(ie)
183 iv3=ale_connectivity%ee_connect%connected(iad2 + 3 - 1)
184 IF(iv3 <= 0) iv3=ie
185 d3x(i) = - xc(i)
186 . +x(1,ixs(2,iv3))+x(1,ixs(3,iv3))+x(1,ixs(4,iv3))+x(1,ixs(5,iv3))
187 . +x(1,ixs(6,iv3))+x(1,ixs(7,iv3))+x(1,ixs(8,iv3))+x(1,ixs(9,iv3))
188 d3y(i) = - yc(i)
189 . +x(2,ixs(2,iv3))+x(2,ixs(3,iv3))+x(2,ixs(4,iv3))+x(2,ixs(5,iv3))
190 . +x(2,ixs(6,iv3))+x(2,ixs(7,iv3))+x(2,ixs(8,iv3))+x(2,ixs(9,iv3))
191 d3z(i) = - zc(i)
192 . +x(3,ixs(2,iv3))+x(3,ixs(3,iv3))+x(3,ixs(4,iv3))+x(3,ixs(5,iv3))
193 . +x(3,ixs(6,iv3))+x(3,ixs(7,iv3))+x(3,ixs(8,iv3))+x(3,ixs(9,iv3))
194 ENDDO
195 DO i=lft,llt
196 ie =nft+i
197 iad2 = ale_connectivity%ee_connect%iad_connect(ie)
198 iv4=ale_connectivity%ee_connect%connected(iad2 + 4 - 1)
199 IF(iv4 <= 0) iv4=ie
200 d4x(i) = - xc(i)
201 . +x(1,ixs(2,iv4))+x(1,ixs(3,iv4))+x(1,ixs(4,iv4))+x(1,ixs(5,iv4))
202 . +x(1,ixs(6,iv4))+x(1,ixs(7,iv4))+x(1,ixs(8,iv4))+x(1,ixs(9,iv4))
203 d4y(i) = - yc(i)
204 . +x(2,ixs(2,iv4))+x(2,ixs(3,iv4))+x(2,ixs(4,iv4))+x(2,ixs(5,iv4))
205 . +x(2,ixs(6,iv4))+x(2,ixs(7,iv4))+x(2,ixs(8,iv4))+x(2,ixs(9,iv4))
206 d4z(i) = - zc(i)
207 . +x(3,ixs(2,iv4))+x(3,ixs(3,iv4))+x(3,ixs(4,iv4))+x(3,ixs(5,iv4))
208 . +x(3,ixs(6,iv4))+x(3,ixs(7,iv4))+x(3,ixs(8,iv4))+x(3,ixs(9,iv4))
209 ENDDO
210 DO i=lft,llt
211 ie =nft+i
212 iad2 = ale_connectivity%ee_connect%iad_connect(ie)
213 iv5=ale_connectivity%ee_connect%connected(iad2 + 5 - 1)
214 IF(iv5 <= 0) iv5=ie
215 d5x(i) = - xc(i)
216 . +x(1,ixs(2,iv5))+x(1,ixs(3,iv5))+x(1,ixs(4,iv5))+x(1,ixs(5,iv5))
217 . +x(1,ixs(6,iv5))+x(1,ixs(7,iv5))+x(1,ixs(8,iv5))+x(1,ixs(9,iv5))
218 d5y(i) = - yc(i)
219 . +x(2,ixs(2,iv5))+x(2,ixs(3,iv5))+x(2,ixs(4,iv5))+x(2,ixs(5,iv5))
220 . +x(2,ixs(6,iv5))+x(2,ixs(7,iv5))+x(2,ixs(8,iv5))+x(2,ixs(9,iv5))
221 d5z(i) = - zc(i)
222 . +x(3,ixs(2,iv5))+x(3,ixs(3,iv5))+x(3,ixs(4,iv5))+x(3,ixs(5,iv5))
223 . +x(3,ixs(6,iv5))+x(3,ixs(7,iv5))+x(3,ixs(8,iv5))+x(3,ixs(9,iv5))
224 ENDDO
225 DO i=lft,llt
226 ie =nft+i
227 iad2 = ale_connectivity%ee_connect%iad_connect(ie)
228 iv6=ale_connectivity%ee_connect%connected(iad2 + 6 - 1)
229 IF(iv6 <= 0) iv6=ie
230 d6x(i) = - xc(i)
231 . +x(1,ixs(2,iv6))+x(1,ixs(3,iv6))+x(1,ixs(4,iv6))+x(1,ixs(5,iv6))
232 . +x(1,ixs(6,iv6))+x(1,ixs(7,iv6))+x(1,ixs(8,iv6))+x(1,ixs(9,iv6))
233 d6y(i) = - yc(i)
234 . +x(2,ixs(2,iv6))+x(2,ixs(3,iv6))+x(2,ixs(4,iv6))+x(2,ixs(5,iv6))
235 . +x(2,ixs(6,iv6))+x(2,ixs(7,iv6))+x(2,ixs(8,iv6))+x(2,ixs(9,iv6))
236 d6z(i) = - zc(i)
237 . +x(3,ixs(2,iv6))+x(3,ixs(3,iv6))+x(3,ixs(4,iv6))+x(3,ixs(5,iv6))
238 . +x(3,ixs(6,iv6))+x(3,ixs(7,iv6))+x(3,ixs(8,iv6))+x(3,ixs(9,iv6))
239 ENDDO
240
241 DO i=lft,llt
242 dd1(i) = d1x(i)**2+d1y(i)**2+d1z(i)**2
243 dd2(i) = d2x(i)**2+d2y(i)**2+d2z(i)**2
244 dd3(i) = d3x(i)**2+d3y(i)**2+d3z(i)**2
245 dd4(i) = d4x(i)**2+d4y(i)**2+d4z(i)**2
246 dd5(i) = d5x(i)**2+d5y(i)**2+d5z(i)**2
247 dd6(i) = d6x(i)**2+d6y(i)**2+d6z(i)**2
248 ENDDO
249
250
251
252 DO i=lft,llt
253 grad(1,i)=four*
max(zero,d1x(i)*n1x(i) + d1y(i)*n1y(i)+d1z(i)*n1z(i)) /
max(em15,dd1(i))
254 grad(2,i)=four*
max(zero,d2x(i)*n2x(i) + d2y(i)*n2y(i)+d2z(i)*n2z(i)) /
max(em15,dd2(i))
255 grad(3,i)=four*
max(zero,d3x(i)*n3x(i) + d3y(i)*n3y(i)+d3z(i)*n3z(i)) /
max(em15,dd3(i))
256 grad(4,i)=four*
max(zero,d4x(i)*n4x(i) + d4y(i)*n4y(i)+d4z(i)*n4z(i)) /
max(em15,dd4(i))
257 grad(5,i)=four*
max(zero,d5x(i)*n5x(i) + d5y(i)*n5y(i)+d5z(i)*n5z(i)) /
max(em15,dd5(i))
258 grad(6,i)=four*
max(zero,d6x(i)*n6x(i) + d6y(i)*n6y(i)+d6z(i)*n6z(i)) /
max(em15,dd6(i))
259 ENDDO
260
261 RETURN