1!----------------------------------------------------------------------- 2! 3! Header with common subroutine 4! 5 include 't320-basis-f.h' 6! 7! Header with QFunctions 8! 9 include 't522-operator-f.h' 10!----------------------------------------------------------------------- 11 program test 12 implicit none 13 include 'ceed/fortran.h' 14 15 integer ceed,err,i,j,k 16 integer stridesqdtet(3),stridesqdhex(3) 17 integer erestrictxtet,erestrictutet,erestrictqditet,& 18& erestrictxhex,erestrictuhex,erestrictqdihex 19 integer bxtet,butet,bxhex,buhex 20 integer qf_setuptet,qf_difftet,qf_setuphex,qf_diffhex 21 integer op_setuptet,op_difftet,op_setuphex,op_diffhex,op_setup,op_diff 22 integer qdatatet,qdatahex,x,u,v 23 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 24 integer row,col,offset 25 parameter(nelemtet=6) 26 parameter(ptet=6) 27 parameter(qtet=4) 28 parameter(nelemhex=6) 29 parameter(phex=3) 30 parameter(qhex=4) 31 parameter(d=2) 32 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 33 parameter(nx=3) 34 parameter(ny=3) 35 parameter(nxtet=3) 36 parameter(nytet=1) 37 parameter(nxhex=3) 38 parameter(ndofs=(nx*2+1)*(ny*2+1)) 39 parameter(nqptstet=nelemtet*qtet) 40 parameter(nqptshex=nelemhex*qhex*qhex) 41 parameter(nqpts=nqptstet+nqptshex) 42 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 43 real*8 arrx(d*ndofs) 44 integer*8 voffset,xoffset 45 46 real*8 qref(d*qtet) 47 real*8 qweight(qtet) 48 real*8 interp(ptet*qtet) 49 real*8 grad(d*ptet*qtet) 50 51 real*8 hv(ndofs) 52 real*8 total 53 54 character arg*32 55 56 external setup,diff 57 58 call getarg(1,arg) 59 60 call ceedinit(trim(arg)//char(0),ceed,err) 61 62! DoF Coordinates 63 do i=0,ny*2 64 do j=0,nx*2 65 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 66 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 67 enddo 68 enddo 69 70 call ceedvectorcreate(ceed,d*ndofs,x,err) 71 xoffset=0 72 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 73 74! Qdata Vectors 75 call ceedvectorcreate(ceed,nqptstet*d*(d+1)/2,qdatatet,err) 76 call ceedvectorcreate(ceed,nqptshex*d*(d+1)/2,qdatahex,err) 77 78! Tet Elements 79 do i=0,2 80 col=mod(i,nx) 81 row=i/nx 82 offset=col*2+row*(nx*2+1)*2 83 84 indxtet(i*2*ptet+1)=2+offset 85 indxtet(i*2*ptet+2)=9+offset 86 indxtet(i*2*ptet+3)=16+offset 87 indxtet(i*2*ptet+4)=1+offset 88 indxtet(i*2*ptet+5)=8+offset 89 indxtet(i*2*ptet+6)=0+offset 90 91 indxtet(i*2*ptet+7)=14+offset 92 indxtet(i*2*ptet+8)=7+offset 93 indxtet(i*2*ptet+9)=0+offset 94 indxtet(i*2*ptet+10)=15+offset 95 indxtet(i*2*ptet+11)=8+offset 96 indxtet(i*2*ptet+12)=16+offset 97 enddo 98 99! -- Restrictions 100 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,& 101 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 102 103 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,& 104 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 105 stridesqdtet=[1,qtet,qtet*d*(d+1)/2] 106 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,d*(d+1)/2,& 107 & d*(d+1)/2*nqptstet,stridesqdtet,erestrictqditet,err) 108 109! -- Bases 110 call buildmats(qref,qweight,interp,grad) 111 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 112 & qweight,bxtet,err) 113 call buildmats(qref,qweight,interp,grad) 114 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 115 & qweight,butet,err) 116 117! -- QFunctions 118 call ceedqfunctioncreateinterior(ceed,1,setup,& 119 &SOURCE_DIR& 120 &//'t522-operator.h:setup'//char(0),qf_setuptet,err) 121 call ceedqfunctionaddinput(qf_setuptet,'weight',1,ceed_eval_weight,err) 122 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 123 call ceedqfunctionaddoutput(qf_setuptet,'rho',d*(d+1)/2,ceed_eval_none,& 124 & err) 125 126 call ceedqfunctioncreateinterior(ceed,1,diff,& 127 &SOURCE_DIR& 128 &//'t522-operator.h:diff'//char(0),qf_difftet,err) 129 call ceedqfunctionaddinput(qf_difftet,'rho',d*(d+1)/2,ceed_eval_none,err) 130 call ceedqfunctionaddinput(qf_difftet,'u',d,ceed_eval_grad,err) 131 call ceedqfunctionaddoutput(qf_difftet,'v',d,ceed_eval_grad,err) 132 133! -- Operators 134! ---- Setup Tet 135 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 136 & ceed_qfunction_none,op_setuptet,err) 137 call ceedoperatorsetfield(op_setuptet,'weight',& 138 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 139 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 140 & bxtet,ceed_vector_active,err) 141 call ceedoperatorsetfield(op_setuptet,'rho',erestrictqditet,& 142 ceed_basis_none,qdatatet,err) 143! ---- diff Tet 144 call ceedoperatorcreate(ceed,qf_difftet,ceed_qfunction_none,& 145 & ceed_qfunction_none,op_difftet,err) 146 call ceedoperatorsetfield(op_difftet,'rho',erestrictqditet,& 147 ceed_basis_none,qdatatet,err) 148 call ceedoperatorsetfield(op_difftet,'u',erestrictutet,& 149 & butet,ceed_vector_active,err) 150 call ceedoperatorsetfield(op_difftet,'v',erestrictutet,& 151 & butet,ceed_vector_active,err) 152 153! Hex Elements 154 do i=0,nelemhex-1 155 col=mod(i,nx) 156 row=i/nx 157 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 158 do j=0,phex-1 159 do k=0,phex-1 160 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 161 enddo 162 enddo 163 enddo 164 165! -- Restrictions 166 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,& 167 & d*ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 168 169 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,ndofs,& 170 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 171 stridesqdhex=[1,qhex*qhex,qhex*qhex*d*(d+1)/2] 172 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 173 & d*(d+1)/2,d*(d+1)/2*nqptshex,stridesqdhex,erestrictqdihex,err) 174 175! -- Bases 176 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 177 & bxhex,err) 178 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 179 & buhex,err) 180 181! -- QFunctions 182 call ceedqfunctioncreateinterior(ceed,1,setup,& 183 &SOURCE_DIR& 184 &//'t522-operator.h:setup'//char(0),qf_setuphex,err) 185 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 186 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 187 call ceedqfunctionaddoutput(qf_setuphex,'rho',d*(d+1)/2,ceed_eval_none,& 188 & err) 189 190 call ceedqfunctioncreateinterior(ceed,1,diff,& 191 &SOURCE_DIR& 192 &//'t522-operator.h:diff'//char(0),qf_diffhex,err) 193 call ceedqfunctionaddinput(qf_diffhex,'rho',d*(d+1)/2,ceed_eval_none,err) 194 call ceedqfunctionaddinput(qf_diffhex,'u',d,ceed_eval_grad,err) 195 call ceedqfunctionaddoutput(qf_diffhex,'v',d,ceed_eval_grad,err) 196 197! -- Operators 198! ---- Setup Hex 199 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 200 & ceed_qfunction_none,op_setuphex,err) 201 call ceedoperatorsetfield(op_setuphex,'_weight',& 202 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 203 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 204 & bxhex,ceed_vector_active,err) 205 call ceedoperatorsetfield(op_setuphex,'rho',erestrictqdihex,& 206 ceed_basis_none,qdatahex,err) 207! ---- diff Hex 208 call ceedoperatorcreate(ceed,qf_diffhex,ceed_qfunction_none,& 209 & ceed_qfunction_none,op_diffhex,err) 210 call ceedoperatorsetfield(op_diffhex,'rho',erestrictqdihex,& 211 ceed_basis_none,qdatahex,err) 212 call ceedoperatorsetfield(op_diffhex,'u',erestrictuhex,& 213 & buhex,ceed_vector_active,err) 214 call ceedoperatorsetfield(op_diffhex,'v',erestrictuhex,& 215 & buhex,ceed_vector_active,err) 216 217! Composite Operators 218 call ceedcompositeoperatorcreate(ceed,op_setup,err) 219 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 220 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 221 222 call ceedcompositeoperatorcreate(ceed,op_diff,err) 223 call ceedcompositeoperatoraddsub(op_diff,op_difftet,err) 224 call ceedcompositeoperatoraddsub(op_diff,op_diffhex,err) 225 226! Apply Setup Operator 227 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 228 & ceed_request_immediate,err) 229 230! Apply diff Operator 231 call ceedvectorcreate(ceed,ndofs,u,err) 232 call ceedvectorsetvalue(u,1.d0,err) 233 call ceedvectorcreate(ceed,ndofs,v,err) 234 235 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 236 237! Check Output 238 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 239 do i=1,ndofs 240 if (abs(hv(voffset+i))>1.0d-14) then 241! LCOV_EXCL_START 242 write(*,*) 'Computed: ',total,' != True: 0.0' 243! LCOV_EXCL_STOP 244 endif 245 enddo 246 call ceedvectorrestorearrayread(v,hv,voffset,err) 247 248! Cleanup 249 call ceedqfunctiondestroy(qf_setuptet,err) 250 call ceedqfunctiondestroy(qf_difftet,err) 251 call ceedoperatordestroy(op_setuptet,err) 252 call ceedoperatordestroy(op_difftet,err) 253 call ceedqfunctiondestroy(qf_setuphex,err) 254 call ceedqfunctiondestroy(qf_diffhex,err) 255 call ceedoperatordestroy(op_setuphex,err) 256 call ceedoperatordestroy(op_diffhex,err) 257 call ceedoperatordestroy(op_setup,err) 258 call ceedoperatordestroy(op_diff,err) 259 call ceedelemrestrictiondestroy(erestrictutet,err) 260 call ceedelemrestrictiondestroy(erestrictxtet,err) 261 call ceedelemrestrictiondestroy(erestrictqditet,err) 262 call ceedelemrestrictiondestroy(erestrictuhex,err) 263 call ceedelemrestrictiondestroy(erestrictxhex,err) 264 call ceedelemrestrictiondestroy(erestrictqdihex,err) 265 call ceedbasisdestroy(butet,err) 266 call ceedbasisdestroy(bxtet,err) 267 call ceedbasisdestroy(buhex,err) 268 call ceedbasisdestroy(bxhex,err) 269 call ceedvectordestroy(x,err) 270 call ceedvectordestroy(u,err) 271 call ceedvectordestroy(v,err) 272 call ceedvectordestroy(qdatatet,err) 273 call ceedvectordestroy(qdatahex,err) 274 call ceeddestroy(ceed,err) 275 end 276!----------------------------------------------------------------------- 277