1!----------------------------------------------------------------------- 2! 3! Header with QFunctions 4! 5 include 't531-operator-f.h' 6!----------------------------------------------------------------------- 7 program test 8 implicit none 9 include 'ceedf.h' 10 11 integer ceed,err,i,j,k 12 integer stridesu(3),stridesqd(3) 13 integer erestrictx,erestrictu,erestrictui 14 integer erestrictqi,erestrictlini 15 integer bx,bu 16 integer qf_setup,qf_diff,qf_diff_lin 17 integer op_setup,op_diff,op_diff_lin 18 integer qdata,x,a,u,v 19 integer nelem,p,q,d 20 integer row,col,offset 21 parameter(nelem=6) 22 parameter(p=3) 23 parameter(q=4) 24 parameter(d=2) 25 integer ndofs,nqpts,nx,ny 26 parameter(nx=3) 27 parameter(ny=2) 28 parameter(ndofs=(nx*2+1)*(ny*2+1)) 29 parameter(nqpts=nelem*q*q) 30 integer indx(nelem*p*p) 31 real*8 arrx(d*ndofs),vv(ndofs) 32 integer*8 xoffset,voffset 33 34 character arg*32 35 36 external setup,diff,diff_lin 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*d*(d+1)/2,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 stridesqd=[1,q*q,q*q*d*(d+1)/2] 79 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,d*(d+1)/2,& 80 & d*(d+1)/2*nqpts,stridesqd,erestrictqi,err) 81 82! Bases 83 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 84 & bx,err) 85 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 86 & bu,err) 87 88! QFunction - setup 89 call ceedqfunctioncreateinterior(ceed,1,setup,& 90 &SOURCE_DIR& 91 &//'t531-operator.h:setup'//char(0),qf_setup,err) 92 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 93 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 94 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 95 96! Operator - setup 97 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 98 & ceed_qfunction_none,op_setup,err) 99 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 100 & bx,ceed_vector_active,err) 101 call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,& 102 & bx,ceed_vector_none,err) 103 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 104 & ceed_basis_collocated,ceed_vector_active,err) 105 106! Apply Setup Operator 107 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 108 109! QFunction - apply 110 call ceedqfunctioncreateinterior(ceed,1,diff,& 111 &SOURCE_DIR& 112 &//'t531-operator.h:diff'//char(0),qf_diff,err) 113 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 114 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 115 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 116 117! Operator - apply 118 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 119 & ceed_qfunction_none,op_diff,err) 120 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 121 & bu,ceed_vector_active,err) 122 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 123 & ceed_basis_collocated,qdata,err) 124 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 125 & bu,ceed_vector_active,err) 126 127! Apply original Poisson Operator 128 call ceedvectorcreate(ceed,ndofs,u,err) 129 call ceedvectorsetvalue(u,1.d0,err) 130 call ceedvectorcreate(ceed,ndofs,v,err) 131 call ceedvectorsetvalue(v,0.d0,err) 132 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 133 134! Check Output 135 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 136 do i=1,ndofs 137 if (abs(vv(voffset+i))>1.0d-14) then 138! LCOV_EXCL_START 139 write(*,*) 'Error: Operator computed v[i] = ',vv(voffset+i),' != 0.0' 140! LCOV_EXCL_STOP 141 endif 142 enddo 143 call ceedvectorrestorearrayread(v,vv,voffset,err) 144 145! Assemble QFunction 146 call ceedoperatorlinearassembleqfunction(op_diff,a,erestrictlini,& 147 & ceed_request_immediate,err) 148 149! QFunction - apply linearized 150 call ceedqfunctioncreateinterior(ceed,1,diff_lin,& 151 &SOURCE_DIR& 152 &//'t531-operator.h:diff_lin'//char(0),qf_diff_lin,err) 153 call ceedqfunctionaddinput(qf_diff_lin,'du',d,ceed_eval_grad,err) 154 call ceedqfunctionaddinput(qf_diff_lin,'qdata',d*d,ceed_eval_none,err) 155 call ceedqfunctionaddoutput(qf_diff_lin,'dv',d,ceed_eval_grad,err) 156 157! Operator - apply linearized 158 call ceedoperatorcreate(ceed,qf_diff_lin,ceed_qfunction_none,& 159 & ceed_qfunction_none,op_diff_lin,err) 160 call ceedoperatorsetfield(op_diff_lin,'du',erestrictu,& 161 & bu,ceed_vector_active,err) 162 call ceedoperatorsetfield(op_diff_lin,'qdata',erestrictlini,& 163 & ceed_basis_collocated,a,err) 164 call ceedoperatorsetfield(op_diff_lin,'dv',erestrictu,& 165 & bu,ceed_vector_active,err) 166 167! Apply linearized Poisson Operator 168 call ceedvectorsetvalue(v,0.d0,err) 169 call ceedoperatorapply(op_diff_lin,u,v,ceed_request_immediate,err) 170 171! Check Output 172 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 173 do i=1,ndofs 174 if (abs(vv(voffset+i))>1.0d-14) then 175! LCOV_EXCL_START 176 write(*,*) 'Error: Linearized operator computed v[i] = ',vv(voffset+i),& 177 & ' != 0.0' 178! LCOV_EXCL_STOP 179 endif 180 enddo 181 call ceedvectorrestorearrayread(v,vv,voffset,err) 182 183! Cleanup 184 call ceedqfunctiondestroy(qf_setup,err) 185 call ceedqfunctiondestroy(qf_diff,err) 186 call ceedqfunctiondestroy(qf_diff_lin,err) 187 call ceedoperatordestroy(op_setup,err) 188 call ceedoperatordestroy(op_diff,err) 189 call ceedoperatordestroy(op_diff_lin,err) 190 call ceedelemrestrictiondestroy(erestrictu,err) 191 call ceedelemrestrictiondestroy(erestrictx,err) 192 call ceedelemrestrictiondestroy(erestrictui,err) 193 call ceedelemrestrictiondestroy(erestrictqi,err) 194 call ceedelemrestrictiondestroy(erestrictlini,err) 195 call ceedbasisdestroy(bu,err) 196 call ceedbasisdestroy(bx,err) 197 call ceedvectordestroy(x,err) 198 call ceedvectordestroy(a,err) 199 call ceedvectordestroy(u,err) 200 call ceedvectordestroy(v,err) 201 call ceedvectordestroy(qdata,err) 202 call ceeddestroy(ceed,err) 203 end 204!----------------------------------------------------------------------- 205