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 erestrictx,erestrictu,erestrictxi,erestrictui,erestrictlini 40 integer bx,bu 41 integer qf_setup,qf_mass 42 integer op_setup,op_mass 43 integer qdata,x,a,u,v 44 integer nelem,p,q,d 45 integer row,col,offset 46 parameter(nelem=6) 47 parameter(p=3) 48 parameter(q=4) 49 parameter(d=2) 50 integer ndofs,nqpts,nx,ny 51 parameter(nx=3) 52 parameter(ny=2) 53 parameter(ndofs=(nx*2+1)*(ny*2+1)) 54 parameter(nqpts=nelem*q*q) 55 integer indx(nelem*p*p) 56 real*8 arrx(d*ndofs),aa(nqpts),qq(nqpts),vv(ndofs) 57 integer*8 xoffset,aoffset,qoffset,voffset 58 real*8 total 59 60 character arg*32 61 62 external setup,mass 63 64 call getarg(1,arg) 65 66 call ceedinit(trim(arg)//char(0),ceed,err) 67 68! DoF Coordinates 69 do i=0,nx*2 70 do j=0,ny*2 71 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 72 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 73 enddo 74 enddo 75 call ceedvectorcreate(ceed,d*ndofs,x,err) 76 xoffset=0 77 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 78 79! Qdata Vector 80 call ceedvectorcreate(ceed,nqpts,qdata,err) 81 82! Element Setup 83 do i=0,nelem-1 84 col=mod(i,nx) 85 row=i/nx 86 offset=col*(p-1)+row*(nx*2+1)*(p-1) 87 do j=0,p-1 88 do k=0,p-1 89 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 90 enddo 91 enddo 92 enddo 93 94! Restrictions 95 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 96 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 97 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,p*p,& 98 & nelem*p*p,d,erestrictxi,err) 99 100 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 101 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 102 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 103 & 1,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',erestrictxi,& 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(erestrictxi,err) 220 call ceedelemrestrictiondestroy(erestrictlini,err) 221 call ceedbasisdestroy(bu,err) 222 call ceedbasisdestroy(bx,err) 223 call ceedvectordestroy(x,err) 224 call ceedvectordestroy(a,err) 225 call ceedvectordestroy(u,err) 226 call ceedvectordestroy(v,err) 227 call ceedvectordestroy(qdata,err) 228 call ceeddestroy(ceed,err) 229 end 230!----------------------------------------------------------------------- 231