xref: /libCEED/tests/t523-operator-f.f90 (revision ccedf6b0eedde5bdb26f1d628a64390ea2c01243)
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 stridesxtet(3),stridesutet(3),stridesxhex(3),stridesuhex(3)
49      integer erestrictxtet,erestrictutet,erestrictxitet,erestrictuitet,&
50&             erestrictxhex,erestrictuhex,erestrictxihex,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      stridesxtet=[1,ptet,ptet*d]
123      call ceedelemrestrictioncreatestrided(ceed,nelemtet,ptet,&
124     & nelemtet*ptet,d,stridesxtet,erestrictxitet,err)
125
126      call ceedelemrestrictioncreate(ceed,imode,nelemtet,ptet,ndofs,1,&
127     & ceed_mem_host,ceed_use_pointer,indxtet,erestrictutet,err)
128      stridesutet=[1,qtet,qtet]
129      call ceedelemrestrictioncreatestrided(ceed,nelemtet,qtet,nqptstet,&
130     & 1,stridesutet,erestrictuitet,err)
131
132! -- Bases
133      call buildmats(qref,qweight,interp,grad)
134      call ceedbasiscreateh1(ceed,ceed_triangle,d,ptet,qtet,interp,grad,qref,&
135     & qweight,bxtet,err)
136      call buildmats(qref,qweight,interp,grad)
137      call ceedbasiscreateh1(ceed,ceed_triangle,1,ptet,qtet,interp,grad,qref,&
138     & qweight,butet,err)
139
140! -- QFunctions
141      call ceedqfunctioncreateinterior(ceed,1,setup,&
142     &SOURCE_DIR&
143     &//'t510-operator.h:setup'//char(0),qf_setuptet,err)
144      call ceedqfunctionaddinput(qf_setuptet,'_weight',1,ceed_eval_weight,err)
145      call ceedqfunctionaddinput(qf_setuptet,'dx',d*d,ceed_eval_grad,err)
146      call ceedqfunctionaddoutput(qf_setuptet,'rho',1,ceed_eval_none,err)
147
148      call ceedqfunctioncreateinterior(ceed,1,mass,&
149     &SOURCE_DIR&
150     &//'t510-operator.h:mass'//char(0),qf_masstet,err)
151      call ceedqfunctionaddinput(qf_masstet,'rho',1,ceed_eval_none,err)
152      call ceedqfunctionaddinput(qf_masstet,'u',1,ceed_eval_interp,err)
153      call ceedqfunctionaddoutput(qf_masstet,'v',1,ceed_eval_interp,err)
154
155! -- Operators
156! ---- Setup Tet
157      call ceedoperatorcreate(ceed,qf_setuptet,ceed_qfunction_none,&
158     & ceed_qfunction_none,op_setuptet,err)
159      call ceedoperatorsetfield(op_setuptet,'_weight',erestrictxitet,&
160     & bxtet,ceed_vector_none,err)
161      call ceedoperatorsetfield(op_setuptet,'dx',erestrictxtet,&
162     & bxtet,ceed_vector_active,err)
163      call ceedoperatorsetfield(op_setuptet,'rho',erestrictuitet,&
164     & ceed_basis_collocated,qdatatet,err)
165! ---- Mass Tet
166      call ceedoperatorcreate(ceed,qf_masstet,ceed_qfunction_none,&
167     & ceed_qfunction_none,op_masstet,err)
168      call ceedoperatorsetfield(op_masstet,'rho',erestrictuitet,&
169     & ceed_basis_collocated,qdatatet,err)
170      call ceedoperatorsetfield(op_masstet,'u',erestrictutet,&
171     & butet,ceed_vector_active,err)
172      call ceedoperatorsetfield(op_masstet,'v',erestrictutet,&
173     & butet,ceed_vector_active,err)
174
175! Hex Elements
176      do i=0,nelemhex-1
177        col=mod(i,nx)
178        row=i/nx
179        offset=(nxtet*2+1)*(nytet*2)*(1+row)+col*2
180        do j=0,phex-1
181          do k=0,phex-1
182            indxhex(phex*(phex*i+k)+j+1)=offset+k*(nxhex*2+1)+j
183          enddo
184        enddo
185      enddo
186
187! -- Restrictions
188      call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,d,&
189     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictxhex,err)
190      stridesxhex=[1,phex*phex,phex*phex*d]
191      call ceedelemrestrictioncreatestrided(ceed,nelemhex,phex*phex,&
192     & nelemhex*phex*phex,d,stridesxhex,erestrictxihex,err)
193
194      call ceedelemrestrictioncreate(ceed,imode,nelemhex,phex*phex,ndofs,1,&
195     & ceed_mem_host,ceed_use_pointer,indxhex,erestrictuhex,err)
196      stridesuhex=[1,qhex*qhex,qhex*qhex]
197      call ceedelemrestrictioncreatestrided(ceed,nelemhex,qhex*qhex,&
198     & nqptshex,1,stridesuhex,erestrictuihex,err)
199
200! -- Bases
201      call ceedbasiscreatetensorh1lagrange(ceed,d,d,phex,qhex,ceed_gauss,&
202     & bxhex,err)
203      call ceedbasiscreatetensorh1lagrange(ceed,d,1,phex,qhex,ceed_gauss,&
204     & buhex,err)
205
206! -- QFunctions
207      call ceedqfunctioncreateinterior(ceed,1,setup,&
208     &SOURCE_DIR&
209     &//'t510-operator.h:setup'//char(0),qf_setuphex,err)
210      call ceedqfunctionaddinput(qf_setuphex,'_weight',1,ceed_eval_weight,err)
211      call ceedqfunctionaddinput(qf_setuphex,'dx',d*d,ceed_eval_grad,err)
212      call ceedqfunctionaddoutput(qf_setuphex,'rho',1,ceed_eval_none,err)
213
214      call ceedqfunctioncreateinterior(ceed,1,mass,&
215     &SOURCE_DIR&
216     &//'t510-operator.h:mass'//char(0),qf_masshex,err)
217      call ceedqfunctionaddinput(qf_masshex,'rho',1,ceed_eval_none,err)
218      call ceedqfunctionaddinput(qf_masshex,'u',1,ceed_eval_interp,err)
219      call ceedqfunctionaddoutput(qf_masshex,'v',1,ceed_eval_interp,err)
220
221! -- Operators
222! ---- Setup Hex
223      call ceedoperatorcreate(ceed,qf_setuphex,ceed_qfunction_none,&
224     & ceed_qfunction_none,op_setuphex,err)
225      call ceedoperatorsetfield(op_setuphex,'_weight',erestrictxihex,&
226     & bxhex,ceed_vector_none,err)
227      call ceedoperatorsetfield(op_setuphex,'dx',erestrictxhex,&
228     & bxhex,ceed_vector_active,err)
229      call ceedoperatorsetfield(op_setuphex,'rho',erestrictuihex,&
230     & ceed_basis_collocated,qdatahex,err)
231! ---- Mass Hex
232      call ceedoperatorcreate(ceed,qf_masshex,ceed_qfunction_none,&
233     & ceed_qfunction_none,op_masshex,err)
234      call ceedoperatorsetfield(op_masshex,'rho',erestrictuihex,&
235     & ceed_basis_collocated,qdatahex,err)
236      call ceedoperatorsetfield(op_masshex,'u',erestrictuhex,&
237     & buhex,ceed_vector_active,err)
238      call ceedoperatorsetfield(op_masshex,'v',erestrictuhex,&
239     & buhex,ceed_vector_active,err)
240
241! Composite Operators
242      call ceedcompositeoperatorcreate(ceed,op_setup,err)
243      call ceedcompositeoperatoraddsub(op_setup,op_setuptet,err)
244      call ceedcompositeoperatoraddsub(op_setup,op_setuphex,err)
245
246      call ceedcompositeoperatorcreate(ceed,op_mass,err)
247      call ceedcompositeoperatoraddsub(op_mass,op_masstet,err)
248      call ceedcompositeoperatoraddsub(op_mass,op_masshex,err)
249
250! View
251      call ceedoperatorview(op_setup,err)
252      call ceedoperatorview(op_mass,err)
253
254! Cleanup
255      call ceedqfunctiondestroy(qf_setuptet,err)
256      call ceedqfunctiondestroy(qf_masstet,err)
257      call ceedoperatordestroy(op_setuptet,err)
258      call ceedoperatordestroy(op_masstet,err)
259      call ceedqfunctiondestroy(qf_setuphex,err)
260      call ceedqfunctiondestroy(qf_masshex,err)
261      call ceedoperatordestroy(op_setuphex,err)
262      call ceedoperatordestroy(op_masshex,err)
263      call ceedoperatordestroy(op_setup,err)
264      call ceedoperatordestroy(op_mass,err)
265      call ceedelemrestrictiondestroy(erestrictutet,err)
266      call ceedelemrestrictiondestroy(erestrictxtet,err)
267      call ceedelemrestrictiondestroy(erestrictuitet,err)
268      call ceedelemrestrictiondestroy(erestrictxitet,err)
269      call ceedelemrestrictiondestroy(erestrictuhex,err)
270      call ceedelemrestrictiondestroy(erestrictxhex,err)
271      call ceedelemrestrictiondestroy(erestrictuihex,err)
272      call ceedelemrestrictiondestroy(erestrictxihex,err)
273      call ceedbasisdestroy(butet,err)
274      call ceedbasisdestroy(bxtet,err)
275      call ceedbasisdestroy(buhex,err)
276      call ceedbasisdestroy(bxhex,err)
277      call ceedvectordestroy(x,err)
278      call ceedvectordestroy(qdatatet,err)
279      call ceedvectordestroy(qdatahex,err)
280      call ceeddestroy(ceed,err)
281      end
282!-----------------------------------------------------------------------
283