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 erestrictxtet,erestrictutet,erestrictxitet,erestrictuitet,& 45& erestrictxhex,erestrictuhex,erestrictxihex,erestrictuihex 46 integer bxtet,butet,bxhex,buhex 47 integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 48 integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 49 integer qdatatet,qdatahex,x,u,v 50 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 51 integer row,col,offset 52 parameter(nelemtet=6) 53 parameter(ptet=6) 54 parameter(qtet=4) 55 parameter(nelemhex=6) 56 parameter(phex=3) 57 parameter(qhex=4) 58 parameter(d=2) 59 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 60 parameter(nx=3) 61 parameter(ny=3) 62 parameter(nxtet=3) 63 parameter(nytet=1) 64 parameter(nxhex=3) 65 parameter(ndofs=(nx*2+1)*(ny*2+1)) 66 parameter(nqptstet=nelemtet*qtet) 67 parameter(nqptshex=nelemhex*qhex*qhex) 68 parameter(nqpts=nqptstet+nqptshex) 69 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 70 real*8 arrx(d*ndofs) 71 integer*8 voffset,xoffset 72 73 real*8 qref(d*qtet) 74 real*8 qweight(qtet) 75 real*8 interp(ptet*qtet) 76 real*8 grad(d*ptet*qtet) 77 78 real*8 hv(ndofs) 79 real*8 total 80 81 character arg*32 82 83 external setup,mass 84 85 call getarg(1,arg) 86 87 call ceedinit(trim(arg)//char(0),ceed,err) 88 89! DoF Coordinates 90 do i=0,ny*2 91 do j=0,nx*2 92 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 93 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 94 enddo 95 enddo 96 97 call ceedvectorcreate(ceed,d*ndofs,x,err) 98 xoffset=0 99 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 100 101! Qdata Vectors 102 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 103 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 104 105! Tet Elements 106 do i=0,2 107 col=mod(i,nx) 108 row=i/nx 109 offset=col*2+row*(nx*2+1)*2 110 111 indxtet(i*2*ptet+1)=2+offset 112 indxtet(i*2*ptet+2)=9+offset 113 indxtet(i*2*ptet+3)=16+offset 114 indxtet(i*2*ptet+4)=1+offset 115 indxtet(i*2*ptet+5)=8+offset 116 indxtet(i*2*ptet+6)=0+offset 117 118 indxtet(i*2*ptet+7)=14+offset 119 indxtet(i*2*ptet+8)=7+offset 120 indxtet(i*2*ptet+9)=0+offset 121 indxtet(i*2*ptet+10)=15+offset 122 indxtet(i*2*ptet+11)=8+offset 123 indxtet(i*2*ptet+12)=16+offset 124 enddo 125 126! -- Restrictions 127 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,d,& 128 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 129 call ceedelemrestrictioncreateidentity(ceed,imode,nelemtet,ptet,& 130 & nelemtet*ptet,d,erestrictxitet,err) 131 132 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,& 133 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 134 call ceedelemrestrictioncreateidentity(ceed,imode,nelemtet,qtet,nqptstet,& 135 & 1,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',erestrictxitet,& 165 & 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 call ceedelemrestrictioncreateidentity(ceed,imode,nelemhex,phex*phex,& 196 & nelemhex*phex*phex,d,erestrictxihex,err) 197 198 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,& 199 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 200 call ceedelemrestrictioncreateidentity(ceed,imode,nelemhex,qhex*qhex,& 201 & nqptshex,1,erestrictuihex,err) 202 203! -- Bases 204 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 205 & bxhex,err) 206 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 207 & buhex,err) 208 209! -- QFunctions 210 call ceedqfunctioncreateinterior(ceed,1,setup,& 211 &SOURCE_DIR& 212 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 213 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 214 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 215 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 216 217 call ceedqfunctioncreateinterior(ceed,1,mass,& 218 &SOURCE_DIR& 219 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 220 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 221 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 222 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 223 224! -- Operators 225! ---- Setup Hex 226 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 227 & ceed_qfunction_none,op_setuphex,& 228 & err) 229 call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,& 230 & bxhex,ceed_vector_none,err) 231 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 232 & bxhex,ceed_vector_active,err) 233 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 234 & ceed_basis_collocated,qdatahex,err) 235! ---- Mass Hex 236 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 237 & ceed_qfunction_none,op_masshex,& 238 & err) 239 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 240 & ceed_basis_collocated,qdatahex,err) 241 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 242 & buhex,ceed_vector_active,err) 243 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 244 & buhex,ceed_vector_active,err) 245 246! Composite Operators 247 call ceedcompositeoperatorcreate(ceed,op_setup,err) 248 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 249 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 250 251 call ceedcompositeoperatorcreate(ceed,op_mass,err) 252 call ceedcompositeoperatoraddsub(op_mass,op_masstet,err) 253 call ceedcompositeoperatoraddsub(op_mass,op_masshex,err) 254 255! Apply Setup Operator 256 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 257 & ceed_request_immediate,err) 258 259! Apply Mass Operator 260 call ceedvectorcreate(ceed,ndofs,u,err) 261 call ceedvectorsetvalue(u,1.d0,err) 262 call ceedvectorcreate(ceed,ndofs,v,err) 263 call ceedvectorsetvalue(v,0.d0,err) 264 265 call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err) 266 267! Check Output 268 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 269 total=0. 270 do i=1,ndofs 271 total=total+hv(voffset+i) 272 enddo 273 if (abs(total-1.)>1.0d-10) then 274! LCOV_EXCL_START 275 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 276! LCOV_EXCL_STOP 277 endif 278 call ceedvectorrestorearrayread(v,hv,voffset,err) 279 280 call ceedvectorsetvalue(v,1.d0,err) 281 call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err) 282 283! Check Output 284 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 285 total=-ndofs 286 do i=1,ndofs 287 total=total+hv(voffset+i) 288 enddo 289 if (abs(total-1.)>1.0d-10) then 290! LCOV_EXCL_START 291 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 292! LCOV_EXCL_STOP 293 endif 294 call ceedvectorrestorearrayread(v,hv,voffset,err) 295 296! Cleanup 297 call ceedqfunctiondestroy(qf_setuptet,err) 298 call ceedqfunctiondestroy(qf_masstet,err) 299 call ceedoperatordestroy(op_setuptet,err) 300 call ceedoperatordestroy(op_masstet,err) 301 call ceedqfunctiondestroy(qf_setuphex,err) 302 call ceedqfunctiondestroy(qf_masshex,err) 303 call ceedoperatordestroy(op_setuphex,err) 304 call ceedoperatordestroy(op_masshex,err) 305 call ceedoperatordestroy(op_setup,err) 306 call ceedoperatordestroy(op_mass,err) 307 call ceedelemrestrictiondestroy(erestrictutet,err) 308 call ceedelemrestrictiondestroy(erestrictxtet,err) 309 call ceedelemrestrictiondestroy(erestrictuitet,err) 310 call ceedelemrestrictiondestroy(erestrictxitet,err) 311 call ceedelemrestrictiondestroy(erestrictuhex,err) 312 call ceedelemrestrictiondestroy(erestrictxhex,err) 313 call ceedelemrestrictiondestroy(erestrictuihex,err) 314 call ceedelemrestrictiondestroy(erestrictxihex,err) 315 call ceedbasisdestroy(butet,err) 316 call ceedbasisdestroy(bxtet,err) 317 call ceedbasisdestroy(buhex,err) 318 call ceedbasisdestroy(bxhex,err) 319 call ceedvectordestroy(x,err) 320 call ceedvectordestroy(u,err) 321 call ceedvectordestroy(v,err) 322 call ceedvectordestroy(qdatatet,err) 323 call ceedvectordestroy(qdatahex,err) 324 call ceeddestroy(ceed,err) 325 end 326!----------------------------------------------------------------------- 327