1!----------------------------------------------------------------------- 2! 3! Header with common subroutine 4! 5 include 't320-basis-f.h' 6! 7! Header with QFunctions 8! 9 include 't510-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 stridesutet(3),stridesuhex(3) 17 integer erestrictxtet,erestrictutet,erestrictuitet,& 18& erestrictxhex,erestrictuhex,erestrictuihex 19 integer bxtet,butet,bxhex,buhex 20 integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 21 integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 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,mass 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,qdatatet,err) 76 call ceedvectorcreate(ceed,nqptshex,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 stridesutet=[1,qtet,qtet] 106 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,nqptstet,& 107 & stridesutet,erestrictuitet,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 &//'t510-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',1,ceed_eval_none,err) 124 125 call ceedqfunctioncreateinterior(ceed,1,mass,& 126 &SOURCE_DIR& 127 &//'t510-operator.h:mass'//char(0),qf_masstet,err) 128 call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err) 129 call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err) 130 call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err) 131 132! -- Operators 133! ---- Setup Tet 134 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 135 & ceed_qfunction_none,op_setuptet,err) 136 call ceedoperatorsetfield(op_setuptet,'weight',& 137 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 138 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 139 & bxtet,ceed_vector_active,err) 140 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 141 ceed_basis_none,qdatatet,err) 142! ---- Mass Tet 143 call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 144 & ceed_qfunction_none,op_masstet,err) 145 call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 146 ceed_basis_none,qdatatet,err) 147 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 148 & butet,ceed_vector_active,err) 149 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 150 & butet,ceed_vector_active,err) 151 152! Hex Elements 153 do i=0,nelemhex-1 154 col=mod(i,nx) 155 row=i/nx 156 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 157 do j=0,phex-1 158 do k=0,phex-1 159 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 160 enddo 161 enddo 162 enddo 163 164! -- Restrictions 165 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,d*ndofs,& 166 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 167 168 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,ndofs,& 169 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 170 stridesuhex=[1,qhex*qhex,qhex*qhex] 171 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 172 & 1,nqptshex,stridesuhex,erestrictuihex,err) 173 174! -- Bases 175 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 176 & bxhex,err) 177 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 178 & buhex,err) 179 180! -- QFunctions 181 call ceedqfunctioncreateinterior(ceed,1,setup,& 182 &SOURCE_DIR& 183 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 184 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 185 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 186 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 187 188 call ceedqfunctioncreateinterior(ceed,1,mass,& 189 &SOURCE_DIR& 190 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 191 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 192 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 193 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 194 195! -- Operators 196! ---- Setup Hex 197 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 198 & ceed_qfunction_none,op_setuphex,& 199 & err) 200 call ceedoperatorsetfield(op_setuphex,'_weight',& 201 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 202 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 203 & bxhex,ceed_vector_active,err) 204 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 205 ceed_basis_none,qdatahex,err) 206! ---- Mass Hex 207 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 208 & ceed_qfunction_none,op_masshex,& 209 & err) 210 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 211 ceed_basis_none,qdatahex,err) 212 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 213 & buhex,ceed_vector_active,err) 214 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 215 & buhex,ceed_vector_active,err) 216 217! Composite Operators 218 call ceedoperatorcreatecomposite(ceed,op_setup,err) 219 call ceedoperatorcompositeaddsub(op_setup,op_setuptet,err) 220 call ceedoperatorcompositeaddsub(op_setup,op_setuphex,err) 221 222 call ceedoperatorcreatecomposite(ceed,op_mass,err) 223 call ceedoperatorcompositeaddsub(op_mass,op_masstet,err) 224 call ceedoperatorcompositeaddsub(op_mass,op_masshex,err) 225 226! Apply Setup Operator 227 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 228 & ceed_request_immediate,err) 229 230! Apply Mass Operator 231 call ceedvectorcreate(ceed,ndofs,u,err) 232 call ceedvectorsetvalue(u,1.d0,err) 233 call ceedvectorcreate(ceed,ndofs,v,err) 234 call ceedvectorsetvalue(v,0.d0,err) 235 236 call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err) 237 238! Check Output 239 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 240 total=0. 241 do i=1,ndofs 242 total=total+hv(voffset+i) 243 enddo 244 if (abs(total-1.)>1.0d-10) then 245! LCOV_EXCL_START 246 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 247! LCOV_EXCL_STOP 248 endif 249 call ceedvectorrestorearrayread(v,hv,voffset,err) 250 251 call ceedvectorsetvalue(v,1.d0,err) 252 call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err) 253 254! Check Output 255 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 256 total=-ndofs 257 do i=1,ndofs 258 total=total+hv(voffset+i) 259 enddo 260 if (abs(total-1.)>1.0d-10) then 261! LCOV_EXCL_START 262 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 263! LCOV_EXCL_STOP 264 endif 265 call ceedvectorrestorearrayread(v,hv,voffset,err) 266 267! Cleanup 268 call ceedqfunctiondestroy(qf_setuptet,err) 269 call ceedqfunctiondestroy(qf_masstet,err) 270 call ceedoperatordestroy(op_setuptet,err) 271 call ceedoperatordestroy(op_masstet,err) 272 call ceedqfunctiondestroy(qf_setuphex,err) 273 call ceedqfunctiondestroy(qf_masshex,err) 274 call ceedoperatordestroy(op_setuphex,err) 275 call ceedoperatordestroy(op_masshex,err) 276 call ceedoperatordestroy(op_setup,err) 277 call ceedoperatordestroy(op_mass,err) 278 call ceedelemrestrictiondestroy(erestrictutet,err) 279 call ceedelemrestrictiondestroy(erestrictxtet,err) 280 call ceedelemrestrictiondestroy(erestrictuitet,err) 281 call ceedelemrestrictiondestroy(erestrictuhex,err) 282 call ceedelemrestrictiondestroy(erestrictxhex,err) 283 call ceedelemrestrictiondestroy(erestrictuihex,err) 284 call ceedbasisdestroy(butet,err) 285 call ceedbasisdestroy(bxtet,err) 286 call ceedbasisdestroy(buhex,err) 287 call ceedbasisdestroy(bxhex,err) 288 call ceedvectordestroy(x,err) 289 call ceedvectordestroy(u,err) 290 call ceedvectordestroy(v,err) 291 call ceedvectordestroy(qdatatet,err) 292 call ceedvectordestroy(qdatahex,err) 293 call ceeddestroy(ceed,err) 294 end 295!----------------------------------------------------------------------- 296