OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
soltosph.F File Reference
#include "implicit_f.inc"
#include "sphcom.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine soltosphx8 (nsphdir, ncell, inod, ids, idmax, x, ixs, kxsp, ipartsp, nod2sp, irst)
subroutine soltosphx4 (nsphdir, ncell, inod, ids, idmax, x, ixs, kxsp, ipartsp, nod2sp, irst)
subroutine soltosphv8 (nsphdir, rho, ncell, x, spbuf, ixs, kxsp, ipartsp, irst)
subroutine soltosphv4 (nsphdir, rho, ncell, x, spbuf, ixs)

Function/Subroutine Documentation

◆ soltosphv4()

subroutine soltosphv4 ( integer nsphdir,
rho,
integer ncell,
x,
spbuf,
integer, dimension(nixs) ixs )

Definition at line 588 of file soltosph.F.

591C-----------------------------------------------
592 USE message_mod
593 use element_mod , only : nixs
594C-----------------------------------------------
595C I m p l i c i t T y p e s
596C-----------------------------------------------
597#include "implicit_f.inc"
598C-----------------------------------------------
599C C o m m o n B l o c k s
600C-----------------------------------------------
601#include "sphcom.inc"
602C-----------------------------------------------
603C D u m m y A r g u m e n t s
604C-----------------------------------------------
605 INTEGER NSPHDIR, NCELL, IXS(NIXS)
606 my_real
607 . rho, x(3,*), spbuf(nspbuf,*)
608C-----------------------------------------------
609C L o c a l V a r i a b l e s
610C-----------------------------------------------
611 INTEGER IT,IP,N1, N2, N3, N4, NT, NP
612C
613 my_real
614 . x1,x2,x3,x4,
615 . y1,y2,y3,y4,
616 . z1,z2,z3,z4,
617 . vol, det, wi,
618 . jac1 ,jac2 ,jac3 ,
619 . jac4 ,jac5 ,jac6 ,
620 . jac7 ,jac8 ,jac9 ,
621 . jac_59_68, jac_67_49, jac_48_57
622C-----------------------------------------------
623 np = 0
624 nt = 0
625 DO it=1,nsphdir
626 nt=nt+it
627 np=np+nt
628 END DO
629C-----------------------------------------------
630 n1=ixs(2)
631 x1=x(1,n1)
632 y1=x(2,n1)
633 z1=x(3,n1)
634 n2=ixs(4)
635 x2=x(1,n2)
636 y2=x(2,n2)
637 z2=x(3,n2)
638 n3=ixs(7)
639 x3=x(1,n3)
640 y3=x(2,n3)
641 z3=x(3,n3)
642 n4=ixs(6)
643 x4=x(1,n4)
644 y4=x(2,n4)
645 z4=x(3,n4)
646C
647C Jacobian matrix
648 jac1=x1-x4
649 jac2=y1-y4
650 jac3=z1-z4
651 jac4=x2-x4
652 jac5=y2-y4
653 jac6=z2-z4
654 jac7=x3-x4
655 jac8=y3-y4
656 jac9=z3-z4
657C
658 jac_59_68=jac5*jac9-jac6*jac8
659 jac_67_49=jac6*jac7-jac4*jac9
660 jac_48_57=jac4*jac8-jac5*jac7
661C
662 det=jac1*jac_59_68+jac2*jac_67_49+jac3*jac_48_57
663C-----------------------------------------------
664 DO ip=1,ncell
665C
666 wi = one/(six*np)
667C
668 vol= wi*det
669 IF(det<zero)
670 . CALL ancmsg(msgid=1038,
671 . msgtype=msgerror,
672 . anmode=aninfo,
673 . i1=ixs(nixs))
674C
675C Particle volume will be replaced by particle mass later (spinit3)
676 spbuf(2,ip) =vol*rho
677 spbuf(12,ip)=vol
678C
679 ENDDO
680C-----------------------------------------------
681 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
message_mod
Definition message_mod.F:1257
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:895

◆ soltosphv8()

subroutine soltosphv8 ( integer nsphdir,
rho,
integer ncell,
x,
spbuf,
integer, dimension(nixs) ixs,
integer, dimension(nisp,*) kxsp,
integer, dimension(*) ipartsp,
integer, dimension(3,*) irst )

Definition at line 335 of file soltosph.F.

338C-----------------------------------------------
339C M o d u l e s
340C-----------------------------------------------
341 USE message_mod
342 use element_mod , only : nixs
343C-----------------------------------------------
344C I m p l i c i t T y p e s
345C-----------------------------------------------
346#include "implicit_f.inc"
347C-----------------------------------------------
348C G l o b a l P a r a m e t e r s
349C-----------------------------------------------
350C C o m m o n B l o c k s
351C-----------------------------------------------
352#include "sphcom.inc"
353C-----------------------------------------------
354C D u m m y A r g u m e n t s
355C-----------------------------------------------
356 INTEGER NSPHDIR, NCELL, IXS(NIXS), KXSP(NISP,*),
357 . IPARTSP(*), IRST(3,*)
358 my_real
359 . rho, x(3,*), spbuf(nspbuf,*)
360C-----------------------------------------------
361C L o c a l V a r i a b l e s
362C-----------------------------------------------
363 INTEGER I, J, IR, IS, IT, IP,
364 . N1, N2, N3, N4, N5, N6, N7, N8, NP
365C
366 my_real
367 . x1,x2,x3,x4,x5,x6,x7,x8,
368 . y1,y2,y3,y4,y5,y6,y7,y8,
369 . z1,z2,z3,z4,z5,z6,z7,z8,
370 . x17 , x28 , x35 , x46 ,
371 . y17 , y28 , y35 , y46 ,
372 . z17 , z28 , z35 , z46 ,
373 . vol, hx(4), hy(4), hz(4), det,
374 . jac1 ,jac2 ,jac3 ,
375 . jac4 ,jac5 ,jac6 ,
376 . jac7 ,jac8 ,jac9 ,
377 . cj1 ,cj2 ,cj3 ,
378 . cj4 ,cj5 ,cj6 ,
379 . cj7 ,cj8 ,cj9 ,
380 . jac_59_68, jac_67_49, jac_48_57,
381 . jac_38_29, jac_19_37, jac_27_18,
382 . jac_26_35, jac_34_16, jac_15_24,
383 . x_17_46 , x_28_35 ,
384 . y_17_46 , y_28_35 ,
385 . z_17_46 , z_28_35 ,
386 . ksi, eta, zeta, wi,
387 . xi, yi, zi
388C-----------------------------------------------
389 my_real
390 . w_gauss(9,9),a_gauss(9,9)
391 DATA w_gauss /
392 1 2. ,0. ,0. ,
393 1 0. ,0. ,0. ,
394 1 0. ,0. ,0. ,
395 2 1. ,1. ,0. ,
396 2 0. ,0. ,0. ,
397 2 0. ,0. ,0. ,
398 3 0.555555555555556,0.888888888888889,0.555555555555556,
399 3 0. ,0. ,0. ,
400 3 0. ,0. ,0. ,
401 4 0.347854845137454,0.652145154862546,0.652145154862546,
402 4 0.347854845137454,0. ,0. ,
403 4 0. ,0. ,0. ,
404 5 0.236926885056189,0.478628670499366,0.568888888888889,
405 5 0.478628670499366,0.236926885056189,0. ,
406 5 0. ,0. ,0. ,
407 6 0.171324492379170,0.360761573048139,0.467913934572691,
408 6 0.467913934572691,0.360761573048139,0.171324492379170,
409 6 0. ,0. ,0. ,
410 7 0.129484966168870,0.279705391489277,0.381830050505119,
411 7 0.417959183673469,0.381830050505119,0.279705391489277,
412 7 0.129484966168870,0. ,0. ,
413 8 0.101228536290376,0.222381034453374,0.313706645877887,
414 8 0.362683783378362,0.362683783378362,0.313706645877887,
415 8 0.222381034453374,0.101228536290376,0. ,
416 9 0.081274388361574,0.180648160694857,0.260610696402935,
417 9 0.312347077040003,0.330239355001260,0.312347077040003,
418 9 0.260610696402935,0.180648160694857,0.081274388361574/
419 DATA a_gauss /
420 1 0. ,0. ,0. ,
421 1 0. ,0. ,0. ,
422 1 0. ,0. ,0. ,
423 2 -.5 ,0.5 ,0. ,
424 2 0. ,0. ,0. ,
425 2 0. ,0. ,0. ,
426 3 -.666666666666666,0. ,0.666666666666666,
427 3 0. ,0. ,0. ,
428 3 0. ,0. ,0. ,
429 4 -.75 ,-.25 ,0.25 ,
430 4 0.75 ,0. ,0. ,
431 4 0. ,0. ,0. ,
432 5 -.8 ,-.4 ,0. ,
433 5 0.4 ,0.8 ,0. ,
434 5 0. ,0. ,0. ,
435 6 -.833333333333333,-.5 ,-.166666666666666,
436 6 0.166666666666666,0.5 ,0.833333333333333,
437 6 0. ,0. ,0. ,
438 7 -.857142857142857,-.571428571428571,-.285714285714285,
439 7 0. ,0.285714285714285,0.571428571428571,
440 7 0.857142857142857,0. ,0. ,
441 8 -.875 ,-.625 ,-.375 ,
442 8 -.125 ,0.125 ,0.375,
443 8 0.625 ,0.875 ,0. ,
444 9 -.888888888888888,-.666666666666666,-.444444444444444,
445 9 -.222222222222222,0. ,0.222222222222222,
446 9 0.444444444444444,0.666666666666666,0.888888888888888/
447C-----------------------------------------------
448 np = nsphdir*nsphdir*nsphdir
449C
450 n1=ixs(2)
451 x1=x(1,n1)
452 y1=x(2,n1)
453 z1=x(3,n1)
454 n2=ixs(3)
455 x2=x(1,n2)
456 y2=x(2,n2)
457 z2=x(3,n2)
458 n3=ixs(4)
459 x3=x(1,n3)
460 y3=x(2,n3)
461 z3=x(3,n3)
462 n4=ixs(5)
463 x4=x(1,n4)
464 y4=x(2,n4)
465 z4=x(3,n4)
466 n5=ixs(6)
467 x5=x(1,n5)
468 y5=x(2,n5)
469 z5=x(3,n5)
470 n6=ixs(7)
471 x6=x(1,n6)
472 y6=x(2,n6)
473 z6=x(3,n6)
474 n7=ixs(8)
475 x7=x(1,n7)
476 y7=x(2,n7)
477 z7=x(3,n7)
478 n8=ixs(9)
479 x8=x(1,n8)
480 y8=x(2,n8)
481 z8=x(3,n8)
482C
483 x17=x7-x1
484 x28=x8-x2
485 x35=x5-x3
486 x46=x6-x4
487 y17=y7-y1
488 y28=y8-y2
489 y35=y5-y3
490 y46=y6-y4
491 z17=z7-z1
492 z28=z8-z2
493 z35=z5-z3
494 z46=z6-z4
495C
496C Jacobian matrix
497 cj4=x17+x28-x35-x46
498 cj5=y17+y28-y35-y46
499 cj6=z17+z28-z35-z46
500 x_17_46=x17+x46
501 x_28_35=x28+x35
502 y_17_46=y17+y46
503 y_28_35=y28+y35
504 z_17_46=z17+z46
505 z_28_35=z28+z35
506C
507 cj7=x_17_46+x_28_35
508 cj8=y_17_46+y_28_35
509 cj9=z_17_46+z_28_35
510 cj1=x_17_46-x_28_35
511 cj2=y_17_46-y_28_35
512 cj3=z_17_46-z_28_35
513C Hourglass
514C mode 1
515C 1 1 -1 -1 -1 -1 1 1
516 hx(1)=(x1+x2-x3-x4-x5-x6+x7+x8)
517 hy(1)=(y1+y2-y3-y4-y5-y6+y7+y8)
518 hz(1)=(z1+z2-z3-z4-z5-z6+z7+z8)
519C mode 2
520C 1 -1 -1 1 -1 1 1 -1
521 hx(2)=(x1-x2-x3+x4-x5+x6+x7-x8)
522 hy(2)=(y1-y2-y3+y4-y5+y6+y7-y8)
523 hz(2)=(z1-z2-z3+z4-z5+z6+z7-z8)
524C mode 3
525C 1 -1 1 -1 1 -1 1 -1
526 hx(3)=(x1-x2+x3-x4+x5-x6+x7-x8)
527 hy(3)=(y1-y2+y3-y4+y5-y6+y7-y8)
528 hz(3)=(z1-z2+z3-z4+z5-z6+z7-z8)
529C mode 4
530C -1 1 -1 1 1 -1 1 -1
531 hx(4)=(-x1+x2-x3+x4+x5-x6+x7-x8)
532 hy(4)=(-y1+y2-y3+y4+y5-y6+y7-y8)
533 hz(4)=(-z1+z2-z3+z4+z5-z6+z7-z8)
534C-----------------------------------------------
535 DO ip=1,ncell
536 ir=irst(1,ip)
537 is=irst(2,ip)
538 it=irst(3,ip)
539C
540 ksi = a_gauss(it,nsphdir)
541 eta = a_gauss(ir,nsphdir)
542 zeta = a_gauss(is,nsphdir)
543C
544 wi = eight/np
545C
546C Jacobian matrix
547 jac1=cj1+hx(3)*eta+(hx(2)+hx(4)*eta)*zeta
548 jac2=cj2+hy(3)*eta+(hy(2)+hy(4)*eta)*zeta
549 jac3=cj3+hz(3)*eta+(hz(2)+hz(4)*eta)*zeta
550C
551 jac4=cj4+hx(1)*zeta+(hx(3)+hx(4)*zeta)*ksi
552 jac5=cj5+hy(1)*zeta+(hy(3)+hy(4)*zeta)*ksi
553 jac6=cj6+hz(1)*zeta+(hz(3)+hz(4)*zeta)*ksi
554C
555 jac7=cj7+hx(2)*ksi+(hx(1)+hx(4)*ksi)*eta
556 jac8=cj8+hy(2)*ksi+(hy(1)+hy(4)*ksi)*eta
557 jac9=cj9+hz(2)*ksi+(hz(1)+hz(4)*ksi)*eta
558C
559 jac_59_68=jac5*jac9-jac6*jac8
560 jac_67_49=jac6*jac7-jac4*jac9
561 jac_48_57=jac4*jac8-jac5*jac7
562C
563 det=one_over_512*(jac1*jac_59_68+jac2*jac_67_49+jac3*jac_48_57)
564 vol= wi*det
565 IF(det<zero)
566 . CALL ancmsg(msgid=1038,
567 . msgtype=msgerror,
568 . anmode=aninfo,
569 . i1=ixs(nixs))
570C
571C Particle volume will be replaced by particle mass later (spinit3)
572 spbuf(2,ip) =vol*rho
573 spbuf(12,ip)=vol
574C
575 ENDDO
576C-----------------------------------------------
577 RETURN

◆ soltosphx4()

subroutine soltosphx4 ( integer nsphdir,
integer ncell,
integer inod,
integer ids,
integer idmax,
x,
integer, dimension(nixs) ixs,
integer, dimension(nisp,*) kxsp,
integer, dimension(*) ipartsp,
integer, dimension(*) nod2sp,
integer, dimension(3,*) irst )

Definition at line 186 of file soltosph.F.

190C-----------------------------------------------
191C M o d u l e s
192C-----------------------------------------------
193 USE message_mod
194 use element_mod , only : nixs
195C-----------------------------------------------
196C I m p l i c i t T y p e s
197C-----------------------------------------------
198#include "implicit_f.inc"
199C-----------------------------------------------
200C G l o b a l P a r a m e t e r s
201C-----------------------------------------------
202C C o m m o n B l o c k s
203C-----------------------------------------------
204#include "sphcom.inc"
205C-----------------------------------------------
206C D u m m y A r g u m e n t s
207C-----------------------------------------------
208 INTEGER NSPHDIR, NCELL, INOD, IDS, IDMAX, IXS(NIXS),
209 . KXSP(NISP,*), IPARTSP(*), NOD2SP(*), IRST(3,*)
210 my_real
211 . x(3,*)
212C-----------------------------------------------
213C L o c a l V a r i a b l e s
214C-----------------------------------------------
215 INTEGER IR, IS, IT,N1,N2,N3,N4
216C
217 my_real
218 . x1,x2,x3,x4,y1,y2,y3,y4,
219 . z1,z2,z3,z4,phi1,phi2,phi3,phi4,ksi,
220 . eta,zeta,xi,yi,zi,a_gauss_tetra(9,9)
221C-----------------------------------------------
222C E x t e r n a l F u n c t i o n s
223C-----------------------------------------------
224 my_real
225 . checkvolume_4n
226C-----------------------------------------------
227C A_GAUSS Generated with (2*IR-1)/(2(NSPHDIR+1))
228C-----------------------------------------------
229 DATA a_gauss_tetra /
230 1 0.250000000000000,0.000000000000000,0.000000000000000,
231 1 0.000000000000000,0.000000000000000,0.000000000000000,
232 1 0.000000000000000,0.000000000000000,0.000000000000000,
233 2 0.166666666666667,0.500000000000000,0.000000000000000,
234 2 0.000000000000000,0.000000000000000,0.000000000000000,
235 2 0.000000000000000,0.000000000000000,0.000000000000000,
236 3 0.125000000000000,0.375000000000000,0.625000000000000,
237 3 0.000000000000000,0.000000000000000,0.000000000000000,
238 3 0.000000000000000,0.000000000000000,0.000000000000000,
239 4 0.100000000000000,0.300000000000000,0.500000000000000,
240 4 0.700000000000000,0.000000000000000,0.000000000000000,
241 4 0.000000000000000,0.000000000000000,0.000000000000000,
242 5 0.083333333333333,0.250000000000000,0.416666666666667,
243 5 0.583333333333333,0.750000000000000,0.000000000000000,
244 5 0.000000000000000,0.000000000000000,0.000000000000000,
245 6 0.071428571428571,0.214285714285714,0.357142857142857,
246 6 0.500000000000000,0.642857142857143,0.785714285714286,
247 6 0.000000000000000,0.000000000000000,0.000000000000000,
248 7 0.062500000000000,0.187500000000000,0.312500000000000,
249 7 0.437500000000000,0.562500000000000,0.687500000000000,
250 7 0.812500000000000,0.000000000000000,0.000000000000000,
251 8 0.055555555555556,0.166666666666667,0.277777777777778,
252 8 0.388888888888889,0.500000000000000,0.611111111111111,
253 8 0.722222222222222,0.833333333333333,0.000000000000000,
254 9 0.050000000000000,0.150000000000000,0.250000000000000,
255 9 0.350000000000000,0.450000000000000,0.550000000000000,
256 9 0.650000000000000,0.750000000000000,0.850000000000000/
257C-----------------------------------------------
258C---- KSI - R - N4->N1
259C---- ETA - S - N4->N2
260C---- ZETA- T - N4->N3
261C-----------------------------------------------
262 n1=ixs(2)
263 x1=x(1,n1)
264 y1=x(2,n1)
265 z1=x(3,n1)
266 n2=ixs(4)
267 x2=x(1,n2)
268 y2=x(2,n2)
269 z2=x(3,n2)
270 n3=ixs(7)
271 x3=x(1,n3)
272 y3=x(2,n3)
273 z3=x(3,n3)
274 n4=ixs(6)
275 x4=x(1,n4)
276 y4=x(2,n4)
277 z4=x(3,n4)
278
279C------------------------------------------------
280C Renumber tetrahedron in case of negative volume
281 IF (checkvolume_4n(x ,ixs(1)) < zero) THEN
282 n2=ixs(6)
283 n4=ixs(4)
284 ENDIF
285C-----------------------------------------------
286 DO ir=1,nsphdir
287 DO is=1,nsphdir-ir+1
288 DO it=1,nsphdir-is-ir+2
289C
290 ksi = a_gauss_tetra(ir,nsphdir)
291 eta = a_gauss_tetra(is,nsphdir)
292 zeta = a_gauss_tetra(it,nsphdir)
293C
294 phi1=ksi
295 phi2=eta
296 phi3=zeta
297 phi4=one-ksi-eta-zeta
298C
299 xi=phi1*x1+phi2*x2+phi3*x3+phi4*x4
300 yi=phi1*y1+phi2*y2+phi3*y3+phi4*y4
301 zi=phi1*z1+phi2*z2+phi3*z3+phi4*z4
302C
303 ncell=ncell+1
304 ipartsp(ncell)=ids
305 inod =inod+1
306 kxsp(3,ncell) =inod
307 x(1,inod)=xi
308 x(2,inod)=yi
309 x(3,inod)=zi
310 nod2sp(inod) =ncell
311 kxsp(2,ncell)=-1
312 idmax=idmax+1
313 kxsp(nisp,ncell)=idmax
314 irst(1,ncell-first_sphsol+1)=ir
315 irst(2,ncell-first_sphsol+1)=is
316 irst(3,ncell-first_sphsol+1)=it
317C
318 ENDDO
319 ENDDO
320 ENDDO
321C
322C-----------------------------------------------
323 RETURN
checkvolume_4n
function checkvolume_4n(x, ixs)
Definition checksvolume.F:217

◆ soltosphx8()

subroutine soltosphx8 ( integer nsphdir,
integer ncell,
integer inod,
integer ids,
integer idmax,
x,
integer, dimension(nixs) ixs,
integer, dimension(nisp,*) kxsp,
integer, dimension(*) ipartsp,
integer, dimension(*) nod2sp,
integer, dimension(3,*) irst )

Definition at line 30 of file soltosph.F.

34C-----------------------------------------------
35C M o d u l e s
36C-----------------------------------------------
37 USE message_mod
38 use element_mod , only : nixs
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C G l o b a l P a r a m e t e r s
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "sphcom.inc"
49C-----------------------------------------------
50C D u m m y A r g u m e n t s
51C-----------------------------------------------
52 INTEGER NSPHDIR, NCELL, INOD, IDS, IDMAX, IXS(NIXS),
53 . KXSP(NISP,*), IPARTSP(*), NOD2SP(*), IRST(3,*)
54 my_real
55 . x(3,*)
56C-----------------------------------------------
57C L o c a l V a r i a b l e s
58C-----------------------------------------------
59 INTEGER I, J, IR, IS, IT,
60 . N1, N2, N3, N4, N5, N6, N7, N8
61C 12
62 my_real
63 . x1,x2,x3,x4,x5,x6,x7,x8,
64 . y1,y2,y3,y4,y5,y6,y7,y8,
65 . z1,z2,z3,z4,z5,z6,z7,z8,
66 . phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
67 . ksi, eta, zeta, wi,
68 . xi, yi, zi
69C-----------------------------------------------
70 my_real
71 . a_gauss(9,9)
72 DATA a_gauss /
73 1 0. ,0. ,0. ,
74 1 0. ,0. ,0. ,
75 1 0. ,0. ,0. ,
76 2 -.5 ,0.5 ,0. ,
77 2 0. ,0. ,0. ,
78 2 0. ,0. ,0. ,
79 3 -.666666666666666,0. ,0.666666666666666,
80 3 0. ,0. ,0. ,
81 3 0. ,0. ,0. ,
82 4 -.75 ,-.25 ,0.25 ,
83 4 0.75 ,0. ,0. ,
84 4 0. ,0. ,0. ,
85 5 -.8 ,-.4 ,0. ,
86 5 0.4 ,0.8 ,0. ,
87 5 0. ,0. ,0. ,
88 6 -.833333333333333,-.5 ,-.166666666666666,
89 6 0.166666666666666,0.5 ,0.833333333333333,
90 6 0. ,0. ,0. ,
91 7 -.857142857142857,-.571428571428571,-.285714285714285,
92 7 0. ,0.285714285714285,0.571428571428571,
93 7 0.857142857142857,0. ,0. ,
94 8 -.875 ,-.625 ,-.375 ,
95 8 -.125 ,0.125 ,0.375,
96 8 0.625 ,0.875 ,0. ,
97 9 -.888888888888888,-.666666666666666,-.444444444444444,
98 9 -.222222222222222,0. ,0.222222222222222,
99 9 0.444444444444444,0.666666666666666,0.888888888888888/
100C-----------------------------------------------
101 n1=ixs(2)
102 x1=x(1,n1)
103 y1=x(2,n1)
104 z1=x(3,n1)
105 n2=ixs(3)
106 x2=x(1,n2)
107 y2=x(2,n2)
108 z2=x(3,n2)
109 n3=ixs(4)
110 x3=x(1,n3)
111 y3=x(2,n3)
112 z3=x(3,n3)
113 n4=ixs(5)
114 x4=x(1,n4)
115 y4=x(2,n4)
116 z4=x(3,n4)
117 n5=ixs(6)
118 x5=x(1,n5)
119 y5=x(2,n5)
120 z5=x(3,n5)
121 n6=ixs(7)
122 x6=x(1,n6)
123 y6=x(2,n6)
124 z6=x(3,n6)
125 n7=ixs(8)
126 x7=x(1,n7)
127 y7=x(2,n7)
128 z7=x(3,n7)
129 n8=ixs(9)
130 x8=x(1,n8)
131 y8=x(2,n8)
132 z8=x(3,n8)
133C-----------------------------------------------
134 DO ir=1,nsphdir
135 DO is=1,nsphdir
136 DO it=1,nsphdir
137 ksi = a_gauss(ir,nsphdir)
138 eta = a_gauss(is,nsphdir)
139 zeta = a_gauss(it,nsphdir)
140C
141 phi1=(one-ksi)*(one-eta)*(one-zeta)
142 phi2=(one-ksi)*(one-eta)*(one+zeta)
143 phi3=(one+ksi)*(one-eta)*(one+zeta)
144 phi4=(one+ksi)*(one-eta)*(one-zeta)
145 phi5=(one-ksi)*(one+eta)*(one-zeta)
146 phi6=(one-ksi)*(one+eta)*(one+zeta)
147 phi7=(one+ksi)*(one+eta)*(one+zeta)
148 phi8=(one+ksi)*(one+eta)*(one-zeta)
149 xi=one_over_8*(phi1*x1+phi2*x2+phi3*x3+phi4*x4+
150 . phi5*x5+phi6*x6+phi7*x7+phi8*x8)
151 yi=one_over_8*(phi1*y1+phi2*y2+phi3*y3+phi4*y4+
152 . phi5*y5+phi6*y6+phi7*y7+phi8*y8)
153 zi=one_over_8*(phi1*z1+phi2*z2+phi3*z3+phi4*z4+
154 . phi5*z5+phi6*z6+phi7*z7+phi8*z8)
155C
156 ncell=ncell+1
157 ipartsp(ncell)=ids
158 inod =inod+1
159 kxsp(3,ncell) =inod
160 x(1,inod)=xi
161 x(2,inod)=yi
162 x(3,inod)=zi
163 nod2sp(inod) =ncell
164 kxsp(2,ncell)=-1
165 idmax=idmax+1
166 kxsp(nisp,ncell)=idmax
167 irst(1,ncell-first_sphsol+1)=ir
168 irst(2,ncell-first_sphsol+1)=is
169 irst(3,ncell-first_sphsol+1)=it
170C
171 ENDDO
172 ENDDO
173 ENDDO
174C-----------------------------------------------
175 RETURN