31
32
33
34#include "implicit_f.inc"
35
36
37
38#include "mvsiz_p.inc"
39
40
41
42 INTEGER NELA, NPTR, NPTS, NPTT
43 INTEGER NC1(*), NC2(*), NC3(*), NC4(*), NC5(*), NC6(*), NC7(*), NC8(*)
44 my_real x(3,*), volg(*), volp(mvsiz,*)
45
46
47
48
49 INTEGER I, J
50 INTEGER IR,IS,IT,IP
51
53 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
54 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
55 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz)
57 . aj1(mvsiz),aj2(mvsiz),aj3(mvsiz),
58 . aj4(mvsiz),aj5(mvsiz),aj6(mvsiz),
59 . aj7(mvsiz),aj8(mvsiz),aj9(mvsiz)
61 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz)
63 . x17(mvsiz) , x28(mvsiz) , x35(mvsiz) , x46(mvsiz),
64 . y17(mvsiz) , y28(mvsiz) , y35(mvsiz) , y46(mvsiz),
65 . z17(mvsiz) , z28(mvsiz) , z35(mvsiz) , z46(mvsiz),
66 . a17(mvsiz) , a28(mvsiz) ,
67 . b17(mvsiz) , b28(mvsiz) ,
68 . c17(mvsiz) , c28(mvsiz)
70 . hx(mvsiz,4), hy(mvsiz,4), hz(mvsiz,4),
71 . hx1pg(mvsiz), hx2pg(mvsiz), hx3pg(mvsiz),hx4pg2(mvsiz),
72 . hy1pg(mvsiz), hy2pg(mvsiz), hy3pg(mvsiz),hy4pg2(mvsiz),
73 . hz1pg(mvsiz), hz2pg(mvsiz), hz3pg(mvsiz),hz4pg2(mvsiz),
74 . ajp1(mvsiz,8),ajp2(mvsiz,8),ajp3(mvsiz,8),
75 . ajp4(mvsiz,8),ajp5(mvsiz,8),ajp6(mvsiz,8),
76 . ajp7(mvsiz,8),ajp8(mvsiz,8),ajp9(mvsiz,8)
78 . pg2, pg, wi
80 . w_gauss(9,9)
81 DATA w_gauss /
82 1 2. ,0. ,0. ,
83 1 0. ,0. ,0. ,
84 1 0. ,0. ,0. ,
85 2 1. ,1. ,0. ,
86 2 0. ,0. ,0. ,
87 2 0. ,0. ,0. ,
88 3 0.555555555555556,0.888888888888889,0.555555555555556,
89 3 0. ,0. ,0. ,
90 3 0. ,0. ,0. ,
91 4 0.347854845137454,0.652145154862546,0.652145154862546,
92 4 0.347854845137454,0. ,0. ,
93 4 0. ,0. ,0. ,
94 5 0.236926885056189,0.478628670499366,0.568888888888889,
95 5 0.478628670499366,0.236926885056189,0. ,
96 5 0. ,0. ,0. ,
97 6 0.171324492379170,0.360761573048139,0.467913934572691,
98 6 0.467913934572691,0.360761573048139,0.171324492379170,
99 6 0. ,0. ,0. ,
100 7 0.129484966168870,0.279705391489277,0.381830050505119,
101 7 0.417959183673469,0.381830050505119,0.279705391489277,
102 7 0.129484966168870,0. ,0. ,
103 8 0.101228536290376,0.222381034453374,0.313706645877887,
104 8 0.362683783378362,0.362683783378362,0.313706645877887,
105 8 0.222381034453374,0.101228536290376,0. ,
106 9 0.081274388361574,0.180648160694857,0.260610696402935,
107 9 0.312347077040003,0.330239355001260,0.312347077040003,
108 9 0.260610696402935,0.180648160694857,0.081274388361574/
109
110 parameter(pg=.577350269189626)
111
112
113 DO i=1,nela
114 x1(i)=x(1,nc1(i))
115 y1(i)=x(2,nc1(i))
116 z1(i)=x(3,nc1(i))
117 x2(i)=x(1,nc2(i))
118 y2(i)=x(2,nc2(i))
119 z2(i)=x(3,nc2(i))
120 x3(i)=x(1,nc3(i))
121 y3(i)=x(2,nc3(i))
122 z3(i)=x(3,nc3(i))
123 x4(i)=x(1,nc4(i))
124 y4(i)=x(2,nc4(i))
125 z4(i)=x(3,nc4(i))
126 x5(i)=x(1,nc5(i))
127 y5(i)=x(2,nc5(i))
128 z5(i)=x(3,nc5(i))
129 x6(i)=x(1,nc6(i))
130 y6(i)=x(2,nc6(i))
131 z6(i)=x(3,nc6(i))
132 x7(i)=x(1,nc7(i))
133 y7(i)=x(2,nc7(i))
134 z7(i)=x(3,nc7(i))
135 x8(i)=x(1,nc8(i))
136 y8(i)=x(2,nc8(i))
137 z8(i)=x(3,nc8(i))
138 ENDDO
139
140 DO i=1,nela
141 x17(i)=x7(i)-x1(i)
142 x28(i)=x8(i)-x2(i)
143 x35(i)=x5(i)-x3(i)
144 x46(i)=x6(i)-x4(i)
145
146 y17(i)=y7(i)-y1(i)
147 y28(i)=y8(i)-y2(i)
148 y35(i)=y5(i)-y3(i)
149 y46(i)=y6(i)-y4(i)
150
151 z17(i)=z7(i)-z1(i)
152 z28(i)=z8(i)-z2(i)
153 z35(i)=z5(i)-z3(i)
154 z46(i)=z6(i)-z4(i)
155
156 aj4(i)=x17(i)+x28(i)-x35(i)-x46(i)
157 aj5(i)=y17(i)+y28(i)-y35(i)-y46(i)
158 aj6(i)=z17(i)+z28(i)-z35(i)-z46(i)
159
160 a17(i)=x17(i)+x46(i)
161 a28(i)=x28(i)+x35(i)
162 b17(i)=y17(i)+y46(i)
163 b28(i)=y28(i)+y35(i)
164 c17(i)=z17(i)+z46(i)
165 c28(i)=z28(i)+z35(i)
166
167 aj7(i)=a17(i)+a28(i)
168 aj8(i)=b17(i)+b28(i)
169 aj9(i)=c17(i)+c28(i)
170 aj1(i)=a17(i)-a28(i)
171 aj2(i)=b17(i)-b28(i)
172 aj3(i)=c17(i)-c28(i)
173
174 jac_59_68(i)=aj5(i)*aj9(i)-aj6(i)*aj8(i)
175 jac_67_49(i)=aj6(i)*aj7(i)-aj4(i)*aj9(i)
176 jac_48_57(i)=aj4(i)*aj8(i)-aj5(i)*aj7(i)
177
178 volg(i)=one_over_64*(aj1(i)*jac_59_68(i)+aj2(i)*jac_67_49(i)+aj3(i)*jac_48_57(i))
179 ENDDO
180
181
182
183 DO i=1,nela
184 hx(i,1)=(x1(i)+x2(i)-x3(i)-x4(i)-x5(i)-x6(i)+x7(i)+x8(i))
185 hy(i,1)=(y1(i)+y2(i)-y3(i)-y4(i)-y5(i)-y6(i)+y7(i)+y8(i))
186 hz(i,1)=(z1(i)+z2(i)-z3(i)-z4(i)-z5(i)-z6(i)+z7(i)+z8(i))
187 ENDDO
188
189
190 DO i=1,nela
191 hx(i,2)=(x1(i)-x2(i)-x3(i)+x4(i)-x5(i)+x6(i)+x7(i)-x8(i))
192 hy(i,2)=(y1(i)-y2(i)-y3(i)+y4(i)-y5(i)+y6(i)+y7(i)-y8(i))
193 hz(i,2)=(z1(i)-z2(i)-z3(i)+z4(i)-z5(i)+z6(i)+z7(i)-z8(i))
194 ENDDO
195
196
197 DO i=1,nela
198 hx(i,3)=(x1(i)-x2(i)+x3(i)-x4(i)+x5(i)-x6(i)+x7(i)-x8(i))
199 hy(i,3)=(y1(i)-y2(i)+y3(i)-y4(i)+y5(i)-y6(i)+y7(i)-y8(i))
200 hz(i,3)=(z1(i)-z2(i)+z3(i)-z4(i)+z5(i)-z6(i)+z7(i)-z8(i))
201 ENDDO
202
203
204 DO i=1,nela
205 hx(i,4)=(-x1(i)+x2(i)-x3(i)+x4(i)+x5(i)-x6(i)+x7(i)-x8(i))
206 hy(i,4)=(-y1(i)+y2(i)-y3(i)+y4(i)+y5(i)-y6(i)+y7(i)-y8(i))
207 hz(i,4)=(-z1(i)+z2(i)-z3(i)+z4(i)+z5(i)-z6(i)+z7(i)-z8(i))
208 ENDDO
209
210 pg2=pg*pg
211
212 DO i=1,nela
213 hx1pg(i) =hx(i,1)*pg
214 hx2pg(i) =hx(i,2)*pg
215 hx3pg(i) =hx(i,3)*pg
216 hx4pg2(i)=hx(i,4)*pg2
217 hy1pg(i) =hy(i,1)*pg
218 hy2pg(i) =hy(i,2)*pg
219 hy3pg(i) =hy(i,3)*pg
220 hy4pg2(i)=hy(i,4)*pg2
221 hz1pg(i) =hz(i,1)*pg
222 hz2pg(i) =hz(i,2)*pg
223 hz3pg(i) =hz(i,3)*pg
224 hz4pg2(i)=hz(i,4)*pg2
225 ENDDO
226 ip=1
227
228 DO i=1,nela
229 ajp1(i,ip)=aj1(i)-hx3pg(i)-hx2pg(i)+hx4pg2(i)
230 ajp2(i,ip)=aj2(i)-hy3pg(i)-hy2pg(i)+hy4pg2(i)
231 ajp3(i,ip)=aj3(i)-hz3pg(i)-hz2pg(i)+hz4pg2(i)
232
233 ajp4(i,ip)=aj4(i)-hx1pg(i)-hx3pg(i)+hx4pg2(i)
234 ajp5(i,ip)=aj5(i)-hy1pg(i)-hy3pg(i)+hy4pg2(i)
235 ajp6(i,ip)=aj6(i)-hz1pg(i)-hz3pg(i)+hz4pg2(i)
236
237 ajp7(i,ip)=aj7(i)-hx2pg(i)-hx1pg(i)+hx4pg2(i)
238 ajp8(i,ip)=aj8(i)-hy2pg(i)-hy1pg(i)+hy4pg2(i)
239 ajp9(i,ip)=aj9(i)-hz2pg(i)-hz1pg(i)+hz4pg2(i)
240 ENDDO
241 ip=2
242
243 DO i=1,nela
244 ajp1(i,ip)=aj1(i)-hx3pg(i)-hx2pg(i)+hx4pg2(i)
245 ajp2(i,ip)=aj2(i)-hy3pg(i)-hy2pg(i)+hy4pg2(i)
246 ajp3(i,ip)=aj3(i)-hz3pg(i)-hz2pg(i)+hz4pg2(i)
247
248 ajp4(i,ip)=aj4(i)-hx1pg(i)+hx3pg(i)-hx4pg2(i)
249 ajp5(i,ip)=aj5(i)-hy1pg(i)+hy3pg(i)-hy4pg2(i)
250 ajp6(i,ip)=aj6(i)-hz1pg(i)+hz3pg(i)-hz4pg2(i)
251
252 ajp7(i,ip)=aj7(i)+hx2pg(i)-hx1pg(i)-hx4pg2(i)
253 ajp8(i,ip)=aj8(i)+hy2pg(i)-hy1pg(i)-hy4pg2(i)
254 ajp9(i,ip)=aj9(i)+hz2pg(i)-hz1pg(i)-hz4pg2(i)
255 ENDDO
256 ip=3
257
258 DO i=1,nela
259 ajp1(i,ip)=aj1(i)+hx3pg(i)-hx2pg(i)-hx4pg2(i)
260 ajp2(i,ip)=aj2(i)+hy3pg(i)-hy2pg(i)-hy4pg2(i)
261 ajp3(i,ip)=aj3(i)+hz3pg(i)-hz2pg(i)-hz4pg2(i)
262
263 ajp4(i,ip)=aj4(i)-hx1pg(i)-hx3pg(i)+hx4pg2(i)
264 ajp5(i,ip)=aj5(i)-hy1pg(i)-hy3pg(i)+hy4pg2(i)
265 ajp6(i,ip)=aj6(i)-hz1pg(i)-hz3pg(i)+hz4pg2(i)
266
267 ajp7(i,ip)=aj7(i)-hx2pg(i)+hx1pg(i)-hx4pg2(i)
268 ajp8(i,ip)=aj8(i)-hy2pg(i)+hy1pg(i)-hy4pg2(i)
269 ajp9(i,ip)=aj9(i)-hz2pg(i)+hz1pg(i)-hz4pg2(i)
270 ENDDO
271 ip=4
272
273 DO i=1,nela
274 ajp1(i,ip)=aj1(i)+hx3pg(i)-hx2pg(i)-hx4pg2(i)
275 ajp2(i,ip)=aj2(i)+hy3pg(i)-hy2pg(i)-hy4pg2(i)
276 ajp3(i,ip)=aj3(i)+hz3pg(i)-hz2pg(i)-hz4pg2(i)
277
278 ajp4(i,ip)=aj4(i)-hx1pg(i)+hx3pg(i)-hx4pg2(i)
279 ajp5(i,ip)=aj5(i)-hy1pg(i)+hy3pg(i)-hy4pg2(i)
280 ajp6(i,ip)=aj6(i)-hz1pg(i)+hz3pg(i)-hz4pg2(i)
281
282 ajp7(i,ip)=aj7(i)+hx2pg(i)+hx1pg(i)+hx4pg2(i)
283 ajp8(i,ip)=aj8(i)+hy2pg(i)+hy1pg(i)+hy4pg2(i)
284 ajp9(i,ip)=aj9(i)+hz2pg(i)+hz1pg(i)+hz4pg2(i)
285 ENDDO
286 ip=5
287
288 DO i=1,nela
289 ajp1(i,ip)=aj1(i)-hx3pg(i)+hx2pg(i)-hx4pg2(i)
290 ajp2(i,ip)=aj2(i)-hy3pg(i)+hy2pg(i)-hy4pg2(i)
291 ajp3(i,ip)=aj3(i)-hz3pg(i)+hz2pg(i)-hz4pg2(i)
292
293 ajp4(i,ip)=aj4(i)+hx1pg(i)-hx3pg(i)-hx4pg2(i)
294 ajp5(i,ip)=aj5(i)+hy1pg(i)-hy3pg(i)-hy4pg2(i)
295 ajp6(i,ip)=aj6(i)+hz1pg(i)-hz3pg(i)-hz4pg2(i)
296
297 ajp7(i,ip)=aj7(i)-hx2pg(i)-hx1pg(i)+hx4pg2(i)
298 ajp8(i,ip)=aj8(i)-hy2pg(i)-hy1pg(i)+hy4pg2(i)
299 ajp9(i,ip)=aj9(i)-hz2pg(i)-hz1pg(i)+hz4pg2(i)
300 ENDDO
301 ip=6
302
303 DO i=1,nela
304 ajp1(i,ip)=aj1(i)-hx3pg(i)+hx2pg(i)-hx4pg2(i)
305 ajp2(i,ip)=aj2(i)-hy3pg(i)+hy2pg(i)-hy4pg2(i)
306 ajp3(i,ip)=aj3(i)-hz3pg(i)+hz2pg(i)-hz4pg2(i)
307
308 ajp4(i,ip)=aj4(i)+hx1pg(i)+hx3pg(i)+hx4pg2(i)
309 ajp5(i,ip)=aj5(i)+hy1pg(i)+hy3pg(i)+hy4pg2(i)
310 ajp6(i,ip)=aj6(i)+hz1pg(i)+hz3pg(i)+hz4pg2(i)
311
312 ajp7(i,ip)=aj7(i)+hx2pg(i)-hx1pg(i)-hx4pg2(i)
313 ajp8(i,ip)=aj8(i)+hy2pg(i)-hy1pg(i)-hy4pg2(i)
314 ajp9(i,ip)=aj9(i)+hz2pg(i)-hz1pg(i)-hz4pg2(i)
315 ENDDO
316 ip=7
317
318 DO i=1,nela
319 ajp1(i,ip)=aj1(i)+hx3pg(i)+hx2pg(i)+hx4pg2(i)
320 ajp2(i,ip)=aj2(i)+hy3pg(i)+hy2pg(i)+hy4pg2(i)
321 ajp3(i,ip)=aj3(i)+hz3pg(i)+hz2pg(i)+hz4pg2(i)
322
323 ajp4(i,ip)=aj4(i)+hx1pg(i)-hx3pg(i)-hx4pg2(i)
324 ajp5(i,ip)=aj5(i)+hy1pg(i)-hy3pg(i)-hy4pg2(i)
325 ajp6(i,ip)=aj6(i)+hz1pg(i)-hz3pg(i)-hz4pg2(i)
326
327 ajp7(i,ip)=aj7(i)-hx2pg(i)+hx1pg(i)-hx4pg2(i)
328 ajp8(i,ip)=aj8(i)-hy2pg(i)+hy1pg(i)-hy4pg2(i)
329 ajp9(i,ip)=aj9(i)-hz2pg(i)+hz1pg(i)-hz4pg2(i)
330 ENDDO
331 ip=8
332
333 DO i=1,nela
334 ajp1(i,ip)=aj1(i)+hx3pg(i)+hx2pg(i)+hx4pg2(i)
335 ajp2(i,ip)=aj2(i)+hy3pg(i)+hy2pg(i)+hy4pg2(i)
336 ajp3(i,ip)=aj3(i)+hz3pg(i)+hz2pg(i)+hz4pg2(i)
337
338 ajp4(i,ip)=aj4(i)+hx1pg(i)+hx3pg(i)+hx4pg2(i)
339 ajp5(i,ip)=aj5(i)+hy1pg(i)+hy3pg(i)+hy4pg2(i)
340 ajp6(i,ip)=aj6(i)+hz1pg(i)+hz3pg(i)+hz4pg2(i)
341
342 ajp7(i,ip)=aj7(i)+hx2pg(i)+hx1pg(i)+hx4pg2(i)
343 ajp8(i,ip)=aj8(i)+hy2pg(i)+hy1pg(i)+hy4pg2(i)
344 ajp9(i,ip)=aj9(i)+hz2pg(i)+hz1pg(i)+hz4pg2(i)
345 ENDDO
346
347
348 DO ir=1,nptr
349 DO is=1,npts
350 DO it=1,nptt
351 ip = ir + ( (is-1) + (it-1)*npts )*nptr
352 wi = w_gauss(ir,nptr)*w_gauss(is,npts)*w_gauss(it,nptt)
353
354 DO i=1,nela
355 jac_59_68(i)=ajp5(i,ip)*ajp9(i,ip)-ajp6(i,ip)*ajp8(i,ip)
356 jac_67_49(i)=ajp6(i,ip)*ajp7(i,ip)-ajp4(i,ip)*ajp9(i,ip)
357 jac_48_57(i)=ajp4(i,ip)*ajp8(i,ip)-ajp5(i,ip)*ajp7(i,ip)
358 ENDDO
359
360 DO i=1,nela
361 volp(i,ip)=one_over_512*wi*(ajp1(i,ip)*jac_59_68(i)+ajp2(i,ip)*jac_67_49(i)+ajp3(i,ip)*jac_48_57(i))
362 ENDDO
363 ENDDO
364 ENDDO
365 ENDDO
366
367 RETURN