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 integer q,ierr 9 10 do i=1,q 11 v1(i)=u1(i)*(u2(i+q*0)*u2(i+q*3)-u2(i+q*1)*u2(i+q*2)) 12 enddo 13 14 ierr=0 15 end 16!----------------------------------------------------------------------- 17 subroutine mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 18& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 19 real*8 ctx 20 real*8 u1(1) 21 real*8 u2(1) 22 real*8 v1(1) 23 integer q,ierr 24 25 do i=1,q 26 v1(i)=u2(i)*u1(i) 27 enddo 28 29 ierr=0 30 end 31!----------------------------------------------------------------------- 32 program test 33 34 include 'ceedf.h' 35 36 integer ceed,err,i,j,k 37 integer imode 38 parameter(imode=ceed_noninterlaced) 39 integer stridesx(3),stridesu(3) 40 integer erestrictx,erestrictu,erestrictxi,erestrictui 41 integer bx,bu 42 integer qf_setup,qf_mass 43 integer op_setup,op_mass 44 integer qdata,x,a,u,v 45 integer nelem,p,q,d 46 integer row,col,offset 47 parameter(nelem=6) 48 parameter(p=3) 49 parameter(q=4) 50 parameter(d=2) 51 integer ndofs,nqpts,nx,ny 52 parameter(nx=3) 53 parameter(ny=2) 54 parameter(ndofs=(nx*2+1)*(ny*2+1)) 55 parameter(nqpts=nelem*q*q) 56 integer indx(nelem*p*p) 57 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 58 integer*8 xoffset,aoffset,uoffset,voffset 59 60 character arg*32 61 62 external setup,mass 63 64 call getarg(1,arg) 65 66 call ceedinit(trim(arg)//char(0),ceed,err) 67 68! DoF Coordinates 69 do i=0,nx*2 70 do j=0,ny*2 71 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 72 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 73 enddo 74 enddo 75 call ceedvectorcreate(ceed,d*ndofs,x,err) 76 xoffset=0 77 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 78 79! Qdata Vector 80 call ceedvectorcreate(ceed,nqpts,qdata,err) 81 82! Element Setup 83 do i=0,nelem-1 84 col=mod(i,nx) 85 row=i/nx 86 offset=col*(p-1)+row*(nx*2+1)*(p-1) 87 do j=0,p-1 88 do k=0,p-1 89 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 90 enddo 91 enddo 92 enddo 93 94! Restrictions 95 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 96 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 97 stridesx=[1,p*p,p*p*d] 98 call ceedelemrestrictioncreatestrided(ceed,nelem,p*p,& 99 & nelem*p*p,d,stridesx,erestrictxi,err) 100 101 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 102 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 103 stridesu=[1,q*q,q*q] 104 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,nqpts,& 105 & 1,stridesu,erestrictui,err) 106 107! Bases 108 call ceedbasiscreatetensorh1lagrange(ceed,d,d,p,q,ceed_gauss,& 109 & bx,err) 110 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,& 111 & bu,err) 112 113! QFunctions 114! -- Setup 115 call ceedqfunctioncreateinterior(ceed,1,setup,& 116 &SOURCE_DIR& 117 &//'t510-operator.h:setup'//char(0),qf_setup,err) 118 call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err) 119 call ceedqfunctionaddinput(qf_setup,'dx',d*d,ceed_eval_grad,err) 120 call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err) 121! -- Mass 122 call ceedqfunctioncreateinterior(ceed,1,mass,& 123 &SOURCE_DIR& 124 &//'t510-operator.h:mass'//char(0),qf_mass,err) 125 call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err) 126 call ceedqfunctionaddinput(qf_mass,'u',1,ceed_eval_interp,err) 127 call ceedqfunctionaddoutput(qf_mass,'v',1,ceed_eval_interp,err) 128 129! Operators 130! -- Setup 131 call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,& 132 & ceed_qfunction_none,op_setup,err) 133 call ceedoperatorsetfield(op_setup,'_weight',erestrictxi,& 134 & bx,ceed_vector_none,err) 135 call ceedoperatorsetfield(op_setup,'dx',erestrictx,& 136 & bx,ceed_vector_active,err) 137 call ceedoperatorsetfield(op_setup,'rho',erestrictui,& 138 & ceed_basis_collocated,ceed_vector_active,err) 139! -- Mass 140 call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,& 141 & ceed_qfunction_none,op_mass,err) 142 call ceedoperatorsetfield(op_mass,'rho',erestrictui,& 143 & ceed_basis_collocated,qdata,err) 144 call ceedoperatorsetfield(op_mass,'u',erestrictu,& 145 & bu,ceed_vector_active,err) 146 call ceedoperatorsetfield(op_mass,'v',erestrictu,& 147 & bu,ceed_vector_active,err) 148 149! Apply Setup Operator 150 call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err) 151 152! Assemble Diagonal 153 call ceedoperatorassemblelineardiagonal(op_mass,a,& 154 & ceed_request_immediate,err) 155 156! Manually assemble diagonal 157 call ceedvectorcreate(ceed,ndofs,u,err) 158 call ceedvectorsetvalue(u,0.d0,err) 159 call ceedvectorcreate(ceed,ndofs,v,err) 160 do i=1,ndofs 161 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 162 uu(i+uoffset)=1.d0 163 if (i>1) then 164 uu(i-1+uoffset)=0.d0 165 endif 166 call ceedvectorrestorearray(u,uu,uoffset,err) 167 168 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 169 170 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 171 atrue(i)=vv(voffset+i) 172 call ceedvectorrestorearrayread(v,vv,voffset,err) 173 enddo 174 175! Check Output 176 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 177 do i=1,ndofs 178 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 179! LCOV_EXCL_START 180 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 181 & atrue(i) 182! LCOV_EXCL_STOP 183 endif 184 enddo 185 call ceedvectorrestorearrayread(a,aa,aoffset,err) 186 187! Cleanup 188 call ceedqfunctiondestroy(qf_setup,err) 189 call ceedqfunctiondestroy(qf_mass,err) 190 call ceedoperatordestroy(op_setup,err) 191 call ceedoperatordestroy(op_mass,err) 192 call ceedelemrestrictiondestroy(erestrictu,err) 193 call ceedelemrestrictiondestroy(erestrictx,err) 194 call ceedelemrestrictiondestroy(erestrictui,err) 195 call ceedelemrestrictiondestroy(erestrictxi,err) 196 call ceedbasisdestroy(bu,err) 197 call ceedbasisdestroy(bx,err) 198 call ceedvectordestroy(x,err) 199 call ceedvectordestroy(a,err) 200 call ceedvectordestroy(u,err) 201 call ceedvectordestroy(v,err) 202 call ceedvectordestroy(qdata,err) 203 call ceeddestroy(ceed,err) 204 end 205!----------------------------------------------------------------------- 206