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