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