117
118
119
120#include "implicit_f.inc"
121
122
123
124 INTEGER NSN, NINT, NC, ISKIP, NCF_S
125 INTEGER LLL(*), JLL(*), SLL(*), IADLL(*),
126 . IRECT(4,*),NSV(*),IRTL(*),COMNTAG(*),ICFTAG(*),JCFTAG(*)
128 . ltsm(*),xll(*),x(3,*),ms(*),in(*),v(3,*),vr(3,*),
129 . a(3,*),ar(3,*)
130
131
132
133 INTEGER INOD(5), I,J,K,L,IC,JC,II,JJ,IK,ISK,IAD,NIR,NDL,
134 . NC_INI,NCL,COMFLAG
136 . rx(4),ry(4),rz(4),hloc(6,6),lloc(6),
137 . b1,b2,b3,c1,c2,c3,det,s,hij,
138 . x0,x1,x2,x3,x4,xs,y0,y1,y2,y3,y4,ys,z0,z1,z2,z3,z4,zs,
139 . x12,x22,x32,x42,y12,y22,y32,y42,z12,z22,z32,z42,
140 . xx,yy,zz,xxx,yyy,zzz,xy,yz,zx,xy2,yz2,zx2,fact
141
142 ncl = 6
143 DO ii=1,nsn
144 l=irtl(ii)
145 nir = 4
146 DO jj=1,nir
147 inod(jj) = irect(jj,l)
148 ENDDO
149 IF(inod(4)==inod(3)) nir=3
150 fact = one / nir
151 inod(nir+1) = nsv(ii)
152 ndl = 3*nir+1
153 nc_ini = nc
154 comflag = 0
155 DO jj=1,nir+1
156 IF (comntag(inod(jj))>1) comflag = 1
157 ENDDO
158
159 IF (nir==4) THEN
160 x1=x(1,inod(1))
161 y1=x(2,inod(1))
162 z1=x(3,inod(1))
163 x2=x(1,inod(2))
164 y2=x(2,inod(2))
165 z2=x(3,inod(2))
166 x3=x(1,inod(3))
167 y3=x(2,inod(3))
168 z3=x(3,inod(3))
169 x4=x(1,inod(4))
170 y4=x(2,inod(4))
171 z4=x(3,inod(4))
172 x0=fourth*(x1+x2+x3+x4)
173 y0=fourth*(y1+y2+y3+y4)
174 z0=fourth*(z1+z2+z3+z4)
175 x1=x1-x0
176 y1=y1-y0
177 z1=z1-z0
178 x2=x2-x0
179 y2=y2-y0
180 z2=z2-z0
181 x3=x3-x0
182 y3=y3-y0
183 z3=z3-z0
184 x4=x4-x0
185 y4=y4-y0
186 z4=z4-z0
187
188 rx(1) = x1
189 rx(2) = x2
190 rx(3) = x3
191 rx(4) = x4
192 ry(1) = y1
193 ry(2) = y2
194 ry(3) = y3
195 ry(4) = y4
196 rz(1) = z1
197 rz(2) = z2
198 rz(3) = z3
199 rz(4) = z4
200
201 x12=x1*x1
202 x22=x2*x2
203 x32=x3*x3
204 x42=x4*x4
205 y12=y1*y1
206 y22=y2*y2
207 y32=y3*y3
208 y42=y4*y4
209 z12=z1*z1
210 z22=z2*z2
211 z32=z3*z3
212 z42=z4*z4
213 xx=x12 + x22 + x32 + x42
214 yy=y12 + y22 + y32 + y42
215 zz=z12 + z22 + z32 + z42
216 xy=x1*y1 + x2*y2 + x3*y3 + x4*y4
217 yz=y1*z1 + y2*z2 + y3*z3 + y4*z4
218 zx=z1*x1 + z2*x2 + z3*x3 + z4*x4
219 ELSEIF (nir==3) THEN
220 x1=x(1,inod(1))
221 y1=x(2,inod(1))
222 z1=x(3,inod(1))
223 x2=x(1,inod(2))
224 y2=x(2,inod(2))
225 z2=x(3,inod(2))
226 x3=x(1,inod(3))
227 y3=x(2,inod(3))
228 z3=x(3,inod(3))
229 x0=third*(x1+x2+x3)
230 y0=third*(y1+y2+y3)
231 z0=third*(z1+z2+z3)
232 x1=x1-x0
233 y1=y1-y0
234 z1=z1-z0
235 x2=x2-x0
236 y2=y2-y0
237 z2=z2-z0
238 x3=x3-x0
239 y3=y3-y0
240 z3=z3-z0
241
242 rx(1) = x1
243 rx(2) = x2
244 rx(3) = x3
245 ry(1) = y1
246 ry(2) = y2
247 ry(3) = y3
248 rz(1) = z1
249 rz(2) = z2
250 rz(3) = z3
251
252 x12=x1*x1
253 x22=x2*x2
254 x32=x3*x3
255 y12=y1*y1
256 y22=y2*y2
257 y32=y3*y3
258 z12=z1*z1
259 z22=z2*z2
260 z32=z3*z3
261 xx=x12 + x22 + x32
262 yy=y12 + y22 + y32
263 zz=z12 + z22 + z32
264 xy=x1*y1 + x2*y2 + x3*y3
265 yz=y1*z1 + y2*z2 + y3*z3
266 zx=z1*x1 + z2*x2 + z3*x3
267 ENDIF
268 xs=x(1,inod(nir+1))-x0
269 ys=x(2,inod(nir+1))-y0
270 zs=x(3,inod(nir+1))-z0
271 zzz=xx+yy
272 xxx=yy+zz
273 yyy=zz+xx
274 xy2=xy*xy
275 yz2=yz*yz
276 zx2=zx*zx
277 det=xxx*yyy*zzz - xxx*yz2 - yyy*zx2 - zzz*xy2 - two*xy*yz*zx
278 det=one/det
279 b1=zzz*yyy-yz2
280 b2=xxx*zzz-zx2
281 b3=yyy*xxx-xy2
282 c3=zzz*xy+yz*zx
283 c1=xxx*yz+zx*xy
284 c2=yyy*zx+xy*yz
285
286
287
288 nc = nc + 1
289 iadll(nc+1) = iadll(nc) + ndl
290 iad = iadll(nc) -1
291 DO jj=1,nir
292 ik = iad+jj
293 lll(ik) = inod(jj)
294 jll(ik) = 1
295 sll(ik) = 0
296 xll(ik) = fact
297 . + det*zs*(b2*rz(jj) - c1*ry(jj))
298 . - det*ys*(c1*rz(jj) - b3*ry(jj))
299 ENDDO
300 iad = iad + nir
301 DO jj=1,nir
302 ik = iad+jj
303 lll(ik) = inod(jj)
304 jll(ik) = 2
305 sll(ik) = 0
306 xll(ik) = det*zs*(c1*rx(jj) - c3*rz(jj))
307 . - det*ys*(b3*rx(jj) - c2*rz(jj))
308 ENDDO
309 iad = iad + nir
310 DO jj=1,nir
311 ik = iad+jj
312 lll(ik) = inod(jj)
313 jll(ik) = 3
314 sll(ik) = 0
315 xll(ik) = det*zs*(c3*ry(jj) - b2*rx(jj))
316 . - det*ys*(c2*ry(jj) - c1*rx(jj))
317 ENDDO
318 ik = iad + nir+1
319 lll(ik) = inod(nir + 1)
320 jll(ik) = 1
321 sll(ik) = nint
322 xll(ik) = -one
323
324 nc = nc + 1
325 iadll(nc+1) = iadll(nc) + ndl
326 iad = iadll(nc) -1
327 DO jj=1,nir
328 ik = iad+jj
329 lll(ik) = inod(jj)
330 jll(ik) = 1
331 sll(ik) = 0
332 xll(ik) = det*xs*(c1*rz(jj) - b3*ry(jj))
333 . - det*zs*(c3*rz(jj) - c2*ry(jj))
334 ENDDO
335 iad = iad + nir
336 DO jj=1,nir
337 ik = iad+jj
338 lll(ik) = inod(jj)
339 jll(ik) = 2
340 sll(ik) = 0
341 xll(ik) = fact
342 . + det*xs*(b3*rx(jj) - c2*rz(jj))
343 . - det*zs*(c2*rx(jj) - b1*rz(jj))
344 ENDDO
345 iad = iad + nir
346 DO jj=1,nir
347 ik = iad+jj
348 lll(ik) = inod(jj)
349 jll(ik) = 3
350 sll(ik) = 0
351 xll(ik) = det*xs*(c2*ry(jj) - c1*rx(jj))
352 . - det*zs*(b1*ry(jj) - c3*rx(jj))
353 ENDDO
354 ik = iad + nir+1
355 lll(ik) = inod(nir + 1)
356 jll(ik) = 2
357 sll(ik) = nint
358 xll(ik) = -one
359
360 nc = nc + 1
361 iadll(nc+1) = iadll(nc) + ndl
362 iad = iadll(nc) -1
363 DO jj=1,nir
364 ik = iad+jj
365 lll(ik) = inod(jj)
366 jll(ik) = 1
367 sll(ik) = 0
368 xll(ik) = det*ys*(c3*rz(jj) - c2*ry(jj))
369 . - det*xs*(b2*rz(jj) - c1*ry(jj))
370 ENDDO
371 iad = iad + nir
372 DO jj=1,nir
373 ik = iad+jj
374 lll(ik) = inod(jj)
375 jll(ik) = 2
376 sll(ik) = 0
377 xll(ik) = det*ys*(c2*rx(jj) - b1*rz(jj))
378 . - det*xs*(c1*rx(jj) - c3*rz(jj))
379 ENDDO
380 iad = iad + nir
381 DO jj=1,nir
382 ik = iad+jj
383 lll(ik) = inod(jj)
384 jll(ik) = 3
385 sll(ik) = 0
386 xll(ik) = fact
387 . + det*ys*(b1*ry(jj) - c3*rx(jj))
388 . - det*xs*(c3*ry(jj) - b2*rx(jj))
389 ENDDO
390 ik = iad + nir+1
391 lll(ik) = inod(nir+1)
392 jll(ik) = 3
393 sll(ik) = nint
394 xll(ik) = -one
395
396
397
398 nc = nc + 1
399 iadll(nc+1) = iadll(nc) + ndl
400 iad = iadll(nc) -1
401 DO jj=1,nir
402 ik = iad+jj
403 lll(ik) = inod(jj)
404 jll(ik) = 1
405 sll(ik) = 0
406 xll(ik) = det*(c3*rz(jj) - c2*ry(jj))
407 ENDDO
408 iad = iad + nir
409 DO jj=1,nir
410 ik = iad+jj
411 lll(ik) = inod(jj)
412 jll(ik) = 2
413 sll(ik) = 0
414 xll(ik) = det*(c2*rx(jj) - b1*rz(jj))
415 ENDDO
416 iad = iad + nir
417 DO jj=1,nir
418 ik = iad+jj
419 lll(ik) = inod(jj)
420 jll(ik) = 3
421 sll(ik) = 0
422 xll(ik) = det*(b1*ry(jj) - c3*rx(jj))
423 ENDDO
424 ik = iad + nir + 1
425 lll(ik) = inod(nir+1)
426 jll(ik) = 4
427 sll(ik) = nint
428 xll(ik) = -one
429
430 nc = nc + 1
431 iadll(nc+1) = iadll(nc) + ndl
432 iad = iadll(nc) -1
433 DO jj=1,nir
434 ik = iad+jj
435 lll(ik) = inod(jj)
436 jll(ik) = 1
437 sll(ik) = 0
438 xll(ik) = det*(b2*rz(jj) - c1*ry(jj))
439 ENDDO
440 iad = iad + nir
441 DO jj=1,nir
442 ik = iad+jj
443 lll(ik) = inod(jj)
444 jll(ik) = 2
445 sll(ik) = 0
446 xll(ik) = det*(c1*rx(jj) - c3*rz(jj))
447 ENDDO
448 iad = iad + nir
449 DO jj=1,nir
450 ik = iad+jj
451 lll(ik) = inod(jj)
452 jll(ik) = 3
453 sll(ik) = 0
454 xll(ik) = det*(c3*ry(jj) - b2*rx(jj))
455 ENDDO
456 ik = iad + nir + 1
457 lll(ik) = inod(nir+1)
458 jll(ik) = 5
459 sll(ik) = nint
460 xll(ik) = -one
461
462 nc = nc + 1
463 iadll(nc+1) = iadll(nc) + ndl
464 iad = iadll(nc) -1
465 DO jj=1,nir
466 ik = iad+jj
467 lll(ik) = inod(jj)
468 jll(ik) = 1
469 sll(ik) = 0
470 xll(ik) = det*(c1*rz(jj) - b3*ry(jj))
471 ENDDO
472 iad = iad + nir
473 DO jj=1,nir
474 ik = iad+jj
475 lll(ik) = inod(jj)
476 jll(ik) = 2
477 sll(ik) = 0
478 xll(ik) = det*(b3*rx(jj) - c2*rz(jj))
479 ENDDO
480 iad = iad + nir
481 DO jj=1,nir
482 ik = iad+jj
483 lll(ik) = inod(jj)
484 jll(ik) = 3
485 sll(ik) = 0
486 xll(ik) = det*(c2*ry(jj) - c1*rx(jj))
487 ENDDO
488 ik = iad + nir + 1
489 lll(ik) = inod(nir+1)
490 jll(ik) = 6
491 sll(ik) = nint
492 xll(ik) = -one
493
495 1 iadll ,lll ,jll ,xll ,ltsm ,
496 2 v ,vr ,a ,ar ,ms ,
497 3 in ,nc_ini ,ncl )
498 IF (comflag==0) THEN
499 iskip = iskip + ncl
500 nc = nc_ini
501 ELSE
502 ic = nc_ini - ncf_s
503 DO k=1,ncl
504 ic = ic + 1
505 icftag(ic) = ic + iskip
506 jcftag(ic+iskip) = nc_ini + k
507 ENDDO
508 ENDIF
509 ENDDO
510
511 RETURN
subroutine lag_direct(iadll, lll, jll, xll, ltsm, v, vr, a, ar, ms, in, nc_ini, ncl)