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