1!----------------------------------------------------------------------- 2 subroutine setup_mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,& 3& u14,u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,& 4& v16,ierr) 5 real*8 ctx 6 real*8 u1(1) 7 real*8 u2(1) 8 real*8 v1(1) 9 integer q,ierr 10 11 do i=1,q 12 v1(i)=u2(i)*(u1(i+q*0)*u1(i+q*3)-u1(i+q*1)*u1(i+q*2)) 13 enddo 14 15 ierr=0 16 end 17!----------------------------------------------------------------------- 18 subroutine apply(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,& 19& u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr) 20 real*8 ctx 21 real*8 u1(1) 22 real*8 u2(1) 23 real*8 v1(1) 24 integer q,ierr 25 26 do i=1,q 27! mass 28 v1(i) = u1(i)*u2(i) 29 enddo 30 31 ierr=0 32 end 33!----------------------------------------------------------------------- 34 program test 35 36 include 'ceedf.h' 37 38 integer ceed,err,i 39 integer stridesx(3),stridesu(3),stridesq(3) 40 integer erestrictxi,erestrictui,erestrictqi 41 integer bx,bu 42 integer qf_setup_mass,qf_apply 43 integer op_setup_mass,op_apply,op_inv 44 integer qdata_mass,x,u,v 45 integer nelem,p,q,d 46 parameter(nelem=1) 47 parameter(p=4) 48 parameter(q=5) 49 parameter(d=2) 50 integer ndofs,nqpts 51 parameter(ndofs=p*p) 52 parameter(nqpts=nelem*q*q) 53 real*8 arrx(d*nelem*2*2),uu(ndofs) 54 integer*8 xoffset,uoffset 55 56 character arg*32 57 58 external setup_mass,apply 59 60 call getarg(1,arg) 61 62 call ceedinit(trim(arg)//char(0),ceed,err) 63 64! DoF Coordinates 65 do i=0,1 66 do j=0,1 67 arrx(i+1+j*2+0*4)=i 68 arrx(i+1+j*2+1*4)=j 69 enddo 70 enddo 71 call ceedvectorcreate(ceed,d*nelem*2*2,x,err) 72 xoffset=0 73 call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err) 74 75! Qdata Vector 76 call ceedvectorcreate(ceed,nqpts,qdata_mass,err) 77 78! Restrictions 79 stridesx=[1,2*2,2*2*d] 80 call ceedelemrestrictioncreatestrided(ceed,nelem,2*2,d,d*nelem*2*2,& 81 & stridesx,erestrictxi,err) 82 83 stridesu=[1,p*p,p*p] 84 call ceedelemrestrictioncreatestrided(ceed,nelem,p*p,1,ndofs,& 85 & stridesu,erestrictui,err) 86 87 stridesq=[1,q*q,q*q] 88 call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,1,nqpts,& 89 & stridesq,erestrictqi,err) 90 91! Bases 92 call ceedbasiscreatetensorh1lagrange(ceed,d,d,2,q,ceed_gauss,bx,err) 93 call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err) 94 95! QFunction - setup mass 96 call ceedqfunctioncreateinterior(ceed,1,setup_mass,& 97 &SOURCE_DIR& 98 &//'t540-operator.h:setup_mass'//char(0),qf_setup_mass,err) 99 call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err) 100 call ceedqfunctionaddinput(qf_setup_mass,'_weight',1,ceed_eval_weight,err) 101 call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err) 102 103! Operator - setup mass 104 call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,& 105 & ceed_qfunction_none,op_setup_mass,err) 106 call ceedoperatorsetfield(op_setup_mass,'dx',erestrictxi,& 107 & bx,ceed_vector_active,err) 108 call ceedoperatorsetfield(op_setup_mass,'_weight',& 109 & ceed_elemrestriction_none,bx,ceed_vector_none,err) 110 call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictqi,& 111 & ceed_basis_collocated,ceed_vector_active,err) 112 113! Apply Setup Operators 114 call ceedoperatorapply(op_setup_mass,x,qdata_mass,& 115 & ceed_request_immediate,err) 116 117! QFunction - apply 118 call ceedqfunctioncreateinterior(ceed,1,apply,& 119 &SOURCE_DIR& 120 &//'t540-operator.h:apply'//char(0),qf_apply,err) 121 call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err) 122 call ceedqfunctionaddinput(qf_apply,'qdata_mass',1,ceed_eval_none,err) 123 call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err) 124 125! Operator - apply 126 call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,& 127 & ceed_qfunction_none,op_apply,err) 128 call ceedoperatorsetfield(op_apply,'u',erestrictui,& 129 & bu,ceed_vector_active,err) 130 call ceedoperatorsetfield(op_apply,'qdata_mass',erestrictqi,& 131 & ceed_basis_collocated,qdata_mass,err) 132 call ceedoperatorsetfield(op_apply,'v',erestrictui,& 133 & bu,ceed_vector_active,err) 134 135! Apply original operator 136 call ceedvectorcreate(ceed,ndofs,u,err) 137 call ceedvectorsetvalue(u,1.d0,err) 138 call ceedvectorcreate(ceed,ndofs,v,err) 139 call ceedvectorsetvalue(v,0.d0,err) 140 call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err) 141 142! Create FDM element inverse 143 call ceedoperatorcreatefdmelementinverse(op_apply,op_inv,& 144 & ceed_request_immediate,err) 145 146! Apply FDM element inverse 147 call ceedoperatorapply(op_inv,v,u,ceed_request_immediate,err) 148 149! Check Output 150 call ceedvectorgetarrayread(u,ceed_mem_host,uu,uoffset,err) 151 do i=1,ndofs 152 if (abs(uu(uoffset+i)-1.0)>5.0d-14) then 153! LCOV_EXCL_START 154 write(*,*) '[',i,'] Error in inverse: ',uu(uoffset+i),' != 1.0' 155! LCOV_EXCL_STOP 156 endif 157 enddo 158 call ceedvectorrestorearrayread(u,uu,uoffset,err) 159 160! Cleanup 161 call ceedqfunctiondestroy(qf_setup_mass,err) 162 call ceedqfunctiondestroy(qf_apply,err) 163 call ceedoperatordestroy(op_setup_mass,err) 164 call ceedoperatordestroy(op_apply,err) 165 call ceedoperatordestroy(op_inv,err) 166 call ceedelemrestrictiondestroy(erestrictxi,err) 167 call ceedelemrestrictiondestroy(erestrictui,err) 168 call ceedelemrestrictiondestroy(erestrictqi,err) 169 call ceedbasisdestroy(bu,err) 170 call ceedbasisdestroy(bx,err) 171 call ceedvectordestroy(x,err) 172 call ceedvectordestroy(u,err) 173 call ceedvectordestroy(v,err) 174 call ceedvectordestroy(qdata_mass,err) 175 call ceeddestroy(ceed,err) 176 end 177!----------------------------------------------------------------------- 178