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