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