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