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