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