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 erestrictxtet,erestrictutet,erestrictxitet,erestrictuitet,& 45& erestrictxhex,erestrictuhex,erestrictxihex,erestrictuihex 46 integer bxtet,butet,bxhex,buhex 47 integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 48 integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 49 integer qdatatet,qdatahex,x,u,v 50 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 51 integer row,col,offset 52 parameter(nelemtet=6) 53 parameter(ptet=6) 54 parameter(qtet=4) 55 parameter(nelemhex=6) 56 parameter(phex=3) 57 parameter(qhex=4) 58 parameter(d=2) 59 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 60 parameter(nx=3) 61 parameter(ny=3) 62 parameter(nxtet=3) 63 parameter(nytet=1) 64 parameter(nxhex=3) 65 parameter(ndofs=(nx*2+1)*(ny*2+1)) 66 parameter(nqptstet=nelemtet*qtet) 67 parameter(nqptshex=nelemhex*qhex*qhex) 68 parameter(nqpts=nqptstet+nqptshex) 69 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 70 real*8 arrx(d*ndofs) 71 integer*8 voffset,xoffset 72 73 real*8 qref(d*qtet) 74 real*8 qweight(qtet) 75 real*8 interp(ptet*qtet) 76 real*8 grad(d*ptet*qtet) 77 78 real*8 hv(ndofs) 79 real*8 total 80 81 character arg*32 82 83 external setup,mass 84 85 call getarg(1,arg) 86 87 call ceedinit(trim(arg)//char(0),ceed,err) 88 89! DoF Coordinates 90 do i=0,ny*2 91 do j=0,nx*2 92 arrx(i+j*(ny*2+1)+0*ndofs+1)=1.d0*i/(2*ny) 93 arrx(i+j*(ny*2+1)+1*ndofs+1)=1.d0*j/(2*nx) 94 enddo 95 enddo 96 97 call ceedvectorcreate(ceed,d*ndofs,x,err) 98 xoffset=0 99 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 100 101! Qdata Vectors 102 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 103 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 104 105! Tet Elements 106 do i=0,2 107 col=mod(i,nx) 108 row=i/nx 109 offset=col*2+row*(nx*2+1)*2 110 111 indxtet(i*2*ptet+1)=2+offset 112 indxtet(i*2*ptet+2)=9+offset 113 indxtet(i*2*ptet+3)=16+offset 114 indxtet(i*2*ptet+4)=1+offset 115 indxtet(i*2*ptet+5)=8+offset 116 indxtet(i*2*ptet+6)=0+offset 117 118 indxtet(i*2*ptet+7)=14+offset 119 indxtet(i*2*ptet+8)=7+offset 120 indxtet(i*2*ptet+9)=0+offset 121 indxtet(i*2*ptet+10)=15+offset 122 indxtet(i*2*ptet+11)=8+offset 123 indxtet(i*2*ptet+12)=16+offset 124 enddo 125 126! -- Restrictions 127 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,d,& 128 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 129 call ceedelemrestrictioncreateidentity(ceed,imode,nelemtet,ptet,& 130 & nelemtet*ptet,d,erestrictxitet,err) 131 132 call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,& 133 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 134 call ceedelemrestrictioncreateidentity(ceed,imode,nelemtet,qtet,nqptstet,& 135 & 1,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',erestrictxitet,& 165 & 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 call ceedelemrestrictioncreateidentity(ceed,imode,nelemhex,phex*phex,& 196 & nelemhex*phex*phex,d,erestrictxihex,err) 197 198 call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,& 199 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 200 call ceedelemrestrictioncreateidentity(ceed,imode,nelemhex,qhex*qhex,& 201 & nqptshex,1,erestrictuihex,err) 202 203! -- Bases 204 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 205 & bxhex,err) 206 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 207 & buhex,err) 208 209! -- QFunctions 210 call ceedqfunctioncreateinterior(ceed,1,setup,& 211 &SOURCE_DIR& 212 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 213 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 214 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 215 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 216 217 call ceedqfunctioncreateinterior(ceed,1,mass,& 218 &SOURCE_DIR& 219 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 220 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 221 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 222 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 223 224! -- Operators 225! ---- Setup Hex 226 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 227 & ceed_qfunction_none,op_setuphex,err) 228 call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,& 229 & bxhex,ceed_vector_none,err) 230 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 231 & bxhex,ceed_vector_active,err) 232 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 233 & ceed_basis_collocated,qdatahex,err) 234! ---- Mass Hex 235 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 236 & ceed_qfunction_none,op_masshex,err) 237 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 238 & ceed_basis_collocated,qdatahex,err) 239 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 240 & buhex,ceed_vector_active,err) 241 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 242 & 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 262 call ceedoperatorapply(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! Cleanup 278 call ceedqfunctiondestroy(qf_setuptet,err) 279 call ceedqfunctiondestroy(qf_masstet,err) 280 call ceedoperatordestroy(op_setuptet,err) 281 call ceedoperatordestroy(op_masstet,err) 282 call ceedqfunctiondestroy(qf_setuphex,err) 283 call ceedqfunctiondestroy(qf_masshex,err) 284 call ceedoperatordestroy(op_setuphex,err) 285 call ceedoperatordestroy(op_masshex,err) 286 call ceedoperatordestroy(op_setup,err) 287 call ceedoperatordestroy(op_mass,err) 288 call ceedelemrestrictiondestroy(erestrictutet,err) 289 call ceedelemrestrictiondestroy(erestrictxtet,err) 290 call ceedelemrestrictiondestroy(erestrictuitet,err) 291 call ceedelemrestrictiondestroy(erestrictxitet,err) 292 call ceedelemrestrictiondestroy(erestrictuhex,err) 293 call ceedelemrestrictiondestroy(erestrictxhex,err) 294 call ceedelemrestrictiondestroy(erestrictuihex,err) 295 call ceedelemrestrictiondestroy(erestrictxihex,err) 296 call ceedbasisdestroy(butet,err) 297 call ceedbasisdestroy(bxtet,err) 298 call ceedbasisdestroy(buhex,err) 299 call ceedbasisdestroy(bxhex,err) 300 call ceedvectordestroy(x,err) 301 call ceedvectordestroy(u,err) 302 call ceedvectordestroy(v,err) 303 call ceedvectordestroy(qdatatet,err) 304 call ceedvectordestroy(qdatahex,err) 305 call ceeddestroy(ceed,err) 306 end 307!----------------------------------------------------------------------- 308