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 erestrictx,erestrictu,erestrictxi,erestrictui,erestrictqi 48 integer bx,bu 49 integer qf_setup,qf_diff 50 integer op_setup,op_diff 51 integer qdata,x,a,u,v 52 integer nelem,p,q,d 53 integer row,col,offset 54 parameter(nelem=6) 55 parameter(p=3) 56 parameter(q=4) 57 parameter(d=2) 58 integer ndofs,nqpts,nx,ny 59 parameter(nx=3) 60 parameter(ny=2) 61 parameter(ndofs=(nx*2+1)*(ny*2+1)) 62 parameter(nqpts=nelem*q*q) 63 integer indx(nelem*p*p) 64 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 65 integer*8 xoffset,aoffset,uoffset,voffset 66 67 character arg*32 68 69 external setup,diff,diff_lin 70 71 call getarg(1,arg) 72 73 call ceedinit(trim(arg)//char(0),ceed,err) 74 75! DoF Coordinates 76 do i=0,nx*2 77 do j=0,ny*2 78 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 79 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 80 enddo 81 enddo 82 call ceedvectorcreate(ceed,d*ndofs,x,err) 83 xoffset=0 84 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 85 86! Qdata Vector 87 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata,err) 88 89! Element Setup 90 do i=0,nelem-1 91 col=mod(i,nx) 92 row=i/nx 93 offset=col*(p-1)+row*(nx*2+1)*(p-1) 94 do j=0,p-1 95 do k=0,p-1 96 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 97 enddo 98 enddo 99 enddo 100 101! Restrictions 102 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 103 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 104 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,p*p,& 105 & nelem*p*p,d,erestrictxi,err) 106 107 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 108 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 109 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 110 & 1,erestrictui,err) 111 112 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 113 & d*(d+1)/2,erestrictqi,err) 114 115! Bases 116 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 117 & bx,err) 118 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 119 & bu,err) 120 121! QFunction - setup 122 call ceedqfunctioncreateinterior(ceed,1,setup,& 123 &SOURCE_DIR& 124 &//'t531-operator.h:setup'//char(0),qf_setup,err) 125 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 126 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 127 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 128 129! Operator - setup 130 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 131 & ceed_qfunction_none,op_setup,err) 132 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 133 & bx,ceed_vector_active,err) 134 call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,& 135 & bx,ceed_vector_none,err) 136 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 137 & ceed_basis_collocated,ceed_vector_active,err) 138 139! Apply Setup Operator 140 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 141 142! QFunction - apply 143 call ceedqfunctioncreateinterior(ceed,1,diff,& 144 &SOURCE_DIR& 145 &//'t531-operator.h:diff'//char(0),qf_diff,err) 146 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 147 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 148 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 149 150! Operator - apply 151 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 152 & ceed_qfunction_none,op_diff,err) 153 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 154 & bu,ceed_vector_active,err) 155 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 156 & ceed_basis_collocated,qdata,err) 157 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 158 & bu,ceed_vector_active,err) 159 160! Assemble Diagonal 161 call ceedoperatorassemblelineardiagonal(op_diff,a,& 162 & ceed_request_immediate,err) 163 164! Manually assemble diagonal 165 call ceedvectorcreate(ceed,ndofs,u,err) 166 call ceedvectorsetvalue(u,0.d0,err) 167 call ceedvectorcreate(ceed,ndofs,v,err) 168 do i=1,ndofs 169 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 170 uu(i+uoffset)=1.d0 171 if (i>1) then 172 uu(i-1+uoffset)=0.d0 173 endif 174 call ceedvectorrestorearray(u,uu,uoffset,err) 175 176 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 177 178 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 179 atrue(i)=vv(voffset+i) 180 call ceedvectorrestorearrayread(v,vv,voffset,err) 181 enddo 182 183! Check Output 184 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 185 do i=1,ndofs 186 if (abs(aa(aoffset+i)-atrue(i))>1.0d-13) then 187! LCOV_EXCL_START 188 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 189 & atrue(i) 190! LCOV_EXCL_STOP 191 endif 192 enddo 193 call ceedvectorrestorearrayread(a,aa,aoffset,err) 194 195! Cleanup 196 call ceedqfunctiondestroy(qf_setup,err) 197 call ceedqfunctiondestroy(qf_diff,err) 198 call ceedoperatordestroy(op_setup,err) 199 call ceedoperatordestroy(op_diff,err) 200 call ceedelemrestrictiondestroy(erestrictu,err) 201 call ceedelemrestrictiondestroy(erestrictx,err) 202 call ceedelemrestrictiondestroy(erestrictxi,err) 203 call ceedelemrestrictiondestroy(erestrictui,err) 204 call ceedelemrestrictiondestroy(erestrictqi,err) 205 call ceedbasisdestroy(bu,err) 206 call ceedbasisdestroy(bx,err) 207 call ceedvectordestroy(x,err) 208 call ceedvectordestroy(a,err) 209 call ceedvectordestroy(u,err) 210 call ceedvectordestroy(v,err) 211 call ceedvectordestroy(qdata,err) 212 call ceeddestroy(ceed,err) 213 end 214!----------------------------------------------------------------------- 215