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