subroutine rhslhs( i pri,nelgrp,coeg,nocon,itri,itedge,segdat,isegd, i rhsdat,conddat,convdat,absodat,refchr, o u,d2ph,ul2,uene,umax,cstab) implicit double precision (a-h,o-z) * ********************************************************* * solves the equation system ********************************************************* * * output: * * ul2 l_2-norm of u * uene energy norm of u * umax maximum of u * character*3 refchr dimension pri(*),nelgrp(5,*),itri(4,*),coeg(2,*),d2ph(*), & nocon(*),u(*),itedge(3,*),segdat(2,*),isegd(*) dimension rhsdat(*),conddat(*),convdat(*),absodat(*) * dimension iv(3),ele(3,3),fele(3) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact common /poisson/ nrhs,ncond,nconv,nabso,numseg * nele=nrtri * * call dwmatg(rhsdat,1,nrtri,0,'rd','rhsdat') * call dwmatg(segdat,2,numseg,0,'rd','segdat2') * call dwmatg(isegd,1,numseg,0,'i','isegd') nvtot=3*nvact neva=3 neno=3 nnodv=neva/neno nv=nvact * call dimpr(pri,nvtot,1,8,lauold,1) call dimpr(pri,nvtot,1,4,lajb,1) call dimpr(pri,nvtot,1,4,lajd,1) call dimpr(pri,nvtot,1,8,larhs,1) call dimpr(pri,nno,1,4,lajnt,1) do 10 i=1,nvtot pri(lauold+i-1)=u(i) 10 continue * ********************************************************* * Renumber to minimize bandwidth, new numbers in jnt ********************************************************* * call renum( i pri,nno,nrtri,itri, o pri(lajnt)) * ********************************************************* * Setup skyline pointers ********************************************************* * call dimpr(pri,nno,1,8,lajod,1) call jbjd( i pri(lajnt),itri,nocon,nrtri,nno, i neno,nnodv,pri(lajod), o pri(lajb),pri(lajd),npv,neq,isize) call dimpr(pri,nno,1,8,lajod,-1) * ********************************************************* * Boundary conditions ********************************************************* * call dimpr(pri,npv,1,4,lajp,1) call dimpr(pri,npv,1,8,lapv,1) * call jppv(u,itri,itedge,nrtri, & segdat,isegd,pri(lajp),pri(lapv),npv,pri(lajnt)) * call dimpr(pri,isize,1,8,lauma,1) call dimpr(pri,isize,1,8,lalma,1) * do 100 iel=1,nrtri * ********************************************************* * Compute element matrices and assemble * for triangle iel ******************************************************* * * call dwmatg(coeg,2,nno,0,'rd','coeg') call eltot(coeg,itri,iel,nelgrp,nocon,itedge, & segdat,isegd,numseg,rhsdat,nrhs,conddat,ncond, & convdat,nconv,absodat,nabso, & ele,fele,pri(lajnt),iv,u) * ********************************************************* * Temporary nodes used in the assembly process ********************************************************* * call shift(iv,pri(lajnt),3) * call assemb( i fele,iv,ele,pri(lajd),neva, m pri(larhs),pri(lauma),pri(lalma)) * 100 continue * ********************************************************* * Solve ********************************************************* * call usolv( & pri(lauma),pri(lalma),pri(lajb),pri(lajd), & neq,u,pri(larhs),pri(lapv),pri(lajp),npv,1,0,4) * * call dwmatg(pri(lapv),1,npv,0,'rd','pv') call shiftv(pri(lajnt),nno,u,pri(larhs)) call eml(u,conddat,ncond,nno,nrtri,itri,nelgrp, $ coeg,ul2,uene,umax) * call dwmatg(u,1,nno,0,'rd','u') * cstab=1.d0 if(refchr.eq.'ENE') goto 9000 * ********************************************************* *additional computations for error control ********************************************************* * call dimpr(pri,nvtot,1,8,lau,1) call dimpr(pri,nvtot,1,8,laux,1) call dimpr(pri,nvtot,1,8,lauy,1) call dimpr(pri,nvtot,1,8,laama,1) * do 120 i=1,isize pri(lauma+i-1)=0.d0 pri(lalma+i-1)=0.d0 120 continue * ********************************************************* * compute element matrices for the dual problem ********************************************************* * do 200 iel=1,nrtri * call elstab(coeg,itri,iel,nocon,itedge, & segdat,isegd,numseg,pri(laama),conddat,ncond, & convdat,nconv,absodat,nabso,nelgrp, & ele,fele,pri(lajnt),iv,u) * * call shift(iv,pri(lajnt),3) * call assemb( i fele,iv,ele,pri(lajd),neva, m pri(larhs),pri(lauma),pri(lalma)) 200 continue * * triangulate the dual system matrix * call usolv( & pri(lauma),pri(lalma),pri(lajb),pri(lajd), & neq,pri(lau),pri(larhs),pri(lapv),pri(lajp),npv,1,0,1) * ********************************************************* * for the maximum norm we must solve several times ********************************************************* * if(refchr.eq.'MAX') then do 210 i=1,nno pri(larhs+i-1)=u(i)-pri(lauold+i-1) 210 continue call stabilm( & pri(lauma),pri(lalma),pri(lajb),pri(lajd), & neq,pri(lau),pri(larhs),pri(lapv),pri(lajp), & npv,coeg,itri,nocon,nno,nrtri,pri(lajnt), & pri(laux),pri(lauy),pri(laama),cstab) endif * if(refchr.eq.'RMS') then * ********************************************************* * data for the dual problem is the error ********************************************************* * do 220 i=1,nno pri(larhs+i-1)=u(i)-pri(lauold+i-1) pri(laama+i-1)=0.d0 220 continue * call dwmatg(pri(larhs),1,nno,0,'rd','rhs') * call stabilr( & pri(lauma),pri(lalma),pri(lajb),pri(lajd), & neq,pri(lau),pri(larhs),pri(lapv),pri(lajp), & npv,coeg,itri,nocon,nno,nrtri,pri(lajnt), & pri(laux),pri(lauy),pri(laama),d2ph) endif * call dimpr(pri,nvtot,1,8,laama,-1) call dimpr(pri,nvtot,1,8,lauy,-1) call dimpr(pri,nvtot,1,8,laux,-1) call dimpr(pri,nvtot,1,8,lau,-1) * 9000 continue call dimpr(pri,isize,1,8,lalma,-1) call dimpr(pri,isize,1,8,lauma,-1) call dimpr(pri,npv,1,8,lapv,-1) call dimpr(pri,npv,1,4,lajp,-1) call dimpr(pri,nno,1,4,lajnt,-1) call dimpr(pri,nvtot,1,8,larhs,-1) call dimpr(pri,nvtot,1,4,lajd,-1) call dimpr(pri,nvtot,1,4,lajb,-1) call dimpr(pri,nvtot,1,8,lauold,-1) * return end * * * subroutine eml( i u,conddat,ncond,nno,nrtri,itri,nelgrp,coeg, o ul2,uene,umax) * implicit double precision (a-h,o-z) dimension u(*),itri(4,*),coeg(2,*),nelgrp(5,*) dimension conddat(ncond,*),condu(2,2) dimension iv(3),fi(3),fix(3),fiy(3),xc(3),yc(3) character*80 rhscha,cndcha(3),cnvcha(2),abscha,ex * common /chardata/ idrhs,idcnd,idcnv,idabs, $ rhscha,cndcha,cnvcha,abscha * ********************************************************* * compute the l2-norm of the solution ********************************************************* * ul2=0.d0 uene=0.d0 umax=0.d0 * do 500 iel=1,nrtri do 300 i=1,3 iv(i)=itri(i,iel) umax=dmax1(umax,dabs(u(iv(i)))) xc(i)=coeg(1,iv(i)) yc(i)=coeg(2,iv(i)) 300 continue * ngau=3 do 400 igm=1,ngau call gaussx(ngau,igm,xc,yc,xm,ym) call base(xm,ym,xc,yc,fi,fix,fiy,area) * ********************************************************* * conduction ********************************************************* * do i=1,2 do j=1,2 condu(i,j)=0.d0 enddo enddo * if(idcnd.eq.0) then ex=cndcha(1) call evalexpr(ex, xm, ym, c1, ierr) condu(1,1)=c1 ex=cndcha(2) call evalexpr(ex, xm, ym, c1, ierr) condu(1,2)=c1 ex=cndcha(3) call evalexpr(ex, xm, ym, c1, ierr) condu(2,2)=c1 condu(2,1)=condu(1,2) endif if(idcnd.eq.1) then call shepardv( i xm,ym,ncond,conddat(1,1),conddat(1,2),conddat(1,3),3, o condu(1,1)) condu(2,2)=condu(1,2) condu(1,2)=condu(2,1) endif if(idcnd.eq.2) then ielfa=macrotr(iel,itri,nelgrp) condu(1,1)=conddat(1,ielfa) condu(2,1)=conddat(2,ielfa) condu(1,2)=condu(2,1) condu(2,2)=conddat(3,ielfa) endif * uac=fi(1)*u(iv(1))+fi(2)*u(iv(2))+fi(3)*u(iv(3)) ux=fix(1)*u(iv(1))+fix(2)*u(iv(2))+fix(3)*u(iv(3)) uy=fiy(1)*u(iv(1))+fiy(2)*u(iv(2))+fiy(3)*u(iv(3)) * uxn=condu(1,1)*ux+condu(1,2)*uy uyn=condu(2,1)*ux+condu(2,2)*uy * ul2=ul2+uac*uac*area/ngau uene=uene+(ux*uxn+uy*uyn)*area/ngau 400 continue 500 continue ul2=dsqrt(ul2) uene=dsqrt(uene) * return end * * * subroutine shiftv( i jnt,nno, m u, d usta) * implicit double precision (a-h,o-z) dimension jnt(*),u(*),usta(*) * ********************************************************* * give the solution the right numbering ********************************************************* * do 100 i=1,nno usta(i)=0.e0 itemp=jnt(i) if(itemp.gt.0) usta(i)=u(itemp) 100 continue do 200 i=1,nno u(i)=usta(i) 200 continue * return end * * * subroutine jbjd( i jnod,itri,nocon,nrtri,nno, i neno,nnodv,jod, o jb,jd,npv,neq,isize) implicit double precision (a-h,o-z) * ********************************************************* * compute pointer arrays for a skyline * matrix stored as a one-dimensional array * * jb(i) points to the first nonzeroelement in * row/col i starting from top of row/col * jd(i) points to the diagonal element in row/col i * starting from top of matrix ********************************************************* * dimension jnod(*),itri(4,*) dimension jb(*),jd(*),nocon(*),jod(*) * ********************************************************* * compute neq, npv * compute jod, jod(i) is the minimum node * number coupled to node i * (only used internally) ********************************************************* * neq=0 npv=0 do 10 i=1,nno if(nocon(i).eq.-1) then npv=npv+1 endif neq=neq+nnodv jod(i)=i 10 continue do 50 iel=1,nrtri kb=jnod(itri(1,iel)) * do 30 i=2,neno kb=min0(kb,jnod(itri(i,iel))) 30 continue * do 40 i=1,neno ki=jnod(itri(i,iel)) jod(ki)=min0(jod(ki),kb) 40 continue 50 continue * ********************************************************* * calculate jb jd, isize ********************************************************* * je=jod(1) do 70 i=2,nno ja=jod(i) jg=(i-1)*nnodv jc=(i-2)*nnodv+1 * do 60 j=jc,jg jb(j)=je 60 continue je=(ja-1)*nnodv+1 70 continue jc=jg+1 jg=jg+nnodv * do 130 j=jc,jg 130 jb(j)=je * jd(1)=1 do 140 i=2,neq 140 jd(i)=jd(i-1)+i-jb(i)+1 * isize=jd(neq) * return end * function fslask(r) double precision fslask,r fslask=r return end * subroutine assemb( i fele,iv,ele,jd,neva, m fu,a,b) implicit double precision (a-h,o-z) * ********************************************************* * assembles element data to global data * * in: fele element load vector * ele element stiffness matrix * iv node numbers * jd pointer to the diagonal of the * global matrices * neva number of element variables * * modified: fu global load vector * a the upper right half of * the system matrix * b the lower left half of * the system matrix ********************************************************* * dimension i iv(*),ele(neva,*),jd(*),fele(*), m fu(*),a(*),b(*) * do 1000 i=1,neva li=iv(i) fu(li)=fu(li)+fele(i) * do 2000 j=1,neva lj=iv(j) if(lj-li) 100,200,300 * 100 continue ij=jd(li)-li+lj b(ij)=b(ij)+ele(i,j) goto 2000 * 200 continue ij=jd(lj) a(ij)=a(ij)+ele(i,j) b(ij)=b(ij)+ele(i,j) goto 2000 * 300 continue ij=jd(lj)-lj+li a(ij)=a(ij)+ele(i,j) 2000 continue 1000 continue return end * * * subroutine usolv(a,b,jb,jd,neq,u,f,pv,jp,npv,npvz,ipe,iwsolv) implicit double precision (a-h,o-z) * ************************************************************************* * purpose: solves an asymmetric system of equations * * makes a partial decomposition (static condensation) of the * non-symmetric system equation su=f to equation ipe-1. * makes a back substitution of the partly decomposed system * equation, when the last (neq-ipe+1) values in u are given * known values. * * a - the upper right half of the system matrix * stored with consequtive columns in a one dimen- * sional array in double precision. size jd(neq) * b - the lower left half of the system matrix * stored with consequtive rows in a one dimen- * sional array in double precision. size jd(neq) * jb - contains for each column in the system matrix the * row/col number of the first non-zero element. size neq * jd - contains the position of the diagonal elements of * the system matrix when it is stored in a one * dimensional array. size neq * neq - number of equations * u - contains solution (including prescribed) as output. * f - contains load as input. * contains load and reactions (if npvz.ne.0) as output. * pv - contains the magnitude of the prescribed variables. * jp - contains the numbers of the prescribed variables. * the numbers have to be stored in an increasing * order. size npv * npv - number of prescribed variables (may be zero). * npvz - npvz=0 all prescribed variables are zero, or the * load vector is modified elsewhere. * ipe - ipe=0 gives a normal solution * ipe>0 a partial decomposition is made. ipe = the * first row number in the part of the matrix which * is not triangulated * iwsolv - iwsolv=1 triangulate a with respect to ipe * iwsolv=2 triangulate u with respect to ipe * (the triangulated matrix must be given as input)* * iwsolv=3 backsubstitution with respect to ipe * (a triang. and u triang. or with a part of * the solution as input) * iwsolv=4 triangulate and backsubstitute (1+2+3) * with respect to ipe ************************************************************************* * dimension a(*),b(*),u(*),f(*),jd(*),jb(*),jp(*),pv(*) * do 10 i=1,neq 10 u(i)=f(i) * ip1=ipe-1 if(ip1.lt.0) ip1=neq * ********************************************************* * change jb for variables which have * prescribed values ********************************************************* * if (npv.eq.0) go to 30 do 20 j=1,npv k=jp(j) 20 jb(k)=-iabs(jb(k)) * 30 goto (1000,2000,3000,1000),iwsolv 1000 continue * ********************************************************* * triangulate the matrix ********************************************************* * do 90 i=2,neq jbi=jb(i) je=i-1 if(jbi.lt.0.or.jbi.gt.je) goto 90 jdi=jd(i)-i * do 60 j=jbi,je ij=jdi+j rb=b(ij) jdj=jd(j) jbj=jb(j) if(jbj.lt.0) goto 60 if(j.eq.1) goto 50 jdk=jdj-j kb=jbi if(jbi.lt.jbj) kb=jbj ke=j-1 if(ke.gt.ip1) ke=ip1 if(kb.gt.ke) goto 50 * do 40 k=kb,ke if(jb(k).lt.0) goto 40 ik=jdi+k kj=jdk+k rb=rb-b(ik)*a(kj) 40 continue 50 continue b(ij)=rb/a(jdj) 60 continue * * lb=2 if(jbi.gt.2) lb=jbi * do 80 l=lb,i li=jdi+l jbl=jb(l) if(jbl.lt.0) goto 80 kb=jbi if(kb.lt.jbl) kb=jbl ke=l-1 if(ke.gt.ip1) ke=ip1 if(kb.gt.ke) goto 80 jdl=jd(l)-l ra=a(li) * do 70 k=kb,ke if(jb(k).lt.0) goto 70 lk=jdl+k ki=jdi+k ra=ra-b(lk)*a(ki) 70 continue a(li)=ra 80 continue 90 continue if(iwsolv.eq.4) goto 2000 goto 9000 2000 continue * ********************************************************* * triangulate the load array ********************************************************* * if (npv.eq.0.or.npvz.eq.0) go to 160 * ********************************************************* * change the load array for variables which have * prescribed values ********************************************************* * jbe=1 * do 150 ip=1,npv ieq=jp(ip) if(ieq.gt.1) jbe=jd(ieq-1)+1 xp=pv(ip) if (xp) 100,150,100 100 jend=jd(ieq)-1 if(jend.lt.jbe) goto 130 i1=0 jbi=-jb(ieq) * do 120 i=jbe,jend 110 i1=i1+1 if(jbi.gt.i1) goto 110 u(i1)=u(i1)-a(i)*xp 120 continue 130 continue jbk=1 i2=0 * do 140 j=ieq,neq if(j.gt.1) jbk=jd(j-1)+1 jp1=jd(j)-i2 i2=i2+1 if(jbk.gt.jp1) goto 140 u(j)=u(j)-b(jp1)*xp 140 continue 150 continue * ********************************************************* * triangulate ********************************************************* * 160 do 180 l=2,neq jbl=jb(l) if(jbl.lt.0) goto 180 kb=jbl ke=l-1 if(ke.gt.ip1) ke=ip1 if(kb.gt.ke) goto 180 jdl=jd(l)-l rf=u(l) * do 170 k=kb,ke if(jb(k).lt.0) goto 170 lk=jdl+k rf=rf-b(lk)*u(k) 170 continue u(l)=rf 180 continue * i2=neq if(ipe.le.neq) i2=ip1 * do 190 i=1,i2 if(jb(i).lt.0) goto 190 aii=a(jd(i)) u(i)=u(i)/aii 190 continue if(iwsolv.eq.2.or.iwsolv.eq.4) goto 3000 goto 9000 * ********************************************************* * backsubstitute ********************************************************* * 3000 continue * do 240 j=1,neq i=neq-j+1 i1=jb(i) if (i1) 240,240,200 200 if (i1-i) 210,240,210 210 i2=i-1 if ((ip1+1).lt.i) i2=ip1 if (i1-i2) 220,220,240 220 jr3=jd(i)-i * do 230 k=i1,i2 ai3=a(jr3+k) ki=-k+i jdk=jd(k) 230 u(k)=u(k)-ai3*u(k+ki)/a(jdk) 240 continue * if (npv.eq.0) go to 9000 if (npvz.eq.0) go to 330 * ********************************************************* * calculate reactions ********************************************************* * * do 320 j=1,npv ieq=jp(j) u(ieq)=0.d0 jbieq=-jb(ieq) idi=ieq-1 if(jbieq.gt.idi) goto 280 i1=-1 * do 270 i=jbieq,idi i1=i1+1 xp=u(i) if(jb(i).gt.0) goto 260 * do 250 k=1,npv if(jp(k).eq.i) xp=pv(k) if (jp(k)-i) 250,260,250 250 continue 260 continue ja=jd(ieq)-ieq+iabs(jb(ieq))+i1 u(ieq)=u(ieq)+xp*b(ja) 270 continue 280 continue i2=0 * do 310 i=ieq,neq jbi=iabs(jb(i)) if(jbi.gt.ieq) goto 310 ia=jd(i)-i2 i2=i2+1 xp=u(i) if(jb(i).gt.0) goto 300 * do 290 k=1,npv if(jp(k).eq.i) xp=pv(k) if (jp(k)-i) 290,300,290 290 continue 300 continue u(ieq)=u(ieq)+xp*a(ia) 310 continue 320 continue 330 continue * ********************************************************* * change jb,f,u ********************************************************* * do 350 j=1,npv k=jp(j) if(npvz.eq.0) goto 340 f(k)=-f(k)+u(k) 340 u(k)=pv(j) jb(k)=-jb(k) 350 continue 9000 continue return end * * * subroutine jtcal( i itri,nrtri, o jt) * implicit double precision (a-h,o-z) dimension itri(4,*),jt(*) * do 10 i=1,nrtri do 30 j=1,3 jt(nrtri*(j-1)+i)=itri(j,i) 30 continue 10 continue return end * * * subroutine bands(nno,nrtri,jt,memjt,jmem,idiff,maxin) implicit double precision (a-h,o-z) * dimension jt(*),memjt(*),jmem(*) * ********************************************************* * size of jt: (number of nodes)*(nodes on element) * size of memjt: number of nodes times number of * connecting nodes (maybe 10 maximum) * size of jmem: number of nodes ********************************************************* * idiff=nno call zerop(jmem,nno*4) do 60 j=1,nrtri do 50 i=1,3 * jnti=jt(nrtri*(i-1)+j) * if(jnti.eq.0) goto 60 jsub=(jnti-1)*maxin * do 40 ii=1,3 if(ii.eq.i) goto 40 jjt=jt(nrtri*(ii-1)+j) if(jjt.eq.0) goto 50 mem1=jmem(jnti) if(mem1.eq.0) goto 30 do 20 iii=1,mem1 if(memjt(jsub+iii).eq.jjt) goto 40 20 continue 30 continue jmem(jnti)=jmem(jnti)+1 memjt(jsub+jmem(jnti))=jjt if(iabs(jnti-jjt).gt.idiff) idiff=iabs(jnti-jjt) 40 continue 50 continue 60 continue return end * * * subroutine optnum(nno,memjt,jmem,jnt,idiff,newjt, & joint,maxin) implicit double precision (a-h,o-z) * ********************************************************* * jnt contains the new node numbers (temporary) ********************************************************* * dimension memjt(*),jmem(*),jnt(*),newjt(*),joint(*) * minmax=idiff do 60 ik=1,nno if(jmem(ik).eq.0) goto 60 call zerop(joint,nno*4) call zerop(newjt,nno*4) max=0 i=1 newjt(1)=ik joint(ik)=1 k=1 30 continue k4=jmem(newjt(i)) if(k4.eq.0) goto 45 jsub=(newjt(i)-1)*maxin do 40 jj=1,k4 k5=memjt(jsub+jj) if(joint(k5).gt.0) goto 40 k=k+1 newjt(k)=k5 joint(k5)=k ndiff=iabs(i-k) if(ndiff.ge.minmax) goto 60 if(ndiff.gt.max) max=ndiff 40 continue if(k.eq.nno) goto 50 45 continue i=i+1 if(i.gt.nno) goto 50 goto 30 50 continue minmax=max do 55 j=1,nno jnt(j)=joint(j) 55 continue 60 continue return end * * * subroutine shift(iv,jnt,neno) * implicit double precision (a-h,o-z) integer iv(*),jnt(*) * do 10 i=1,neno iv(i)=jnt(iv(i)) 10 continue * return end * * * subroutine gaussx(in,i,xc,yc,xm,ym) implicit double precision (a-h,o-z) * ********************************************************* * compute coordinate corresponding to gauss point * * in number of gauss points in the triangle * i number of the gauss point for which * coordinate is sought ********************************************************* * dimension xc(*),yc(*) go to (10,200,20),in 10 continue xm=(xc(1)+xc(2)+xc(3))/3.d0 ym=(yc(1)+yc(2)+yc(3))/3.d0 go to 200 20 continue go to (21,22,23),i 21 xm=(xc(1)+xc(2))/2.d0 ym=(yc(1)+yc(2))/2.d0 go to 200 22 xm=(xc(2)+xc(3))/2.d0 ym=(yc(2)+yc(3))/2.d0 go to 200 23 xm=(xc(1)+xc(3))/2.d0 ym=(yc(1)+yc(3))/2.d0 200 continue return end * * * subroutine base( i x,y,xc,yc, o fi,fix,fiy,area) implicit double precision (a-h,o-z) * ********************************************************* * base functions and area measure * for constant strain elements * * x,y coordinates * xc(3),yc(3) node coordinates * fi(3) base functions * fix(3) ... derivatives * area area of element ********************************************************* * dimension xc(*),yc(*),fi(*),fiy(*),fix(*) dimension dy(3),dx(3) * det=0.d0 do 10 j=1,3 j1 = (3*j*j-13*j+16)/2 j2 = (-3*j*j+11*j-4)/2 dx(j) = xc(j1)-xc(j2) dy(j) = yc(j1)-yc(j2) det = det-xc(j)*dy(j) 10 continue * do 20 j = 1,3 fi(j) = 1.d0-((x-xc(j))*dy(j)-(y-yc(j))*dx(j))/det fix(j) = -dy(j)/det fiy(j) = dx(j)/det 20 continue * area=det/2.d0 return end * * * subroutine meshpl(pri,coeg,u,itri,cond,itedge,nelgrp, & arrow,imesh) * ********************************************************* * write for plot ********************************************************* * implicit double precision (a-h,o-z) common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact common /poisson/ nrhs,ncond,nconv,nabso,numseg * dimension pri(*),cond(*),arrow(2,*),nelgrp(5,*) dimension u(*),coeg(2,*),itri(4,*),itedge(3,*) * call dimpr(pri,nno,1,8,laama,1) call arrows(u,cond,ncond,coeg,itri,nno,nrtri, $ nelgrp,arrow,pri(laama)) call dimpr(pri,nno,1,8,laama,-1) * write(77) nrtri,nno write(77) ((itri(i,j),i=1,3),j=1,nrtri) write(77) ((itedge(i,j),i=1,3),j=1,nrtri) write(77) ((coeg(i,j),i=1,2),j=1,nno) write(77) (u(j),j=1,nno) write(77) ((arrow(i,j),i=1,2),j=1,nno) * return end subroutine wrsol(iunf,iout,coeg,u,itri,itedge,imesh) * ********************************************************* * write for plot ********************************************************* * implicit double precision (a-h,o-z) * dimension u(*),coeg(2,*),itri(3,*),itedge(3,*) * write(iout,*) 'STATE' write(iout,*) 24 write(iout,*) 'SOLUTION' write(iout,*) imesh rewind(iunf) * do 10 im=1,imesh read(iunf) nrtri,nno read(iunf) ((itri(i,j),i=1,3),j=1,nrtri) read(iunf) ((itedge(i,j),i=1,3),j=1,nrtri) read(iunf) ((coeg(i,j),i=1,2),j=1,nno) read(iunf) (u(j),j=1,nno) * write(iout,*) nrtri,nno,2,3 write(iout,*) ((itri(i,j),i=1,3),j=1,nrtri) write(iout,*) ((itedge(i,j),i=1,3),j=1,nrtri) write(iout,100) ((coeg(i,j),i=1,2),j=1,nno) write(iout,*) 'Concentration;',1 write(iout,100) (u(j),j=1,nno) * read(iunf) ((coeg(i,j),i=1,2),j=1,nno) write(iout,*) 'Diffusive flux;',2 write(iout,100) ((coeg(i,j),i=1,2),j=1,nno) 10 continue write(iout,*) 'END' close(unit=iunf,status='delete') endfile(unit=iout) close(unit=iout,status='keep') * return 100 format(1x,e12.5) end * * * subroutine gener(nelgrp,nocon,level) * implicit double precision (a-h,o-z) dimension nelgrp(5,*),nocon(*),level(*) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * nelr=mact ielfr=mact+1 inofr=nno+1 lvl=1 nelgrp(1,ielfr)=-(ielfr+1) nocon(inofr)=inofr+1 do 10 i=1,nelr level(i)=1 10 continue * return end * * * subroutine jppv(u,itri,itedge,nrtri, & segda,isegd,jp,pv,npv,jnt) implicit double precision (a-h,o-z) * dimension itri(4,*),itedge(3,*),segda(2,*),isegd(*) dimension u(*),jp(*),pv(*),jnt(*) * npv=0 do 400 iel=1,nrtri do 500 i=1,3 if(itedge(i,iel).gt.0) goto 500 isname= -itedge(i,iel) if(isegd(isname).ge.0) goto 500 * i1=(8-5*i+i*i)/2 i2=(4+3*i-i*i)/2 iv1=itri(i1,iel) iv2=itri(i2,iel) * if(npv.gt.0) then do 10 j=1,npv if(jp(j).eq.jnt(iv1)) goto 20 10 continue npv=npv+1 jp(npv)=jnt(iv1) pv(npv)=segda(1,isname) 20 continue do 30 j=1,npv if(jp(j).eq.jnt(iv2)) goto 40 30 continue npv=npv+1 jp(npv)=jnt(iv2) pv(npv)=segda(1,isname) 40 continue else npv=npv+1 jp(npv)=jnt(iv1) pv(npv)=segda(1,isname) npv=npv+1 jp(npv)=jnt(iv2) pv(npv)=segda(1,isname) endif * 500 continue 400 continue return end * * * subroutine renum(pri,nno,nrtri,itri,jnt) implicit double precision (a-h,o-z) * dimension pri(*),jnt(*),itri(4,*) * maxin=12 call dimpr(pri,nrtri,3,4,lajt,1) call dimpr(pri,nno,maxin,4,lamemj,1) call dimpr(pri,nno,1,4,lajmem,1) call dimpr(pri,nno,1,4,lajoin,1) call dimpr(pri,nno,1,4,lanewj,1) * call jtcal(itri,nrtri,pri(lajt)) * call bands(nno,nrtri,pri(lajt),pri(lamemj), & pri(lajmem),idiff,maxin) * call optnum(nno,pri(lamemj), & pri(lajmem),jnt,idiff,pri(lanewj),pri(lajoin),maxin) * call dimpr(pri,nno,1,4,lanewj,-1) call dimpr(pri,nno,1,4,lajoin,-1) call dimpr(pri,nno,1,4,lajmem,-1) call dimpr(pri,nno,maxin,4,lamemj,-1) call dimpr(pri,nrtri,3,4,lajt,-1) * return end subroutine eltot(coeg,itri,iel,nelgrp,nocon,itedge, & segdat,isegd,numseg,rhsdat,nrhs,conddat,ncond, & convdat,nconv,absodat,nabso, & ele,fele,jnt,iv,u) * implicit double precision (a-h,o-z) dimension coeg(2,*),itri(4,*),nocon(*),itedge(3,*), & ele(3,*),fele(*),jnt(*),iv(*),u(*),nelgrp(5,*) dimension rhsdat(nrhs,*),conddat(ncond,*),segdat(2,*) dimension convdat(nconv,*),absodat(nabso,*),isegd(*) dimension fi(3),fix(3),fiy(3),xc(3),yc(3),fiv(3),condu(2,2) dimension fiex(3),fiey(3),veloc(2) character*80 rhscha,cndcha(3),cnvcha(2),abscha,ex * common /chardata/ idrhs,idcnd,idcnv,idabs, $ rhscha,cndcha,cnvcha,abscha * ********************************************************* * element matrix and load vector ********************************************************* * eps=1.d-4 * if(idcnd.eq.2.or.idabs.eq.2.or.idcnv.eq.2.or. $ idrhs.eq.2) ielfa=macrotr(iel,itri,nelgrp) * xm=0.d0 ym=0.d0 do 20 i=1,3 fele(i)=0.d0 iv(i)=itri(i,iel) xc(i)=coeg(1,iv(i)) yc(i)=coeg(2,iv(i)) xm=xm+xc(i)/3. ym=ym+yc(i)/3. do 20 j=1,3 ele(i,j)=0.d0 20 continue * ngau=1 do 400 igm=1,ngau call gaussx(ngau,igm,xc,yc,xm,ym) call base(xm,ym,xc,yc,fi,fix,fiy,area) wei=area/ngau h1=dsqrt(2.d0*area) uac=fi(1)*u(iv(1))+fi(2)*u(iv(2))+fi(3)*u(iv(3)) do 30 i=1,3 fiex(i)=0.d0 fiey(i)=0.d0 30 continue * ********************************************************* * conduction ********************************************************* * do i=1,2 do j=1,2 condu(i,j)=0.d0 enddo enddo * if(idcnd.eq.0) then ex=cndcha(1) call evalexpr(ex, xm, ym, c1, ierr) condu(1,1)=c1 ex=cndcha(2) call evalexpr(ex, xm, ym, c1, ierr) condu(1,2)=c1 ex=cndcha(3) call evalexpr(ex, xm, ym, c1, ierr) condu(2,2)=c1 condu(2,1)=condu(1,2) endif if(idcnd.eq.1) then call shepardv( i xm,ym,ncond,conddat(1,1),conddat(1,2),conddat(1,3),3, o condu(1,1)) condu(2,2)=condu(1,2) condu(1,2)=condu(2,1) endif if(idcnd.eq.2) then condu(1,1)=conddat(1,ielfa) condu(2,1)=conddat(2,ielfa) condu(1,2)=condu(2,1) condu(2,2)=conddat(3,ielfa) endif * do i=1,3 fiex(i)=fiex(i)+condu(1,1)*fix(i)+condu(1,2)*fiy(i) fiey(i)=fiey(i)+condu(2,1)*fix(i)+condu(2,2)*fiy(i) enddo * ********************************************************* * absorption ********************************************************* * if(idabs.eq.0) call evalexpr(abscha, xm, ym, prod, ierr) if(idabs.eq.1) call shepardv( i xm,ym,nabso,absodat(1,1),absodat(1,2),absodat(1,3),1, o prod) if(idabs.eq.2) prod=absodat(1,ielfa) * ********************************************************* * convection ********************************************************* * if(idcnv.eq.0) then ex=cnvcha(1) call evalexpr(ex, xm, ym, veloc(1), ierr) ex=cnvcha(2) call evalexpr(ex, xm, ym, veloc(2), ierr) endif if(idcnv.eq.1) call shepardv( i xm,ym,nconv,convdat(1,1),convdat(1,2),convdat(1,3),2, o veloc) if(idcnv.eq.2) then veloc(1)=convdat(1,ielfa) veloc(2)=convdat(2,ielfa) endif * ********************************************************* * load ********************************************************* * if(idrhs.eq.0) call evalexpr(rhscha, xm, ym, fxy, ierr) if(idrhs.eq.1) call shepardv( i xm,ym,nrhs,rhsdat(1,1),rhsdat(1,2),rhsdat(1,3),1, o fxy) if(idrhs.eq.2) fxy=rhsdat(1,ielfa) * ********************************************************* * streamline diffusion ********************************************************* * if(ierr.eq.1) then write(6,*) 'IERR=',ierr stop endif h2=0.d0 a11=condu(1,1) a22=condu(2,2) a12=condu(2,1) f1=(a11+a22)/2.d0 gl1=f1+dsqrt(f1*f1-(a11*a22-a12*a12)) if(gl1-h1.lt.-eps*h1) h2=h1*dexp(gl1/(gl1-h1)) * velx=veloc(1) vely=veloc(2) vabs=h2/dmax1(sqrt(velx*velx+vely*vely),eps) * do 350 i=1,3 fiv(i)=fi(i)+(velx*fix(i)+vely*fiy(i))*vabs fele(i)=fele(i)+fiv(i)*fxy*wei do 300 j=1,3 * ele(i,j)=ele(i,j)+((velx*fix(j)+vely*fiy(j))*fiv(i) & +fiy(i)*fiey(j)+fix(i)*fiex(j))*wei * ele(i,j)=ele(i,j)+prod*fi(j)*fiv(i)*wei 300 continue 350 continue 400 continue * ********************************************************* * boundary production term, boundary data ********************************************************* * do 500 i=1,3 if(itedge(i,iel).gt.0) goto 500 isname= -itedge(i,iel) if(isegd(isname).le.0) goto 500 iv1=(8-5*i+i*i)/2 iv2=(4+3*i-i*i)/2 dx=xc(iv1)-xc(iv2) dy=yc(iv1)-yc(iv2) xmid=(xc(iv1)+xc(iv2))/2.d0 ymid=(yc(iv1)+yc(iv2))/2.d0 call base(xmid,ymid,xc,yc,fi,fix,fiy,area) dl=dsqrt(dx*dx+dy*dy) fneu=segdat(1,isname) gamma=segdat(2,isname) * write(6,*) 'gamma=',gamma,' fneu=',fneu if(gamma.lt.0) then gamma=-gamma fneu=-fneu endif ele(iv1,iv1)=ele(iv1,iv1)+gamma*dl/3.d0 ele(iv1,iv2)=ele(iv1,iv2)+gamma*dl/6.d0 ele(iv2,iv1)=ele(iv2,iv1)+gamma*dl/6.d0 ele(iv2,iv2)=ele(iv2,iv2)+gamma*dl/3.d0 fele(iv1)=fele(iv1)+fneu*dl/2.d0 fele(iv2)=fele(iv2)+fneu*dl/2.d0 500 continue 600 continue * return end subroutine shepardv( i x,y,npoin,xpoin,ypoin,fpoin,nfxy, o fxy) implicit double precision (a-h,o-z) dimension xpoin(*),ypoin(*),fpoin(npoin,*),fxy(*) * * use shepard's distance-weighted metod * for interpolating scattered data. * the subroutine returns the interpolated * function fxy evaluated at (x,y). * * xpoin - x-coordiantes for all points * ypoin - y-coordiantes for all points * fpoin - function values for all points * npoin - number of points * do 10 i=1,nfxy fxy(i)=0.d0 10 continue if(npoin.eq.0) return denom=0.d0 do 30 i=1,npoin xi=xpoin(i) yi=ypoin(i) dist=dsqrt((x-xi)*(x-xi)+(y-yi)*(y-yi)) dist=dmax1(dist,1.d-10) do 20 j=1,nfxy fxy(j)=fxy(j)+fpoin(i,j)/dist 20 continue denom=denom+1.d0/dist 30 continue do 40 i=1,nfxy fxy(i)=fxy(i)/denom 40 continue return end subroutine stabilm( & ama,bma,jb,jd,neq,u,rhs,pv,jp,npv, & coeg,itri,nocon,nno,nrtri,jnt,uxnod, & uynod,amass,cstab) implicit double precision (a-h,o-z) dimension ama(*),bma(*),jb(*),jd(*),u(*),coeg(2,*) dimension rhs(*),pv(*),jp(*),nocon(*),itri(4,*) dimension uxnod(*),uynod(*),amass(*),jnt(*) dimension xc(3),yc(3),jn(2) * * find a location where the solution has changed a lot * rmax=0.d0 jnew=0 do 5 i=1,nno if(rmax.lt.dabs(rhs(i)).and.nocon(i).eq.0) then jold=jnew rmax=dabs(rhs(i)) jnew=i endif 5 continue jn(1)=jold jn(2)=jnew * * compute stability constants * cstab=0.d0 call zerop(pv,npv*8) call zerop(rhs,nno*8) iold=1 * do 10 ij=1,2 i=jn(ij) * rhs(iold)=0.d0 iold=jnt(i) rhs(iold)=1.d0 * call usolv( & ama,bma,jb,jd, & neq,u,rhs,pv,jp,npv,0,0,2) * call shiftv(jnt,nno,u,amass) * * compute the integral of the second derivatives * do 11 j=1,nno uxnod(j)=0.d0 uynod(j)=0.d0 amass(j)=0.d0 11 continue * gradu=0.d0 do 20 iel=1,nrtri do 15 j=1,3 ivj=itri(j,iel) xc(j)=coeg(1,ivj) yc(j)=coeg(2,ivj) 15 continue * call grad(ux,uy,coeg,u,itri(1,iel)) da = tarea(xc,yc) gradu=gradu+(dabs(ux)+dabs(uy))*da do 16 j=1,3 ivj=itri(j,iel) uxnod(ivj)=uxnod(ivj)+da*ux/3.d0 uynod(ivj)=uynod(ivj)+da*uy/3.d0 amass(ivj)=amass(ivj)+da/3.d0 16 continue 20 continue do 21 j=1,nno if(amass(j).eq.0.d0) goto 21 uxnod(j)=uxnod(j)/amass(j) uynod(j)=uynod(j)/amass(j) 21 continue * * second derivatives * d2phi=0.d0 * do 30 iel=1,nrtri do 25 j=1,3 ivj=itri(j,iel) xc(j)=coeg(1,ivj) yc(j)=coeg(2,ivj) 25 continue call grad(uxx,uxy,coeg,uxnod,itri(1,iel)) call grad(uyx,uyy,coeg,uynod,itri(1,iel)) uxy=(uxy+uyx)/2.d0 d2max=dmax1(dabs(uxx),dabs(uxy),dabs(uyy)) d2phi=d2phi + d2max*tarea(xc,yc) 30 continue cstab=dmax1(d2phi,cstab) 10 continue * return end subroutine stabilr( & ama,bma,jb,jd,neq,u,rhs,pv,jp,npv, & coeg,itri,nocon,nno,nrtri,jnt,uxnod, & uynod,amass,d2ph) implicit double precision (a-h,o-z) dimension ama(*),bma(*),jb(*),jd(*),u(*),coeg(2,*) dimension rhs(*),pv(*),jp(*),nocon(*),itri(4,*) dimension uxnod(*),uynod(*),amass(*),jnt(*),d2ph(*) dimension xc(3),yc(3),iv(3) * * compute second derivatives for l_2-control * call zerop(pv,npv*8) errnorm=0.d0 do 10 iel=1,nrtri rhsel=0.d0 do 5 j=1,3 ivj=itri(j,iel) iv(j)=jnt(ivj) xc(j)=coeg(1,ivj) yc(j)=coeg(2,ivj) rhsel=rhsel+rhs(ivj)/3.d0 5 continue da = tarea(xc,yc) amass(iv(1))=amass(iv(1))+da*rhsel amass(iv(2))=amass(iv(2))+da*rhsel amass(iv(3))=amass(iv(3))+da*rhsel errnorm=errnorm+rhsel*rhsel*da 10 continue * errnorm=dsqrt(errnorm) do 11 j=1,nno rhs(j)=amass(j)/errnorm uxnod(j)=0.d0 uynod(j)=0.d0 amass(j)=0.d0 11 continue * call usolv( & ama,bma,jb,jd, & neq,u,rhs,pv,jp,npv,0,0,2) * call shiftv(jnt,nno,u,rhs) * call dwmatg(u,1,nno,0,'rd','u') * * compute the second derivatives * gradu=0.d0 do 20 iel=1,nrtri do 15 j=1,3 ivj=itri(j,iel) xc(j)=coeg(1,ivj) yc(j)=coeg(2,ivj) 15 continue * call grad(ux,uy,coeg,u,itri(1,iel)) da = tarea(xc,yc) gradu=dmax1(gradu,dabs(ux),dabs(uy)) do 16 j=1,3 ivj=itri(j,iel) uxnod(ivj)=uxnod(ivj)+da*ux/3.d0 uynod(ivj)=uynod(ivj)+da*uy/3.d0 amass(ivj)=amass(ivj)+da/3.d0 16 continue 20 continue do 21 j=1,nno if(amass(j).eq.0.d0) goto 21 uxnod(j)=uxnod(j)/amass(j) uynod(j)=uynod(j)/amass(j) 21 continue * * second derivatives * do 30 iel=1,nrtri do 25 j=1,3 ivj=itri(j,iel) xc(j)=coeg(1,ivj) yc(j)=coeg(2,ivj) 25 continue call grad(uxx,uxy,coeg,uxnod,itri(1,iel)) call grad(uyx,uyy,coeg,uynod,itri(1,iel)) uxy=(uyx+uxy)/2.d0 d2max=dmax1(dabs(uxx),dabs(uxy),dabs(uyy)) d2ph(iel)=d2max 30 continue return end subroutine elstab(coeg,itri,iel,nocon,itedge, & segdat,isegd,numseg,rhs,conddat,ncond, & convdat,nconv,absodat,nabso,nelgrp, & ele,fele,jnt,iv,u) * implicit double precision (a-h,o-z) dimension coeg(2,*),itri(4,*),nocon(*),itedge(3,*), & ele(3,*),fele(*),jnt(*),iv(*),u(*),nelgrp(5,*) dimension rhs(*),conddat(ncond,*),segdat(2,*) dimension convdat(nconv,*),absodat(nabso,*),isegd(*) dimension fi(3),fix(3),fiy(3),xc(3),yc(3),fiv(3),condu(2,2) dimension fiex(3),fiey(3),veloc(2) character*80 rhscha,cndcha(3),cnvcha(2),abscha,ex * common /chardata/ idrhs,idcnd,idcnv,idabs, $ rhscha,cndcha,cnvcha,abscha * * ********************************************************* * element matrix and load vector ********************************************************* * eps=1.d-4 * if(idcnd.eq.2.or.idabs.eq.2.or.idcnv.eq.2) $ ielfa=macrotr(iel,itri,nelgrp) * do 20 i=1,3 fele(i)=0.d0 iv(i)=itri(i,iel) xc(i)=coeg(1,iv(i)) yc(i)=coeg(2,iv(i)) xm=xm+xc(i)/3. ym=ym+yc(i)/3. do 20 j=1,3 ele(i,j)=0.d0 20 continue * ngau=1 do 400 igm=1,ngau call gaussx(ngau,igm,xc,yc,xm,ym) call base(xm,ym,xc,yc,fi,fix,fiy,area) wei=area/ngau h1=dsqrt(2.d0*area) uac=fi(1)*u(iv(1))+fi(2)*u(iv(2))+fi(3)*u(iv(3)) do 30 i=1,3 fiex(i)=0.d0 fiey(i)=0.d0 30 continue * ********************************************************* * conduction ********************************************************* * do i=1,2 do j=1,2 condu(i,j)=0.d0 enddo enddo * if(idcnd.eq.0) then ex=cndcha(1) call evalexpr(ex, xm, ym, c1, ierr) condu(1,1)=c1 ex=cndcha(2) call evalexpr(ex, xm, ym, c1, ierr) condu(1,2)=c1 ex=cndcha(3) call evalexpr(ex, xm, ym, c1, ierr) condu(2,2)=c1 condu(2,1)=condu(1,2) endif if(idcnd.eq.1) then call shepardv( i xm,ym,ncond,conddat(1,1),conddat(1,2),conddat(1,3),3, o condu(1,1)) condu(2,2)=condu(1,2) condu(1,2)=condu(2,1) endif if(idcnd.eq.2) then condu(1,1)=conddat(1,ielfa) condu(2,1)=conddat(2,ielfa) condu(1,2)=condu(2,1) condu(2,2)=conddat(3,ielfa) endif do 40 i=1,3 fiex(i)=fiex(i)+condu(1,1)*fix(i)+condu(1,2)*fiy(i) fiey(i)=fiey(i)+condu(2,1)*fix(i)+condu(2,2)*fiy(i) 40 continue * ********************************************************* * absorption ********************************************************* * if(idabs.eq.0) call evalexpr(abscha, xm, ym, prod, ierr) if(idabs.eq.1) call shepardv( i xm,ym,nabso,absodat(1,1),absodat(1,2),absodat(1,3),1, o prod) if(idabs.eq.2) prod=absodat(1,ielfa) * ********************************************************* * convection ********************************************************* * if(idcnv.eq.0) then ex=cnvcha(1) call evalexpr(ex, xm, ym, veloc(1), ierr) ex=cnvcha(2) call evalexpr(ex, xm, ym, veloc(2), ierr) endif if(idcnv.eq.1) call shepardv( i xm,ym,nconv,convdat(1,1),convdat(1,2),convdat(1,3),2, o veloc) if(idcnv.eq.2) then veloc(1)=convdat(1,ielfa) veloc(2)=convdat(2,ielfa) endif * ********************************************************* * load ********************************************************* * fxy=fi(1)*rhs(iv(1))+fi(2)*rhs(iv(2))+fi(3)*rhs(iv(3)) * ********************************************************* * streamline diffusion ********************************************************* * h2=0.d0 a11=condu(1,1) a22=condu(2,2) a12=condu(2,1) f1=(a11+a22)/2.d0 gl1=f1+dsqrt(f1*f1-(a11*a22-a12*a12)) if(gl1-h1.lt.-eps*h1) h2=h1*dexp(gl1/(gl1-h1)) * velx=-veloc(1) vely=-veloc(2) vabs=h2/dmax1(sqrt(velx*velx+vely*vely),eps) * do 350 i=1,3 fiv(i)=fi(i)+(velx*fix(i)+vely*fiy(i))*vabs fele(i)=fele(i)+fiv(i)*fxy*wei do 300 j=1,3 ele(i,j)=ele(i,j)+((velx*fix(j)+vely*fiy(j))*fiv(i) & +fiy(i)*fiey(j)+fix(i)*fiex(j))*wei ele(i,j)=ele(i,j)+prod*fi(j)*fiv(i)*wei 300 continue 350 continue 400 continue * ********************************************************* * boundary production term, boundary data ********************************************************* * do 500 i=1,3 if(itedge(i,iel).gt.0) goto 500 isname= -itedge(i,iel) if(isegd(isname).le.0) goto 500 iv1=(8-5*i+i*i)/2 iv2=(4+3*i-i*i)/2 dx=xc(iv1)-xc(iv2) dy=yc(iv1)-yc(iv2) dl=dsqrt(dx*dx+dy*dy) gamma=segdat(2,isname) ele(iv1,iv1)=ele(iv1,iv1)+gamma*dl/3.d0 ele(iv1,iv2)=ele(iv1,iv2)+gamma*dl/6.d0 ele(iv2,iv1)=ele(iv2,iv1)+gamma*dl/6.d0 ele(iv2,iv2)=ele(iv2,iv2)+gamma*dl/3.d0 500 continue 600 continue * return end integer function macrotr(iel,itri,nelgrp) dimension nelgrp(5,*),itri(4,*) iel2=itri(4,iel) 10 continue iel3=nelgrp(5,iel2) if(iel3.eq.0) then macrotr=iel2 return else iel2=iel3 goto 10 endif return end