xref: /libCEED/tests/t523-operator-f.f90 (revision ded9b81dfe1c5b0a720aeaeab78a4be02eef9bb5)
1!-----------------------------------------------------------------------
2!
3! Header with common subroutine
4!
5      include 't320-basis-f.h'
6!
7! Header with QFunctions
8!
9      include 't510-operator-f.h'
10!-----------------------------------------------------------------------
11      program test
12      implicit none
13      include 'ceedf.h'
14
15      integer ceed,err,i,j,k
16      integer stridesutet(3),stridesuhex(3)
17      integer erestrictxtet,erestrictutet,erestrictuitet,&
18&             erestrictxhex,erestrictuhex,erestrictuihex
19      integer bxtet,butet,bxhex,buhex
20      integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex
21      integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass
22      integer qdatatet,qdatahex,x
23      integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d
24      integer row,col,offset
25      parameter(nelemtet=6)
26      parameter(ptet=6)
27      parameter(qtet=4)
28      parameter(nelemhex=6)
29      parameter(phex=3)
30      parameter(qhex=4)
31      parameter(d=2)
32      integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex
33      parameter(nx=3)
34      parameter(ny=3)
35      parameter(nxtet=3)
36      parameter(nytet=1)
37      parameter(nxhex=3)
38      parameter(ndofs=(nx*2+1)*(ny*2+1))
39      parameter(nqptstet=nelemtet*qtet)
40      parameter(nqptshex=nelemhex*qhex*qhex)
41      parameter(nqpts=nqptstet+nqptshex)
42      integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex)
43
44      real*8 qref(d*qtet)
45      real*8 qweight(qtet)
46      real*8 interp(ptet*qtet)
47      real*8 grad(d*ptet*qtet)
48
49      character arg*32
50
51! LCOV_EXCL_START
52      external setup,mass
53! LCOV_EXCL_STOP
54
55      call getarg(1,arg)
56
57      call ceedinit(trim(arg)//char(0),ceed,err)
58
59! DoF Coordinates
60      call ceedvectorcreate(ceed,d*ndofs,x,err)
61
62! Qdata Vectors
63      call ceedvectorcreate(ceed,nqptstet,qdatatet,err)
64      call ceedvectorcreate(ceed,nqptshex,qdatahex,err)
65
66! Tet Elements
67      do i=0,2
68        col=mod(i,nx)
69        row=i/nx
70        offset=col*2+row*(nx*2+1)*2
71
72        indxtet(i*2*ptet+1)=2+offset
73        indxtet(i*2*ptet+2)=9+offset
74        indxtet(i*2*ptet+3)=16+offset
75        indxtet(i*2*ptet+4)=1+offset
76        indxtet(i*2*ptet+5)=8+offset
77        indxtet(i*2*ptet+6)=0+offset
78
79        indxtet(i*2*ptet+7)=14+offset
80        indxtet(i*2*ptet+8)=7+offset
81        indxtet(i*2*ptet+9)=0+offset
82        indxtet(i*2*ptet+10)=15+offset
83        indxtet(i*2*ptet+11)=8+offset
84        indxtet(i*2*ptet+12)=16+offset
85      enddo
86
87! -- Restrictions
88      call ceedelemrestrictioncreate(ceed,nelemtet,ptet,d,ndofs,d*ndofs,&
89     & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err)
90
91      call ceedelemrestrictioncreate(ceed,nelemtet,ptet,1,1,ndofs,&
92     & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err)
93      stridesutet=[1,qtet,qtet]
94      call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,1,nqptstet,&
95     & stridesutet,erestrictuitet,err)
96
97! -- Bases
98      call buildmats(qref,qweight,interp,grad)
99      call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,&
100     & qweight,bxtet,err)
101      call buildmats(qref,qweight,interp,grad)
102      call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,&
103     & qweight,butet,err)
104
105! -- QFunctions
106      call ceedqfunctioncreateinterior(ceed,1,setup,&
107     &SOURCE_DIR&
108     &//'t510-operator.h:setup'//char(0),qf_setuptet,err)
109      call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err)
110      call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err)
111      call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err)
112
113      call ceedqfunctioncreateinterior(ceed,1,mass,&
114     &SOURCE_DIR&
115     &//'t510-operator.h:mass'//char(0),qf_masstet,err)
116      call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err)
117      call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err)
118      call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err)
119
120! -- Operators
121! ---- Setup Tet
122      call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,&
123     & ceed_qfunction_none,op_setuptet,err)
124      call ceedoperatorsetfield(op_setuptet,'_weight',&
125     & ceed_elemrestriction_none,bxtet,ceed_vector_none,err)
126      call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,&
127     & bxtet,ceed_vector_active,err)
128      call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,&
129     & ceed_basis_collocated,qdatatet,err)
130! ---- Mass Tet
131      call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,&
132     & ceed_qfunction_none,op_masstet,err)
133      call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,&
134     & ceed_basis_collocated,qdatatet,err)
135      call ceedoperatorsetfield(op_masstet,'u',erestrictutet,&
136     & butet,ceed_vector_active,err)
137      call ceedoperatorsetfield(op_masstet,'v',erestrictutet,&
138     & butet,ceed_vector_active,err)
139
140! Hex Elements
141      do i=0,nelemhex-1
142        col=mod(i,nx)
143        row=i/nx
144        offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2
145        do j=0,phex-1
146          do k=0,phex-1
147            indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j
148          enddo
149        enddo
150      enddo
151
152! -- Restrictions
153      call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,d*ndofs,&
154     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err)
155
156      call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,ndofs,&
157     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err)
158      stridesuhex=[1,qhex*qhex,qhex*qhex]
159      call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,&
160     & 1,nqptshex,stridesuhex,erestrictuihex,err)
161
162! -- Bases
163      call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,&
164     & bxhex,err)
165      call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,&
166     & buhex,err)
167
168! -- QFunctions
169      call ceedqfunctioncreateinterior(ceed,1,setup,&
170     &SOURCE_DIR&
171     &//'t510-operator.h:setup'//char(0),qf_setuphex,err)
172      call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err)
173      call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err)
174      call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err)
175
176      call ceedqfunctioncreateinterior(ceed,1,mass,&
177     &SOURCE_DIR&
178     &//'t510-operator.h:mass'//char(0),qf_masshex,err)
179      call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err)
180      call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err)
181      call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err)
182
183! -- Operators
184! ---- Setup Hex
185      call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,&
186     & ceed_qfunction_none,op_setuphex,err)
187      call ceedoperatorsetfield(op_setuphex,'_weight',&
188     & ceed_elemrestriction_none,bxhex,ceed_vector_none,err)
189      call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,&
190     & bxhex,ceed_vector_active,err)
191      call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,&
192     & ceed_basis_collocated,qdatahex,err)
193! ---- Mass Hex
194      call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,&
195     & ceed_qfunction_none,op_masshex,err)
196      call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,&
197     & ceed_basis_collocated,qdatahex,err)
198      call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,&
199     & buhex,ceed_vector_active,err)
200      call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,&
201     & buhex,ceed_vector_active,err)
202
203! Composite Operators
204      call ceedcompositeoperatorcreate(ceed,op_setup,err)
205      call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err)
206      call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err)
207
208      call ceedcompositeoperatorcreate(ceed,op_mass,err)
209      call ceedcompositeoperatoraddsub(op_mass,op_masstet,err)
210      call ceedcompositeoperatoraddsub(op_mass,op_masshex,err)
211
212! View
213      call ceedoperatorview(op_setup,err)
214      call ceedoperatorview(op_mass,err)
215
216! Cleanup
217      call ceedqfunctiondestroy(qf_setuptet,err)
218      call ceedqfunctiondestroy(qf_masstet,err)
219      call ceedoperatordestroy(op_setuptet,err)
220      call ceedoperatordestroy(op_masstet,err)
221      call ceedqfunctiondestroy(qf_setuphex,err)
222      call ceedqfunctiondestroy(qf_masshex,err)
223      call ceedoperatordestroy(op_setuphex,err)
224      call ceedoperatordestroy(op_masshex,err)
225      call ceedoperatordestroy(op_setup,err)
226      call ceedoperatordestroy(op_mass,err)
227      call ceedelemrestrictiondestroy(erestrictutet,err)
228      call ceedelemrestrictiondestroy(erestrictxtet,err)
229      call ceedelemrestrictiondestroy(erestrictuitet,err)
230      call ceedelemrestrictiondestroy(erestrictuhex,err)
231      call ceedelemrestrictiondestroy(erestrictxhex,err)
232      call ceedelemrestrictiondestroy(erestrictuihex,err)
233      call ceedbasisdestroy(butet,err)
234      call ceedbasisdestroy(bxtet,err)
235      call ceedbasisdestroy(buhex,err)
236      call ceedbasisdestroy(bxhex,err)
237      call ceedvectordestroy(x,err)
238      call ceedvectordestroy(qdatatet,err)
239      call ceedvectordestroy(qdatahex,err)
240      call ceeddestroy(ceed,err)
241      end
242!-----------------------------------------------------------------------
243