1!----------------------------------------------------------------------- 2! 3! Header with common subroutine 4! 5 include 't320-basis-f.h' 6! 7! Header with QFunctions 8! 9 include 't510-operator-f.h' 10!----------------------------------------------------------------------- 11 program test 12 implicit none 13 include 'ceed/fortran.h' 14 15 integer ceed,err,i,j,k 16 integer stridesutet(3),stridesuhex(3) 17 integer erestrictxtet,erestrictutet,erestrictuitet,& 18& erestrictxhex,erestrictuhex,erestrictuihex 19 integer bxtet,butet,bxhex,buhex 20 integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex 21 integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass 22 integer qdatatet,qdatahex,x 23 integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d 24 integer row,col,offset 25 parameter(nelemtet=6) 26 parameter(ptet=6) 27 parameter(qtet=4) 28 parameter(nelemhex=6) 29 parameter(phex=3) 30 parameter(qhex=4) 31 parameter(d=2) 32 integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex 33 parameter(nx=3) 34 parameter(ny=3) 35 parameter(nxtet=3) 36 parameter(nytet=1) 37 parameter(nxhex=3) 38 parameter(ndofs=(nx*2+1)*(ny*2+1)) 39 parameter(nqptstet=nelemtet*qtet) 40 parameter(nqptshex=nelemhex*qhex*qhex) 41 parameter(nqpts=nqptstet+nqptshex) 42 integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex) 43 44 real*8 qref(d*qtet) 45 real*8 qweight(qtet) 46 real*8 interp(ptet*qtet) 47 real*8 grad(d*ptet*qtet) 48 49 character arg*32 50 51! LCOV_EXCL_START 52 external setup,mass 53! LCOV_EXCL_STOP 54 55 call getarg(1,arg) 56 57 call ceedinit(trim(arg)//char(0),ceed,err) 58 59! DoF Coordinates 60 call ceedvectorcreate(ceed,d*ndofs,x,err) 61 62! Qdata Vectors 63 call ceedvectorcreate(ceed,nqptstet,qdatatet,err) 64 call ceedvectorcreate(ceed,nqptshex,qdatahex,err) 65 66! Tet Elements 67 do i=0,2 68 col=mod(i,nx) 69 row=i/nx 70 offset=col*2+row*(nx*2+1)*2 71 72 indxtet(i*2*ptet+1)=2+offset 73 indxtet(i*2*ptet+2)=9+offset 74 indxtet(i*2*ptet+3)=16+offset 75 indxtet(i*2*ptet+4)=1+offset 76 indxtet(i*2*ptet+5)=8+offset 77 indxtet(i*2*ptet+6)=0+offset 78 79 indxtet(i*2*ptet+7)=14+offset 80 indxtet(i*2*ptet+8)=7+offset 81 indxtet(i*2*ptet+9)=0+offset 82 indxtet(i*2*ptet+10)=15+offset 83 indxtet(i*2*ptet+11)=8+offset 84 indxtet(i*2*ptet+12)=16+offset 85 enddo 86 87! -- Restrictions 88 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,& 89 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err) 90 91 call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,& 92 & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err) 93 stridesutet=[1,qtet,qtet] 94 call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,nqptstet,& 95 & stridesutet,erestrictuitet,err) 96 97! -- Bases 98 call buildmats(qref,qweight,interp,grad) 99 call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,& 100 & qweight,bxtet,err) 101 call buildmats(qref,qweight,interp,grad) 102 call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,& 103 & qweight,butet,err) 104 105! -- QFunctions 106 call ceedqfunctioncreateinterior(ceed,1,setup,& 107 &SOURCE_DIR& 108 &//'t510-operator.h:setup'//char(0),qf_setuptet,err) 109 call ceedqfunctionaddinput(qf_setuptet,'weight',1,ceed_eval_weight,err) 110 call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err) 111 call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err) 112 113 call ceedqfunctioncreateinterior(ceed,1,mass,& 114 &SOURCE_DIR& 115 &//'t510-operator.h:mass'//char(0),qf_masstet,err) 116 call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err) 117 call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err) 118 call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err) 119 120! -- Operators 121! ---- Setup Tet 122 call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,& 123 & ceed_qfunction_none,op_setuptet,err) 124 call ceedoperatorsetname(op_setuptet,'triangle elements',err) 125 call ceedoperatorsetfield(op_setuptet,'weight',& 126 & ceed_elemrestriction_none,bxtet,ceed_vector_none,err) 127 call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,& 128 & bxtet,ceed_vector_active,err) 129 call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,& 130 ceed_basis_none,qdatatet,err) 131! ---- Mass Tet 132 call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,& 133 & ceed_qfunction_none,op_masstet,err) 134 call ceedoperatorsetname(op_masstet,'triangle elements',err) 135 call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,& 136 ceed_basis_none,qdatatet,err) 137 call ceedoperatorsetfield(op_masstet,'u',erestrictutet,& 138 & butet,ceed_vector_active,err) 139 call ceedoperatorsetfield(op_masstet,'v',erestrictutet,& 140 & butet,ceed_vector_active,err) 141 142! Hex Elements 143 do i=0,nelemhex-1 144 col=mod(i,nx) 145 row=i/nx 146 offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2 147 do j=0,phex-1 148 do k=0,phex-1 149 indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j 150 enddo 151 enddo 152 enddo 153 154! -- Restrictions 155 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,d*ndofs,& 156 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err) 157 158 call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,ndofs,& 159 & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err) 160 stridesuhex=[1,qhex*qhex,qhex*qhex] 161 call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,& 162 & 1,nqptshex,stridesuhex,erestrictuihex,err) 163 164! -- Bases 165 call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,& 166 & bxhex,err) 167 call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,& 168 & buhex,err) 169 170! -- QFunctions 171 call ceedqfunctioncreateinterior(ceed,1,setup,& 172 &SOURCE_DIR& 173 &//'t510-operator.h:setup'//char(0),qf_setuphex,err) 174 call ceedqfunctionaddinput(qf_setuphex,'weight',1,ceed_eval_weight,err) 175 call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err) 176 call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err) 177 178 call ceedqfunctioncreateinterior(ceed,1,mass,& 179 &SOURCE_DIR& 180 &//'t510-operator.h:mass'//char(0),qf_masshex,err) 181 call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err) 182 call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err) 183 call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err) 184 185! -- Operators 186! ---- Setup Hex 187 call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,& 188 & ceed_qfunction_none,op_setuphex,err) 189 call ceedoperatorsetname(op_setuphex,'quadrilateral elements',err) 190 call ceedoperatorsetfield(op_setuphex,'weight',& 191 & ceed_elemrestriction_none,bxhex,ceed_vector_none,err) 192 call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,& 193 & bxhex,ceed_vector_active,err) 194 call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,& 195 ceed_basis_none,qdatahex,err) 196! ---- Mass Hex 197 call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,& 198 & ceed_qfunction_none,op_masshex,err) 199 call ceedoperatorsetname(op_masshex,'quadrilateral elements',err) 200 call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,& 201 ceed_basis_none,qdatahex,err) 202 call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,& 203 & buhex,ceed_vector_active,err) 204 call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,& 205 & buhex,ceed_vector_active,err) 206 207! Composite Operators 208 call ceedoperatorcreatecomposite(ceed,op_setup,err) 209 call ceedoperatorsetname(op_setup,'setup',err) 210 call ceedoperatorcompositeaddsub(op_setup,op_setuptet,err) 211 call ceedoperatorcompositeaddsub(op_setup,op_setuphex,err) 212 213 call ceedoperatorcreatecomposite(ceed,op_mass,err) 214 call ceedoperatorsetname(op_mass,'mass',err) 215 call ceedoperatorcompositeaddsub(op_mass,op_masstet,err) 216 call ceedoperatorcompositeaddsub(op_mass,op_masshex,err) 217 218! View 219 call ceedoperatorview(op_setup,err) 220 call ceedoperatorview(op_mass,err) 221 222! Cleanup 223 call ceedqfunctiondestroy(qf_setuptet,err) 224 call ceedqfunctiondestroy(qf_masstet,err) 225 call ceedoperatordestroy(op_setuptet,err) 226 call ceedoperatordestroy(op_masstet,err) 227 call ceedqfunctiondestroy(qf_setuphex,err) 228 call ceedqfunctiondestroy(qf_masshex,err) 229 call ceedoperatordestroy(op_setuphex,err) 230 call ceedoperatordestroy(op_masshex,err) 231 call ceedoperatordestroy(op_setup,err) 232 call ceedoperatordestroy(op_mass,err) 233 call ceedelemrestrictiondestroy(erestrictutet,err) 234 call ceedelemrestrictiondestroy(erestrictxtet,err) 235 call ceedelemrestrictiondestroy(erestrictuitet,err) 236 call ceedelemrestrictiondestroy(erestrictuhex,err) 237 call ceedelemrestrictiondestroy(erestrictxhex,err) 238 call ceedelemrestrictiondestroy(erestrictuihex,err) 239 call ceedbasisdestroy(butet,err) 240 call ceedbasisdestroy(bxtet,err) 241 call ceedbasisdestroy(buhex,err) 242 call ceedbasisdestroy(bxhex,err) 243 call ceedvectordestroy(x,err) 244 call ceedvectordestroy(qdatatet,err) 245 call ceedvectordestroy(qdatahex,err) 246 call ceeddestroy(ceed,err) 247 end 248!----------------------------------------------------------------------- 249