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 stridesu(3) 40 integer erestrictx,erestrictu,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 99 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 100 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 101 stridesu=[1,q*q,q*q] 102 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 103 & 1,stridesu,erestrictui,err) 104 105! Bases 106 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 107 & bx,err) 108 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 109 & bu,err) 110 111! QFunctions 112! -- Setup 113 call ceedqfunctioncreateinterior(ceed,1,setup,& 114 &SOURCE_DIR& 115 &//'t510-operator.h:setup'//char(0),qf_setup,err) 116 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 117 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 118 call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err) 119! -- Mass 120 call ceedqfunctioncreateinterior(ceed,1,mass,& 121 &SOURCE_DIR& 122 &//'t510-operator.h:mass'//char(0),qf_mass,err) 123 call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err) 124 call ceedqfunctionaddinput(qf_mass,'u',1,ceed_eval_interp,err) 125 call ceedqfunctionaddoutput(qf_mass,'v',1,ceed_eval_interp,err) 126 127! Operators 128! -- Setup 129 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 130 & ceed_qfunction_none,op_setup,err) 131 call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,& 132 & bx,ceed_vector_none,err) 133 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 134 & bx,ceed_vector_active,err) 135 call ceedoperatorsetfield(op_setup,'rho',erestrictui,& 136 & ceed_basis_collocated,ceed_vector_active,err) 137! -- Mass 138 call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,& 139 & ceed_qfunction_none,op_mass,err) 140 call ceedoperatorsetfield(op_mass,'rho',erestrictui,& 141 & ceed_basis_collocated,qdata,err) 142 call ceedoperatorsetfield(op_mass,'u',erestrictu,& 143 & bu,ceed_vector_active,err) 144 call ceedoperatorsetfield(op_mass,'v',erestrictu,& 145 & bu,ceed_vector_active,err) 146 147! Apply Setup Operator 148 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 149 150! Assemble QFunction 151 call ceedoperatorassemblelinearqfunction(op_mass,a,erestrictlini,& 152 & ceed_request_immediate,err) 153 154! Check Output 155 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 156 call ceedvectorgetarrayread(qdata,ceed_mem_host,qq,qoffset,err) 157 do i=1,nqpts 158 if (abs(qq(qoffset+i)-aa(aoffset+i))>1.0d-9) then 159! LCOV_EXCL_START 160 write(*,*) 'Error: A[',i,'] = ',aa(aoffset+i),' != ',& 161 & qq(qoffset+i) 162! LCOV_EXCL_STOP 163 endif 164 enddo 165 call ceedvectorrestorearrayread(a,aa,aoffset,err) 166 call ceedvectorrestorearrayread(qdata,qq,qoffset,err) 167 168! Apply original Mass Operator 169 call ceedvectorcreate(ceed,ndofs,u,err) 170 call ceedvectorsetvalue(u,1.d0,err) 171 call ceedvectorcreate(ceed,ndofs,v,err) 172 call ceedvectorsetvalue(v,0.d0,err) 173 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 174 175! Check Output 176 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 177 total=0. 178 do i=1,ndofs 179 total=total+vv(voffset+i) 180 enddo 181 if (abs(total-1.)>1.0d-14) then 182! LCOV_EXCL_START 183 write(*,*) 'Error: True operator computed area = ',total,' != 1.0' 184! LCOV_EXCL_STOP 185 endif 186 call ceedvectorrestorearrayread(v,vv,voffset,err) 187 188! Switch to new qdata 189 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 190 call ceedvectorsetarray(qdata,ceed_mem_host,ceed_copy_values,aa,& 191 & aoffset,err) 192 call ceedvectorrestorearrayread(a,aa,aoffset,err) 193 194! Apply new Mass Operator 195 call ceedvectorsetvalue(v,0.d0,err) 196 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 197 198! Check Output 199 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 200 total=0. 201 do i=1,ndofs 202 total=total+vv(voffset+i) 203 enddo 204 if (abs(total-1.)>1.0d-10) then 205! LCOV_EXCL_START 206 write(*,*) 'Error: Assembled operator computed area = ',total,' != 1.0' 207! LCOV_EXCL_STOP 208 endif 209 call ceedvectorrestorearrayread(v,vv,voffset,err) 210 211! Cleanup 212 call ceedqfunctiondestroy(qf_setup,err) 213 call ceedqfunctiondestroy(qf_mass,err) 214 call ceedoperatordestroy(op_setup,err) 215 call ceedoperatordestroy(op_mass,err) 216 call ceedelemrestrictiondestroy(erestrictu,err) 217 call ceedelemrestrictiondestroy(erestrictx,err) 218 call ceedelemrestrictiondestroy(erestrictui,err) 219 call ceedelemrestrictiondestroy(erestrictlini,err) 220 call ceedbasisdestroy(bu,err) 221 call ceedbasisdestroy(bx,err) 222 call ceedvectordestroy(x,err) 223 call ceedvectordestroy(a,err) 224 call ceedvectordestroy(u,err) 225 call ceedvectordestroy(v,err) 226 call ceedvectordestroy(qdata,err) 227 call ceeddestroy(ceed,err) 228 end 229!----------------------------------------------------------------------- 230