xref: /libCEED/tests/t523-operator-f.f90 (revision 823118015a84d7808a28beecbdbc99e4bf09a8a7)
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 erestrictxtet,erestrictutet,erestrictxitet,erestrictuitet,&
49&             erestrictxhex,erestrictuhex,erestrictxihex,erestrictuihex
50      integer bxtet,butet,bxhex,buhex
51      integer qf_setuptet,qf_masstet,qf_setuphex,qf_masshex
52      integer op_setuptet,op_masstet,op_setuphex,op_masshex,op_setup,op_mass
53      integer qdatatet,qdatahex,x
54      integer nelemtet,nelemhex,ptet,phex,qtet,qhex,d
55      integer row,col,offset
56      parameter(nelemtet=6)
57      parameter(ptet=6)
58      parameter(qtet=4)
59      parameter(nelemhex=6)
60      parameter(phex=3)
61      parameter(qhex=4)
62      parameter(d=2)
63      integer ndofs,nqptstet,nqptshex,nqpts,nx,ny,nxtet,nytet,nxhex
64      parameter(nx=3)
65      parameter(ny=3)
66      parameter(nxtet=3)
67      parameter(nytet=1)
68      parameter(nxhex=3)
69      parameter(ndofs=(nx*2+1)*(ny*2+1))
70      parameter(nqptstet=nelemtet*qtet)
71      parameter(nqptshex=nelemhex*qhex*qhex)
72      parameter(nqpts=nqptstet+nqptshex)
73      integer indxtet(nelemtet*ptet),indxhex(nelemhex*phex*phex)
74
75      real*8 qref(d*qtet)
76      real*8 qweight(qtet)
77      real*8 interp(ptet*qtet)
78      real*8 grad(d*ptet*qtet)
79
80      character arg*32
81
82! LCOV_EXCL_START
83      external setup,mass
84! LCOV_EXCL_STOP
85
86      call getarg(1,arg)
87
88      call ceedinit(trim(arg)//char(0),ceed,err)
89
90! DoF Coordinates
91      call ceedvectorcreate(ceed,d*ndofs,x,err)
92
93! Qdata Vectors
94      call ceedvectorcreate(ceed,nqptstet,qdatatet,err)
95      call ceedvectorcreate(ceed,nqptshex,qdatahex,err)
96
97! Tet Elements
98      do i=0,2
99        col=mod(i,nx)
100        row=i/nx
101        offset=col*2+row*(nx*2+1)*2
102
103        indxtet(i*2*ptet+1)=2+offset
104        indxtet(i*2*ptet+2)=9+offset
105        indxtet(i*2*ptet+3)=16+offset
106        indxtet(i*2*ptet+4)=1+offset
107        indxtet(i*2*ptet+5)=8+offset
108        indxtet(i*2*ptet+6)=0+offset
109
110        indxtet(i*2*ptet+7)=14+offset
111        indxtet(i*2*ptet+8)=7+offset
112        indxtet(i*2*ptet+9)=0+offset
113        indxtet(i*2*ptet+10)=15+offset
114        indxtet(i*2*ptet+11)=8+offset
115        indxtet(i*2*ptet+12)=16+offset
116      enddo
117
118! -- Restrictions
119      call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,d,&
120     & ceed_mem_host,ceed_use_pointer,indxtet,erestrictxtet,err)
121      call ceedelemrestrictioncreateidentity(ceed,imode,nelemtet,ptet,&
122     & nelemtet*ptet,d,erestrictxitet,err)
123
124      call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,&
125     & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err)
126      call ceedelemrestrictioncreateidentity(ceed,imode,nelemtet,qtet,nqptstet,&
127     & 1,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',erestrictxitet,&
157     & 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      call ceedelemrestrictioncreateidentity(ceed,imode,nelemhex,phex*phex,&
188     & nelemhex*phex*phex,d,erestrictxihex,err)
189
190      call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,&
191     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err)
192      call ceedelemrestrictioncreateidentity(ceed,imode,nelemhex,qhex*qhex,&
193     & nqptshex,1,erestrictuihex,err)
194
195! -- Bases
196      call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,&
197     & bxhex,err)
198      call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,&
199     & buhex,err)
200
201! -- QFunctions
202      call ceedqfunctioncreateinterior(ceed,1,setup,&
203     &SOURCE_DIR&
204     &//'t510-operator.h:setup'//char(0),qf_setuphex,err)
205      call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err)
206      call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err)
207      call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err)
208
209      call ceedqfunctioncreateinterior(ceed,1,mass,&
210     &SOURCE_DIR&
211     &//'t510-operator.h:mass'//char(0),qf_masshex,err)
212      call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err)
213      call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err)
214      call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err)
215
216! -- Operators
217! ---- Setup Hex
218      call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,&
219     & ceed_qfunction_none,op_setuphex,err)
220      call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,&
221     & bxhex,ceed_vector_none,err)
222      call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,&
223     & bxhex,ceed_vector_active,err)
224      call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,&
225     & ceed_basis_collocated,qdatahex,err)
226! ---- Mass Hex
227      call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,&
228     & ceed_qfunction_none,op_masshex,err)
229      call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,&
230     & ceed_basis_collocated,qdatahex,err)
231      call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,&
232     & buhex,ceed_vector_active,err)
233      call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,&
234     & buhex,ceed_vector_active,err)
235
236! Composite Operators
237      call ceedcompositeoperatorcreate(ceed,op_setup,err)
238      call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err)
239      call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err)
240
241      call ceedcompositeoperatorcreate(ceed,op_mass,err)
242      call ceedcompositeoperatoraddsub(op_mass,op_masstet,err)
243      call ceedcompositeoperatoraddsub(op_mass,op_masshex,err)
244
245! View
246      call ceedoperatorview(op_setup,err)
247      call ceedoperatorview(op_mass,err)
248
249! Cleanup
250      call ceedqfunctiondestroy(qf_setuptet,err)
251      call ceedqfunctiondestroy(qf_masstet,err)
252      call ceedoperatordestroy(op_setuptet,err)
253      call ceedoperatordestroy(op_masstet,err)
254      call ceedqfunctiondestroy(qf_setuphex,err)
255      call ceedqfunctiondestroy(qf_masshex,err)
256      call ceedoperatordestroy(op_setuphex,err)
257      call ceedoperatordestroy(op_masshex,err)
258      call ceedoperatordestroy(op_setup,err)
259      call ceedoperatordestroy(op_mass,err)
260      call ceedelemrestrictiondestroy(erestrictutet,err)
261      call ceedelemrestrictiondestroy(erestrictxtet,err)
262      call ceedelemrestrictiondestroy(erestrictuitet,err)
263      call ceedelemrestrictiondestroy(erestrictxitet,err)
264      call ceedelemrestrictiondestroy(erestrictuhex,err)
265      call ceedelemrestrictiondestroy(erestrictxhex,err)
266      call ceedelemrestrictiondestroy(erestrictuihex,err)
267      call ceedelemrestrictiondestroy(erestrictxihex,err)
268      call ceedbasisdestroy(butet,err)
269      call ceedbasisdestroy(bxtet,err)
270      call ceedbasisdestroy(buhex,err)
271      call ceedbasisdestroy(bxhex,err)
272      call ceedvectordestroy(x,err)
273      call ceedvectordestroy(qdatatet,err)
274      call ceedvectordestroy(qdatahex,err)
275      call ceeddestroy(ceed,err)
276      end
277!-----------------------------------------------------------------------
278