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 real*8 total 79 80 character arg*32 81 82 external setup,mass 83 84 call getarg(1,arg) 85 86 call ceedinit(trim(arg)//char(0),ceed,err) 87 88! DoF Coordinates 89 do i=0,ny*2 90 do j=0,nx*2 91 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 92 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 93 enddo 94 enddo 95 96 call ceedvectorcreate(ceed,d*ndofs,x,err) 97 xoffset=0 98 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 99 100! Qdata Vectors 101 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 102 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 103 104! Tet Elements 105 do i=0,2 106 col=mod(i,nx) 107 row=i/nx 108 offset=col*2+row*(nx*2+1)*2 109 110 indxtet(i*2*ptet+1)=2+offset 111 indxtet(i*2*ptet+2)=9+offset 112 indxtet(i*2*ptet+3)=16+offset 113 indxtet(i*2*ptet+4)=1+offset 114 indxtet(i*2*ptet+5)=8+offset 115 indxtet(i*2*ptet+6)=0+offset 116 117 indxtet(i*2*ptet+7)=14+offset 118 indxtet(i*2*ptet+8)=7+offset 119 indxtet(i*2*ptet+9)=0+offset 120 indxtet(i*2*ptet+10)=15+offset 121 indxtet(i*2*ptet+11)=8+offset 122 indxtet(i*2*ptet+12)=16+offset 123 enddo 124 125! -- Restrictions 126 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,& 127 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 128 129 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,& 130 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 131 stridesutet=[1,qtet,qtet] 132 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,nqptstet,& 133 & stridesutet,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',& 163 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 164 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 165 & bxtet,ceed_vector_active,err) 166 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 167 & 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_basis_collocated,qdatatet,err) 173 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 174 & butet,ceed_vector_active,err) 175 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 176 & 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,d,ndofs,d*ndofs,& 192 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 193 194 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,ndofs,& 195 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 196 stridesuhex=[1,qhex*qhex,qhex*qhex] 197 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 198 & 1,nqptshex,stridesuhex,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,& 225 & err) 226 call ceedoperatorsetfield(op_setuphex,'_weight',& 227 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 228 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 229 & bxhex,ceed_vector_active,err) 230 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 231 & ceed_basis_collocated,qdatahex,err) 232! ---- Mass Hex 233 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 234 & ceed_qfunction_none,op_masshex,& 235 & err) 236 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 237 & ceed_basis_collocated,qdatahex,err) 238 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 239 & buhex,ceed_vector_active,err) 240 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 241 & buhex,ceed_vector_active,err) 242 243! Composite Operators 244 call ceedcompositeoperatorcreate(ceed,op_setup,err) 245 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 246 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 247 248 call ceedcompositeoperatorcreate(ceed,op_mass,err) 249 call ceedcompositeoperatoraddsub(op_mass,op_masstet,err) 250 call ceedcompositeoperatoraddsub(op_mass,op_masshex,err) 251 252! Apply Setup Operator 253 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 254 & ceed_request_immediate,err) 255 256! Apply Mass Operator 257 call ceedvectorcreate(ceed,ndofs,u,err) 258 call ceedvectorsetvalue(u,1.d0,err) 259 call ceedvectorcreate(ceed,ndofs,v,err) 260 call ceedvectorsetvalue(v,0.d0,err) 261 262 call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err) 263 264! Check Output 265 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 266 total=0. 267 do i=1,ndofs 268 total=total+hv(voffset+i) 269 enddo 270 if (abs(total-1.)>1.0d-10) then 271! LCOV_EXCL_START 272 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 273! LCOV_EXCL_STOP 274 endif 275 call ceedvectorrestorearrayread(v,hv,voffset,err) 276 277 call ceedvectorsetvalue(v,1.d0,err) 278 call ceedoperatorapplyadd(op_mass,u,v,ceed_request_immediate,err) 279 280! Check Output 281 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 282 total=-ndofs 283 do i=1,ndofs 284 total=total+hv(voffset+i) 285 enddo 286 if (abs(total-1.)>1.0d-10) then 287! LCOV_EXCL_START 288 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 289! LCOV_EXCL_STOP 290 endif 291 call ceedvectorrestorearrayread(v,hv,voffset,err) 292 293! Cleanup 294 call ceedqfunctiondestroy(qf_setuptet,err) 295 call ceedqfunctiondestroy(qf_masstet,err) 296 call ceedoperatordestroy(op_setuptet,err) 297 call ceedoperatordestroy(op_masstet,err) 298 call ceedqfunctiondestroy(qf_setuphex,err) 299 call ceedqfunctiondestroy(qf_masshex,err) 300 call ceedoperatordestroy(op_setuphex,err) 301 call ceedoperatordestroy(op_masshex,err) 302 call ceedoperatordestroy(op_setup,err) 303 call ceedoperatordestroy(op_mass,err) 304 call ceedelemrestrictiondestroy(erestrictutet,err) 305 call ceedelemrestrictiondestroy(erestrictxtet,err) 306 call ceedelemrestrictiondestroy(erestrictuitet,err) 307 call ceedelemrestrictiondestroy(erestrictuhex,err) 308 call ceedelemrestrictiondestroy(erestrictxhex,err) 309 call ceedelemrestrictiondestroy(erestrictuihex,err) 310 call ceedbasisdestroy(butet,err) 311 call ceedbasisdestroy(bxtet,err) 312 call ceedbasisdestroy(buhex,err) 313 call ceedbasisdestroy(bxhex,err) 314 call ceedvectordestroy(x,err) 315 call ceedvectordestroy(u,err) 316 call ceedvectordestroy(v,err) 317 call ceedvectordestroy(qdatatet,err) 318 call ceedvectordestroy(qdatahex,err) 319 call ceeddestroy(ceed,err) 320 end 321!----------------------------------------------------------------------- 322