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 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 69 real*8 qref(d*qtet) 70 real*8 qweight(qtet) 71 real*8 interp(ptet*qtet) 72 real*8 grad(d*ptet*qtet) 73 74 character arg*32 75 76 external setup,mass 77 78 call getarg(1,arg) 79 80 call ceedinit(trim(arg)//char(0),ceed,err) 81 82! DoF Coordinates 83 call ceedvectorcreate(ceed,d*ndofs,x,err) 84 85! Qdata Vectors 86 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 87 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 88 89! Tet Elements 90 do i=0,2 91 col=mod(i,nx) 92 row=i/nx 93 offset=col*2+row*(nx*2+1)*2 94 95 indxtet(i*2*ptet+1)=2+offset 96 indxtet(i*2*ptet+2)=9+offset 97 indxtet(i*2*ptet+3)=16+offset 98 indxtet(i*2*ptet+4)=1+offset 99 indxtet(i*2*ptet+5)=8+offset 100 indxtet(i*2*ptet+6)=0+offset 101 102 indxtet(i*2*ptet+7)=14+offset 103 indxtet(i*2*ptet+8)=7+offset 104 indxtet(i*2*ptet+9)=0+offset 105 indxtet(i*2*ptet+10)=15+offset 106 indxtet(i*2*ptet+11)=8+offset 107 indxtet(i*2*ptet+12)=16+offset 108 enddo 109 110! -- Restrictions 111 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,ndofs,d,ceed_mem_host,& 112 & ceed_use_pointer,indxtet,erestrictxtet,err) 113 call ceedelemrestrictioncreateidentity(ceed,nelemtet,ptet,nelemtet*ptet,& 114 & d,erestrictxitet,err) 115 116 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,ndofs,1,ceed_mem_host,& 117 & ceed_use_pointer,indxtet,erestrictutet,err) 118 call ceedelemrestrictioncreateidentity(ceed,nelemtet,qtet,nqptstet,1,& 119 & erestrictuitet,err) 120 121! -- Bases 122 call buildmats(qref,qweight,interp,grad) 123 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 124 & qweight,bxtet,err) 125 call buildmats(qref,qweight,interp,grad) 126 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 127 & qweight,butet,err) 128 129! -- QFunctions 130 call ceedqfunctioncreateinterior(ceed,1,setup,& 131 &SOURCE_DIR& 132 &//'t520-operator.h:setup'//char(0),qf_setuptet,err) 133 call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err) 134 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 135 call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err) 136 137 call ceedqfunctioncreateinterior(ceed,1,mass,& 138 &SOURCE_DIR& 139 &//'t520-operator.h:mass'//char(0),qf_masstet,err) 140 call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err) 141 call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err) 142 call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err) 143 144! -- Operators 145! ---- Setup Tet 146 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 147 & ceed_qfunction_none,op_setuptet,err) 148 call ceedoperatorsetfield(op_setuptet,'_weight',erestrictxitet,& 149 & ceed_notranspose,bxtet,ceed_vector_none,err) 150 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 151 & ceed_notranspose,bxtet,ceed_vector_active,err) 152 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 153 & ceed_notranspose,ceed_basis_collocated,qdatatet,err) 154! ---- Mass Tet 155 call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 156 & ceed_qfunction_none,op_masstet,err) 157 call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 158 & ceed_notranspose,ceed_basis_collocated,qdatatet,err) 159 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 160 & ceed_notranspose,butet,ceed_vector_active,err) 161 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 162 & ceed_notranspose,butet,ceed_vector_active,err) 163 164! Hex Elements 165 do i=0,nelemhex-1 166 col=mod(i,nx) 167 row=i/nx 168 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 169 do j=0,phex-1 170 do k=0,phex-1 171 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 172 enddo 173 enddo 174 enddo 175 176! -- Restrictions 177 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,ndofs,d,& 178 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 179 call ceedelemrestrictioncreateidentity(ceed,nelemhex,phex*phex,& 180 & nelemhex*phex*phex,d,erestrictxihex,err) 181 182 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,ndofs,1,& 183 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 184 call ceedelemrestrictioncreateidentity(ceed,nelemhex,qhex*qhex,nqptshex,& 185 & 1,erestrictuihex,err) 186 187! -- Bases 188 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 189 & bxhex,err) 190 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 191 & buhex,err) 192 193! -- QFunctions 194 call ceedqfunctioncreateinterior(ceed,1,setup,& 195 &SOURCE_DIR& 196 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 197 call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err) 198 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 199 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 200 201 call ceedqfunctioncreateinterior(ceed,1,mass,& 202 &SOURCE_DIR& 203 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 204 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 205 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 206 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 207 208! -- Operators 209! ---- Setup Hex 210 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 211 & ceed_qfunction_none,op_setuphex,err) 212 call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,& 213 & ceed_notranspose,bxhex,ceed_vector_none,err) 214 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 215 & ceed_notranspose,bxhex,ceed_vector_active,err) 216 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 217 & ceed_notranspose,ceed_basis_collocated,qdatahex,err) 218! ---- Mass Hex 219 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 220 & ceed_qfunction_none,op_masshex,err) 221 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 222 & ceed_notranspose,ceed_basis_collocated,qdatahex,err) 223 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 224 & ceed_notranspose,buhex,ceed_vector_active,err) 225 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 226 & ceed_notranspose,buhex,ceed_vector_active,err) 227 228! Composite Operators 229 call ceedcompositeoperatorcreate(ceed,op_setup,err) 230 call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err) 231 call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err) 232 233 call ceedcompositeoperatorcreate(ceed,op_mass,err) 234 call ceedcompositeoperatoraddsub(op_mass,op_masstet,err) 235 call ceedcompositeoperatoraddsub(op_mass,op_masshex,err) 236 237! View 238 call ceedoperatorview(op_setup,err) 239 call ceedoperatorview(op_mass,err) 240 241! Cleanup 242 call ceedqfunctiondestroy(qf_setuptet,err) 243 call ceedqfunctiondestroy(qf_masstet,err) 244 call ceedoperatordestroy(op_setuptet,err) 245 call ceedoperatordestroy(op_masstet,err) 246 call ceedqfunctiondestroy(qf_setuphex,err) 247 call ceedqfunctiondestroy(qf_masshex,err) 248 call ceedoperatordestroy(op_setuphex,err) 249 call ceedoperatordestroy(op_masshex,err) 250 call ceedoperatordestroy(op_setup,err) 251 call ceedoperatordestroy(op_mass,err) 252 call ceedelemrestrictiondestroy(erestrictutet,err) 253 call ceedelemrestrictiondestroy(erestrictxtet,err) 254 call ceedelemrestrictiondestroy(erestrictuitet,err) 255 call ceedelemrestrictiondestroy(erestrictxitet,err) 256 call ceedelemrestrictiondestroy(erestrictuhex,err) 257 call ceedelemrestrictiondestroy(erestrictxhex,err) 258 call ceedelemrestrictiondestroy(erestrictuihex,err) 259 call ceedelemrestrictiondestroy(erestrictxihex,err) 260 call ceedbasisdestroy(butet,err) 261 call ceedbasisdestroy(bxtet,err) 262 call ceedbasisdestroy(buhex,err) 263 call ceedbasisdestroy(bxhex,err) 264 call ceedvectordestroy(x,err) 265 call ceedvectordestroy(qdatatet,err) 266 call ceedvectordestroy(qdatahex,err) 267 call ceeddestroy(ceed,err) 268 end 269!----------------------------------------------------------------------- 270