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