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 erestrictx,erestrictu,erestrictxi,erestrictui 40 integer bx,bu 41 integer qf_setup,qf_mass 42 integer op_setup,op_mass 43 integer qdata,x,a,u,v 44 integer nelem,p,q,d 45 integer row,col,offset 46 parameter(nelem=6) 47 parameter(p=3) 48 parameter(q=4) 49 parameter(d=2) 50 integer ndofs,nqpts,nx,ny 51 parameter(nx=3) 52 parameter(ny=2) 53 parameter(ndofs=(nx*2+1)*(ny*2+1)) 54 parameter(nqpts=nelem*q*q) 55 integer indx(nelem*p*p) 56 real*8 arrx(d*ndofs),aa(nqpts),uu(ndofs),vv(ndofs),atrue(ndofs) 57 integer*8 xoffset,aoffset,uoffset,voffset 58 59 character arg*32 60 61 external setup,mass 62 63 call getarg(1,arg) 64 65 call ceedinit(trim(arg)//char(0),ceed,err) 66 67! DoF Coordinates 68 do i=0,nx*2 69 do j=0,ny*2 70 arrx(i+j*(nx*2+1)+0*ndofs+1)=1.d0*i/(2*nx) 71 arrx(i+j*(nx*2+1)+1*ndofs+1)=1.d0*j/(2*ny) 72 enddo 73 enddo 74 call ceedvectorcreate(ceed,d*ndofs,x,err) 75 xoffset=0 76 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 77 78! Qdata Vector 79 call ceedvectorcreate(ceed,nqpts,qdata,err) 80 81! Element Setup 82 do i=0,nelem-1 83 col=mod(i,nx) 84 row=i/nx 85 offset=col*(p-1)+row*(nx*2+1)*(p-1) 86 do j=0,p-1 87 do k=0,p-1 88 indx(p*(p*i+k)+j+1)=offset+k*(nx*2+1)+j 89 enddo 90 enddo 91 enddo 92 93! Restrictions 94 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,d,& 95 & ceed_mem_host,ceed_use_pointer,indx,erestrictx,err) 96 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,p*p,& 97 & nelem*p*p,d,erestrictxi,err) 98 99 call ceedelemrestrictioncreate(ceed,imode,nelem,p*p,ndofs,1,& 100 & ceed_mem_host,ceed_use_pointer,indx,erestrictu,err) 101 call ceedelemrestrictioncreateidentity(ceed,imode,nelem,q*q,nqpts,& 102 & 1,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',erestrictxi,& 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 ceedelemrestrictiondestroy(erestrictxi,err) 193 call ceedbasisdestroy(bu,err) 194 call ceedbasisdestroy(bx,err) 195 call ceedvectordestroy(x,err) 196 call ceedvectordestroy(a,err) 197 call ceedvectordestroy(u,err) 198 call ceedvectordestroy(v,err) 199 call ceedvectordestroy(qdata,err) 200 call ceeddestroy(ceed,err) 201 end 202!----------------------------------------------------------------------- 203