412
415
416
417
418#include "implicit_f.inc"
419
420
421
422#include "param_c.inc"
423
424
425
426 INTEGER ISECT,INTR,,IG
429 CHARACTER(LEN=NCHARTITLE)::IDTITL
430
431
432
433 INTEGER I,J,IP,NC,NS,,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.
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,-.339981043584856,0.339981043584856,
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 = 2,nip
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)
697
698 IF (intr == 17) THEN
699 nip = 17
700 r1 = 0.5477225575*l(1)
701 r2 = 0.8062257748*l(1)
702 r3 = l(1)
703 d2 = a_lobatto(4,5)*l(1)
704 phi = zero
705 dphi= pi * fourth
706 ip = 1
707 geo(ipy+ip) = zero
708 geo(ipz+ip) = zero
709 geo(ipa+ip) = pi*r1*r1
710 ai = pi * (r2*r2 - r1*r1)/eight
711 DO ip = 2,nip-1,2
712 geo(ipy+ip) = d2*cos(phi)
713 geo(ipz+ip) = d2*sin(phi)
714 geo(ipa+ip) = ai
715 phi = phi + dphi
716 ENDDO
717 phi = zero
718 ai = pi * (r3*r3 - r2*r2)/eight
719 DO ip = 3,nip,2
720 geo(ipy+ip) = l(1)*cos(phi)
721 geo(ipz+ip) = l(1)*sin(phi)
722 geo(ipa+ip) = ai
723 phi = phi + dphi
724 ENDDO
725 ELSEIF (intr == 25) THEN
726 nip = 25
727 r1 = 0.46291005*l(1)
728 r2 = 0.69006559*l(1)
729 r3 = 0.859124693*l(1)
730 r4 = l(1)
731 d2 = a_lobatto(5,7)*l(1)
732 d3 = a_lobatto(6,7)*l(1)
733 d4 = a_lobatto(7,7)*l(1)
734 ip = 1
735 geo(ipy+ip) = zero
736 geo(ipz+ip) = zero
737 geo(ipa+ip) = pi*r1*r1
738 phi = zero
739 dphi= pi * fourth
740 ai = pi * (r2*r2 - r1*r1)/eight
741 DO ip = 2,nip-2,3
742 geo(ipy+ip) = d2*cos(phi)
743 geo(ipz+ip) = d2*sin(phi)
744 geo(ipa+ip) = ai
745 phi = phi + dphi
746 ENDDO
747 phi = zero
748 ai = pi * (r3*r3 - r2*r2)/eight
749 DO ip = 3,nip-1,3
750 geo(ipy+ip) = d3*cos(phi)
751 geo(ipz+ip) = d3*sin(phi)
752 geo(ipa+ip) = ai
753 phi = phi + dphi
754 ENDDO
755 phi = zero
756 ai = pi * (r4*r4 - r3*r3)/eight
757 DO ip = 4,nip,3
758 geo(ipy+ip) = d4*cos(phi)
759 geo(ipz+ip) = d4*sin(phi)
760 geo(ipa+ip) = ai
761 phi = phi + dphi
762 ENDDO
763
764
765 ENDIF
766
767
768 CASE (5)
769
771 IF (intr /= 1.AND. intr /= 9 .AND. intr /= 17 ) THEN
773 . msgtype=msgerror,
774 . anmode=aninfo_blind_1,
775 . i1=ig,
776 . c1=idtitl)
777
778 ELSEIF (intr == 1) THEN
779 nip = 1
780 geo(ipy+1) = zero
781 geo(ipz+1) = zero
783 ELSEIF (intr == 9) THEN
784 nip = intr
785 r2 = 0.57346235*l(1)
786 ip = 1
787 geo(ipy+ip) = zero
788 geo(ipz+ip) = zero
789 geo(ipa+ip) = pi*r2*r2
790 r1 = 0.7*l(1)
791 dphi = pi*half
792 phi = zero
793 ai = pi*(l(1)*l(1) - r2*r2)/eight
794 DO ip = 2,nip-1,2
795 geo(ipy+ip) = l(1)*sin(phi)
796 geo(ipz+ip) = l(1)*cos(phi)
797 geo(ipa+ip) = ai
798 phi = phi + dphi
799 ENDDO
800 phi = pi*fourth
801 DO ip = 3,nip,2
802 geo(ipy+ip) = r1*cos(phi)
803 geo(ipz+ip) = r1*sin(phi)
804 geo(ipa+ip) = ai
805 phi = phi + dphi
806 ENDDO
807 ELSEIF (intr == 17) THEN
808 nip = intr
809 r1 = 0.4472136*l(1)
810 r2 = 0.774597 *l(1)
811 r3 = l(1)
812 d2 = half*l(1)
813 ai = pi * (r3*r3 - r2*r2)/eight
814 phi = zero
815 dphi= pi * fourth
816 ip = 1
817 geo(ipy+ip) = zero
818 geo(ipz+ip) = zero
819 geo(ipa+ip) = pi*r1*r1
820 DO ip = 2,nip-1,2
821 geo(ipy+ip) = l(1)*cos(phi)
822 geo(ipz+ip) = l(1)*sin(phi)
823 geo(ipa+ip) = ai
824 phi = phi + dphi
825 ENDDO
826 phi = pi / eight
827 ai = pi * (r2*r2 - r1*r1)/eight
828 DO ip = 3,nip,2
829 geo(ipy+ip) = d2*cos(phi)
830 geo(ipz+ip) = d2*sin(phi)
831 geo(ipa+ip) = ai
832 phi = phi + dphi
833 ENDDO
834 ENDIF
835
836 CASE (6)
837
838 nc = 3+intr
839 ns = 4*nc
840 nip = nc*ns
843 dphi = pi*two/ns
844 r1 = zero
845 r2 = l(1) / sqrt(em20+nc)
846 ip = 0
847 DO i = 1,nc
848 len
849 r1 = r2
850 r2 = l(1)*sqrt((i+one)/nc)
851 phi = zero
852 DO j = 1,ns
853 ip = ip+1
854 geo(ipy+ip) = len(i)*sin(phi)
855 geo(ipz+ip) = len(i)*cos(phi)
856 geo(ipa+ip) = ai
857 phi = phi+dphi
859 ENDDO
860 ENDDO
861 CASE DEFAULT
862 END SELECT
863
864 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)