31
32
33
34#include "implicit_f.inc"
35
36
37
38 INTEGER N0,N1,N2,NC, LLL(*),JLL(*),SLL(*),IADLL(*)
40 . x(3,*),xll(*),sk(9),l1(3),l2(3),
alpha
41
42
43
44 INTEGER IK,IC,IAD
46 . xt1(3),yt1(3),zt1(3),xt2(3),yt2(3),zt2(3),
47 . x0,y0,z0,x1,y1,z1,x2,y2,z2,xr,yr,zr,xm,ym,zm,
norm
48
49
50
51
52
53
54
55
56
57
58
59
60
61 xt1(1) = sk(1)*l1(1)+sk(4)*l1(2)+sk(7)*l1(3)
62 xt1(2) = sk(2)*l1(1)+sk(5)*l1(2)+sk(8)*l1(3)
63 xt1(3) = sk(3)*l1(1)+sk(6)*l1(2)+sk(9)*l1(3)
64 norm = one/sqrt(xt1(1)*xt1(1)+xt1(2)*xt1(2)+xt1(3)*xt1(3))
68 IF (abs(xt1(1))>zep99) THEN
69 zt1(1) =-xt1(3)
70 zt1(2) = zero
71 zt1(3) = xt1(1)
72 ELSE
73 zt1(1) = zero
74 zt1(2) =-xt1(3)
75 zt1(3) = xt1(2)
76 ENDIF
77 norm = one/sqrt(zt1(1)*zt1(1)+zt1(2)*zt1(2)+zt1(3)*zt1(3))
81 yt1(1) = zt1(2)*xt1(3)-zt1(3)*xt1(2)
82 yt1(2) = zt1(3)*xt1(1)-zt1(1)*xt1(3)
83 yt1(3) = zt1(1)*xt1(2)-zt1(2)*xt1(1)
84
85 xt2(1) = sk(1)*l2(1)+sk(4)*l2(2)+sk(7)*l2(3)
86 xt2(2) = sk(2)*l2(1)+sk(5)*l2(2)+sk(8)*l2(3)
87 xt2(3) = sk(3)*l2(1)+sk(6)*l2(2)+sk(9)*l2(3)
88 norm = one/sqrt(xt2(1)*xt2(1)+xt2(2)*xt2(2)+xt2(3)*xt2(3))
92 IF (abs(xt2(1))>zep99) THEN
93 zt2(1) =-xt2(3)
94 zt2(2) = zero
95 zt2(3) = xt2(1)
96 ELSE
97 zt2(1) = zero
98 zt2(2) =-xt2(3)
99 zt2(3) = xt2(2)
100 ENDIF
101 norm = one/sqrt(zt1(1)*zt1(1)+zt1(2)*zt1(2)+zt1(3)*zt1(3))
105 yt2(1) = zt2(2)*xt2(3)-zt2(3)*xt2(2)
106 yt2(2) = zt2(3)*xt2(1)-zt2(1)*xt2(3)
107 yt2(3) = zt2(1)*xt2(2)-zt2(2)*xt2(1)
108
109 x0=x(1,n0)
110 y0=x(2,n0)
111 z0=x(3,n0)
112 x1=x(1,n1) - x0
113 y1=x(2,n1) - y0
114 z1=x(3,n1) - z0
115 x2=x(1,n2) - x0
116 y2=x(2,n2) - y0
117 z2=x(3,n2) - z0
118
119
120 nc = nc + 1
121 iadll(nc+1) = iadll(nc) + 4
122 ik = iadll(nc)
123 lll(ik) = n1
124 jll(ik) = 1
125 sll(ik) = 0
126 xll(ik) =-one
127 ik = ik+1
128 lll(ik) = n0
129 jll(ik) = 1
130 sll(ik) = 0
131 xll(ik) = one
132 ik = ik+1
133 lll(ik) = n0
134 jll(ik) = 5
135 sll(ik) = 0
136 xll(ik) = z1
137 ik = ik+1
138 lll(ik) = n0
139 jll(ik) = 6
140 sll(ik) = 0
141 xll(ik) =-y1
142
143
144 nc = nc + 1
145 iadll(nc+1) = iadll(nc) + 4
146 ik = iadll(nc)
147 lll(ik) = n1
148 jll(ik) = 2
149 sll(ik) = 0
150 xll(ik) =-one
151 ik = ik+1
152 lll(ik) = n0
153 jll(ik) = 2
154 sll(ik) = 0
155 xll(ik) = one
156 ik = ik+1
157 lll(ik) = n0
158 jll(ik) = 6
159 sll(ik) = 0
160 xll(ik) = x1
161 ik = ik+1
162 lll(ik) = n0
163 jll(ik) = 4
164 sll(ik) = 0
165 xll(ik) =-z1
166
167
168 nc = nc + 1
169 iadll(nc+1) = iadll(nc) + 4
170 ik = iadll(nc)
171 lll(ik) = n1
172 jll(ik) = 3
173 sll(ik) = 0
174 xll(ik) =-one
175 ik = ik+1
176 lll(ik) = n0
177 jll(ik) = 3
178 sll(ik) = 0
179 xll(ik) = one
180 ik = ik+1
181 lll(ik) = n0
182 jll(ik) = 4
183 sll(ik) = 0
184 xll(ik) = y1
185 ik = ik+1
186 lll(ik) = n0
187 jll(ik) = 5
188 sll(ik) = 0
189 xll(ik) =-x1
190
191
192
193
194 nc = nc + 1
195 iadll(nc+1) = iadll(nc) + 2
196 ik = iadll(nc)
197 lll(ik) = n2
198 jll(ik) = 4
199 sll(ik) = 0
200 xll(ik) = one
201 ik = ik+1
202 lll(ik) = n0
203 jll(ik) = 4
204 sll(ik) = 0
205 xll(ik) =-one
206
207
208 nc = nc + 1
209 iadll(nc+1) = iadll(nc) + 2
210 ik = iadll(nc)
211 lll(ik) = n2
212 jll(ik) = 5
213 sll(ik) = 0
214 xll(ik) = one
215 ik = ik+1
216 lll(ik) = n0
217 jll(ik) = 5
218 sll(ik) = 0
219 xll(ik) =-one
220
221
222 nc = nc + 1
223 iadll(nc+1) = iadll(nc) + 2
224 ik = iadll(nc)
225 lll(ik) = n2
226 jll(ik) = 6
227 sll(ik) = 0
228 xll(ik) = one
229 ik = ik+1
230 lll(ik) = n0
231 jll(ik) = 6
232 sll(ik) = 0
233 xll(ik) =-one
234
235
236
237
238 nc = nc + 1
239 iadll(nc+1) = iadll(nc) + 12
240 ik = iadll(nc)
241 lll(ik) = n1
242 jll(ik) = 4
243 sll(ik) = 0
244 xll(ik) =
alpha*xt1(1)
245 ik = ik+1
246 lll(ik) = n1
247 jll(ik) = 5
248 sll(ik) = 0
249 xll(ik) =
alpha*xt1(2)
250 ik = ik+1
251 lll(ik) = n1
252 jll(ik) = 6
253 sll(ik) = 0
254 xll(ik) =
alpha*xt1(3)
255
256 ik = ik+1
257 lll(ik) = n0
258 jll(ik) = 4
259 sll(ik) = 0
260 xll(ik) =-
alpha*xt1(1)
261 ik = ik+1
262 lll(ik) = n0
263 jll(ik) = 5
264 sll(ik) = 0
265 xll(ik) =-
alpha*xt1(2)
266 ik = ik+1
267 lll(ik) = n0
268 jll(ik) = 6
269 sll(ik) = 0
270 xll(ik) =-
alpha*xt1(3)
271
272 ik = ik+1
273 lll(ik) = n2
274 jll(ik) = 1
275 sll(ik) = 0
276 xll(ik) = xt2(1)
277 ik = ik+1
278 lll(ik) = n2
279 jll(ik) = 2
280 sll(ik) = 0
281 xll(ik) = xt2(2)
282 ik = ik+1
283 lll(ik) = n2
284 jll(ik) = 3
285 sll(ik) = 0
286 xll(ik) = xt2(3)
287
288 ik = ik+1
289 lll(ik) = n0
290 jll(ik) = 1
291 sll(ik) = 0
292 xll(ik) =-xt2(1)
293 ik = ik+1
294 lll(ik) = n0
295 jll(ik) = 2
296 sll(ik) = 0
297 xll(ik) =-xt2(2)
298 ik = ik+1
299 lll(ik) = n0
300 jll(ik) = 3
301 sll(ik) = 0
302 xll(ik) =-xt2(3)
303
304
305 nc = nc + 1
306 iadll(nc+1) = iadll(nc) + 9
307 ik = iadll(nc)
308 lll(ik) = n2
309 jll(ik) = 1
310 sll(ik) = 0
311 xll(ik) =-yt2(1)
312 ik = ik+1
313 lll(ik) = n2
314 jll(ik) = 2
315 sll(ik) = 0
316 xll(ik) =-yt2(2)
317 ik = ik+1
318 lll(ik) = n2
319 jll(ik) = 3
320 sll(ik) = 0
321 xll(ik) =-yt2(3)
322 ik = ik+1
323 lll(ik) = n0
324 jll(ik) = 1
325 sll(ik) = 0
326 xll(ik) = yt2(1)
327 ik = ik+1
328 lll(ik) = n0
329 jll(ik) = 2
330 sll(ik) = 0
331 xll(ik) = yt2(2)
332 ik = ik+1
333 lll(ik) = n0
334 jll(ik) = 3
335 sll(ik) = 0
336 xll(ik) = yt2(3)
337
338 ik = ik+1
339 lll(ik) = n2
340 jll(ik) = 4
341 sll(ik) = 0
342 xll(ik) = y2*yt2(3) - z2*yt2(2)
343 ik = ik+1
344 lll(ik) = n2
345 jll(ik) = 5
346 sll(ik) = 0
347 xll(ik) = z2*yt2(1) - x2*yt2(3)
348 ik = ik+1
349 lll(ik) = n2
350 jll(ik) = 6
351 sll(ik) = 0
352 xll(ik) = x2*yt2(2) - y2*yt2(1)
353
354
355 nc = nc + 1
356 iadll(nc+1) = iadll(nc) + 9
357 ik = iadll(nc)
358 lll(ik) = n2
359 jll(ik) = 1
360 sll(ik) = 0
361 xll(ik) =-zt2(1)
362 ik = ik+1
363 lll(ik) = n2
364 jll(ik) = 2
365 sll(ik) = 0
366 xll(ik) =-zt2(2)
367 ik = ik+1
368 lll(ik) = n2
369 jll(ik) = 3
370 sll(ik) = 0
371 xll(ik) =-zt2(3)
372 ik = ik+1
373 lll(ik) = n0
374 jll(ik) = 1
375 sll(ik) = 0
376 xll(ik) = zt2(1)
377 ik = ik+1
378 lll(ik) = n0
379 jll(ik) = 2
380 sll(ik) = 0
381 xll(ik) = zt2(2)
382 ik = ik+1
383 lll(ik) = n0
384 jll(ik) = 3
385 sll(ik) = 0
386 xll(ik) = zt2(3)
387
388 ik = ik+1
389 lll(ik) = n2
390 jll(ik) = 4
391 sll(ik) = 0
392 xll(ik) = y2*zt2(3) - z2*zt2(2)
393 ik = ik+1
394 lll(ik) = n2
395 jll(ik) = 5
396 sll(ik) = 0
397 xll(ik) = z2*zt2(1) - x2*zt2(3)
398 ik = ik+1
399 lll(ik) = n2
400 jll(ik) = 6
401 sll(ik) = 0
402 xll(ik) = x1*zt2(2) - y2*zt2(1)
403
404 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB