1!----------------------------------------------------------------------- 2 subroutine setup(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 3& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 4 real*8 ctx 5 real*8 u1(1) 6 real*8 u2(1) 7 real*8 v1(1) 8 integer q,ierr 9 10 do i=1,q 11 v1(i)=u1(i)*(u2(i+q*0)*u2(i+q*3)-u2(i+q*1)*u2(i+q*2)) 12 enddo 13 14 ierr=0 15 end 16!----------------------------------------------------------------------- 17 subroutine mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 18& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 19 real*8 ctx 20 real*8 u1(1) 21 real*8 u2(1) 22 real*8 v1(1) 23 integer q,ierr 24 25 do i=1,q 26 v1(i)=u2(i)*u1(i) 27 enddo 28 29 ierr=0 30 end 31!----------------------------------------------------------------------- 32 program test 33 34 include 'ceedf.h' 35 36 integer ceed,err,i,j,k 37 integer erestrictx,erestrictu,erestrictxi,erestrictui,erestrictlini 38 integer bx,bu 39 integer qf_setup,qf_mass 40 integer op_setup,op_mass 41 integer qdata,x,a,u,v 42 integer nelem,p,q,d 43 integer row,col,offset 44 parameter(nelem=6) 45 parameter(p=3) 46 parameter(q=4) 47 parameter(d=2) 48 integer ndofs,nqpts,nx,ny 49 parameter(nx=3) 50 parameter(ny=2) 51 parameter(ndofs=(nx*2+1)*(ny*2+1)) 52 parameter(nqpts=nelem*q*q) 53 integer indx(nelem*p*p) 54 real*8 arrx(d*ndofs),aa(nqpts),qq(nqpts),vv(ndofs) 55 integer*8 xoffset,aoffset,qoffset,voffset 56 real*8 total 57 58 character arg*32 59 60 external setup,mass 61 62 call getarg(1,arg) 63 64 call ceedinit(trim(arg)//char(0),ceed,err) 65 66! DoF Coordinates 67 do i=0,nx*2 68 do j=0,ny*2 69 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 70 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 71 enddo 72 enddo 73 call ceedvectorcreate(ceed,d*ndofs,x,err) 74 xoffset=0 75 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 76 77! Qdata Vector 78 call ceedvectorcreate(ceed,nqpts,qdata,err) 79 80! Element Setup 81 do i=0,nelem-1 82 col=mod(i,nx) 83 row=i/nx 84 offset=col*(p-1)+row*(nx*2+1)*(p-1) 85 do j=0,p-1 86 do k=0,p-1 87 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 88 enddo 89 enddo 90 enddo 91 92! Restrictions 93 call ceedelemrestrictioncreate(ceed,nelem,p*p,ndofs,d,& 94 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 95 call ceedelemrestrictioncreateidentity(ceed,nelem,p*p,& 96 & nelem*p*p,d,erestrictxi,err) 97 98 call ceedelemrestrictioncreate(ceed,nelem,p*p,ndofs,1,& 99 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 100 call ceedelemrestrictioncreateidentity(ceed,nelem,q*q,nqpts,& 101 & 1,erestrictui,err) 102 103! Bases 104 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 105 & bx,err) 106 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 107 & bu,err) 108 109! QFunctions 110! -- Setup 111 call ceedqfunctioncreateinterior(ceed,1,setup,& 112 &SOURCE_DIR& 113 &//'t510-operator.h:setup'//char(0),qf_setup,err) 114 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 115 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 116 call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err) 117! -- Mass 118 call ceedqfunctioncreateinterior(ceed,1,mass,& 119 &SOURCE_DIR& 120 &//'t510-operator.h:mass'//char(0),qf_mass,err) 121 call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err) 122 call ceedqfunctionaddinput(qf_mass,'u',1,ceed_eval_interp,err) 123 call ceedqfunctionaddoutput(qf_mass,'v',1,ceed_eval_interp,err) 124 125! Operators 126! -- Setup 127 call ceedoperatorcreate(ceed,qf_setup,ceed_null,ceed_null,op_setup,& 128 & err) 129 call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,& 130 & ceed_notranspose,bx,ceed_vector_none,err) 131 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 132 & ceed_notranspose,bx,ceed_vector_active,err) 133 call ceedoperatorsetfield(op_setup,'rho',erestrictui,& 134 & ceed_notranspose,ceed_basis_collocated,ceed_vector_active,err) 135! -- Mass 136 call ceedoperatorcreate(ceed,qf_mass,ceed_null,ceed_null,op_mass,& 137 & err) 138 call ceedoperatorsetfield(op_mass,'rho',erestrictui,& 139 & ceed_notranspose,ceed_basis_collocated,qdata,err) 140 call ceedoperatorsetfield(op_mass,'u',erestrictu,& 141 & ceed_notranspose,bu,ceed_vector_active,err) 142 call ceedoperatorsetfield(op_mass,'v',erestrictu,& 143 & ceed_notranspose,bu,ceed_vector_active,err) 144 145! Apply Setup Operator 146 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 147 148! Assemble QFunction 149 call ceedoperatorassemblelinearqfunction(op_mass,a,erestrictlini,& 150 & ceed_request_immediate,err) 151 152! Check Output 153 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 154 call ceedvectorgetarrayread(qdata,ceed_mem_host,qq,qoffset,err) 155 do i=1,nqpts 156 if (abs(qq(qoffset+i)-aa(aoffset+i))>1.0d-9) then 157! LCOV_EXCL_START 158 write(*,*) 'Error: A[',i,'] = ',aa(aoffset+i),' != ',& 159 & qq(qoffset+i) 160! LCOV_EXCL_STOP 161 endif 162 enddo 163 call ceedvectorrestorearrayread(a,aa,aoffset,err) 164 call ceedvectorrestorearrayread(qdata,qq,qoffset,err) 165 166! Apply original Mass Operator 167 call ceedvectorcreate(ceed,ndofs,u,err) 168 call ceedvectorsetvalue(u,1.d0,err) 169 call ceedvectorcreate(ceed,ndofs,v,err) 170 call ceedvectorsetvalue(v,0.d0,err) 171 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 172 173! Check Output 174 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 175 total=0. 176 do i=1,ndofs 177 total=total+vv(voffset+i) 178 enddo 179 if (abs(total-1.)>1.0d-14) then 180! LCOV_EXCL_START 181 write(*,*) 'Error: True operator computed area = ',total,' != 1.0' 182! LCOV_EXCL_STOP 183 endif 184 call ceedvectorrestorearrayread(v,vv,voffset,err) 185 186! Switch to new qdata 187 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 188 call ceedvectorsetarray(qdata,ceed_mem_host,ceed_copy_values,aa,& 189 & aoffset,err) 190 call ceedvectorrestorearrayread(a,aa,aoffset,err) 191 192! Apply new Mass Operator 193 call ceedvectorsetvalue(v,0.d0,err) 194 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 195 196! Check Output 197 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 198 total=0. 199 do i=1,ndofs 200 total=total+vv(voffset+i) 201 enddo 202 if (abs(total-1.)>1.0d-10) then 203! LCOV_EXCL_START 204 write(*,*) 'Error: Assembled operator computed area = ',total,' != 1.0' 205! LCOV_EXCL_STOP 206 endif 207 call ceedvectorrestorearrayread(v,vv,voffset,err) 208 209! Cleanup 210 call ceedqfunctiondestroy(qf_setup,err) 211 call ceedqfunctiondestroy(qf_mass,err) 212 call ceedoperatordestroy(op_setup,err) 213 call ceedoperatordestroy(op_mass,err) 214 call ceedelemrestrictiondestroy(erestrictu,err) 215 call ceedelemrestrictiondestroy(erestrictx,err) 216 call ceedelemrestrictiondestroy(erestrictui,err) 217 call ceedelemrestrictiondestroy(erestrictxi,err) 218 call ceedelemrestrictiondestroy(erestrictlini,err) 219 call ceedbasisdestroy(bu,err) 220 call ceedbasisdestroy(bx,err) 221 call ceedvectordestroy(x,err) 222 call ceedvectordestroy(a,err) 223 call ceedvectordestroy(u,err) 224 call ceedvectordestroy(v,err) 225 call ceedvectordestroy(qdata,err) 226 call ceeddestroy(ceed,err) 227 end 228!----------------------------------------------------------------------- 229