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