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