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