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