412
415
416
417
418#include "implicit_f.inc"
419
420
421
422#include "param_c.inc"
423
424
425
426 INTEGER ISECT,INTR,NIP,IG
429 CHARACTER(LEN=NCHARTITLE)::IDTITL
430
431
432
433 INTEGER I,J,IP,NC,NS,IPY,IPZ,IPA
434 my_real ai,yi,zi,wi,phi,dphi,r1,r2,r3,r4,d2, d3, d4
435 my_real w_gauss(9,9),a_gauss(9,9),w_lobatto(9,9),a_lobatto(9,9),len(10)
436
437 DATA w_gauss /
438 1 2. ,0. ,0. ,
439 1 0. ,0. ,0. ,
440 1 0. ,0. ,0. ,
441 2 1. ,1. ,0. ,
442 2 0. ,0. ,0.
443 2 0. ,0. ,0. ,
444 3 0.555555555555556,0.888888888888889,0.555555555555556,
445 3 0. ,0. ,0. ,
446 3 0. ,0. ,0. ,
447 4 0.347854845137454,0.652145154862546,0.652145154862546,
448 4 0.347854845137454,0. ,0. ,
449 4 0. ,0. ,0. ,
450 5 0.236926885056189,0.478628670499366,0.568888888888889,
451 5 0.478628670499366,0.236926885056189,0. ,
452 5 0. ,0. ,0. ,
453 6 0.171324492379170,0.360761573048139,0.467913934572691,
454 6 0.467913934572691,0.360761573048139,0.171324492379170,
455 6 0. ,0. ,0. ,
456 7 0.129484966168870,0.279705391489277,0.381830050505119,
457 7 0.417959183673469,0.381830050505119,0.279705391489277,
458 7 0.129484966168870,0. ,0. ,
459 8 0.101228536290376,0.222381034453374,0.313706645877887,
460 8 0.362683783378362,0.362683783378362,0.313706645877887,
461 8 0.222381034453374,0.101228536290376,0. ,
462 9 0.081274388361574,0.180648160694857,0.260610696402935,
463 9 0.312347077040003,0.330239355001260,0.312347077040003,
464 9 0.260610696402935,0.180648160694857,0.081274388361574/
465 DATA a_gauss /
466 1 0. ,0. ,0. ,
467 1 0. ,0. ,0. ,
468 1 0. ,0. ,0. ,
469 2 -.577350269189626,0.577350269189626,0. ,
470 2 0. ,0. ,0. ,
471 2 0. ,0. ,0. ,
472 3 -.774596669241483,0. ,0.774596669241483,
473 3 0. ,0. ,0. ,
474 3 0. ,0. ,0. ,
475 4 -.861136311594053,-.33998
476 4 0.861136311594053,0. ,0. ,
477 4 0. ,0. ,0. ,
478 5 -.906179845938664,-.538469310105683,0. ,
479 5 0.538469310105683,0.906179845938664,0. ,
480 5 0. ,0. ,0. ,
481 6 -.932469514203152,-.661209386466265,-.238619186083197,
482 6 0.238619186083197,0.661209386466265,0.932469514203152,
483 6 0. ,0. ,0. ,
484 7 -.949107912342759,-.741531185599394,-.405845151377397,
485 7 0. ,0.405845151377397,0.741531185599394,
486 7 0.949107912342759,0. ,0. ,
487 8 -.960289856497536,-.796666477413627,-.525532409916329,
488 8 -.183434642495650,0.183434642495650,0.525532409916329,
489 8 0.796666477413627,0.960289856497536,0. ,
490 9 -.968160239507626,-.836031107326636,-.613371432700590,
491 9 -.324253423403809,0. ,0.324253423403809,
492 9 0.613371432700590,0.836031107326636,0.968160239507626/
493
494 DATA w_lobatto /
495 1 2. ,0. ,0. ,
496 1 0. ,0. ,0. ,
497 1 0. ,0. ,0. ,
498 2 1. ,1. ,0. ,
499 2 0. ,0. ,0. ,
500 2 0. ,0. ,0. ,
501 3 0.333333333333333,1.333333333333333,0.333333333333333,
502 3 0. ,0. ,0. ,
503 3 0. ,0. ,0. ,
504 4 0.166666666666667,0.833333333333333,0.833333333333333,
505 4 0.166666666666667,0. ,0. ,
506 4 0. ,0. ,0. ,
507 5 0.100000000000000,0.544444444444444,0.711111111111111,
508 5 0.544444444444444,0.100000000000000,0. ,
509 5 0. ,0. ,0. ,
510 6 0.066666666666667,0.378474956297847,0.554858377035486,
511 6 0.554858377035486,0.378474956297847,0.066666666666667,
512 6 0. ,0. ,0. ,
513 7 0.047619047619048,0.276826047361566,0.431745381209863,
514 7 0.487619047619048,0.431745381209863,0.276826047361566,
515 7 0.047619047619048,0. ,0. ,
516 8 0.035714285714286,0.210704227143506,0.341122692483504,
517 8 0.412458794658704,0.412458794658704,0.341122692483504,
518 8 0.210704227143506,0.035714285714286,0. ,
519 9 0.027777777777778,0.165495361560806,0.274538712500162,
520 9 0.346428510973046,0.371519274376417,0.346428510973046,
521 9 0.274538712500162,0.165495361560806,0.027777777777778/
522 DATA a_lobatto /
523 1 0. ,0. ,0. ,
524 1 0. ,0. ,0. ,
525 1 0. ,0. ,0. ,
526 2 -1.00000000000000,1.000000000000000,0. ,
527 2 0. ,0. ,0. ,
528 2 0. ,0. ,0. ,
529 3 -1.00000000000000,0. ,1.000000000000000,
530 3 0. ,0. ,0. ,
531 3 0. ,0. ,0. ,
532 4 -1.00000000000000,-.447213595499958,0.447213595499958,
533 4 1.000000000000000,0. ,0. ,
534 4 0. ,0. ,0. ,
535 5 -1.00000000000000,-.654653670707977,0. ,
536 5 0.654653670707977,1.000000000000000,0. ,
537 5 0. ,0. ,0. ,
538 6 -1.00000000000000,-.765055323929465,-.285231516480645,
539 6 0.285231516480645,0.765055323929465,1.000000000000000,
540 6 0. ,0. ,0. ,
541 7 -1.00000000000000,-.830223896278567,-.468848793470714,
542 7 0. ,0.468848793470714,0.830223896278567,
543 7 1.000000000000000,0. ,0. ,
544 8 -1.00000000000000,-.871740148509607,-.591700181433142,
545 8 -.209299217902479,0.209299217902479,0.591700181433142,
546 8 0.871740148509607,1.000000000000000,0. ,
547 9 -1.00000000000000,-.899757995411460,-.677186279510737,
548 9 -.363117463826178,0. ,0.363117463826178,
549 9 0.677186279510737,0.899757995411460,1.000000000000000/
550
551
552 ipy = 200
553 ipz = 300
554 ipa = 400
555
556 SELECT CASE (isect)
557
558 CASE (1)
559
561
562 IF (intr < 1 .OR. intr > 9) THEN
563
564 ELSEIF (intr == 1) THEN
565 nip = intr
566 geo(ipy+1) = zero
567 geo(ipz+1) = zero
569 ELSE
570 nip = intr*intr
571 r1 = l(1)*half
572 r2 = l(2)*half
574 ip = 0
575 DO i = 1,intr
576 DO j = 1,intr
577 ip = ip+1
578 geo(ipy+ip)=a_gauss(i,intr)*r1
579 geo(ipz+ip)=a_gauss(j,intr)*r2
580 geo(ipa+ip)=w_gauss(i,intr)*w_gauss(j,intr)*ai
581 ENDDO
582 ENDDO
583 ENDIF
584
585 CASE (2)
586
588
589 IF (intr < 1 .OR. intr > 9) THEN
590
591 ELSEIF (intr == 1) THEN
592 nip = 1
593 geo(ipy+1) = zero
594 geo(ipz+1) = zero
596 ELSEIF (intr == 2) THEN
597 nip = intr*intr
599 r1 = l(1)/sqr3
600 r1 = l(1)*sqr2*half
601 dphi = two*pi/nip
602 phi = dphi*half
603 DO ip = 1,nip
604 geo(ipy+ip) = r1*sin(phi)
605 geo(ipz+ip) = r1*cos(phi)
606 geo(ipa+ip) = ai
607 phi = phi + dphi
608 ENDDO
609 ELSEIF (intr == 3) THEN
610 nip = intr*intr
612 ip = 1
613 geo(ipy+ip) = zero
614 geo(ipz+ip) = zero
615 geo(ipa+ip) = ai*four
616 r1 = l(1)*sqr3*half
617 dphi = pi*fourth
618 phi = zero
619 DO ip
620 geo(ipy+ip) = r1*sin(phi)
621 geo(ipz+ip) = r1*cos(phi)
622 geo(ipa+ip) = ai
623 phi = phi + dphi
624 ENDDO
625 ELSEIF (intr == 4) THEN
626 nip = 7
627 r1 = l(1)*sqr2/sqr3
628 dphi = pi*third
629 ip = 1
630 geo(ipy+ip) = zero
631 geo(ipz+ip) = zero
632 geo(ipa+ip) =
area*fourth
634 phi = zero
635 DO ip = 2,nip
636 geo(ipy+ip) = r1*sin(phi)
637 geo(ipz+ip) = r1*cos(phi)
638 geo(ipa+ip) = ai
639 phi = phi + dphi
640 ENDDO
641 ELSEIF (intr == 5) THEN
642 nip = 21
643 ip = 1
644 geo(ipy+ip) = zero
645 geo(ipz+ip) = zero
646 geo(ipa+ip) =
area/nine
647 ai =
area*(sixteen + sqr6)/360.
648 r1 = sqrt((six-sqr6)/ten)*l(1)
649 phi = pi/five
650 DO ip = 2,11
651 geo(ipy+ip) = r1*cos(phi*ip)
652 geo(ipz+ip) = r1*sin(phi*ip)
653 geo(ipa+ip) = ai
654 ENDDO
655 ai =
area*(sixteen - sqr6)/360.
656 r1 = sqrt((six+sqr6)/ten)*l(1)
657 DO ip = 12,21
658 geo(ipy+ip) = r1*cos(phi*ip)
659 geo(ipz+ip) = r1*sin(phi*ip)
660 geo(ipa+ip) = ai
661 ENDDO
662 ENDIF
663
664 CASE (3)
665
667
668 IF (intr < 1 .OR. intr > 9) THEN
670 . msgtype=msgerror,
671 . anmode=aninfo_blind_1,
672 . i1=ig,
673 . c1=idtitl)
674 ELSEIF (intr == 1) THEN
675 nip = intr
676 geo(ipy+1) = zero
677 geo(ipz+1) = zero
679 ELSE
680 nip = intr*intr
681 r1 = l(1)*half
682 r2 = l(2)*half
684 ip = 0
685 DO i = 1,intr
686 DO j = 1,intr
687 ip = ip+1
688 geo(ipy+ip)=a_lobatto(i,intr)*r1
689 geo(ipz+ip)=a_lobatto(j,intr)*r2
690 geo(ipa+ip)=w_lobatto(i,intr)*w_lobatto(j,intr)*ai
691 ENDDO
692 ENDDO
693 ENDIF
694
695 CASE (4)
696
697 IF (intr == 17) THEN
698 nip = 17
699 r1 = 0.5477225575*l(1)
700 r2 = 0.8062257748*l(1)
701 r3 = l(1)
702 d2 = a_lobatto(4,5)*l(1)
703 phi = zero
704 dphi= pi * fourth
705 ip = 1
706 geo(ipy+ip) = zero
707 geo(ipz+ip) = zero
708 geo(ipa+ip) = pi*r1*r1
709 ai = pi * (r2*r2 - r1*r1)/eight
710 DO ip = 2,nip-1,2
711 geo(ipy+ip) = d2*cos(phi)
712 geo(ipz+ip) = d2*sin(phi)
713 geo(ipa+ip) = ai
714 phi = phi + dphi
715 ENDDO
716 phi = zero
717 ai = pi * (r3*r3 - r2*r2)/eight
718 DO ip = 3,nip,2
719 geo(ipy+ip) = l(1)*cos(phi)
720 geo(ipz+ip) = l(1)*sin(phi)
721 geo(ipa+ip) = ai
722 phi = phi + dphi
723 ENDDO
724 ELSEIF (intr == 25) THEN
725 nip = 25
726 r1 = 0.46291005*l(1)
727 r2 = 0.69006559*l(1)
728 r3 = 0.859124693*l(1)
729 r4 = l(1)
730 d2 = a_lobatto(5,7)*l(1)
731 d3 = a_lobatto(6,7)*l(1)
732 d4 = a_lobatto(7,7)*l(1)
733 ip = 1
734 geo(ipy+ip) = zero
735 geo(ipz+ip) = zero
736 geo(ipa+ip) = pi*r1*r1
737 phi = zero
738 dphi= pi * fourth
739 ai = pi * (r2*r2 - r1*r1)/eight
740 DO ip = 2,nip-2,3
741 geo(ipy+ip) = d2*cos(phi)
742 geo(ipz+ip) = d2*sin(phi)
743 geo(ipa+ip) = ai
744 phi = phi + dphi
745 ENDDO
746 phi = zero
747 ai = pi * (r3*r3 - r2*r2)/eight
748 DO ip = 3,nip-1,3
749 geo(ipy+ip) = d3*cos(phi)
750 geo(ipz+ip) = d3*sin(phi)
751 geo(ipa+ip) = ai
752 phi = phi + dphi
753 ENDDO
754 phi = zero
755 ai = pi * (r4*r4 - r3*r3)/eight
756 DO ip = 4,nip,3
757 geo(ipy+ip) = d4*cos(phi)
758
759 geo(ipa+ip) = ai
760 phi = phi + dphi
761 ENDDO
762
763
764 ENDIF
765
766
767 CASE (5)
768
770 IF (intr /= 1.AND. intr /= 9 .AND. intr /= 17 ) THEN
772 . msgtype=msgerror,
773 . anmode=aninfo_blind_1,
774 . i1=ig,
775 . c1=idtitl)
776
777 ELSEIF (intr == 1) THEN
778 nip = 1
779 geo(ipy+1) = zero
780 geo(ipz+1) = zero
782 ELSEIF (intr == 9) THEN
783 nip = intr
784 r2 = 0.57346235*l(1)
785 ip = 1
786 geo(ipy+ip) = zero
787 geo(ipz+ip) = zero
788 geo(ipa+ip) = pi*r2*r2
789 r1 = 0.7*l(1)
790 dphi = pi*half
791 phi = zero
792 ai = pi*(l(1)*l(1) - r2*r2)/eight
793 DO ip = 2,nip-1,2
794 geo(ipy+ip) = l(1)*sin(phi)
795 geo(ipz+ip) = l(1)*cos(phi)
796 geo(ipa+ip) = ai
797 phi = phi + dphi
798 ENDDO
799 phi = pi*fourth
800 DO ip = 3,nip,2
801 geo(ipy+ip) = r1*cos(phi)
802 geo
803 geo(ipa+ip) = ai
804 phi = phi + dphi
805 ENDDO
806 ELSEIF (intr == 17) THEN
807 nip = intr
808 r1 = 0.4472136*l(1)
809 r2 = 0.774597 *l(1)
810 r3 = l(1)
811 d2 = half*l(1)
812 ai = pi * (r3*r3 - r2*r2)/eight
813 phi = zero
814 dphi= pi * fourth
815 ip = 1
816 geo(ipy+ip) = zero
817 geo(ipz+ip) = zero
818 geo(ipa+ip) = pi*r1*r1
819 DO ip = 2,nip-1,2
820 geo(ipy+ip) = l(1)*cos(phi)
821 geo(ipz+ip) = l(1)*sin(phi)
822 geo(ipa+ip) = ai
823 phi = phi + dphi
824 ENDDO
825 phi = pi / eight
826 ai = pi * (r2*r2 - r1*r1)/eight
827 DO ip = 3,nip,2
828 geo(ipy+ip) = d2*cos(phi)
829 geo(ipz+ip) = d2*sin(phi)
830 geo(ipa+ip) = ai
831 phi = phi + dphi
832 ENDDO
833 ENDIF
834
835 CASE (6)
836
837 nc = 3+intr
838 ns = 4*nc
839 nip = nc*ns
842 dphi = pi*two/ns
843 r1 = zero
844 r2 = l(1) / sqrt(em20+nc)
845 ip = 0
846 DO i = 1,nc
847 len(i) = (r2 + r1*(sqr3-one)) / sqr3
848 r1 = r2
849 r2 = l(1)*sqrt((i+one)/nc)
850 phi = zero
851 DO j = 1,ns
852 ip = ip+1
853 geo(ipy+ip) = len(i)*sin(phi)
854 geo(ipz+ip) = len(i)*cos(phi)
855 geo(ipa+ip) = ai
856 phi = phi+dphi
858 ENDDO
859 ENDDO
860 CASE DEFAULT
861 END SELECT
862
863 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
integer, parameter nchartitle
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)