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