1*52bfb9bbSJeremy L Thompson !----------------------------------------------------------------------- 2*52bfb9bbSJeremy L Thompson subroutine buildmats(qref,qweight,interp,grad) 3*52bfb9bbSJeremy L Thompson integer p,q,d 4*52bfb9bbSJeremy L Thompson parameter(p=6) 5*52bfb9bbSJeremy L Thompson parameter(q=4) 6*52bfb9bbSJeremy L Thompson parameter(d=2) 7*52bfb9bbSJeremy L Thompson 8*52bfb9bbSJeremy L Thompson real*8 qref(d*q) 9*52bfb9bbSJeremy L Thompson real*8 qweight(q) 10*52bfb9bbSJeremy L Thompson real*8 interp(p*q) 11*52bfb9bbSJeremy L Thompson real*8 grad(d*p*q) 12*52bfb9bbSJeremy L Thompson real*8 x1,x2 13*52bfb9bbSJeremy L Thompson 14*52bfb9bbSJeremy L Thompson qref=(/2.d-1,6.d-1,1.d0/3.d0,2.d-1,2.d-1,2.d-1,1.d0/3.d0,6.d-1/) 15*52bfb9bbSJeremy L Thompson qweight=(/25.d0/96.d0,25.d0/96.d0,-27.d0/96.d0,25.d0/96.d0/) 16*52bfb9bbSJeremy L Thompson 17*52bfb9bbSJeremy L Thompson do i=0,q-1 18*52bfb9bbSJeremy L Thompson x1 = qref(0*q+i+1) 19*52bfb9bbSJeremy L Thompson x2 = qref(1*q+i+1); 20*52bfb9bbSJeremy L Thompson ! Interp 21*52bfb9bbSJeremy L Thompson interp(i*P+1)=2.*(x1+x2-1.)*(x1+x2-1./2.); 22*52bfb9bbSJeremy L Thompson interp(i*P+2)=-4.*x1*(x1+x2-1.); 23*52bfb9bbSJeremy L Thompson interp(i*P+3)=2.*x1*(x1-1./2.); 24*52bfb9bbSJeremy L Thompson interp(i*P+4)=-4.*x2*(x1+x2-1.); 25*52bfb9bbSJeremy L Thompson interp(i*P+5)=4.*x1*x2; 26*52bfb9bbSJeremy L Thompson interp(i*P+6)=2.*x2*(x2-1./2.); 27*52bfb9bbSJeremy L Thompson ! Grad 28*52bfb9bbSJeremy L Thompson grad((i+0)*P+1)=2.*(1.*(x1+x2-1./2.)+(x1+x2-1.)*1.); 29*52bfb9bbSJeremy L Thompson grad((i+Q)*P+1)=2.*(1.*(x1+x2-1./2.)+(x1+x2-1.)*1.); 30*52bfb9bbSJeremy L Thompson grad((i+0)*P+2)=-4.*(1.*(x1+x2-1.)+x1*1.); 31*52bfb9bbSJeremy L Thompson grad((i+Q)*P+2)=-4.*(x1*1.); 32*52bfb9bbSJeremy L Thompson grad((i+0)*P+3)=2.*(1.*(x1-1./2.)+x1*1.); 33*52bfb9bbSJeremy L Thompson grad((i+Q)*P+3)=2.*0.; 34*52bfb9bbSJeremy L Thompson grad((i+0)*P+4)=-4.*(x2*1.); 35*52bfb9bbSJeremy L Thompson grad((i+Q)*P+4)=-4.*(1.*(x1+x2-1.)+x2*1.); 36*52bfb9bbSJeremy L Thompson grad((i+0)*P+5)=4.*(1.*x2); 37*52bfb9bbSJeremy L Thompson grad((i+Q)*P+5)=4.*(x1*1.); 38*52bfb9bbSJeremy L Thompson grad((i+0)*P+6)=2.*0.; 39*52bfb9bbSJeremy L Thompson grad((i+Q)*P+6)=2.*(1.*(x2-1./2.)+x2*1.); 40*52bfb9bbSJeremy L Thompson enddo 41*52bfb9bbSJeremy L Thompson 42*52bfb9bbSJeremy L Thompson end 43*52bfb9bbSJeremy L Thompson !----------------------------------------------------------------------- 44