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 real*8 w 9 integer q,ierr 10 11 do i=1,q 12 w=u2(i)/(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 13 v1(i+q*0)=w*(u1(i+q*2)*u1(i+q*2)+u1(i+q*3)*u1(i+q*3)) 14 v1(i+q*1)=-w*(u1(i+q*0)*u1(i+q*2)+u1(i+q*2)*u1(i+q*3)) 15 v1(i+q*2)=w*(u1(i+q*0)*u1(i+q*0)+u1(i+q*1)*u1(i+q*1)) 16 enddo 17 18 ierr=0 19 end 20!----------------------------------------------------------------------- 21 subroutine diff(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 22& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 23 real*8 ctx 24 real*8 u1(1) 25 real*8 u2(1) 26 real*8 v1(1) 27 real*8 du0,du1 28 integer q,ierr 29 30 do i=1,q 31 du0=u1(i+q*0) 32 du1=u1(i+q*1) 33 v1(i+q*0)=u2(i+q*0)*du0+u2(i+q*2)*du1 34 v1(i+q*1)=u2(i+q*2)*du0+u2(i+q*1)*du1 35 enddo 36 37 ierr=0 38 end 39!----------------------------------------------------------------------- 40 program test 41 42 include 'ceedf.h' 43 44 integer ceed,err,i,j,k 45 integer imode 46 parameter(imode=ceed_noninterlaced) 47 integer stridesx(3),stridesu(3),stridesqd(3) 48 integer erestrictx,erestrictu,erestrictxi,erestrictui,erestrictqi 49 integer bx,bu 50 integer qf_setup,qf_diff 51 integer op_setup,op_diff 52 integer qdata,x,a,u,v 53 integer nelem,p,q,d 54 integer row,col,offset 55 parameter(nelem=6) 56 parameter(p=3) 57 parameter(q=4) 58 parameter(d=2) 59 integer ndofs,nqpts,nx,ny 60 parameter(nx=3) 61 parameter(ny=2) 62 parameter(ndofs=(nx*2+1)*(ny*2+1)) 63 parameter(nqpts=nelem*q*q) 64 integer indx(nelem*p*p) 65 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 66 integer*8 xoffset,aoffset,uoffset,voffset 67 68 character arg*32 69 70 external setup,diff,diff_lin 71 72 call getarg(1,arg) 73 74 call ceedinit(trim(arg)//char(0),ceed,err) 75 76! DoF Coordinates 77 do i=0,nx*2 78 do j=0,ny*2 79 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 80 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 81 enddo 82 enddo 83 call ceedvectorcreate(ceed,d*ndofs,x,err) 84 xoffset=0 85 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 86 87! Qdata Vector 88 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata,err) 89 90! Element Setup 91 do i=0,nelem-1 92 col=mod(i,nx) 93 row=i/nx 94 offset=col*(p-1)+row*(nx*2+1)*(p-1) 95 do j=0,p-1 96 do k=0,p-1 97 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 98 enddo 99 enddo 100 enddo 101 102! Restrictions 103 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 104 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 105 stridesx=[1,p*p,p*p*d] 106 call ceedelemrestrictioncreatestrided(ceed,nelem,p*p,& 107 & nelem*p*p,d,stridesx,erestrictxi,err) 108 109 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 110 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 111 stridesu=[1,q*q,q*q] 112 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 113 & 1,stridesu,erestrictui,err) 114 115 stridesqd=[1,q*q,q*q*d*(d+1)/2] 116 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 117 & d*(d+1)/2,stridesqd,erestrictqi,err) 118 119! Bases 120 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 121 & bx,err) 122 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 123 & bu,err) 124 125! QFunction - setup 126 call ceedqfunctioncreateinterior(ceed,1,setup,& 127 &SOURCE_DIR& 128 &//'t531-operator.h:setup'//char(0),qf_setup,err) 129 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 130 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 131 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 132 133! Operator - setup 134 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 135 & ceed_qfunction_none,op_setup,err) 136 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 137 & bx,ceed_vector_active,err) 138 call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,& 139 & bx,ceed_vector_none,err) 140 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 141 & ceed_basis_collocated,ceed_vector_active,err) 142 143! Apply Setup Operator 144 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 145 146! QFunction - apply 147 call ceedqfunctioncreateinterior(ceed,1,diff,& 148 &SOURCE_DIR& 149 &//'t531-operator.h:diff'//char(0),qf_diff,err) 150 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 151 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 152 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 153 154! Operator - apply 155 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 156 & ceed_qfunction_none,op_diff,err) 157 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 158 & bu,ceed_vector_active,err) 159 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 160 & ceed_basis_collocated,qdata,err) 161 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 162 & bu,ceed_vector_active,err) 163 164! Assemble Diagonal 165 call ceedoperatorassemblelineardiagonal(op_diff,a,& 166 & ceed_request_immediate,err) 167 168! Manually assemble diagonal 169 call ceedvectorcreate(ceed,ndofs,u,err) 170 call ceedvectorsetvalue(u,0.d0,err) 171 call ceedvectorcreate(ceed,ndofs,v,err) 172 do i=1,ndofs 173 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 174 uu(i+uoffset)=1.d0 175 if (i>1) then 176 uu(i-1+uoffset)=0.d0 177 endif 178 call ceedvectorrestorearray(u,uu,uoffset,err) 179 180 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 181 182 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 183 atrue(i)=vv(voffset+i) 184 call ceedvectorrestorearrayread(v,vv,voffset,err) 185 enddo 186 187! Check Output 188 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 189 do i=1,ndofs 190 if (abs(aa(aoffset+i)-atrue(i))>1.0d-13) then 191! LCOV_EXCL_START 192 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 193 & atrue(i) 194! LCOV_EXCL_STOP 195 endif 196 enddo 197 call ceedvectorrestorearrayread(a,aa,aoffset,err) 198 199! Cleanup 200 call ceedqfunctiondestroy(qf_setup,err) 201 call ceedqfunctiondestroy(qf_diff,err) 202 call ceedoperatordestroy(op_setup,err) 203 call ceedoperatordestroy(op_diff,err) 204 call ceedelemrestrictiondestroy(erestrictu,err) 205 call ceedelemrestrictiondestroy(erestrictx,err) 206 call ceedelemrestrictiondestroy(erestrictxi,err) 207 call ceedelemrestrictiondestroy(erestrictui,err) 208 call ceedelemrestrictiondestroy(erestrictqi,err) 209 call ceedbasisdestroy(bu,err) 210 call ceedbasisdestroy(bx,err) 211 call ceedvectordestroy(x,err) 212 call ceedvectordestroy(a,err) 213 call ceedvectordestroy(u,err) 214 call ceedvectordestroy(v,err) 215 call ceedvectordestroy(qdata,err) 216 call ceeddestroy(ceed,err) 217 end 218!----------------------------------------------------------------------- 219