c...Poisson equation "main engine" c...for use with FELIPE version 3 c... - sample datafile is pohex3.dat or podemo.dat c program poiss c parameter(melem=500,mpoin=300,mnode=3, 1 msets=4,mcomp=2,mband=60) c...these parameters fix the dimensions of the arrays: c...melem = max. no. of elements c...mpoin = max. no. of nodes c...mnode = max. no. of nodes in a single element c...msets = max. no. of mat. prop. sets c...mcomp = max. no. of mat. properties (coeffts ax and ay) c...mband = halfbandwidth+1 of global matrix gstif using band storage c c character*1 char,dbdcmp(8) character*80 title,infile,outfil,prtfil c dimension lnods(melem,mnode),lnodn(melem),lsidn(melem), 1 coord(mpoin,2),nkode(mpoin),aslod(mpoin),ndof(mpoin), 2 lpros(melem),props(msets,mcomp),asdis(mpoin), 3 gstif(mpoin,mband),dbdcrd(8),dbdval(8),lmtyp(melem) 4 ,radcft(9,2),nrad(mpoin),astmp(mpoin) c...lnods(l,n) = global no. of n'th node of element l c...lnodn(l) = no. of nodes in element l ( = 3 for linear triangle) c...lsidn = no. of sides in element l ( = 3 for linear triangle) c...coord(np,1), coord(np,2) = x,y coords of node np c...nkode(np) = fixity code of node np c...aslod is global nodal load (i.e. right hand side) vector c...ndof(np) = no. of dof.s at node np (not used in this program) c...lpros(l) = mat. prop. set no. of element l c...props(m,i) = value of i'th mat. property in mat. prop. set m c...asdis is global nodal unknown u vector c...astmp is a scratch vector used in the equation solution substitutions c...dbdcmp,dbdval,dbdcrd are component, coord, value for Dirichlet c... boundary value sets. c...gstif holds assembled global stiffness matrix in compact form c...lmtyp(l) = material type of elt l (not used in this program) c...radcft(i,1) is the radiation coefft for boundary set i c...radcft(i,2) is the prescribed flow gradient for bdry set i c...nrad(n) is the radiating bdry set of node n (0 if not on bdry) c c c...ifull = 0 for demo version: max of 50 nodes or elements ifull = 0 c write(6,*)'Program POISS: solves 2D Poisson equation.' write(6,*)'Handles linear triangle elements.' write(6,*)'Dimensioned for',melem,' elements, and',mpoin,' nodes.' if (ifull.eq.0)write(6,*)'Demonstration datafile is podemo.dat' write(6,*) c c...read an existing .dat file c 1 continue c c...get the input file name in infile call askfil(infile,iflag) if (iflag.eq.1) then call yesno(char,'Do you want to exit?') if (char.eq.'y') stop goto 1 endif c write(6,*)'Opening file ',infile,'...' lenin = lenstr(infile) open(15,file=infile(1:lenin),status='old',err=999) write(6,*) lm4 = lenin - 4 prtfil = infile(1:lm4)//'.prt' write(6,*)'Creating file ',prtfil(1:lenin),' for results...' write(6,*) open(16,file=prtfil(1:lenin),status='unknown',err=999) outfil = infile(1:lm4)//'.out' write(6,*)'Creating file ',outfil(1:lenin),' for output...' write(6,*) open(17,file=outfil(1:lenin),status='unknown',err=999) c c...read the input mesh data, and echo it to the output files call input(15,melem,mpoin,mnode,msets,mcomp,nelem, 1 npoin,lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props 2 ,ndofn,npros,title,ndof) call print(16,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) call output(17,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) call yesno(char,'Do you want the mesh details listed?') if (char.eq.'y') then call print(6,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) endif c iflag = 0 do 3 l=1,nelem nnode = lnodn(l) nside = lsidn(l) nodsid = 10*nnode + nside if (nodsid.ne.33) then iflag = 1 write(6,*)'ERROR: Cannot handle element ',l write(6,*)'which has ',nside,' sides and ',nnode,' nodes' endif 3 continue if (iflag.eq.1) then call prompt stop endif c c...read data for radiating and Dirichlet boundaries, if present read(15,*) read(15,'(2i5)')nbrad,nbdir c if (nbrad.gt.0) then c...read radiating boundary node sets - store info in radcft, nrad call radby(mpoin,nrad,radcft,nrset) else nrset = 0 endif c if (nbdir.gt.0) then c...read Dirichlet boundary values, if present, from end of .dat file call bndry(ndbdry,dbdcmp,dbdcrd,dbdval) endif c c if (ifull.eq.0.and.(nelem.gt.51.or.npoin.gt.51)) then write(6,*)'Sorry, demonstration version will not process ' write(6,*)'meshes with more than 51 elements or 51 nodes.' call prompt stop endif c c...open a scratch file to write the element stiffness matrices on open(unit=18,status='scratch',form='unformatted') write(6,*)'Scratch file opened ...' write(6,*) c call prompt write(6,*) write(6,*)'Starting subroutine STIFF...' call stiff(melem,mpoin,mnode,msets,mcomp,ndofn, 1 nelem,npoin,lnods,lnodn,coord,lpros,props,aslod, 2 radcft,nrad,nrset) c call prompt write(6,*) write(6,*)'Starting subroutine ASSMB...' call assmb(melem,mpoin,mnode,mband,nelem,npoin,nhbw,lnods, 1 lnodn,gstif) c call prompt write(6,*) write(6,*)'Starting subroutine SOLVE...' call solve(mpoin,mband,npoin,nhbw,nkode,aslod,asdis,astmp, 1 gstif,coord,ndbdry,dbdcmp,dbdcrd,dbdval) c call prompt write(6,*) write(6,*)'Starting subroutine RESULT...' call result(melem,mpoin,mnode,msets,mcomp, 1 nelem,npoin,lnods,lnodn,coord,nkode,lpros,props,aslod,asdis) c write(6,*) write(6,*)'POISS completed. Results are in file ',prtfil(1:lenin) write(6,*) write(6,*) call prompt stop 999 write(6,*) 1'Cannot open file with this name. Terminating program.' call prompt stop end c c c...functions defining r.h.s. vector, and Dirichlet boundary values c function rhsfn(x,y) c...r.h.s. function f(x,y) only used in Poisson equation c...default here is f(x,y) = 1. Fill in your own choice!!! rhsfn = 1.0 return end c function bdryfn(x,y) c...Dirichlet boundary function (used if no values specified) c...The example here is g(x,y) = x*x - y*y. Fill in your own choice!! bdryfn = x*x - y*y return end c c c subroutine stiff(melem,mpoin,mnode,msets,mcomp,ndofn, 1 nelem,npoin,lnods,lnodn,coord,lpros,props,aslod, 2 radcft,nrad,nrset) c...calculates the element stiffness matrices, and writes them c...to the scratch file. dimension props(msets,mcomp),lnods(melem,mnode), 1aslod(mpoin),estif(8,8),lpros(melem),lnodn(melem), 2coord(mpoin,2),radcft(9,2),nrad(mpoin) c rewind 18 c c...initialize the r.h.s. vector do 10 n=1,npoin aslod(n) = 0.0 10 continue c...start loop over all elements do 80 l=1,nelem c c...extract coefficients ax, ay from props (ax=ay=1 for Laplace) lset = lpros(l) ax = props(lset,1) ay = props(lset,2) c c...zero stiffness matrix do 20 j=1,8 do 20 i=1,8 estif(i,j) = 0.0 20 continue c c...we use one-point Gauss integration over the triangle. c...put local coords of gauss-pt in s,t: s = 1./3. t = 1./3. weight = 0.5 c c...get x,y coords of the nodes, and the Gauss point: np1 = lnods(l,1) x1 = coord(np1,1) y1 = coord(np1,2) np2 = lnods(l,2) x2 = coord(np2,1) y2 = coord(np2,2) np3 = lnods(l,3) x3 = coord(np3,1) y3 = coord(np3,2) x21 = x2 - x1 x31 = x3 - x1 x32 = x3 - x2 y21 = y2 - y1 y31 = y3 - y1 y32 = y3 - y2 xg = x1 + x21*s + x31*t yg = y1 + y21*s + y31*t c c...find determinant of the Jacobian matrix detj = x21*y31 - x31*y21 c c...now write down the basis functions and their Cartesian derivatives phi1 = 1. - s - t phi2 = s phi3 = t dphi1x = -y32/detj dphi1y = x32/detj dphi2x = y31/detj dphi2y = -x31/detj dphi3x = -y21/detj dphi3y = x21/detj c c...the integration factor: factor = weight*detj c c...now form the element stiffness matrix estif(1,1) = (ax*dphi1x*dphi1x + ay*dphi1y*dphi1y)*factor estif(1,2) = (ax*dphi1x*dphi2x + ay*dphi1y*dphi2y)*factor estif(1,3) = (ax*dphi1x*dphi3x + ay*dphi1y*dphi3y)*factor estif(2,2) = (ax*dphi2x*dphi2x + ay*dphi2y*dphi2y)*factor estif(2,3) = (ax*dphi2x*dphi3x + ay*dphi2y*dphi3y)*factor estif(3,3) = (ax*dphi3x*dphi3x + ay*dphi3y*dphi3y)*factor c c...add in contribution to right-hand-side vector (not Laplace) f = rhsfn(xg,yg) if (f.ne.0.d0) then aslod(np1) = aslod(np1) + f*phi1*factor aslod(np2) = aslod(np2) + f*phi2*factor aslod(np3) = aslod(np3) + f*phi3*factor endif c c...add in contributions from radiating boundaries if (nrset.gt.0) then c...check each side in turn nr1 = nrad(np1) nr2 = nrad(np2) nr3 = nrad(np3) if(nr1.gt.0.and.nr1.eq.nr2) then c...side joining nodes 1 and 2 is on radiating boundary nr1 c...one-point integration along side: evaluate at mid-pt c...and multiply by length of side. In all cases the two c...relevant phi functions are 0.5, so c...phi_i*phi_j = 0.5*0.5 = 0.25 c...also subtract integral of qbar*phi_i from the c...nodal component in rhs vector aslod alpha = radcft(nr1,1) qbar = radcft(nr1,2) dist = sqrt(x21*x21 + y21*y21) estif(1,1) = estif(1,1) + alpha*0.25*dist estif(1,2) = estif(1,2) + alpha*0.25*dist estif(2,2) = estif(2,2) + alpha*0.25*dist aslod(np1) = aslod(np1) - qbar*0.5*dist aslod(np2) = aslod(np2) - qbar*0.5*dist endif if (nr1.gt.0.and.nr1.eq.nr3) then c...side joining nodes 1 and 3 is on radiating boundary nr1 alpha = radcft(nr1,1) qbar = radcft(nr1,2) dist = sqrt(x31*x31 + y31*y31) estif(1,1) = estif(1,1) + alpha*0.25*dist estif(1,3) = estif(1,3) + alpha*0.25*dist estif(3,3) = estif(3,3) + alpha*0.25*dist aslod(np1) = aslod(np1) - qbar*0.5*dist aslod(np3) = aslod(np3) - qbar*0.5*dist endif if (nr2.gt.0.and.nr2.eq.nr3) then c...side joining nodes 2 and 3 is on radiating boundary nr1 alpha = radcft(nr2,1) qbar = radcft(nr2,2) dist = sqrt(x32*x32 + y32*y32) estif(2,2) = estif(2,2) + alpha*0.25*dist estif(2,3) = estif(2,3) + alpha*0.25*dist estif(3,3) = estif(3,3) + alpha*0.25*dist aslod(np2) = aslod(np2) - qbar*0.5*dist aslod(np3) = aslod(np3) - qbar*0.5*dist endif endif c c...fill in lower triangle, by symmetry estif(2,1) = estif(1,2) estif(3,1) = estif(1,3) estif(3,2) = estif(2,3) c c...write estif to file write (18) estif write(6,*)'Stiffness matrix for elt ',l,' written to file.' 80 continue c...end of element loop c write(6,*)'Elt stiffness matrices written to file' endfile 18 return end