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