c This is an example 'main' program for calling MGGHAT c program main c c In 'commons', parameter statements are used to set the dimension for c the arrays. Also, all of the program parameters are passed through it. c include 'commons' c c solve the problem c call mgghat stop end c c -------- pde c subroutine pde(x,y,p,q,r,f) real x,y,p,q,r,f c c return the values of the pde coefficents at (x,y) c pde is c c -( p(x,y) * u ) -( q(x,y) * u ) + r(x,y) * u = f(x,y) c x x y y c c NOTE: BE CAREFUL TO GET THE SIGNS RIGHT c e.g. p=q=1. means f=-(uxx+uyy) c ipower=5 p=1. q=1. r = 0. f=-ipower*(ipower-1.)*(x**(ipower-2)+y**(ipower-2)) c return end c c -------- bcond c subroutine bcond(x,y,ipiece,c,g,itype) real x,y,c,g integer ipiece,itype c c returns boundary condition coefficients at (x,y) c c boundary condition is either c c u + c(x,y)*u = g(x,y) or u = g(x,y) c n c c (the natural b.c. is the Neumann b.c. because p=q=1) c c ipiece indicates the boundary piece from which the boundary condition c is determined. itype must be set to the type of boundary condition for c that piece (the type cannot change within a piece). c itype = 1 Dirichlet (second condition above) c 2 Neuman (first condition with c = 0) c 3 Mixed (first condition with c != 0) c c (these pieces assume the assignment in the example inittr for c a rectangular domain) c if (ipiece.eq.1) then c left side; Dirichlet b.c.; U=g itype = 1 c = 0. g = true(x,y) elseif (ipiece.eq.2) then c bottom; Mixed b.c.; Un + U = g itype = 3 c = 1. g = true(x,y) - truey(x,y) elseif (ipiece.eq.3) then c right side; Neuman; Un = g itype = 2 c = 0. g = truex(x,y) else c top; dirichlet b.c.; U=g itype = 1 c = 0. g = true(x,y) endif c return end c c -------- true c real function true(x,y) real x,y ipower=5 true = x**ipower + y**ipower return end c c -------- truex c real function truex(x,y) real x,y ipower = 5 truex = ipower*x**(ipower-1) return end c c -------- truey c real function truey(x,y) real x,y ipower = 5 truey = ipower*y**(ipower-1) return end c c -------- inittr c subroutine inittr include 'commons' c c rectangular domain (ax,bx) X (ay,by) with an ngridx X ngridy grid c ax = 0. bx = 1. ay = 0. by = 1. ngridx = 2 ngridy = 2 c if (outlev.ge.3) write(ioutpt,101) ax,bx,ay,by 101 format(' begin initializing triangulation'/ . ' triangulation for rectangle (',f8.2, . ',',f8.2,') X (',f8.2,',',f8.2,')') c c set initial triangulation c c this is a user provided routine to define the domain and c the initial triangulation c c the user provides: c c nvert - number of vertices c ntri - number of triangles c xvert(1..nvert),yvert(1..nvert) - x and y coordinates c of the vertices c vertex(1..3,1..ntri) - the 3 vertices of each triangle c neigh(i,1..ntri) - which piece of the boundary contains the triangle c side opposite vertex i. Need not be set if c that side is not on the domain boundary. c c the peak of each triangle is the third vertex, i.e., c peak(triangle)=vertex(3,triangle) c c assumptions on the triangulation are: c c 1) each vertex is in the same position (first, second, or third) c in every triangle that contains it, e.g., if vertex 2 is c the first vertex of triangle 3 (vertex(1,3)=2). vertex 2 must c also be the first vertex of any other triangle containing it. c thus we can have vertex(1,4)=2, and cannot have vertex(2,4)=2 c c 2) vertices which are first or second vertices can be in at most c 8 triangles. vertices which are third vertices (peaks) can c be in at most 4 triangles c c This version of the routine triangulates the rectangle (ax,bx)X(ay,by). c ngridx and ngridy specify the number of intervals in each dimension. c The boundary pieces (ipiece for subroutine bcond) are: c 1 - left c 2 - bottom c 3 - right c 4 - top c c for the case of a 4X4 grid, the vertex and triangle numbers are: c c 5---10---15---20---25 c |\ | /|\ | /| c | \ 8|15/ | \24|31/ | c | 7\ | /16|23\ | /32| c | \|/ | \|/ | c 4----9---14---19---24 c | /|\ | /|\ | c | 5/ | \14|21/ | \30| c | /6 |13\ | /22|29\ | c |/ | \|/ | \| c 3----8---13---18---23 c |\ | /|\ | /| c | \4 |11/ | \20|27/ | c | 3\ | /12|19\ | /28| c | \|/ | \|/ | c 2----7---12---17---22 c | /|\ | /|\ | c | 1/ | \10|17/ | \26| c | /2 | 9\ | /18|25\ | c |/ | \|/ | \| c 1----6---11---16---21 c dy = by-ay dx = bx-ax dy = dy/ngridy dx = dx/ngridx c c set number of triangles and vertices c nvert = (ngridx+1)*(ngridy+1) ntri = 2*ngridx*ngridy c c set coordinates of vertices c k=0 do 11 i=1,ngridx+1 xtemp=ax+(i-1)*dx if (i.eq.ngridx+1) xtemp=bx do 10 j=1,ngridy+1 k=k+1 xvert(k)=xtemp yvert(k)=ay+(j-1)*dy if (j.eq.ngridy+1) yvert(k)=by 10 continue 11 continue c c set vertices and boundary pieces of triangles c do 21 i=1,ngridx,2 do 20 j=1,ngridy,2 ivbase=(i-1)*(ngridy+1)+j itbase=2*((i-1)*ngridy+j)-1 vertex(1,itbase) = ivbase+ngridy+2 vertex(2,itbase) = ivbase vertex(3,itbase) = ivbase+1 if (i.eq.1) neigh(1,itbase) = -1 if (j.eq.ngridy) neigh(2,itbase) = -4 vertex(1,itbase+1) = ivbase+ngridy+2 vertex(2,itbase+1) = ivbase vertex(3,itbase+1) = ivbase+ngridy+1 if (j.eq.1) neigh(1,itbase+1) = -2 if (i.eq.ngridx) neigh(2,itbase+1) = -3 if (j.ne.ngridy) then vertex(1,itbase+2) = ivbase+ngridy+2 vertex(2,itbase+2) = ivbase+2 vertex(3,itbase+2) = ivbase+1 if (i.eq.1) neigh(1,itbase+2) = -1 vertex(1,itbase+3) = ivbase+ngridy+2 vertex(2,itbase+3) = ivbase+2 vertex(3,itbase+3) = ivbase+ngridy+3 if (j.eq.ngridy-1) neigh(1,itbase+3) = -4 if (i.eq.ngridx) neigh(2,itbase+3) = -3 endif if (i.ne.ngridx) then vertex(1,itbase+2*ngridy) = ivbase+ngridy+2 vertex(2,itbase+2*ngridy) = ivbase+2*ngridy+2 vertex(3,itbase+2*ngridy) = ivbase+ngridy+1 if (j.eq.1) neigh(1,itbase+2*ngridy) = -2 vertex(1,itbase+2*ngridy+1) = ivbase+ngridy+2 vertex(2,itbase+2*ngridy+1) = ivbase+2*ngridy+2 vertex(3,itbase+2*ngridy+1) = ivbase+2*ngridy+3 if (i.eq.ngridx-1) neigh(1,itbase+2*ngridy+1) = -3 if (j.eq.ngridy) neigh(2,itbase+2*ngridy+1) = -4 if (j.ne.ngridy) then vertex(1,itbase+2*ngridy+2) = ivbase+ngridy+2 vertex(2,itbase+2*ngridy+2) = ivbase+2*ngridy+4 vertex(3,itbase+2*ngridy+2) = ivbase+ngridy+3 if (j.eq.ngridy-1) neigh(1,itbase+2*ngridy+2) = -4 vertex(1,itbase+2*ngridy+3) = ivbase+ngridy+2 vertex(2,itbase+2*ngridy+3) = ivbase+2*ngridy+4 vertex(3,itbase+2*ngridy+3) = ivbase+2*ngridy+3 if (i.eq.ngridx-1) neigh(1,itbase+2*ngridy+3) = -3 endif endif 20 continue 21 continue c if (outlev.ge.3) then write(ioutpt,102) ngridx,ngridy,ntri,nvert 102 format(' grid lines ',i3,' X',i3/ 1 ' triangles ',i5/ 2 ' vertices ',i5/ 3 ' initial triangulation complete') endif return end