152d6035fSJeremy L Thompson!----------------------------------------------------------------------- 252d6035fSJeremy L Thompson! 352d6035fSJeremy L Thompson! Header with common subroutine 452d6035fSJeremy L Thompson! 552bfb9bbSJeremy L Thompson include 't320-basis-f.h' 6752c3701SJeremy L Thompson! 7752c3701SJeremy L Thompson! Header with QFunctions 8752c3701SJeremy L Thompson! 9752c3701SJeremy L Thompson include 't510-operator-f.h' 1052d6035fSJeremy L Thompson!----------------------------------------------------------------------- 1152d6035fSJeremy L Thompson program test 121f9a83abSJed Brown implicit none 13ec3da8bcSJed Brown include 'ceed/fortran.h' 1452d6035fSJeremy L Thompson 1552d6035fSJeremy L Thompson integer ceed,err,i,j,k 1615910d16Sjeremylt integer stridesutet(3),stridesuhex(3) 1715910d16Sjeremylt integer erestrictxtet,erestrictutet,erestrictuitet,& 1815910d16Sjeremylt& erestrictxhex,erestrictuhex,erestrictuihex 1952d6035fSJeremy L Thompson integer bxtet,butet,bxhex,buhex 2052d6035fSJeremy L Thompson integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 2152d6035fSJeremy L Thompson integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 2252d6035fSJeremy L Thompson integer qdatatet,qdatahex,x,u,v 2352d6035fSJeremy L Thompson integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 2452d6035fSJeremy L Thompson integer row,col,offset 2552d6035fSJeremy L Thompson parameter(nelemtet=6) 2652d6035fSJeremy L Thompson parameter(ptet=6) 2752d6035fSJeremy L Thompson parameter(qtet=4) 2852d6035fSJeremy L Thompson parameter(nelemhex=6) 2952d6035fSJeremy L Thompson parameter(phex=3) 3052d6035fSJeremy L Thompson parameter(qhex=4) 3152d6035fSJeremy L Thompson parameter(d=2) 3252d6035fSJeremy L Thompson integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 3352d6035fSJeremy L Thompson parameter(nx=3) 3452d6035fSJeremy L Thompson parameter(ny=3) 3552d6035fSJeremy L Thompson parameter(nxtet=3) 3652d6035fSJeremy L Thompson parameter(nytet=1) 3752d6035fSJeremy L Thompson parameter(nxhex=3) 3852d6035fSJeremy L Thompson parameter(ndofs=(nx*2+1)*(ny*2+1)) 3952d6035fSJeremy L Thompson parameter(nqptstet=nelemtet*qtet) 4052d6035fSJeremy L Thompson parameter(nqptshex=nelemhex*qhex*qhex) 4152d6035fSJeremy L Thompson parameter(nqpts=nqptstet+nqptshex) 4252d6035fSJeremy L Thompson integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 4352d6035fSJeremy L Thompson real*8 arrx(d*ndofs) 44c8b9fe72Sjeremylt integer*8 voffset,xoffset 4552d6035fSJeremy L Thompson 4652d6035fSJeremy L Thompson real*8 qref(d*qtet) 4752d6035fSJeremy L Thompson real*8 qweight(qtet) 4852d6035fSJeremy L Thompson real*8 interp(ptet*qtet) 4952d6035fSJeremy L Thompson real*8 grad(d*ptet*qtet) 5052d6035fSJeremy L Thompson 5152d6035fSJeremy L Thompson real*8 hv(ndofs) 5252d6035fSJeremy L Thompson 5352d6035fSJeremy L Thompson character arg*32 5452d6035fSJeremy L Thompson 5552d6035fSJeremy L Thompson external setup,mass 5652d6035fSJeremy L Thompson 5752d6035fSJeremy L Thompson call getarg(1,arg) 5852d6035fSJeremy L Thompson 5952d6035fSJeremy L Thompson call ceedinit(trim(arg)//char(0),ceed,err) 6052d6035fSJeremy L Thompson 6152d6035fSJeremy L Thompson! DoF Coordinates 6252d6035fSJeremy L Thompson do i=0,ny*2 6352d6035fSJeremy L Thompson do j=0,nx*2 6452d6035fSJeremy L Thompson arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 6552d6035fSJeremy L Thompson arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 6652d6035fSJeremy L Thompson enddo 6752d6035fSJeremy L Thompson enddo 6852d6035fSJeremy L Thompson 6952d6035fSJeremy L Thompson call ceedvectorcreate(ceed,d*ndofs,x,err) 70c8b9fe72Sjeremylt xoffset=0 71c8b9fe72Sjeremylt call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 7252d6035fSJeremy L Thompson 7352d6035fSJeremy L Thompson! Qdata Vectors 7452d6035fSJeremy L Thompson call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 7552d6035fSJeremy L Thompson call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 7652d6035fSJeremy L Thompson 7752d6035fSJeremy L Thompson! Tet Elements 7852d6035fSJeremy L Thompson do i=0,2 7952d6035fSJeremy L Thompson col=mod(i,nx) 8052d6035fSJeremy L Thompson row=i/nx 8152d6035fSJeremy L Thompson offset=col*2+row*(nx*2+1)*2 8252d6035fSJeremy L Thompson 8352d6035fSJeremy L Thompson indxtet(i*2*ptet+1)=2+offset 8452d6035fSJeremy L Thompson indxtet(i*2*ptet+2)=9+offset 8552d6035fSJeremy L Thompson indxtet(i*2*ptet+3)=16+offset 8652d6035fSJeremy L Thompson indxtet(i*2*ptet+4)=1+offset 8752d6035fSJeremy L Thompson indxtet(i*2*ptet+5)=8+offset 8852d6035fSJeremy L Thompson indxtet(i*2*ptet+6)=0+offset 8952d6035fSJeremy L Thompson 9052d6035fSJeremy L Thompson indxtet(i*2*ptet+7)=14+offset 9152d6035fSJeremy L Thompson indxtet(i*2*ptet+8)=7+offset 9252d6035fSJeremy L Thompson indxtet(i*2*ptet+9)=0+offset 9352d6035fSJeremy L Thompson indxtet(i*2*ptet+10)=15+offset 9452d6035fSJeremy L Thompson indxtet(i*2*ptet+11)=8+offset 9552d6035fSJeremy L Thompson indxtet(i*2*ptet+12)=16+offset 9652d6035fSJeremy L Thompson enddo 9752d6035fSJeremy L Thompson 9852d6035fSJeremy L Thompson! -- Restrictions 99d979a051Sjeremylt call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,& 100a8d32208Sjeremylt & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 101d979a051Sjeremylt call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,& 102a8d32208Sjeremylt & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 1037509a596Sjeremylt stridesutet=[1,qtet,qtet] 104d979a051Sjeremylt call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,& 105d979a051Sjeremylt & nqptstet,stridesutet,erestrictuitet,err) 10652d6035fSJeremy L Thompson 10752d6035fSJeremy L Thompson! -- Bases 10852d6035fSJeremy L Thompson call buildmats(qref,qweight,interp,grad) 10952d6035fSJeremy L Thompson call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 11052d6035fSJeremy L Thompson & qweight,bxtet,err) 11152d6035fSJeremy L Thompson call buildmats(qref,qweight,interp,grad) 11252d6035fSJeremy L Thompson call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 11352d6035fSJeremy L Thompson & qweight,butet,err) 11452d6035fSJeremy L Thompson 11552d6035fSJeremy L Thompson! -- QFunctions 11652d6035fSJeremy L Thompson call ceedqfunctioncreateinterior(ceed,1,setup,& 1172d50dd3dSjeremylt &SOURCE_DIR& 118872c4ebbSjeremylt &//'t510-operator.h:setup'//char(0),qf_setuptet,err) 119a61c78d6SJeremy L Thompson call ceedqfunctionaddinput(qf_setuptet,'weight',1,ceed_eval_weight,err) 1204d537eeaSYohann call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 12152d6035fSJeremy L Thompson call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err) 12252d6035fSJeremy L Thompson 12352d6035fSJeremy L Thompson call ceedqfunctioncreateinterior(ceed,1,mass,& 1242d50dd3dSjeremylt &SOURCE_DIR& 125872c4ebbSjeremylt &//'t510-operator.h:mass'//char(0),qf_masstet,err) 12652d6035fSJeremy L Thompson call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err) 12752d6035fSJeremy L Thompson call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err) 12852d6035fSJeremy L Thompson call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err) 12952d6035fSJeremy L Thompson 13052d6035fSJeremy L Thompson! -- Operators 13152d6035fSJeremy L Thompson! ---- Setup Tet 132442e7f0bSjeremylt call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 133442e7f0bSjeremylt & ceed_qfunction_none,op_setuptet,err) 134a61c78d6SJeremy L Thompson call ceedoperatorsetfield(op_setuptet,'weight',& 13515910d16Sjeremylt & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 13652d6035fSJeremy L Thompson call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 137a8d32208Sjeremylt & bxtet,ceed_vector_active,err) 13852d6035fSJeremy L Thompson call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 139356036faSJeremy L Thompson ceed_basis_none,qdatatet,err) 14052d6035fSJeremy L Thompson! ---- Mass Tet 141442e7f0bSjeremylt call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 142442e7f0bSjeremylt & ceed_qfunction_none,op_masstet,err) 14352d6035fSJeremy L Thompson call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 144356036faSJeremy L Thompson ceed_basis_none,qdatatet,err) 14552d6035fSJeremy L Thompson call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 146a8d32208Sjeremylt & butet,ceed_vector_active,err) 14752d6035fSJeremy L Thompson call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 148a8d32208Sjeremylt & butet,ceed_vector_active,err) 14952d6035fSJeremy L Thompson 15052d6035fSJeremy L Thompson! Hex Elements 15152d6035fSJeremy L Thompson do i=0,nelemhex-1 15252d6035fSJeremy L Thompson col=mod(i,nx) 15352d6035fSJeremy L Thompson row=i/nx 15452d6035fSJeremy L Thompson offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 15552d6035fSJeremy L Thompson do j=0,phex-1 15652d6035fSJeremy L Thompson do k=0,phex-1 15752d6035fSJeremy L Thompson indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 15852d6035fSJeremy L Thompson enddo 15952d6035fSJeremy L Thompson enddo 16052d6035fSJeremy L Thompson enddo 16152d6035fSJeremy L Thompson 16252d6035fSJeremy L Thompson! -- Restrictions 163d979a051Sjeremylt call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,& 164d979a051Sjeremylt & d*ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 16552d6035fSJeremy L Thompson 166d979a051Sjeremylt call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,& 167d979a051Sjeremylt & ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 1687509a596Sjeremylt stridesuhex=[1,qhex*qhex,qhex*qhex] 1697509a596Sjeremylt call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 170d979a051Sjeremylt & 1,nqptshex,stridesuhex,erestrictuihex,err) 17152d6035fSJeremy L Thompson 17252d6035fSJeremy L Thompson! -- Bases 17352d6035fSJeremy L Thompson call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 17452d6035fSJeremy L Thompson & bxhex,err) 17552d6035fSJeremy L Thompson call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 17652d6035fSJeremy L Thompson & buhex,err) 17752d6035fSJeremy L Thompson 17852d6035fSJeremy L Thompson! -- QFunctions 17952d6035fSJeremy L Thompson call ceedqfunctioncreateinterior(ceed,1,setup,& 1802d50dd3dSjeremylt &SOURCE_DIR& 181a05f9790Sjeremylt &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 18252d6035fSJeremy L Thompson call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 1834d537eeaSYohann call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 18452d6035fSJeremy L Thompson call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 18552d6035fSJeremy L Thompson 18652d6035fSJeremy L Thompson call ceedqfunctioncreateinterior(ceed,1,mass,& 1872d50dd3dSjeremylt &SOURCE_DIR& 188a05f9790Sjeremylt &//'t510-operator.h:mass'//char(0),qf_masshex,err) 18952d6035fSJeremy L Thompson call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 19052d6035fSJeremy L Thompson call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 19152d6035fSJeremy L Thompson call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 19252d6035fSJeremy L Thompson 19352d6035fSJeremy L Thompson! -- Operators 19452d6035fSJeremy L Thompson! ---- Setup Hex 195442e7f0bSjeremylt call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 196442e7f0bSjeremylt & ceed_qfunction_none,op_setuphex,err) 19715910d16Sjeremylt call ceedoperatorsetfield(op_setuphex,'_weight',& 19815910d16Sjeremylt & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 19952d6035fSJeremy L Thompson call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 200a8d32208Sjeremylt & bxhex,ceed_vector_active,err) 20152d6035fSJeremy L Thompson call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 202356036faSJeremy L Thompson ceed_basis_none,qdatahex,err) 20352d6035fSJeremy L Thompson! ---- Mass Hex 204442e7f0bSjeremylt call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 205442e7f0bSjeremylt & ceed_qfunction_none,op_masshex,err) 20652d6035fSJeremy L Thompson call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 207356036faSJeremy L Thompson ceed_basis_none,qdatahex,err) 20852d6035fSJeremy L Thompson call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 209a8d32208Sjeremylt & buhex,ceed_vector_active,err) 21052d6035fSJeremy L Thompson call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 211a8d32208Sjeremylt & buhex,ceed_vector_active,err) 21252d6035fSJeremy L Thompson 21352d6035fSJeremy L Thompson! Composite Operators 214*ed094490SJeremy L Thompson call ceedoperatorcreatecomposite(ceed,op_setup,err) 215*ed094490SJeremy L Thompson call ceedoperatorcompositeaddsub(op_setup,op_setuptet,err) 216*ed094490SJeremy L Thompson call ceedoperatorcompositeaddsub(op_setup,op_setuphex,err) 21752d6035fSJeremy L Thompson 218*ed094490SJeremy L Thompson call ceedoperatorcreatecomposite(ceed,op_mass,err) 219*ed094490SJeremy L Thompson call ceedoperatorcompositeaddsub(op_mass,op_masstet,err) 220*ed094490SJeremy L Thompson call ceedoperatorcompositeaddsub(op_mass,op_masshex,err) 22152d6035fSJeremy L Thompson 22252d6035fSJeremy L Thompson! Apply Setup Operator 223e97ff134Sjeremylt call ceedoperatorapply(op_setup,x,ceed_vector_none,& 224e97ff134Sjeremylt & ceed_request_immediate,err) 22552d6035fSJeremy L Thompson 22652d6035fSJeremy L Thompson! Apply Mass Operator 22752d6035fSJeremy L Thompson call ceedvectorcreate(ceed,ndofs,u,err) 22852d6035fSJeremy L Thompson call ceedvectorsetvalue(u,0.d0,err) 22952d6035fSJeremy L Thompson call ceedvectorcreate(ceed,ndofs,v,err) 23052d6035fSJeremy L Thompson 23152d6035fSJeremy L Thompson call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 23252d6035fSJeremy L Thompson 23352d6035fSJeremy L Thompson! Check Output 23452d6035fSJeremy L Thompson call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 23552d6035fSJeremy L Thompson do i=1,ndofs 23652d6035fSJeremy L Thompson if (abs(hv(voffset+i))>1.0d-10) then 237a2546046Sjeremylt! LCOV_EXCL_START 23852d6035fSJeremy L Thompson write(*,*) '[',i,'] v ',hv(voffset+i),' != 0.0' 239de996c55Sjeremylt! LCOV_EXCL_STOP 24052d6035fSJeremy L Thompson endif 24152d6035fSJeremy L Thompson enddo 24252d6035fSJeremy L Thompson call ceedvectorrestorearrayread(v,hv,voffset,err) 24352d6035fSJeremy L Thompson 24452d6035fSJeremy L Thompson! Cleanup 24552d6035fSJeremy L Thompson call ceedqfunctiondestroy(qf_setuptet,err) 24652d6035fSJeremy L Thompson call ceedqfunctiondestroy(qf_masstet,err) 24752d6035fSJeremy L Thompson call ceedoperatordestroy(op_setuptet,err) 24852d6035fSJeremy L Thompson call ceedoperatordestroy(op_masstet,err) 24952d6035fSJeremy L Thompson call ceedqfunctiondestroy(qf_setuphex,err) 25052d6035fSJeremy L Thompson call ceedqfunctiondestroy(qf_masshex,err) 25152d6035fSJeremy L Thompson call ceedoperatordestroy(op_setuphex,err) 25252d6035fSJeremy L Thompson call ceedoperatordestroy(op_masshex,err) 25352d6035fSJeremy L Thompson call ceedoperatordestroy(op_setup,err) 25452d6035fSJeremy L Thompson call ceedoperatordestroy(op_mass,err) 25552d6035fSJeremy L Thompson call ceedelemrestrictiondestroy(erestrictutet,err) 25652d6035fSJeremy L Thompson call ceedelemrestrictiondestroy(erestrictxtet,err) 25752d6035fSJeremy L Thompson call ceedelemrestrictiondestroy(erestrictuitet,err) 25852d6035fSJeremy L Thompson call ceedelemrestrictiondestroy(erestrictuhex,err) 25952d6035fSJeremy L Thompson call ceedelemrestrictiondestroy(erestrictxhex,err) 26052d6035fSJeremy L Thompson call ceedelemrestrictiondestroy(erestrictuihex,err) 26152d6035fSJeremy L Thompson call ceedbasisdestroy(butet,err) 26252d6035fSJeremy L Thompson call ceedbasisdestroy(bxtet,err) 26352d6035fSJeremy L Thompson call ceedbasisdestroy(buhex,err) 26452d6035fSJeremy L Thompson call ceedbasisdestroy(bxhex,err) 26552d6035fSJeremy L Thompson call ceedvectordestroy(x,err) 26652d6035fSJeremy L Thompson call ceedvectordestroy(u,err) 26752d6035fSJeremy L Thompson call ceedvectordestroy(v,err) 26852d6035fSJeremy L Thompson call ceedvectordestroy(qdatatet,err) 26952d6035fSJeremy L Thompson call ceedvectordestroy(qdatahex,err) 27052d6035fSJeremy L Thompson call ceeddestroy(ceed,err) 27152d6035fSJeremy L Thompson end 27252d6035fSJeremy L Thompson!----------------------------------------------------------------------- 273