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