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