42
43
44
45#include "implicit_f.inc"
46
47
48
49#include "mvsiz_p.inc"
50
51
52
53 INTEGER, INTENT(IN) :: NEL
54 INTEGER IXS(NIXS,*), IPARTS(*), TAGPRT_SMS(*), TAGELEM_SMS(*)
55
57 . deltax(*),volg(*),
58 . x1(*), x2(*), x3(*), x4(*), x5(*), x6(*), x7(*), x8(*),
59 . y1(*), y2(*), y3(*), y4(*), y5(*), y6(*), y7(*), y8(*),
60 . z1(*), z2(*), z3(*), z4(*), z5(*), z6(*), z7(*), z8(*)
61
62
63
64#include "scr17_c.inc"
65#include "sms_c.inc"
66
67
68
69 INTEGER I, J, IT, IPT, IDEGE(MVSIZ), NINDX, NJNDX,
70 . INDEX(MVSIZ), JNDEX(MVSIZ)
71
73 . aj11, aj12, aj13, aj21,
74 . aj22, aj23, aj31, aj32,
75 . aj33, ai11, ai12, ai13,
76 . ai21, ai22, ai23, ai31,
77 . ai32, ai33
78
80 . x12(mvsiz), x34(mvsiz), x56(mvsiz),
81 . x78(mvsiz), y12(mvsiz), y34(mvsiz), y56(mvsiz), y78(mvsiz),
82 . z12(mvsiz), z34(mvsiz), z56(mvsiz), z78(mvsiz), x14(mvsiz),
83 . x23(mvsiz), x58(mvsiz), x67(mvsiz), y14(mvsiz), y23(mvsiz),
84 . y58(mvsiz), y67(mvsiz), z14(mvsiz), z23(mvsiz), z58(mvsiz),
85 . z67(mvsiz), x15(mvsiz), x26(mvsiz), x37(mvsiz), x48(mvsiz),
86 . y15(mvsiz), y26(mvsiz), y37(mvsiz), y48(mvsiz), z15(mvsiz),
87 . z26(mvsiz), z37(mvsiz), z48(mvsiz), h(8), vlinv,
88 . xx1,yy1,zz1,xx2,yy2,zz2,xx3,yy3,zz3,smax(mvsiz),vmin,
89 . p1(8), p2(8), p3(8),vlinc(mvsiz,8)
91 . aream(mvsiz), fac, v_g
92
93
94 nindx=0
95 njndx=0
96
97 IF (isms_selec==1) THEN
98
99 DO i=1,nel
100 njndx=njndx+1
101 jndex(njndx)=i
102 END DO
103 ELSEIF (isms_selec==2) THEN
104
105 DO i=1,nel
106 IF(tagprt_sms(iparts(i))==0)THEN
107 nindx=nindx+1
108 index(nindx)=i
109 ELSE
110 njndx=njndx+1
111 jndex(njndx)=i
112 END IF
113 END DO
114 ELSEIF (isms_selec==3) THEN
115
116 DO i=1,nel
117 IF(tagelem_sms(i)==0)THEN
118 nindx=nindx+1
119 index(nindx)=i
120 ELSE
121 njndx=njndx+1
122 jndex(njndx)=i
123 END IF
124 END DO
125 ELSEIF (isms_selec==4) THEN
126
127 DO i=1,nel
128 IF ((tagelem_sms(i)==0).AND.(tagprt_sms(iparts(i))==0)) THEN
129 nindx=nindx+1
130 index(nindx)=i
131 ELSE
132 njndx=njndx+1
133 jndex(njndx)=i
134 END IF
135 END DO
136 ENDIF
137
138 IF(idts6 > 0)THEN
141 ELSE
143 END IF
144
145
146
147 DO j=1,nindx
148 i=index(j)
149
150 aream(i) =em20
151
152 CALL idege8(x1(i),x2(i),x3(i),x4(i),y1(i),y2(i),y3(i),y4(i),
153 . z1(i),z2(i),z3(i),z4(i),aream(i),fac
154 CALL idege8(x5(i),x6(i),x7(i),x8(i),y5(i),y6(i),y7(i),y8(i),
155 . z5(i),z6(i),z7(i),z8(i),aream(i),fac,it)
156 CALL idege8(x1(i),x2(i),x6(i),x5(i),y1(i),y2(i),y6(i),y5(i),
157 . z1(i),z2(i),z6(i),z5(i),aream(i),fac,it)
158 CALL idege8(x2(i),x3(i),x7(i),x6(i),y2(i),y3(i),y7(i),y6(i),
159 . z2(i),z3(i),z7(i),z6(i),aream(i),fac,it)
160 CALL idege8(x3(i),x4(i),x8(i),x7(i),y3(i),y4(i),y8(i),y7(i),
161 . z3(i),z4(i),z8(i),z7(i),aream(i),fac,it)
162 CALL idege8(x4(i),x1(i),x5(i),x8(i),y4(i),y1(i),y5(i),y8(i),
163 . z4(i),z1(i),z5(i),z8(i),aream(i),fac,it)
164
166
167 IF (
idege(i) > 2)
THEN
168 fac=one_over_9
169 ELSEIF (
idege(i) > 1)
THEN
170 fac=fourth
171 ELSE
172 fac=one
173 END IF
174
175 IF (it ==0 ) aream(i)=fac*aream(i)
176
177 IF (
idege(i) > 3 )
THEN
178 aream(i)=aream(i)*fac
180 . x1(i), x2(i), x3(i), x4(i), x5(i), x6(i), x7(i), x8(i),
181 . y1(i), y2(i), y3(i), y4(i), y5(i), y6(i), y7(i), y8
182 . z1(i), z2(i), z3(i), z4(i), z5(i), z6(i), z7(i), z8(i))
183 ELSE
184 v_g=volg(i)
185 END IF
186 ENDIF
187
188 deltax(i)=four*volg(i)/sqrt(aream(i))
189 END DO
190
191
192
193 DO j=1,njndx
194 i=jndex(j)
196
197 aream(i) =em20
198
199 IF (
idege(i) > 2)
THEN
200 fac=one_over_9
201 ELSEIF (
idege(i) > 1)
THEN
202 fac=fourth
203 ELSE
204 fac=one
205 END IF
206 it = 0
207 CALL idege8(x1(i),x2(i),x3(i),x4(i),y1(i),y2(i),y3(i),y4(i),
208 . z1(i),z2(i),z3(i),z4(i),aream(i),fac,it)
209 CALL idege8(x5(i),x6(i),x7(i),x8(i),y5(i),y6(i),y7(i),y8(i),
210 . z5(i),z6(i),z7(i),z8(i),aream(i),fac,it)
211 CALL idege8(x1(i),x2(i),x6(i),x5(i),y1(i)
212 . z1(i),z2(i),z6(i),z5(i),aream(i),fac,it)
213 CALL idege8(x2(i),x3(i),x7(i),x6(i),y2(i),y3(i),y7(i),y6(i),
214 . z2(i),z3(i),z7(i),z6(i),aream(i),fac,it)
215 CALL idege8(x3(i),x4(i),x8(i),x7(i),y3(i),y4(i),y8(i),y7(i),
216 . z3(i),z4(i),z8(i),z7(i),aream(i),fac,it)
217 CALL idege8(x4(i),x1(i),x5(i),x8(i),y4(i),y1(i),y5(i),y8(i),
218 . z4(i),z1(i),z5(i),z8(i),aream(i),fac,it)
219
220 IF (it ==0 ) aream(i)=fac*aream(i)
221 IF (
idege(i) > 3 )
THEN
222 aream(i)=aream(i)*fac
224 . x1(i), x2(i), x3(i), x4(i), x5(i), x6(i), x7(i), x8(i),
225 . y1(i), y2(i), y3(i), y4(i), y5(i), y6(i), y7(i), y8(i),
226 . z1(i), z2(i), z3(i), z4(i), z5(i), z6(i)
227 ELSE
228 v_g=volg(i)
229 END IF
230 deltax(i) = four*v_g/sqrt(aream(i))
231 END IF
232 END DO
233
234 DO i=1,nel
235 x12(i)=x1(i)-x2(i)
236 y12(i)=y1(i)-y2(i)
237 z12(i)=z1(i)-z2(i)
238 x34(i)=x3(i)-x4(i)
239 y34(i)=y3(i)-y4(i)
240 z34(i)=z3(i)-z4(i)
241 x56(i)=x5(i)-x6(i)
242 y56(i)=y5(i)-y6(i)
243 z56(i)=z5(i)-z6(i)
244 x78(i)=x7(i)-x8(i)
245 y78(i)=y7(i)-y8(i)
246 z78(i)=z7(i)-z8(i)
247 x14(i)=x1(i)-x4(i)
248 y14(i)=y1(i)-y4(i)
249 z14(i)=z1(i)-z4(i)
250 x23(i)=x2(i)-x3(i)
251 y23(i)=y2(i)-y3(i)
252 z23(i)=z2(i)-z3(i)
253 x58(i)=x5(i)-x8(i)
254 y58(i)=y5(i)-y8(i)
255 z58(i)=z5(i)-z8(i)
256 x67(i)=x6(i)-x7(i)
257 y67(i)=y6(i)-y7(i)
258 z67(i)=z6(i)-z7(i)
259 x15(i)=x1(i)-x5(i)
260 y15(i)=y1(i)-y5(i)
261 z15(i)=z1(i)-z5(i)
262 x26(i)=x2(i)-x6(i)
263 y26(i)=y2(i)-y6(i)
264 z26(i)=z2(i)-z6(i)
265 x37(i)=x3(i)-x7(i)
266 y37(i)=y3(i)-y7(i)
267 z37(i)=z3(i)-z7(i)
268 x48(i)=x4(i)-x8(i)
269 y48(i)=y4(i)-y8(i)
270 z48(i)=z4(i)-z8(i)
271 ENDDO
272
273 DO ipt=1,8
274 CALL basisf (h,p1,p2,p3,ipt)
275
276 DO i=1,nel
277 aj11=p1(1)*x12(i)+p1(3)*x34(i)+p1(5)*x56(i)+p1(7)*x78(i)
278 aj12=p1(1)*y12(i)+p1(3)*y34(i)+p1(5)*y56(i)+p1(7)*y78(i)
279 aj13=p1(1)*z12(i)+p1(3)*z34(i)+p1(5)*z56(i)+p1(7)*z78(i)
280 aj21=p2(1)*x14(i)+p2(2)*x23(i)+p2(5)*x58(i)+p2(6)*x67(i)
281 aj22=p2(1)*y14(i)+p2(2)*y23(i)+p2(5)*y58(i)+p2(6)*y67(i)
282 aj23=p2(1)*z14(i)+p2(2)*z23(i)+p2(5)*z58(i)+p2(6)*z67(i)
283 aj31=p3(1)*x15(i)+p3(2)*x26(i)+p3(3)*x37(i)+p3(4)*x48(i)
284 aj32=p3(1)*y15(i)+p3(2)*y26(i)+p3(3)*y37(i)+p3(4)*y48(i)
285 aj33=p3(1)*z15(i)+p3(2)*z26(i)+p3(3)*z37(i)+p3(4)*z48(i)
286 ai11= aj22*aj33-aj23*aj32
287 ai21=-aj21*aj33+aj23*aj31
288 ai31= aj21*aj32-aj22*aj31
289 vlinc(i,ipt)=aj11*ai11+aj12*ai21+aj13*ai31
290 ENDDO
291 END DO
292
293 DO i=1,nel
294
295 xx1 = x1(i) + x2(i) + x3(i) + x4(i)
296 . - x5(i) - x6(i) - x7(i) - x8(i)
297 yy1 = y1(i) + y2(i) + y3(i) + y4(i)
298 . - y5(i) - y6(i) - y7(i) - y8(i)
299 zz1 = z1(i) + z2(i) + z3(i) + z4(i)
300 . - z5(i) - z6(i) - z7(i) - z8(i)
301 xx2 = x1(i) + x2(i) + x5(i) + x6(i)
302 . - x3(i) - x4(i) - x7(i) - x8(i)
303 yy2 = y1(i) + y2(i) + y5(i) + y6(i)
304 . - y3(i) - y4(i) - y7(i) - y8(i)
305 zz2 = z1(i) + z2(i) + z5(i) + z6(i)
306 . - z3(i) - z4(i) - z7(i) - z8(i)
307 xx3 = x1(i) + x4(i) + x5(i) + x8(i)
308 . - x3(i) - x2(i) - x7(i) - x6(i)
309 yy3 = y1(i) + y4(i) + y5(i) + y8(i)
310 . - y3(i) - y2(i) - y7(i) - y6(i)
311 zz3 = z1(i) + z4(i) + z5(i) + z8(i)
312 . - z3(i) - z2(i) - z7(i) - z6(i)
313
314 smax(i) = (yy1 * zz2 - yy2 * zz1)**2
315 . + (zz1 * xx2 - zz2 * xx1)**2
316 . + (xx1 * yy2 - xx2 * yy1)**2
317 smax(i) =
max(smax(i),(yy1 * zz3 - yy3 * zz1)**2
318 . + (zz1 * xx3 - zz3 * xx1)**2
319 . + (xx1 * yy3 - xx3 * yy1)**2)
320 smax(i) =
max(smax(i),(yy3 * zz2 - yy2 * zz3)**2
321 . + (zz3 * xx2 - zz2 * xx3)**2
322 . + (xx3 * yy2 - xx2 * yy3)**2)
323 ENDDO
324
325
326 IF (idts6>0) THEN
327 DO i=1,nel
329 vmin =
min(vlinc(i,1),vlinc(i,2),vlinc(i,3),vlinc(i,4),
330 . vlinc(i,5),vlinc(i,6),vlinc(i,7),vlinc(i,8))
331 deltax(i)=hundred28*vmin/sqrt(smax
332 ENDIF
333 ENDDO
334 ELSE
335 DO i=1,nel
336 vmin =
min(vlinc(i,1),vlinc(i,2),vlinc(i,3),vlinc(i,4),
337 . vlinc(i,5),vlinc(i,6),vlinc(i,7),vlinc(i,8))
338 deltax(i)=hundred28*vmin/sqrt(smax(i))
339 ENDDO
340 ENDIF
341
342 RETURN
subroutine degenes8(ixs, idege, nel)
subroutine deges4v(det, x1, x2, x3, x4, x5, x6, x7, x8, y1, y2, y3, y4, y5, y6, y7, y8, z1, z2, z3, z4, z5, z6, z7, z8)
subroutine idege8(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, amax, fac, it)
subroutine idege(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, a, amax, fac, it4, it, indx, n_indx)
subroutine basisf(h, p1, p2, p3, ipt)