xref: /libCEED/tests/t523-operator-f.f90 (revision 5fb68f377259d3910de46d787b7c5d1587fd01e1)
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 'ceed/fortran.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 ceedoperatorsetname(op_setuptet,'triangle elements',err)
125      call ceedoperatorsetfield(op_setuptet,'weight',&
126     & ceed_elemrestriction_none,bxtet,ceed_vector_none,err)
127      call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,&
128     & bxtet,ceed_vector_active,err)
129      call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,&
130     & ceed_basis_collocated,qdatatet,err)
131! ---- Mass Tet
132      call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,&
133     & ceed_qfunction_none,op_masstet,err)
134      call ceedoperatorsetname(op_masstet,'triangle elements',err)
135      call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,&
136     & ceed_basis_collocated,qdatatet,err)
137      call ceedoperatorsetfield(op_masstet,'u',erestrictutet,&
138     & butet,ceed_vector_active,err)
139      call ceedoperatorsetfield(op_masstet,'v',erestrictutet,&
140     & butet,ceed_vector_active,err)
141
142! Hex Elements
143      do i=0,nelemhex-1
144        col=mod(i,nx)
145        row=i/nx
146        offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2
147        do j=0,phex-1
148          do k=0,phex-1
149            indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j
150          enddo
151        enddo
152      enddo
153
154! -- Restrictions
155      call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,d,ndofs,d*ndofs,&
156     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err)
157
158      call ceedelemrestrictioncreate(ceed,nelemhex,phex*phex,1,1,ndofs,&
159     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err)
160      stridesuhex=[1,qhex*qhex,qhex*qhex]
161      call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,&
162     & 1,nqptshex,stridesuhex,erestrictuihex,err)
163
164! -- Bases
165      call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,&
166     & bxhex,err)
167      call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,&
168     & buhex,err)
169
170! -- QFunctions
171      call ceedqfunctioncreateinterior(ceed,1,setup,&
172     &SOURCE_DIR&
173     &//'t510-operator.h:setup'//char(0),qf_setuphex,err)
174      call ceedqfunctionaddinput(qf_setuphex,'weight',1,ceed_eval_weight,err)
175      call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err)
176      call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err)
177
178      call ceedqfunctioncreateinterior(ceed,1,mass,&
179     &SOURCE_DIR&
180     &//'t510-operator.h:mass'//char(0),qf_masshex,err)
181      call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err)
182      call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err)
183      call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err)
184
185! -- Operators
186! ---- Setup Hex
187      call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,&
188     & ceed_qfunction_none,op_setuphex,err)
189      call ceedoperatorsetname(op_setuphex,'quadrilateral elements',err)
190      call ceedoperatorsetfield(op_setuphex,'weight',&
191     & ceed_elemrestriction_none,bxhex,ceed_vector_none,err)
192      call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,&
193     & bxhex,ceed_vector_active,err)
194      call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,&
195     & ceed_basis_collocated,qdatahex,err)
196! ---- Mass Hex
197      call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,&
198     & ceed_qfunction_none,op_masshex,err)
199      call ceedoperatorsetname(op_masshex,'quadrilateral elements',err)
200      call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,&
201     & ceed_basis_collocated,qdatahex,err)
202      call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,&
203     & buhex,ceed_vector_active,err)
204      call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,&
205     & buhex,ceed_vector_active,err)
206
207! Composite Operators
208      call ceedcompositeoperatorcreate(ceed,op_setup,err)
209      call ceedoperatorsetname(op_setup,'setup',err)
210      call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err)
211      call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err)
212
213      call ceedcompositeoperatorcreate(ceed,op_mass,err)
214      call ceedoperatorsetname(op_mass,'mass',err)
215      call ceedcompositeoperatoraddsub(op_mass,op_masstet,err)
216      call ceedcompositeoperatoraddsub(op_mass,op_masshex,err)
217
218! View
219      call ceedoperatorview(op_setup,err)
220      call ceedoperatorview(op_mass,err)
221
222! Cleanup
223      call ceedqfunctiondestroy(qf_setuptet,err)
224      call ceedqfunctiondestroy(qf_masstet,err)
225      call ceedoperatordestroy(op_setuptet,err)
226      call ceedoperatordestroy(op_masstet,err)
227      call ceedqfunctiondestroy(qf_setuphex,err)
228      call ceedqfunctiondestroy(qf_masshex,err)
229      call ceedoperatordestroy(op_setuphex,err)
230      call ceedoperatordestroy(op_masshex,err)
231      call ceedoperatordestroy(op_setup,err)
232      call ceedoperatordestroy(op_mass,err)
233      call ceedelemrestrictiondestroy(erestrictutet,err)
234      call ceedelemrestrictiondestroy(erestrictxtet,err)
235      call ceedelemrestrictiondestroy(erestrictuitet,err)
236      call ceedelemrestrictiondestroy(erestrictuhex,err)
237      call ceedelemrestrictiondestroy(erestrictxhex,err)
238      call ceedelemrestrictiondestroy(erestrictuihex,err)
239      call ceedbasisdestroy(butet,err)
240      call ceedbasisdestroy(bxtet,err)
241      call ceedbasisdestroy(buhex,err)
242      call ceedbasisdestroy(bxhex,err)
243      call ceedvectordestroy(x,err)
244      call ceedvectordestroy(qdatatet,err)
245      call ceedvectordestroy(qdatahex,err)
246      call ceeddestroy(ceed,err)
247      end
248!-----------------------------------------------------------------------
249