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 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 ceedoperatorassemblelineardiagonal(op_diff,a,& 161 & ceed_request_immediate,err) 162 163! Manually assemble diagonal 164 call ceedvectorcreate(ceed,ndofs,u,err) 165 call ceedvectorsetvalue(u,0.d0,err) 166 call ceedvectorcreate(ceed,ndofs,v,err) 167 do i=1,ndofs 168 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 169 uu(i+uoffset)=1.d0 170 if (i>1) then 171 uu(i-1+uoffset)=0.d0 172 endif 173 call ceedvectorrestorearray(u,uu,uoffset,err) 174 175 call ceedoperatorapply(op_diff,u,v,ceed_request_immediate,err) 176 177 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 178 atrue(i)=vv(voffset+i) 179 call ceedvectorrestorearrayread(v,vv,voffset,err) 180 enddo 181 182! Check Output 183 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 184 do i=1,ndofs 185 if (abs(aa(aoffset+i)-atrue(i))>1.0d-13) then 186! LCOV_EXCL_START 187 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 188 & atrue(i) 189! LCOV_EXCL_STOP 190 endif 191 enddo 192 call ceedvectorrestorearrayread(a,aa,aoffset,err) 193 194! Cleanup 195 call ceedqfunctiondestroy(qf_setup,err) 196 call ceedqfunctiondestroy(qf_diff,err) 197 call ceedoperatordestroy(op_setup,err) 198 call ceedoperatordestroy(op_diff,err) 199 call ceedelemrestrictiondestroy(erestrictu,err) 200 call ceedelemrestrictiondestroy(erestrictx,err) 201 call ceedelemrestrictiondestroy(erestrictui,err) 202 call ceedelemrestrictiondestroy(erestrictqi,err) 203 call ceedbasisdestroy(bu,err) 204 call ceedbasisdestroy(bx,err) 205 call ceedvectordestroy(x,err) 206 call ceedvectordestroy(a,err) 207 call ceedvectordestroy(u,err) 208 call ceedvectordestroy(v,err) 209 call ceedvectordestroy(qdata,err) 210 call ceeddestroy(ceed,err) 211 end 212!----------------------------------------------------------------------- 213