1*8980d4a7Sjeremylt!----------------------------------------------------------------------- 2*8980d4a7Sjeremylt! 3*8980d4a7Sjeremylt! Header with common subroutine 4*8980d4a7Sjeremylt! 5*8980d4a7Sjeremylt include 't310-basis-f.h' 6*8980d4a7Sjeremylt!----------------------------------------------------------------------- 7*8980d4a7Sjeremylt subroutine feval(x1,x2,val) 8*8980d4a7Sjeremylt real*8 x1,x2,val 9*8980d4a7Sjeremylt 10*8980d4a7Sjeremylt val=x1*x1+x2*x2+x1*x2+1 11*8980d4a7Sjeremylt 12*8980d4a7Sjeremylt end 13*8980d4a7Sjeremylt!----------------------------------------------------------------------- 14*8980d4a7Sjeremylt subroutine dfeval(x1,x2,val) 15*8980d4a7Sjeremylt real*8 x1,x2,val 16*8980d4a7Sjeremylt 17*8980d4a7Sjeremylt val=2*x1+x2 18*8980d4a7Sjeremylt 19*8980d4a7Sjeremylt end 20*8980d4a7Sjeremylt!----------------------------------------------------------------------- 21*8980d4a7Sjeremylt program test 22*8980d4a7Sjeremylt 23*8980d4a7Sjeremylt include 'ceedf.h' 24*8980d4a7Sjeremylt 25*8980d4a7Sjeremylt integer ceed,err 26*8980d4a7Sjeremylt integer input,output 27*8980d4a7Sjeremylt integer p,q,d 28*8980d4a7Sjeremylt parameter(p=6) 29*8980d4a7Sjeremylt parameter(q=4) 30*8980d4a7Sjeremylt parameter(d=2) 31*8980d4a7Sjeremylt 32*8980d4a7Sjeremylt real*8 qref(d*q) 33*8980d4a7Sjeremylt real*8 qweight(q) 34*8980d4a7Sjeremylt real*8 interp(p*q) 35*8980d4a7Sjeremylt real*8 grad(d*p*q) 36*8980d4a7Sjeremylt real*8 xq(d*q) 37*8980d4a7Sjeremylt real*8 xr(d*p) 38*8980d4a7Sjeremylt real*8 iinput(p) 39*8980d4a7Sjeremylt real*8 ooutput(d*q) 40*8980d4a7Sjeremylt real*8 val,diff 41*8980d4a7Sjeremylt real*8 x1,x2 42*8980d4a7Sjeremylt integer*8 offset 43*8980d4a7Sjeremylt 44*8980d4a7Sjeremylt integer b 45*8980d4a7Sjeremylt 46*8980d4a7Sjeremylt character arg*32 47*8980d4a7Sjeremylt 48*8980d4a7Sjeremylt xq=(/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/) 49*8980d4a7Sjeremylt xr=(/0.d0,5.d-1,1.d0,0.d0,5.d-1,0.d0,0.d0,0.d0, 0.d0,5.d-1,5.d-1,1.d0/) 50*8980d4a7Sjeremylt 51*8980d4a7Sjeremylt call getarg(1,arg) 52*8980d4a7Sjeremylt 53*8980d4a7Sjeremylt call buildmats(qref,qweight,interp,grad) 54*8980d4a7Sjeremylt 55*8980d4a7Sjeremylt call ceedinit(trim(arg)//char(0),ceed,err) 56*8980d4a7Sjeremylt 57*8980d4a7Sjeremylt call ceedbasiscreateh1(ceed,ceed_triangle,1,p,q,interp,grad,qref,qweight,& 58*8980d4a7Sjeremylt & b,err) 59*8980d4a7Sjeremylt 60*8980d4a7Sjeremylt do i=1,p 61*8980d4a7Sjeremylt x1=xr(0*p+i) 62*8980d4a7Sjeremylt x2=xr(1*p+i) 63*8980d4a7Sjeremylt call feval(x1,x2,val) 64*8980d4a7Sjeremylt iinput(i)=val 65*8980d4a7Sjeremylt enddo 66*8980d4a7Sjeremylt 67*8980d4a7Sjeremylt call ceedvectorcreate(ceed,p,input,err) 68*8980d4a7Sjeremylt call ceedvectorsetarray(input,ceed_mem_host,ceed_use_pointer,iinput,err) 69*8980d4a7Sjeremylt call ceedvectorcreate(ceed,q*d,output,err) 70*8980d4a7Sjeremylt call ceedvectorsetvalue(output,0.d0,err) 71*8980d4a7Sjeremylt 72*8980d4a7Sjeremylt call ceedbasisapply(b,1,ceed_notranspose,ceed_eval_grad,input,output,err) 73*8980d4a7Sjeremylt 74*8980d4a7Sjeremylt call ceedvectorgetarrayread(output,ceed_mem_host,ooutput,offset,err) 75*8980d4a7Sjeremylt do i=1,q 76*8980d4a7Sjeremylt x1=xq(0*q+i) 77*8980d4a7Sjeremylt x2=xq(1*q+i) 78*8980d4a7Sjeremylt call dfeval(x1,x2,val) 79*8980d4a7Sjeremylt diff=val-ooutput(0*q+i+offset) 80*8980d4a7Sjeremylt if (abs(diff)>1.0d-10) then 81*8980d4a7Sjeremylt write(*,'(A,I1,A,F12.8,A,F12.8)') '[',i,'] ',ooutput(i+offset),& 82*8980d4a7Sjeremylt & ' != ',val 83*8980d4a7Sjeremylt endif 84*8980d4a7Sjeremylt call dfeval(x2,x1,val) 85*8980d4a7Sjeremylt diff=val-ooutput(1*q+i+offset) 86*8980d4a7Sjeremylt if (abs(diff)>1.0d-10) then 87*8980d4a7Sjeremylt write(*,'(A,I1,A,F12.8,A,F12.8)') '[',i,'] ',ooutput(i+offset),& 88*8980d4a7Sjeremylt & ' != ',val 89*8980d4a7Sjeremylt endif 90*8980d4a7Sjeremylt enddo 91*8980d4a7Sjeremylt call ceedvectorrestorearrayread(output,ooutput,offset,err) 92*8980d4a7Sjeremylt 93*8980d4a7Sjeremylt call ceedvectordestroy(input,err) 94*8980d4a7Sjeremylt call ceedvectordestroy(output,err) 95*8980d4a7Sjeremylt call ceedbasisdestroy(b,err) 96*8980d4a7Sjeremylt call ceeddestroy(ceed,err) 97*8980d4a7Sjeremylt 98*8980d4a7Sjeremylt end 99*8980d4a7Sjeremylt!----------------------------------------------------------------------- 100