1!----------------------------------------------------------------------- 2! 3! Header with QFunctions 4! 5 include 't534-operator-f.h' 6!----------------------------------------------------------------------- 7 program test 8 implicit none 9 include 'ceed/fortran.h' 10 11 integer ceed,err,i,j,k 12 integer stridesu(3),stridesqd(3) 13 integer erestrictx,erestrictu,erestrictui,erestrictqi 14 integer bx,bu 15 integer qf_setup,qf_diff 16 integer op_setup,op_diff 17 integer qdata,x,a,u,v 18 integer nelem,p,q,d 19 integer row,col,offset 20 parameter(nelem=6) 21 parameter(p=3) 22 parameter(q=4) 23 parameter(d=2) 24 integer ndofs,nqpts,nx,ny 25 parameter(nx=3) 26 parameter(ny=2) 27 parameter(ndofs=(nx*2+1)*(ny*2+1)) 28 parameter(nqpts=nelem*q*q) 29 integer indx(nelem*p*p) 30 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 31 integer*8 xoffset,aoffset,uoffset,voffset 32 33 character arg*32 34 35 external setup,diff,diff_lin 36 37 call getarg(1,arg) 38 39 call ceedinit(trim(arg)//char(0),ceed,err) 40 41! DoF Coordinates 42 do i=0,nx*2 43 do j=0,ny*2 44 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 45 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 46 enddo 47 enddo 48 call ceedvectorcreate(ceed,d*ndofs,x,err) 49 xoffset=0 50 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 51 52! Qdata Vector 53 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata,err) 54 55! Element Setup 56 do i=0,nelem-1 57 col=mod(i,nx) 58 row=i/nx 59 offset=col*(p-1)+row*(nx*2+1)*(p-1) 60 do j=0,p-1 61 do k=0,p-1 62 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 63 enddo 64 enddo 65 enddo 66 67! Restrictions 68 call ceedelemrestrictioncreate(ceed,nelem,p*p,d,ndofs,d*ndofs,& 69 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 70 71 call ceedelemrestrictioncreate(ceed,nelem,p*p,1,1,ndofs,& 72 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 73 stridesu=[1,q*q,q*q] 74 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,1,nqpts,& 75 & stridesu,erestrictui,err) 76 77 stridesqd=[1,q*q,q*q*d*(d+1)/2] 78 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,d*(d+1)/2,& 79 & d*(d+1)/2*nqpts,stridesqd,erestrictqi,err) 80 81! Bases 82 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 83 & bx,err) 84 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 85 & bu,err) 86 87! QFunction - setup 88 call ceedqfunctioncreateinterior(ceed,1,setup,& 89 &SOURCE_DIR& 90 &//'t531-operator.h:setup'//char(0),qf_setup,err) 91 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 92 call ceedqfunctionaddinput(qf_setup,'weight',1,ceed_eval_weight,err) 93 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 94 95! Operator - setup 96 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 97 & ceed_qfunction_none,op_setup,err) 98 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 99 & bx,ceed_vector_active,err) 100 call ceedoperatorsetfield(op_setup,'weight',ceed_elemrestriction_none,& 101 & bx,ceed_vector_none,err) 102 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 103 ceed_basis_none,ceed_vector_active,err) 104 105! Apply Setup Operator 106 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 107 108! QFunction - apply 109 call ceedqfunctioncreateinterior(ceed,1,diff,& 110 &SOURCE_DIR& 111 &//'t531-operator.h:diff'//char(0),qf_diff,err) 112 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 113 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 114 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 115 116! Operator - apply 117 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 118 & ceed_qfunction_none,op_diff,err) 119 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 120 & bu,ceed_vector_active,err) 121 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 122 ceed_basis_none,qdata,err) 123 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 124 & bu,ceed_vector_active,err) 125 126! Assemble Diagonal 127 call ceedvectorcreate(ceed,ndofs,a,err) 128 call ceedoperatorlinearassemblediagonal(op_diff,a,& 129 & ceed_request_immediate,err) 130 131! Manually assemble diagonal 132 call ceedvectorcreate(ceed,ndofs,u,err) 133 call ceedvectorsetvalue(u,0.d0,err) 134 call ceedvectorcreate(ceed,ndofs,v,err) 135 do i=1,ndofs 136 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 137 uu(i+uoffset)=1.d0 138 if (i>1) then 139 uu(i-1+uoffset)=0.d0 140 endif 141 call ceedvectorrestorearray(u,uu,uoffset,err) 142 143 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 144 145 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 146 atrue(i)=vv(voffset+i) 147 call ceedvectorrestorearrayread(v,vv,voffset,err) 148 enddo 149 150! Check Output 151 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 152 do i=1,ndofs 153 if (abs(aa(aoffset+i)-atrue(i))>1.0d-13) then 154! LCOV_EXCL_START 155 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 156 & atrue(i) 157! LCOV_EXCL_STOP 158 endif 159 enddo 160 call ceedvectorrestorearrayread(a,aa,aoffset,err) 161 162! Cleanup 163 call ceedqfunctiondestroy(qf_setup,err) 164 call ceedqfunctiondestroy(qf_diff,err) 165 call ceedoperatordestroy(op_setup,err) 166 call ceedoperatordestroy(op_diff,err) 167 call ceedelemrestrictiondestroy(erestrictu,err) 168 call ceedelemrestrictiondestroy(erestrictx,err) 169 call ceedelemrestrictiondestroy(erestrictui,err) 170 call ceedelemrestrictiondestroy(erestrictqi,err) 171 call ceedbasisdestroy(bu,err) 172 call ceedbasisdestroy(bx,err) 173 call ceedvectordestroy(x,err) 174 call ceedvectordestroy(a,err) 175 call ceedvectordestroy(u,err) 176 call ceedvectordestroy(v,err) 177 call ceedvectordestroy(qdata,err) 178 call ceeddestroy(ceed,err) 179 end 180!----------------------------------------------------------------------- 181