1!----------------------------------------------------------------------- 2! 3! Header with common subroutine 4! 5 include 't320-basis-f.h' 6! 7! Header with QFunctions 8! 9 include 't510-operator-f.h' 10!----------------------------------------------------------------------- 11 program test 12 implicit none 13 include 'ceed/fortran.h' 14 15 integer ceed,err,i,j,k 16 integer stridesutet(3),stridesuhex(3) 17 integer erestrictxtet,erestrictutet,erestrictuitet,& 18& erestrictxhex,erestrictuhex,erestrictuihex 19 integer bxtet,butet,bxhex,buhex 20 integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 21 integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 22 integer qdatatet,qdatahex,x,u,v 23 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 24 integer row,col,offset 25 parameter(nelemtet=6) 26 parameter(ptet=6) 27 parameter(qtet=4) 28 parameter(nelemhex=6) 29 parameter(phex=3) 30 parameter(qhex=4) 31 parameter(d=2) 32 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 33 parameter(nx=3) 34 parameter(ny=3) 35 parameter(nxtet=3) 36 parameter(nytet=1) 37 parameter(nxhex=3) 38 parameter(ndofs=(nx*2+1)*(ny*2+1)) 39 parameter(nqptstet=nelemtet*qtet) 40 parameter(nqptshex=nelemhex*qhex*qhex) 41 parameter(nqpts=nqptstet+nqptshex) 42 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 43 real*8 arrx(d*ndofs) 44 integer*8 voffset,xoffset 45 46 real*8 qref(d*qtet) 47 real*8 qweight(qtet) 48 real*8 interp(ptet*qtet) 49 real*8 grad(d*ptet*qtet) 50 51 real*8 hv(ndofs) 52 53 character arg*32 54 55 external setup,mass 56 57 call getarg(1,arg) 58 59 call ceedinit(trim(arg)//char(0),ceed,err) 60 61! DoF Coordinates 62 do i=0,ny*2 63 do j=0,nx*2 64 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 65 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 66 enddo 67 enddo 68 69 call ceedvectorcreate(ceed,d*ndofs,x,err) 70 xoffset=0 71 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 72 73! Qdata Vectors 74 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 75 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 76 77! Tet Elements 78 do i=0,2 79 col=mod(i,nx) 80 row=i/nx 81 offset=col*2+row*(nx*2+1)*2 82 83 indxtet(i*2*ptet+1)=2+offset 84 indxtet(i*2*ptet+2)=9+offset 85 indxtet(i*2*ptet+3)=16+offset 86 indxtet(i*2*ptet+4)=1+offset 87 indxtet(i*2*ptet+5)=8+offset 88 indxtet(i*2*ptet+6)=0+offset 89 90 indxtet(i*2*ptet+7)=14+offset 91 indxtet(i*2*ptet+8)=7+offset 92 indxtet(i*2*ptet+9)=0+offset 93 indxtet(i*2*ptet+10)=15+offset 94 indxtet(i*2*ptet+11)=8+offset 95 indxtet(i*2*ptet+12)=16+offset 96 enddo 97 98! -- Restrictions 99 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,& 100 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 101 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,& 102 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 103 stridesutet=[1,qtet,qtet] 104 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,& 105 & nqptstet,stridesutet,erestrictuitet,err) 106 107! -- Bases 108 call buildmats(qref,qweight,interp,grad) 109 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 110 & qweight,bxtet,err) 111 call buildmats(qref,qweight,interp,grad) 112 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 113 & qweight,butet,err) 114 115! -- QFunctions 116 call ceedqfunctioncreateinterior(ceed,1,setup,& 117 &SOURCE_DIR& 118 &//'t510-operator.h:setup'//char(0),qf_setuptet,err) 119 call ceedqfunctionaddinput(qf_setuptet,'weight',1,ceed_eval_weight,err) 120 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 121 call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err) 122 123 call ceedqfunctioncreateinterior(ceed,1,mass,& 124 &SOURCE_DIR& 125 &//'t510-operator.h:mass'//char(0),qf_masstet,err) 126 call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err) 127 call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err) 128 call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err) 129 130! -- Operators 131! ---- Setup Tet 132 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 133 & ceed_qfunction_none,op_setuptet,err) 134 call ceedoperatorsetfield(op_setuptet,'weight',& 135 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 136 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 137 & bxtet,ceed_vector_active,err) 138 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 139 ceed_basis_none,qdatatet,err) 140! ---- Mass Tet 141 call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 142 & ceed_qfunction_none,op_masstet,err) 143 call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 144 ceed_basis_none,qdatatet,err) 145 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 146 & butet,ceed_vector_active,err) 147 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 148 & butet,ceed_vector_active,err) 149 150! Hex Elements 151 do i=0,nelemhex-1 152 col=mod(i,nx) 153 row=i/nx 154 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 155 do j=0,phex-1 156 do k=0,phex-1 157 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 158 enddo 159 enddo 160 enddo 161 162! -- Restrictions 163 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,& 164 & d*ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 165 166 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,& 167 & ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 168 stridesuhex=[1,qhex*qhex,qhex*qhex] 169 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 170 & 1,nqptshex,stridesuhex,erestrictuihex,err) 171 172! -- Bases 173 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 174 & bxhex,err) 175 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 176 & buhex,err) 177 178! -- QFunctions 179 call ceedqfunctioncreateinterior(ceed,1,setup,& 180 &SOURCE_DIR& 181 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 182 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 183 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 184 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 185 186 call ceedqfunctioncreateinterior(ceed,1,mass,& 187 &SOURCE_DIR& 188 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 189 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 190 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 191 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 192 193! -- Operators 194! ---- Setup Hex 195 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 196 & ceed_qfunction_none,op_setuphex,err) 197 call ceedoperatorsetfield(op_setuphex,'_weight',& 198 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 199 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 200 & bxhex,ceed_vector_active,err) 201 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 202 ceed_basis_none,qdatahex,err) 203! ---- Mass Hex 204 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 205 & ceed_qfunction_none,op_masshex,err) 206 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 207 ceed_basis_none,qdatahex,err) 208 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 209 & buhex,ceed_vector_active,err) 210 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 211 & buhex,ceed_vector_active,err) 212 213! Composite Operators 214 call ceedcompositeoperatorcreate(ceed,op_setup,err) 215 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 216 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 217 218 call ceedcompositeoperatorcreate(ceed,op_mass,err) 219 call ceedcompositeoperatoraddsub(op_mass,op_masstet,err) 220 call ceedcompositeoperatoraddsub(op_mass,op_masshex,err) 221 222! Apply Setup Operator 223 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 224 & ceed_request_immediate,err) 225 226! Apply Mass Operator 227 call ceedvectorcreate(ceed,ndofs,u,err) 228 call ceedvectorsetvalue(u,0.d0,err) 229 call ceedvectorcreate(ceed,ndofs,v,err) 230 231 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 232 233! Check Output 234 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 235 do i=1,ndofs 236 if (abs(hv(voffset+i))>1.0d-10) then 237! LCOV_EXCL_START 238 write(*,*) '[',i,'] v ',hv(voffset+i),' != 0.0' 239! LCOV_EXCL_STOP 240 endif 241 enddo 242 call ceedvectorrestorearrayread(v,hv,voffset,err) 243 244! Cleanup 245 call ceedqfunctiondestroy(qf_setuptet,err) 246 call ceedqfunctiondestroy(qf_masstet,err) 247 call ceedoperatordestroy(op_setuptet,err) 248 call ceedoperatordestroy(op_masstet,err) 249 call ceedqfunctiondestroy(qf_setuphex,err) 250 call ceedqfunctiondestroy(qf_masshex,err) 251 call ceedoperatordestroy(op_setuphex,err) 252 call ceedoperatordestroy(op_masshex,err) 253 call ceedoperatordestroy(op_setup,err) 254 call ceedoperatordestroy(op_mass,err) 255 call ceedelemrestrictiondestroy(erestrictutet,err) 256 call ceedelemrestrictiondestroy(erestrictxtet,err) 257 call ceedelemrestrictiondestroy(erestrictuitet,err) 258 call ceedelemrestrictiondestroy(erestrictuhex,err) 259 call ceedelemrestrictiondestroy(erestrictxhex,err) 260 call ceedelemrestrictiondestroy(erestrictuihex,err) 261 call ceedbasisdestroy(butet,err) 262 call ceedbasisdestroy(bxtet,err) 263 call ceedbasisdestroy(buhex,err) 264 call ceedbasisdestroy(bxhex,err) 265 call ceedvectordestroy(x,err) 266 call ceedvectordestroy(u,err) 267 call ceedvectordestroy(v,err) 268 call ceedvectordestroy(qdatatet,err) 269 call ceedvectordestroy(qdatahex,err) 270 call ceeddestroy(ceed,err) 271 end 272!----------------------------------------------------------------------- 273