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 &//'t521-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 &//'t521-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,err) 226 call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,& 227 & ceed_notranspose,bxhex,ceed_vector_none,err) 228 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 229 & ceed_notranspose,bxhex,ceed_vector_active,err) 230 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 231 & ceed_notranspose,ceed_basis_collocated,qdatahex,err) 232! ---- Mass Hex 233 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 234 & ceed_qfunction_none,op_masshex,err) 235 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 236 & ceed_notranspose,ceed_basis_collocated,qdatahex,err) 237 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 238 & ceed_notranspose,buhex,ceed_vector_active,err) 239 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 240 & ceed_notranspose,buhex,ceed_vector_active,err) 241 242! Composite Operators 243 call ceedcompositeoperatorcreate(ceed,op_setup,err) 244 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 245 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 246 247 call ceedcompositeoperatorcreate(ceed,op_mass,err) 248 call ceedcompositeoperatoraddsub(op_mass,op_masstet,err) 249 call ceedcompositeoperatoraddsub(op_mass,op_masshex,err) 250 251! Apply Setup Operator 252 call ceedoperatorapply(op_setup,x,ceed_vector_none,& 253 & ceed_request_immediate,err) 254 255! Apply Mass Operator 256 call ceedvectorcreate(ceed,ndofs,u,err) 257 call ceedvectorsetvalue(u,1.d0,err) 258 call ceedvectorcreate(ceed,ndofs,v,err) 259 260 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 261 262! Check Output 263 call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err) 264 total=0. 265 do i=1,ndofs 266 total=total+hv(voffset+i) 267 enddo 268 if (abs(total-1.)>1.0d-10) then 269! LCOV_EXCL_START 270 write(*,*) 'Computed Area: ',total,' != True Area: 1.0' 271! LCOV_EXCL_STOP 272 endif 273 call ceedvectorrestorearrayread(v,hv,voffset,err) 274 275! Cleanup 276 call ceedqfunctiondestroy(qf_setuptet,err) 277 call ceedqfunctiondestroy(qf_masstet,err) 278 call ceedoperatordestroy(op_setuptet,err) 279 call ceedoperatordestroy(op_masstet,err) 280 call ceedqfunctiondestroy(qf_setuphex,err) 281 call ceedqfunctiondestroy(qf_masshex,err) 282 call ceedoperatordestroy(op_setuphex,err) 283 call ceedoperatordestroy(op_masshex,err) 284 call ceedoperatordestroy(op_setup,err) 285 call ceedoperatordestroy(op_mass,err) 286 call ceedelemrestrictiondestroy(erestrictutet,err) 287 call ceedelemrestrictiondestroy(erestrictxtet,err) 288 call ceedelemrestrictiondestroy(erestrictuitet,err) 289 call ceedelemrestrictiondestroy(erestrictxitet,err) 290 call ceedelemrestrictiondestroy(erestrictuhex,err) 291 call ceedelemrestrictiondestroy(erestrictxhex,err) 292 call ceedelemrestrictiondestroy(erestrictuihex,err) 293 call ceedelemrestrictiondestroy(erestrictxihex,err) 294 call ceedbasisdestroy(butet,err) 295 call ceedbasisdestroy(bxtet,err) 296 call ceedbasisdestroy(buhex,err) 297 call ceedbasisdestroy(bxhex,err) 298 call ceedvectordestroy(x,err) 299 call ceedvectordestroy(u,err) 300 call ceedvectordestroy(v,err) 301 call ceedvectordestroy(qdatatet,err) 302 call ceedvectordestroy(qdatahex,err) 303 call ceeddestroy(ceed,err) 304 end 305!----------------------------------------------------------------------- 306