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 erestrictxtet,erestrictutet,erestrictxitet,erestrictqditet,& 48& erestrictxhex,erestrictuhex,erestrictxihex,erestrictqdihex 49 integer bxtet,butet,bxhex,buhex 50 integer qf_setuptet,qf_difftet,qf_setuphex,qf_diffhex 51 integer op_setuptet,op_difftet,op_setuphex,op_diffhex,op_setup,op_diff 52 integer qdatatet,qdatahex,x,u,v 53 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 54 integer row,col,offset 55 parameter(nelemtet=6) 56 parameter(ptet=6) 57 parameter(qtet=4) 58 parameter(nelemhex=6) 59 parameter(phex=3) 60 parameter(qhex=4) 61 parameter(d=2) 62 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 63 parameter(nx=3) 64 parameter(ny=3) 65 parameter(nxtet=3) 66 parameter(nytet=1) 67 parameter(nxhex=3) 68 parameter(ndofs=(nx*2+1)*(ny*2+1)) 69 parameter(nqptstet=nelemtet*qtet) 70 parameter(nqptshex=nelemhex*qhex*qhex) 71 parameter(nqpts=nqptstet+nqptshex) 72 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 73 real*8 arrx(d*ndofs) 74 integer*8 voffset,xoffset 75 76 real*8 qref(d*qtet) 77 real*8 qweight(qtet) 78 real*8 interp(ptet*qtet) 79 real*8 grad(d*ptet*qtet) 80 81 real*8 hv(ndofs) 82 real*8 total 83 84 character arg*32 85 86 external setup,diff 87 88 call getarg(1,arg) 89 90 call ceedinit(trim(arg)//char(0),ceed,err) 91 92! DoF Coordinates 93 do i=0,ny*2 94 do j=0,nx*2 95 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 96 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 97 enddo 98 enddo 99 100 call ceedvectorcreate(ceed,d*ndofs,x,err) 101 xoffset=0 102 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 103 104! Qdata Vectors 105 call ceedvectorcreate(ceed,nqptstet*d*(d+1)/2,qdatatet,err) 106 call ceedvectorcreate(ceed,nqptshex*d*(d+1)/2,qdatahex,err) 107 108! Tet Elements 109 do i=0,2 110 col=mod(i,nx) 111 row=i/nx 112 offset=col*2+row*(nx*2+1)*2 113 114 indxtet(i*2*ptet+1)=2+offset 115 indxtet(i*2*ptet+2)=9+offset 116 indxtet(i*2*ptet+3)=16+offset 117 indxtet(i*2*ptet+4)=1+offset 118 indxtet(i*2*ptet+5)=8+offset 119 indxtet(i*2*ptet+6)=0+offset 120 121 indxtet(i*2*ptet+7)=14+offset 122 indxtet(i*2*ptet+8)=7+offset 123 indxtet(i*2*ptet+9)=0+offset 124 indxtet(i*2*ptet+10)=15+offset 125 indxtet(i*2*ptet+11)=8+offset 126 indxtet(i*2*ptet+12)=16+offset 127 enddo 128 129! -- Restrictions 130 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,ndofs,d,ceed_mem_host,& 131 & ceed_use_pointer,indxtet,erestrictxtet,err) 132 call ceedelemrestrictioncreateidentity(ceed,nelemtet,ptet,nelemtet*ptet,& 133 & d,erestrictxitet,err) 134 135 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,ndofs,1,ceed_mem_host,& 136 & ceed_use_pointer,indxtet,erestrictutet,err) 137 call ceedelemrestrictioncreateidentity(ceed,nelemtet,qtet,nqptstet,& 138 & d*(d+1)/2,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_null,ceed_null,op_setuptet,& 167 & err) 168 call ceedoperatorsetfield(op_setuptet,'_weight',erestrictxitet,& 169 & ceed_notranspose,bxtet,ceed_vector_none,err) 170 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 171 & ceed_notranspose,bxtet,ceed_vector_active,err) 172 call ceedoperatorsetfield(op_setuptet,'rho',erestrictqditet,& 173 & ceed_notranspose,ceed_basis_collocated,qdatatet,err) 174! ---- diff Tet 175 call ceedoperatorcreate(ceed,qf_difftet,ceed_null,ceed_null,op_difftet,& 176 & err) 177 call ceedoperatorsetfield(op_difftet,'rho',erestrictqditet,& 178 & ceed_notranspose,ceed_basis_collocated,qdatatet,err) 179 call ceedoperatorsetfield(op_difftet,'u',erestrictutet,& 180 & ceed_notranspose,butet,ceed_vector_active,err) 181 call ceedoperatorsetfield(op_difftet,'v',erestrictutet,& 182 & ceed_notranspose,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,ndofs,d,& 198 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 199 call ceedelemrestrictioncreateidentity(ceed,nelemhex,phex*phex,& 200 & nelemhex*phex*phex,d,erestrictxihex,err) 201 202 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,ndofs,1,& 203 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 204 call ceedelemrestrictioncreateidentity(ceed,nelemhex,qhex*qhex,nqptshex,& 205 & d*(d+1)/2,erestrictqdihex,err) 206 207! -- Bases 208 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 209 & bxhex,err) 210 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 211 & buhex,err) 212 213! -- QFunctions 214 call ceedqfunctioncreateinterior(ceed,1,setup,& 215 &SOURCE_DIR& 216 &//'t521-operator.h:setup'//char(0),qf_setuphex,err) 217 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 218 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 219 call ceedqfunctionaddoutput(qf_setuphex,'rho',d*(d+1)/2,ceed_eval_none,& 220 & err) 221 222 call ceedqfunctioncreateinterior(ceed,1,diff,& 223 &SOURCE_DIR& 224 &//'t521-operator.h:diff'//char(0),qf_diffhex,err) 225 call ceedqfunctionaddinput(qf_diffhex,'rho',d*(d+1)/2,ceed_eval_none,err) 226 call ceedqfunctionaddinput(qf_diffhex,'u',d,ceed_eval_grad,err) 227 call ceedqfunctionaddoutput(qf_diffhex,'v',d,ceed_eval_grad,err) 228 229! -- Operators 230! ---- Setup Hex 231 call ceedoperatorcreate(ceed,qf_setuphex,ceed_null,ceed_null,op_setuphex,& 232 & err) 233 call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,& 234 & ceed_notranspose,bxhex,ceed_vector_none,err) 235 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 236 & ceed_notranspose,bxhex,ceed_vector_active,err) 237 call ceedoperatorsetfield(op_setuphex,'rho',erestrictqdihex,& 238 & ceed_notranspose,ceed_basis_collocated,qdatahex,err) 239! ---- diff Hex 240 call ceedoperatorcreate(ceed,qf_diffhex,ceed_null,ceed_null,op_diffhex,& 241 & err) 242 call ceedoperatorsetfield(op_diffhex,'rho',erestrictqdihex,& 243 & ceed_notranspose,ceed_basis_collocated,qdatahex,err) 244 call ceedoperatorsetfield(op_diffhex,'u',erestrictuhex,& 245 & ceed_notranspose,buhex,ceed_vector_active,err) 246 call ceedoperatorsetfield(op_diffhex,'v',erestrictuhex,& 247 & ceed_notranspose,buhex,ceed_vector_active,err) 248 249! Composite Operators 250 call ceedcompositeoperatorcreate(ceed,op_setup,err) 251 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 252 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 253 254 call ceedcompositeoperatorcreate(ceed,op_diff,err) 255 call ceedcompositeoperatoraddsub(op_diff,op_difftet,err) 256 call ceedcompositeoperatoraddsub(op_diff,op_diffhex,err) 257 258! Apply Setup Operator 259 call ceedoperatorapply(op_setup,x,ceed_null,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(erestrictxitet,err) 294 call ceedelemrestrictiondestroy(erestrictuhex,err) 295 call ceedelemrestrictiondestroy(erestrictxhex,err) 296 call ceedelemrestrictiondestroy(erestrictqdihex,err) 297 call ceedelemrestrictiondestroy(erestrictxihex,err) 298 call ceedbasisdestroy(butet,err) 299 call ceedbasisdestroy(bxtet,err) 300 call ceedbasisdestroy(buhex,err) 301 call ceedbasisdestroy(bxhex,err) 302 call ceedvectordestroy(x,err) 303 call ceedvectordestroy(u,err) 304 call ceedvectordestroy(v,err) 305 call ceedvectordestroy(qdatatet,err) 306 call ceedvectordestroy(qdatahex,err) 307 call ceeddestroy(ceed,err) 308 end 309!----------------------------------------------------------------------- 310