C PAVPOI SOURCE BP208322 15/06/22 21:21:10 8543 SUBROUTINE PAVPOI(IWRK,IELE,NDEC,NNO) **************************************************************************** * appelee par prom projection par minimisation d une integrale * FONCTION : * distribue une collection de point d integration dans * les elements de reference et les rend a prom * ENTREES : * NDEC nomdre de points dans une direction ( isotrope pour l instant * * iele numero du support geomtrique de l element * NNO nombre de points de l element ( utile pour polygone uniquement) * SORTIES : * Tableau XX dans le segment iwrk et NDD ***************************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC CCREEL *- PARAMETER ( O1=1.D0,O2=2.D0) SEGMENT MTR REAL*8 G(3,6),ppd(6) ENDSEGMENT SEGMENT WRK5 REAL*8 W(200),A(NBN1,NBN1),B(NBN1,NBN1),XPU(3),RHS(NBN1) REAL*8 SHPXXX(6,NBNN,200),XX(3,200),VAL(NBN1),POIDS(200) INTEGER NDD ENDSEGMENT C WRK5=IWRK ip = 0 C write(6,*)' xx(/1) xx(/2)',xx(/1) ,xx(/2) C goto (99,2, 2,4,4,4,4,8,8,8,8,100,100,8,8,4,4,99,99, & 99,99,99,23,23,23,23,99,99,99,99,99,32,99,99,99,99,99, & 99),iele C ELEMENT DE REFERENCE : SEG2 - SEG3 C ------------------------ 2 CONTINUE delx=O2/ndec x=-O1-delx/O2 DO i=1,ndec xx(1,i)=x+delx*i poids(i)=delx ENDDO ip=ndec GOTO 100 4 CONTINUE C ELEMENT DE REFERENCE TRIANGLES + BASE PRISME ppp = O1/(O2*ndec*ndec) r= O1/(3.D0*ndec) p=O1/ndec x0 =-p+r y0 =x0 x1=x0+r y1=x1 ppp = O1/(ndec*ndec*O2) C do k=1,ndec do i=1,ndec-k+1 ip=ip+1 xx(1,ip) = x0+p*i xx(2,ip) = y0+p*k poids(ip)= ppp if(i.eq.ndec-k+1) goto 50 ip=ip+1 xx(1,ip) = x1+p*i xx(2,ip) = y1+p*k poids(ip)= ppp 50 continue enddo enddo if(iele.eq.16.or.iele.eq.17) goto 16 GOTO 100 8 CONTINUE C ELEMENTS DE REFERENCE QUADRANGLE + BASE CUBE delx = O2/ndec dely=delx ppp = O2*O2/(ndec*ndec) x= -O1 - delx/O2 y= -O1 + dely/O2 do i=1,ndec do j=1,ndec ip=ip+1 xx(1,ip)=x+delx*j xx(2,ip)=y poids(ip) = ppp enddo y = y + dely enddo if(iele.eq.14.or.iele.eq.15) goto 14 GOTO 100 C 14 CONTINUE C ELEMENTS DE REFERENCE CUBES delz = delx ppp = ppp*delz z0= -O1+delz/O2 C do i=1,ndec*ndec xx(3,i)=z0 poids(i)=ppp enddo C ip= ndec*ndec do k=1,ndec-1 do i=1,ndec*ndec ip=ip+1 xx(1,ip)=xx(1,i) xx(2,ip)=xx(2,i) xx(3,ip)=xx(3,i)+delz*k poids(ip)=ppp enddo enddo GOTO 100 C 16 CONTINUE C ELEMENTS DE REFRENCE PRISMES delz = O2/ndec z0= -O1+delz/O2 ppp = ppp*delz C do i=1,ndec*ndec xx(3,i)=z0 poids(i)=ppp enddo C ip= ndec*ndec do k=1,ndec-1 do i=1,ndec*ndec ip=ip+1 xx(1,ip)=xx(1,i) xx(2,ip)=xx(2,i) xx(3,ip)=xx(3,i)+delz*k poids(ip)=ppp enddo enddo GOTO 100 23 CONTINUE C ELEMENTS DE REFERENCE TETRAEDRES p =O1/ndec unq= O1/(O2*O2) und=O1/O2 ppp = O1/(6.D0*ndec*ndec*ndec) segini mtr g(1,1) = unq g(2,1) = unq g(3,1) = unq g(1,2) = O1-unq g(2,2) = unq g(3,2) = und g(1,3) = unq g(2,3) = O1-unq g(3,3) = und g(1,4) = und g(2,4) = und g(3,4) = O1-unq g(1,5) = und g(2,5) = und g(3,5) = unq g(1,6) = O1-unq g(2,6) = O1-unq g(3,6) = O1-unq do i=1,6 do j=1,3 g(j,i)=g(j,i)/ndec-p enddo enddo do k=1,ndec tz = p*k do j=1,ndec-k+1 ty = p*j do i=1,ndec-j+1-k+1 tx = p*i C do 500 l=1,6 x=g(1,l)+tx y=g(2,l)+ty z=g(3,l)+tz if((x+y+z).GT.O1) goto 500 ip = ip+1 xx(1,ip)=x xx(2,ip)=y xx(3,ip)=z poids(ip) = ppp 500 continue enddo enddo enddo segsup mtr IF (iele.eq.25.or.iele.eq.26 ) goto 25 GOTO 100 25 CONTINUE C ELEMENTS DE REFERENCE PYRAMIDES ( decoupee en 4 tetraedres) ipo=ip do ik=1,ipo xx(1,ip+ik) = -xx(1,ik) xx(2,ip+ik) = xx(2,ik) xx(3,ip+ik) = xx(3,ik) enddo ipo=ipo*2 ip=ipo do ik=1,ipo xx(1,ip+ik) = xx(1,ik) xx(2,ip+ik) = -xx(2,ik) xx(3,ip+ik) = xx(3,ik) enddo ip = ipo*2 C les poids do ik=1,ip poids(ik)=ppp enddo goto 100 32 continue C ----------- les element polygonaux alpa2= xpi/nno alpa = alpa2*O2 px=O1/ndec py =px * sin(alpa) ppp = sin(alpa)/(ndec*ndec)/o2 rx= O2/3.D0*cos(alpa2)*cos(alpa2)/ndec ry =O2/3.D0*cos(alpa2)*sin(alpa2)/ndec C x0 =-px+rx y0 =-py+ry x1=x0+rx y1=y0+ry C premier triangle do k=1,ndec do i=1,ndec-k+1 ip=ip+1 xx(1,ip) = x0+px*i xx(2,ip) = y0+py*k poids(ip)= ppp if(i.eq.ndec-k+1) goto 57 ip=ip+1 xx(1,ip) = x1+px*i xx(2,ip) = y1+py*k poids(ip)= ppp 57 continue enddo x0 = x0 + px*cos(alpa) enddo ipo=ip C tes triangles complementaires do ik=1,nno-1 alpha= O2*xpi/nno*ik cc = cos(alpha) ss = sin(alpha) do j=1,ipo ip =ip+1 xx(1,ip)=cc*xx(1,j)-ss*xx(2,j) xx(2,ip)=ss*xx(1,j)+cc*xx(2,j) poids(ip)=ppp enddo enddo C write(6,3333) ((xx(i,j),i=1,2),j=1,ip) 3333 format(2E12.5) C goto 100 99 CONTINUE MOTERR(1:4)=NOMTP(IELE) MOTERR(9:14)='PAVPOI' CALL ERREUR(86) C ELEMENTS NON IMPLEMENTES 100 continue ndd =ip C write(6,*) 'pavpoi ip poids som poids ' ,ip,ppp,(ip*ppp) return END