xref: /libCEED/tests/t540-operator-f.f90 (revision d9b786505a4dfcb66b2fcd9e3b61dd507168515d)
1!-----------------------------------------------------------------------
2!
3! Header with QFunctions
4!
5      include 't540-operator-f.h'
6!-----------------------------------------------------------------------
7      program test
8      implicit none
9      include 'ceed/fortran.h'
10
11      integer ceed,err,i,j
12      integer stridesx(3),stridesu(3),stridesq(3)
13      integer erestrictxi,erestrictui,erestrictqi
14      integer bx,bu
15      integer qf_setup_mass,qf_apply
16      integer op_setup_mass,op_apply,op_inv
17      integer qdata_mass,x,u,v
18      integer nelem,p,q,d
19      parameter(nelem=1)
20      parameter(p=4)
21      parameter(q=5)
22      parameter(d=2)
23      integer ndofs,nqpts
24      parameter(ndofs=p*p)
25      parameter(nqpts=nelem*q*q)
26      real*8 arrx(d*nelem*2*2),uu(ndofs)
27      integer*8 xoffset,uoffset
28
29      character arg*32
30
31      external setup_mass,apply
32
33      call getarg(1,arg)
34
35      call ceedinit(trim(arg)//char(0),ceed,err)
36
37! DoF Coordinates
38      do i=0,1
39        do j=0,1
40          arrx(i+1+j*2+0*4)=i
41          arrx(i+1+j*2+1*4)=j
42        enddo
43      enddo
44      call ceedvectorcreate(ceed,d*nelem*2*2,x,err)
45      xoffset=0
46      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err)
47
48! Qdata Vector
49      call ceedvectorcreate(ceed,nqpts,qdata_mass,err)
50
51! Restrictions
52      stridesx=[1,2*2,2*2*d]
53      call ceedelemrestrictioncreatestrided(ceed,nelem,2*2,d,d*nelem*2*2,&
54     & stridesx,erestrictxi,err)
55
56      stridesu=[1,p*p,p*p]
57      call ceedelemrestrictioncreatestrided(ceed,nelem,p*p,1,ndofs,&
58     & stridesu,erestrictui,err)
59
60      stridesq=[1,q*q,q*q]
61      call ceedelemrestrictioncreatestrided(ceed,nelem,q*q,1,nqpts,&
62     & stridesq,erestrictqi,err)
63
64! Bases
65      call ceedbasiscreatetensorh1lagrange(ceed,d,d,2,q,ceed_gauss,bx,err)
66      call ceedbasiscreatetensorh1lagrange(ceed,d,1,p,q,ceed_gauss,bu,err)
67
68! QFunction - setup mass
69      call ceedqfunctioncreateinterior(ceed,1,setup_mass,&
70     &SOURCE_DIR&
71     &//'t540-operator.h:setup_mass'//char(0),qf_setup_mass,err)
72      call ceedqfunctionaddinput(qf_setup_mass,'dx',d*d,ceed_eval_grad,err)
73      call ceedqfunctionaddinput(qf_setup_mass,'weight',1,ceed_eval_weight,err)
74      call ceedqfunctionaddoutput(qf_setup_mass,'qdata',1,ceed_eval_none,err)
75
76! Operator - setup mass
77      call ceedoperatorcreate(ceed,qf_setup_mass,ceed_qfunction_none,&
78     & ceed_qfunction_none,op_setup_mass,err)
79      call ceedoperatorsetfield(op_setup_mass,'dx',erestrictxi,&
80     & bx,ceed_vector_active,err)
81      call ceedoperatorsetfield(op_setup_mass,'weight',&
82     & ceed_elemrestriction_none,bx,ceed_vector_none,err)
83      call ceedoperatorsetfield(op_setup_mass,'qdata',erestrictqi,&
84     ceed_basis_none,ceed_vector_active,err)
85
86! Apply Setup Operators
87      call ceedoperatorapply(op_setup_mass,x,qdata_mass,&
88     & ceed_request_immediate,err)
89
90! QFunction - apply
91      call ceedqfunctioncreateinterior(ceed,1,apply,&
92     &SOURCE_DIR&
93     &//'t540-operator.h:apply'//char(0),qf_apply,err)
94      call ceedqfunctionaddinput(qf_apply,'u',1,ceed_eval_interp,err)
95      call ceedqfunctionaddinput(qf_apply,'mass qdata',1,ceed_eval_none,err)
96      call ceedqfunctionaddoutput(qf_apply,'v',1,ceed_eval_interp,err)
97
98! Operator - apply
99      call ceedoperatorcreate(ceed,qf_apply,ceed_qfunction_none,&
100     & ceed_qfunction_none,op_apply,err)
101      call ceedoperatorsetfield(op_apply,'u',erestrictui,&
102     & bu,ceed_vector_active,err)
103      call ceedoperatorsetfield(op_apply,'mass qdata',erestrictqi,&
104     ceed_basis_none,qdata_mass,err)
105      call ceedoperatorsetfield(op_apply,'v',erestrictui,&
106     & bu,ceed_vector_active,err)
107
108! Apply original operator
109      call ceedvectorcreate(ceed,ndofs,u,err)
110      call ceedvectorsetvalue(u,1.d0,err)
111      call ceedvectorcreate(ceed,ndofs,v,err)
112      call ceedvectorsetvalue(v,0.d0,err)
113      call ceedoperatorapply(op_apply,u,v,ceed_request_immediate,err)
114
115! Create FDM element inverse
116      call ceedoperatorcreatefdmelementinverse(op_apply,op_inv,&
117     & ceed_request_immediate,err)
118
119! Apply FDM element inverse
120      call ceedoperatorapply(op_inv,v,u,ceed_request_immediate,err)
121
122! Check Output
123      call ceedvectorgetarrayread(u,ceed_mem_host,uu,uoffset,err)
124      do i=1,ndofs
125        if (abs(uu(uoffset+i)-1.0)>5.0d-14) then
126! LCOV_EXCL_START
127          write(*,*) '[',i,'] Error in inverse: ',uu(uoffset+i),' != 1.0'
128! LCOV_EXCL_STOP
129        endif
130      enddo
131      call ceedvectorrestorearrayread(u,uu,uoffset,err)
132
133! Cleanup
134      call ceedqfunctiondestroy(qf_setup_mass,err)
135      call ceedqfunctiondestroy(qf_apply,err)
136      call ceedoperatordestroy(op_setup_mass,err)
137      call ceedoperatordestroy(op_apply,err)
138      call ceedoperatordestroy(op_inv,err)
139      call ceedelemrestrictiondestroy(erestrictxi,err)
140      call ceedelemrestrictiondestroy(erestrictui,err)
141      call ceedelemrestrictiondestroy(erestrictqi,err)
142      call ceedbasisdestroy(bu,err)
143      call ceedbasisdestroy(bx,err)
144      call ceedvectordestroy(x,err)
145      call ceedvectordestroy(u,err)
146      call ceedvectordestroy(v,err)
147      call ceedvectordestroy(qdata_mass,err)
148      call ceeddestroy(ceed,err)
149      end
150!-----------------------------------------------------------------------
151