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 implicit none 39 include 'ceedf.h' 40 41 integer ceed,err,i,j,k 42 integer stridesutet(3),stridesuhex(3) 43 integer erestrictxtet,erestrictutet,erestrictuitet,& 44& erestrictxhex,erestrictuhex,erestrictuihex 45 integer bxtet,butet,bxhex,buhex 46 integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 47 integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 48 integer qdatatet,qdatahex,x,u,v 49 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 50 integer row,col,offset 51 parameter(nelemtet=6) 52 parameter(ptet=6) 53 parameter(qtet=4) 54 parameter(nelemhex=6) 55 parameter(phex=3) 56 parameter(qhex=4) 57 parameter(d=2) 58 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 59 parameter(nx=3) 60 parameter(ny=3) 61 parameter(nxtet=3) 62 parameter(nytet=1) 63 parameter(nxhex=3) 64 parameter(ndofs=(nx*2+1)*(ny*2+1)) 65 parameter(nqptstet=nelemtet*qtet) 66 parameter(nqptshex=nelemhex*qhex*qhex) 67 parameter(nqpts=nqptstet+nqptshex) 68 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 69 real*8 arrx(d*ndofs) 70 integer*8 voffset,xoffset 71 72 real*8 qref(d*qtet) 73 real*8 qweight(qtet) 74 real*8 interp(ptet*qtet) 75 real*8 grad(d*ptet*qtet) 76 77 real*8 hv(ndofs) 78 79 character arg*32 80 81 external setup,mass 82 83 call getarg(1,arg) 84 85 call ceedinit(trim(arg)//char(0),ceed,err) 86 87! DoF Coordinates 88 do i=0,ny*2 89 do j=0,nx*2 90 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 91 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 92 enddo 93 enddo 94 95 call ceedvectorcreate(ceed,d*ndofs,x,err) 96 xoffset=0 97 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 98 99! Qdata Vectors 100 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 101 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 102 103! Tet Elements 104 do i=0,2 105 col=mod(i,nx) 106 row=i/nx 107 offset=col*2+row*(nx*2+1)*2 108 109 indxtet(i*2*ptet+1)=2+offset 110 indxtet(i*2*ptet+2)=9+offset 111 indxtet(i*2*ptet+3)=16+offset 112 indxtet(i*2*ptet+4)=1+offset 113 indxtet(i*2*ptet+5)=8+offset 114 indxtet(i*2*ptet+6)=0+offset 115 116 indxtet(i*2*ptet+7)=14+offset 117 indxtet(i*2*ptet+8)=7+offset 118 indxtet(i*2*ptet+9)=0+offset 119 indxtet(i*2*ptet+10)=15+offset 120 indxtet(i*2*ptet+11)=8+offset 121 indxtet(i*2*ptet+12)=16+offset 122 enddo 123 124! -- Restrictions 125 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,& 126 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 127 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,& 128 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 129 stridesutet=[1,qtet,qtet] 130 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,& 131 & nqptstet,stridesutet,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 &SOURCE_DIR& 144 &//'t510-operator.h:setup'//char(0),qf_setuptet,err) 145 call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err) 146 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 147 call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err) 148 149 call ceedqfunctioncreateinterior(ceed,1,mass,& 150 &SOURCE_DIR& 151 &//'t510-operator.h: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_qfunction_none,& 159 & ceed_qfunction_none,op_setuptet,err) 160 call ceedoperatorsetfield(op_setuptet,'_weight',& 161 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 162 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 163 & bxtet,ceed_vector_active,err) 164 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 165 & ceed_basis_collocated,qdatatet,err) 166! ---- Mass Tet 167 call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 168 & ceed_qfunction_none,op_masstet,err) 169 call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 170 & ceed_basis_collocated,qdatatet,err) 171 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 172 & butet,ceed_vector_active,err) 173 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 174 & 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,d,ndofs,& 190 & d*ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 191 192 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,& 193 & ndofs,ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 194 stridesuhex=[1,qhex*qhex,qhex*qhex] 195 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 196 & 1,nqptshex,stridesuhex,erestrictuihex,err) 197 198! -- Bases 199 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 200 & bxhex,err) 201 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 202 & buhex,err) 203 204! -- QFunctions 205 call ceedqfunctioncreateinterior(ceed,1,setup,& 206 &SOURCE_DIR& 207 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 208 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 209 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 210 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 211 212 call ceedqfunctioncreateinterior(ceed,1,mass,& 213 &SOURCE_DIR& 214 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 215 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 216 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 217 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 218 219! -- Operators 220! ---- Setup Hex 221 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 222 & ceed_qfunction_none,op_setuphex,err) 223 call ceedoperatorsetfield(op_setuphex,'_weight',& 224 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 225 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 226 & bxhex,ceed_vector_active,err) 227 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 228 & ceed_basis_collocated,qdatahex,err) 229! ---- Mass Hex 230 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 231 & ceed_qfunction_none,op_masshex,err) 232 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 233 & ceed_basis_collocated,qdatahex,err) 234 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 235 & buhex,ceed_vector_active,err) 236 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 237 & buhex,ceed_vector_active,err) 238 239! Composite Operators 240 call ceedcompositeoperatorcreate(ceed,op_setup,err) 241 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 242 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 243 244 call ceedcompositeoperatorcreate(ceed,op_mass,err) 245 call ceedcompositeoperatoraddsub(op_mass,op_masstet,err) 246 call ceedcompositeoperatoraddsub(op_mass,op_masshex,err) 247 248! Apply Setup Operator 249 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 250 & 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! LCOV_EXCL_START 264 write(*,*) '[',i,'] v ',hv(voffset+i),' != 0.0' 265! LCOV_EXCL_STOP 266 endif 267 enddo 268 call ceedvectorrestorearrayread(v,hv,voffset,err) 269 270! Cleanup 271 call ceedqfunctiondestroy(qf_setuptet,err) 272 call ceedqfunctiondestroy(qf_masstet,err) 273 call ceedoperatordestroy(op_setuptet,err) 274 call ceedoperatordestroy(op_masstet,err) 275 call ceedqfunctiondestroy(qf_setuphex,err) 276 call ceedqfunctiondestroy(qf_masshex,err) 277 call ceedoperatordestroy(op_setuphex,err) 278 call ceedoperatordestroy(op_masshex,err) 279 call ceedoperatordestroy(op_setup,err) 280 call ceedoperatordestroy(op_mass,err) 281 call ceedelemrestrictiondestroy(erestrictutet,err) 282 call ceedelemrestrictiondestroy(erestrictxtet,err) 283 call ceedelemrestrictiondestroy(erestrictuitet,err) 284 call ceedelemrestrictiondestroy(erestrictuhex,err) 285 call ceedelemrestrictiondestroy(erestrictxhex,err) 286 call ceedelemrestrictiondestroy(erestrictuihex,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