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 stridesu(3),stridesqd(3) 48 integer erestrictx,erestrictu,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 106 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 107 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 108 stridesu=[1,q*q,q*q] 109 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 110 & 1,stridesu,erestrictui,err) 111 112 stridesqd=[1,q*q,q*q*d*(d+1)/2] 113 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 114 & d*(d+1)/2,stridesqd,erestrictqi,err) 115 116! Bases 117 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 118 & bx,err) 119 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 120 & bu,err) 121 122! QFunction - setup 123 call ceedqfunctioncreateinterior(ceed,1,setup,& 124 &SOURCE_DIR& 125 &//'t531-operator.h:setup'//char(0),qf_setup,err) 126 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 127 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 128 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 129 130! Operator - setup 131 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 132 & ceed_qfunction_none,op_setup,err) 133 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 134 & bx,ceed_vector_active,err) 135 call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,& 136 & bx,ceed_vector_none,err) 137 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 138 & ceed_basis_collocated,ceed_vector_active,err) 139 140! Apply Setup Operator 141 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 142 143! QFunction - apply 144 call ceedqfunctioncreateinterior(ceed,1,diff,& 145 &SOURCE_DIR& 146 &//'t531-operator.h:diff'//char(0),qf_diff,err) 147 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 148 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 149 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 150 151! Operator - apply 152 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 153 & ceed_qfunction_none,op_diff,err) 154 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 155 & bu,ceed_vector_active,err) 156 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 157 & ceed_basis_collocated,qdata,err) 158 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 159 & bu,ceed_vector_active,err) 160 161! Assemble Diagonal 162 call ceedoperatorassemblelineardiagonal(op_diff,a,& 163 & ceed_request_immediate,err) 164 165! Manually assemble diagonal 166 call ceedvectorcreate(ceed,ndofs,u,err) 167 call ceedvectorsetvalue(u,0.d0,err) 168 call ceedvectorcreate(ceed,ndofs,v,err) 169 do i=1,ndofs 170 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 171 uu(i+uoffset)=1.d0 172 if (i>1) then 173 uu(i-1+uoffset)=0.d0 174 endif 175 call ceedvectorrestorearray(u,uu,uoffset,err) 176 177 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 178 179 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 180 atrue(i)=vv(voffset+i) 181 call ceedvectorrestorearrayread(v,vv,voffset,err) 182 enddo 183 184! Check Output 185 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 186 do i=1,ndofs 187 if (abs(aa(aoffset+i)-atrue(i))>1.0d-13) then 188! LCOV_EXCL_START 189 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 190 & atrue(i) 191! LCOV_EXCL_STOP 192 endif 193 enddo 194 call ceedvectorrestorearrayread(a,aa,aoffset,err) 195 196! Cleanup 197 call ceedqfunctiondestroy(qf_setup,err) 198 call ceedqfunctiondestroy(qf_diff,err) 199 call ceedoperatordestroy(op_setup,err) 200 call ceedoperatordestroy(op_diff,err) 201 call ceedelemrestrictiondestroy(erestrictu,err) 202 call ceedelemrestrictiondestroy(erestrictx,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