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 implicit none 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 ceedvectorcreate(ceed,ndofs,a,err) 149 call ceedoperatorlinearassemblediagonal(op_mass,a,& 150 & ceed_request_immediate,err) 151 152! Manually assemble diagonal 153 call ceedvectorcreate(ceed,ndofs,u,err) 154 call ceedvectorsetvalue(u,0.d0,err) 155 call ceedvectorcreate(ceed,ndofs,v,err) 156 do i=1,ndofs 157 call ceedvectorgetarray(u,ceed_mem_host,uu,uoffset,err) 158 uu(i+uoffset)=1.d0 159 if (i>1) then 160 uu(i-1+uoffset)=0.d0 161 endif 162 call ceedvectorrestorearray(u,uu,uoffset,err) 163 164 call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err) 165 166 call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err) 167 atrue(i)=vv(voffset+i) 168 call ceedvectorrestorearrayread(v,vv,voffset,err) 169 enddo 170 171! Check Output 172 call ceedvectorgetarrayread(a,ceed_mem_host,aa,aoffset,err) 173 do i=1,ndofs 174 if (abs(aa(aoffset+i)-atrue(i))>1.0d-14) then 175! LCOV_EXCL_START 176 write(*,*) '[',i,'] Error in assembly: ',aa(aoffset+i),' != ',& 177 & atrue(i) 178! LCOV_EXCL_STOP 179 endif 180 enddo 181 call ceedvectorrestorearrayread(a,aa,aoffset,err) 182 183! Cleanup 184 call ceedqfunctiondestroy(qf_setup,err) 185 call ceedqfunctiondestroy(qf_mass,err) 186 call ceedoperatordestroy(op_setup,err) 187 call ceedoperatordestroy(op_mass,err) 188 call ceedelemrestrictiondestroy(erestrictu,err) 189 call ceedelemrestrictiondestroy(erestrictx,err) 190 call ceedelemrestrictiondestroy(erestrictui,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