xref: /libCEED/tests/t523-operator-f.f90 (revision 57af422ce4fd2bed083e8288205852f8262d25d3)
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_none,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_none,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_none,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_none,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 ceedoperatorcreatecomposite(ceed,op_setup,err)
209      call ceedoperatorsetname(op_setup,'setup',err)
210      call ceedoperatorcompositeaddsub(op_setup,op_setuptet,err)
211      call ceedoperatorcompositeaddsub(op_setup,op_setuphex,err)
212
213      call ceedoperatorcreatecomposite(ceed,op_mass,err)
214      call ceedoperatorsetname(op_mass,'mass',err)
215      call ceedoperatorsetnumviewtabs(op_mass,1,err)
216      call ceedoperatorcompositeaddsub(op_mass,op_masstet,err)
217      call ceedoperatorcompositeaddsub(op_mass,op_masshex,err)
218
219! View
220      call ceedoperatorview(op_setup,err)
221      call ceedoperatorview(op_mass,err)
222
223! Cleanup
224      call ceedqfunctiondestroy(qf_setuptet,err)
225      call ceedqfunctiondestroy(qf_masstet,err)
226      call ceedoperatordestroy(op_setuptet,err)
227      call ceedoperatordestroy(op_masstet,err)
228      call ceedqfunctiondestroy(qf_setuphex,err)
229      call ceedqfunctiondestroy(qf_masshex,err)
230      call ceedoperatordestroy(op_setuphex,err)
231      call ceedoperatordestroy(op_masshex,err)
232      call ceedoperatordestroy(op_setup,err)
233      call ceedoperatordestroy(op_mass,err)
234      call ceedelemrestrictiondestroy(erestrictutet,err)
235      call ceedelemrestrictiondestroy(erestrictxtet,err)
236      call ceedelemrestrictiondestroy(erestrictuitet,err)
237      call ceedelemrestrictiondestroy(erestrictuhex,err)
238      call ceedelemrestrictiondestroy(erestrictxhex,err)
239      call ceedelemrestrictiondestroy(erestrictuihex,err)
240      call ceedbasisdestroy(butet,err)
241      call ceedbasisdestroy(bxtet,err)
242      call ceedbasisdestroy(buhex,err)
243      call ceedbasisdestroy(bxhex,err)
244      call ceedvectordestroy(x,err)
245      call ceedvectordestroy(qdatatet,err)
246      call ceedvectordestroy(qdatahex,err)
247      call ceeddestroy(ceed,err)
248      end
249!-----------------------------------------------------------------------
250