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