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 stridesutet(3),stridesuhex(3) 45 integer erestrictxtet,erestrictutet,erestrictuitet,& 46& erestrictxhex,erestrictuhex,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 131 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,& 132 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 133 stridesutet=[1,qtet,qtet] 134 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,nqptstet,& 135 & 1,stridesutet,erestrictuitet,err) 136 137! -- Bases 138 call buildmats(qref,qweight,interp,grad) 139 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 140 & qweight,bxtet,err) 141 call buildmats(qref,qweight,interp,grad) 142 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 143 & qweight,butet,err) 144 145! -- QFunctions 146 call ceedqfunctioncreateinterior(ceed,1,setup,& 147 &SOURCE_DIR& 148 &//'t510-operator.h:setup'//char(0),qf_setuptet,err) 149 call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err) 150 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 151 call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err) 152 153 call ceedqfunctioncreateinterior(ceed,1,mass,& 154 &SOURCE_DIR& 155 &//'t510-operator.h:mass'//char(0),qf_masstet,err) 156 call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err) 157 call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err) 158 call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err) 159 160! -- Operators 161! ---- Setup Tet 162 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 163 & ceed_qfunction_none,op_setuptet,err) 164 call ceedoperatorsetfield(op_setuptet,'_weight',& 165 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 166 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 167 & bxtet,ceed_vector_active,err) 168 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 169 & ceed_basis_collocated,qdatatet,err) 170! ---- Mass Tet 171 call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 172 & ceed_qfunction_none,op_masstet,err) 173 call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 174 & ceed_basis_collocated,qdatatet,err) 175 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 176 & butet,ceed_vector_active,err) 177 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 178 & butet,ceed_vector_active,err) 179 180! Hex Elements 181 do i=0,nelemhex-1 182 col=mod(i,nx) 183 row=i/nx 184 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 185 do j=0,phex-1 186 do k=0,phex-1 187 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 188 enddo 189 enddo 190 enddo 191 192! -- Restrictions 193 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,d,& 194 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 195 196 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,& 197 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 198 stridesuhex=[1,qhex*qhex,qhex*qhex] 199 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 200 & nqptshex,1,stridesuhex,erestrictuihex,err) 201 202! -- Bases 203 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 204 & bxhex,err) 205 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 206 & buhex,err) 207 208! -- QFunctions 209 call ceedqfunctioncreateinterior(ceed,1,setup,& 210 &SOURCE_DIR& 211 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 212 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 213 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 214 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 215 216 call ceedqfunctioncreateinterior(ceed,1,mass,& 217 &SOURCE_DIR& 218 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 219 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 220 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 221 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 222 223! -- Operators 224! ---- Setup Hex 225 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 226 & ceed_qfunction_none,op_setuphex,err) 227 call ceedoperatorsetfield(op_setuphex,'_weight',& 228 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 229 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 230 & bxhex,ceed_vector_active,err) 231 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 232 & ceed_basis_collocated,qdatahex,err) 233! ---- Mass Hex 234 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 235 & ceed_qfunction_none,op_masshex,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 261 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 262 263! Check Output 264 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 265 total=0. 266 do i=1,ndofs 267 total=total+hv(voffset+i) 268 enddo 269 if (abs(total-1.)>1.0d-10) then 270! LCOV_EXCL_START 271 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 272! LCOV_EXCL_STOP 273 endif 274 call ceedvectorrestorearrayread(v,hv,voffset,err) 275 276! Cleanup 277 call ceedqfunctiondestroy(qf_setuptet,err) 278 call ceedqfunctiondestroy(qf_masstet,err) 279 call ceedoperatordestroy(op_setuptet,err) 280 call ceedoperatordestroy(op_masstet,err) 281 call ceedqfunctiondestroy(qf_setuphex,err) 282 call ceedqfunctiondestroy(qf_masshex,err) 283 call ceedoperatordestroy(op_setuphex,err) 284 call ceedoperatordestroy(op_masshex,err) 285 call ceedoperatordestroy(op_setup,err) 286 call ceedoperatordestroy(op_mass,err) 287 call ceedelemrestrictiondestroy(erestrictutet,err) 288 call ceedelemrestrictiondestroy(erestrictxtet,err) 289 call ceedelemrestrictiondestroy(erestrictuitet,err) 290 call ceedelemrestrictiondestroy(erestrictuhex,err) 291 call ceedelemrestrictiondestroy(erestrictxhex,err) 292 call ceedelemrestrictiondestroy(erestrictuihex,err) 293 call ceedbasisdestroy(butet,err) 294 call ceedbasisdestroy(bxtet,err) 295 call ceedbasisdestroy(buhex,err) 296 call ceedbasisdestroy(bxhex,err) 297 call ceedvectordestroy(x,err) 298 call ceedvectordestroy(u,err) 299 call ceedvectordestroy(v,err) 300 call ceedvectordestroy(qdatatet,err) 301 call ceedvectordestroy(qdatahex,err) 302 call ceeddestroy(ceed,err) 303 end 304!----------------------------------------------------------------------- 305