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 implicit none 42 include 'ceedf.h' 43 44 integer ceed,err,i,j,k 45 integer stridesu(3),stridesqd(3) 46 integer erestrictx,erestrictu,erestrictui,erestrictqi 47 integer bx,bu 48 integer qf_setup,qf_diff 49 integer op_setup,op_diff 50 integer qdata,x,a,u,v 51 integer nelem,p,q,d 52 integer row,col,offset 53 parameter(nelem=6) 54 parameter(p=3) 55 parameter(q=4) 56 parameter(d=2) 57 integer ndofs,nqpts,nx,ny 58 parameter(nx=3) 59 parameter(ny=2) 60 parameter(ndofs=(nx*2+1)*(ny*2+1)) 61 parameter(nqpts=nelem*q*q) 62 integer indx(nelem*p*p) 63 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 64 integer*8 xoffset,aoffset,uoffset,voffset 65 66 character arg*32 67 68 external setup,diff,diff_lin 69 70 call getarg(1,arg) 71 72 call ceedinit(trim(arg)//char(0),ceed,err) 73 74! DoF Coordinates 75 do i=0,nx*2 76 do j=0,ny*2 77 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 78 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 79 enddo 80 enddo 81 call ceedvectorcreate(ceed,d*ndofs,x,err) 82 xoffset=0 83 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 84 85! Qdata Vector 86 call ceedvectorcreate(ceed,nqpts*d*(d+1)/2,qdata,err) 87 88! Element Setup 89 do i=0,nelem-1 90 col=mod(i,nx) 91 row=i/nx 92 offset=col*(p-1)+row*(nx*2+1)*(p-1) 93 do j=0,p-1 94 do k=0,p-1 95 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 96 enddo 97 enddo 98 enddo 99 100! Restrictions 101 call ceedelemrestrictioncreate(ceed,nelem,p*p,d,ndofs,d*ndofs,& 102 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 103 104 call ceedelemrestrictioncreate(ceed,nelem,p*p,1,1,ndofs,& 105 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 106 stridesu=[1,q*q,q*q] 107 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,1,nqpts,& 108 & stridesu,erestrictui,err) 109 110 stridesqd=[1,q*q,q*q*d*(d+1)/2] 111 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,d*(d+1)/2,& 112 & d*(d+1)/2*nqpts,stridesqd,erestrictqi,err) 113 114! Bases 115 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 116 & bx,err) 117 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 118 & bu,err) 119 120! QFunction - setup 121 call ceedqfunctioncreateinterior(ceed,1,setup,& 122 &SOURCE_DIR& 123 &//'t531-operator.h:setup'//char(0),qf_setup,err) 124 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 125 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 126 call ceedqfunctionaddoutput(qf_setup,'qdata',d*(d+1)/2,ceed_eval_none,err) 127 128! Operator - setup 129 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 130 & ceed_qfunction_none,op_setup,err) 131 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 132 & bx,ceed_vector_active,err) 133 call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,& 134 & bx,ceed_vector_none,err) 135 call ceedoperatorsetfield(op_setup,'qdata',erestrictqi,& 136 & ceed_basis_collocated,ceed_vector_active,err) 137 138! Apply Setup Operator 139 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 140 141! QFunction - apply 142 call ceedqfunctioncreateinterior(ceed,1,diff,& 143 &SOURCE_DIR& 144 &//'t531-operator.h:diff'//char(0),qf_diff,err) 145 call ceedqfunctionaddinput(qf_diff,'du',d,ceed_eval_grad,err) 146 call ceedqfunctionaddinput(qf_diff,'qdata',d*(d+1)/2,ceed_eval_none,err) 147 call ceedqfunctionaddoutput(qf_diff,'dv',d,ceed_eval_grad,err) 148 149! Operator - apply 150 call ceedoperatorcreate(ceed,qf_diff,ceed_qfunction_none,& 151 & ceed_qfunction_none,op_diff,err) 152 call ceedoperatorsetfield(op_diff,'du',erestrictu,& 153 & bu,ceed_vector_active,err) 154 call ceedoperatorsetfield(op_diff,'qdata',erestrictqi,& 155 & ceed_basis_collocated,qdata,err) 156 call ceedoperatorsetfield(op_diff,'dv',erestrictu,& 157 & bu,ceed_vector_active,err) 158 159! Assemble Diagonal 160 call ceedvectorcreate(ceed,ndofs,a,err) 161 call ceedoperatorlinearassemblediagonal(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(erestrictui,err) 203 call ceedelemrestrictiondestroy(erestrictqi,err) 204 call ceedbasisdestroy(bu,err) 205 call ceedbasisdestroy(bx,err) 206 call ceedvectordestroy(x,err) 207 call ceedvectordestroy(a,err) 208 call ceedvectordestroy(u,err) 209 call ceedvectordestroy(v,err) 210 call ceedvectordestroy(qdata,err) 211 call ceeddestroy(ceed,err) 212 end 213!----------------------------------------------------------------------- 214