xref: /libCEED/tests/t502-operator-f.f90 (revision d80fc06a6882f9ed5d7208de2efea3af5b025f69)
1!-----------------------------------------------------------------------
2      subroutine setup(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,&
3&           u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr)
4      real*8 ctx
5      real*8 u1(1)
6      real*8 u2(1)
7      real*8 v1(1)
8      integer q,ierr
9
10      do i=1,q
11        v1(i)=u1(i)*u2(i)
12      enddo
13
14      ierr=0
15      end
16!-----------------------------------------------------------------------
17      subroutine mass(ctx,q,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,&
18&           u15,u16,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,ierr)
19      real*8 ctx
20      real*8 u1(1)
21      real*8 u2(1)
22      real*8 v1(1)
23      integer q,ierr
24
25      do i=1,q
26        v1(i)=u2(i)*u1(i)
27        v1(q+i)=u2(q+i)*u1(i)
28      enddo
29
30      ierr=0
31      end
32!-----------------------------------------------------------------------
33      program test
34
35      include 'ceedf.h'
36
37      integer ceed,err,i,j
38      integer imode
39      parameter(imode=ceed_noninterlaced)
40      integer imodeu
41      parameter(imodeu=ceed_interlaced)
42      integer stridesu(3)
43      integer erestrictx,erestrictu,erestrictui
44      integer bx,bu
45      integer qf_setup,qf_mass
46      integer op_setup,op_mass
47      integer qdata,x,u,v
48      integer nelem,p,q
49      parameter(nelem=15)
50      parameter(p=5)
51      parameter(q=8)
52      integer nx,nu
53      parameter(nx=nelem+1)
54      parameter(nu=nelem*(p-1)+1)
55      integer indx(nelem*2)
56      integer indu(nelem*p)
57      real*8 arrx(nx)
58      integer*8 voffset,xoffset
59
60      real*8 hu(nu*2),hv(nu*2)
61      real*8 total1,total2
62
63      character arg*32
64
65      external setup,mass
66
67      call getarg(1,arg)
68      call ceedinit(trim(arg)//char(0),ceed,err)
69
70      do i=0,nx-1
71        arrx(i+1)=i/(nx-1.d0)
72      enddo
73      do i=0,nelem-1
74        indx(2*i+1)=i
75        indx(2*i+2)=i+1
76      enddo
77
78      call ceedelemrestrictioncreate(ceed,imode,nelem,2,nx,1,ceed_mem_host,&
79     & ceed_use_pointer,indx,erestrictx,err)
80
81      do i=0,nelem-1
82        do j=0,p-1
83          indu(p*i+j+1)=i*(p-1)+j
84        enddo
85      enddo
86
87      call ceedelemrestrictioncreate(ceed,imodeu,nelem,p,nu,2,ceed_mem_host,&
88     & ceed_use_pointer,indu,erestrictu,err)
89      stridesu=[1,q,q]
90      call ceedelemrestrictioncreatestrided(ceed,nelem,q,q*nelem,1,stridesu,&
91     & erestrictui,err)
92
93      call ceedbasiscreatetensorh1lagrange(ceed,1,1,2,q,ceed_gauss,bx,err)
94      call ceedbasiscreatetensorh1lagrange(ceed,1,2,p,q,ceed_gauss,bu,err)
95
96      call ceedqfunctioncreateinterior(ceed,1,setup,&
97     &SOURCE_DIR&
98     &//'t502-operator.h:setup'//char(0),qf_setup,err)
99      call ceedqfunctionaddinput(qf_setup,'_weight',1,ceed_eval_weight,err)
100      call ceedqfunctionaddinput(qf_setup,'dx',1,ceed_eval_grad,err)
101      call ceedqfunctionaddoutput(qf_setup,'rho',1,ceed_eval_none,err)
102
103      call ceedqfunctioncreateinterior(ceed,1,mass,&
104     &SOURCE_DIR&
105     &//'t502-operator.h:mass'//char(0),qf_mass,err)
106      call ceedqfunctionaddinput(qf_mass,'rho',1,ceed_eval_none,err)
107      call ceedqfunctionaddinput(qf_mass,'u',2,ceed_eval_interp,err)
108      call ceedqfunctionaddoutput(qf_mass,'v',2,ceed_eval_interp,err)
109
110      call ceedoperatorcreate(ceed,qf_setup,ceed_qfunction_none,&
111     & ceed_qfunction_none,op_setup,err)
112      call ceedoperatorcreate(ceed,qf_mass,ceed_qfunction_none,&
113     & ceed_qfunction_none,op_mass,err)
114
115      call ceedvectorcreate(ceed,nx,x,err)
116      xoffset=0
117      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,arrx,xoffset,err)
118      call ceedvectorcreate(ceed,nelem*q,qdata,err)
119
120      call ceedoperatorsetfield(op_setup,'_weight',ceed_elemrestriction_none,&
121     & bx,ceed_vector_none,err)
122      call ceedoperatorsetfield(op_setup,'dx',erestrictx,bx,&
123     & ceed_vector_active,err)
124      call ceedoperatorsetfield(op_setup,'rho',erestrictui,&
125     & ceed_basis_collocated,ceed_vector_active,err)
126      call ceedoperatorsetfield(op_mass,'rho',erestrictui,&
127     & ceed_basis_collocated,qdata,err)
128      call ceedoperatorsetfield(op_mass,'u',erestrictu,bu,&
129     & ceed_vector_active,err)
130      call ceedoperatorsetfield(op_mass,'v',erestrictu,bu,&
131     & ceed_vector_active,err)
132
133      call ceedoperatorapply(op_setup,x,qdata,ceed_request_immediate,err)
134
135      call ceedvectorcreate(ceed,2*nu,u,err)
136      call ceedvectorgetarray(u,ceed_mem_host,hu,voffset,err)
137      do i=1,nu
138        hu(voffset+2*i-1)=1.
139        hu(voffset+2*i)=2.
140      enddo
141      call ceedvectorrestorearray(u,hu,voffset,err)
142      call ceedvectorcreate(ceed,2*nu,v,err)
143      call ceedoperatorapply(op_mass,u,v,ceed_request_immediate,err)
144
145      call ceedvectorgetarrayread(v,ceed_mem_host,hv,voffset,err)
146      total1=0.
147      total2=0.
148      do i=1,nu
149        total1=total1+hv(voffset+2*i-1)
150        total2=total2+hv(voffset+2*i)
151      enddo
152      if (abs(total1-1.)>1.0d-10) then
153! LCOV_EXCL_START
154        write(*,*) 'Computed Area: ',total1,' != True Area: 1.0'
155! LCOV_EXCL_STOP
156      endif
157      if (abs(total2-2.)>1.0d-10) then
158! LCOV_EXCL_START
159        write(*,*) 'Computed Area: ',total2,' != True Area: 2.0'
160! LCOV_EXCL_STOP
161      endif
162      call ceedvectorrestorearrayread(v,hv,voffset,err)
163
164      call ceedvectordestroy(qdata,err)
165      call ceedvectordestroy(x,err)
166      call ceedvectordestroy(u,err)
167      call ceedvectordestroy(v,err)
168      call ceedoperatordestroy(op_mass,err)
169      call ceedoperatordestroy(op_setup,err)
170      call ceedqfunctiondestroy(qf_mass,err)
171      call ceedqfunctiondestroy(qf_setup,err)
172      call ceedbasisdestroy(bu,err)
173      call ceedbasisdestroy(bx,err)
174      call ceedelemrestrictiondestroy(erestrictu,err)
175      call ceedelemrestrictiondestroy(erestrictx,err)
176      call ceedelemrestrictiondestroy(erestrictui,err)
177      call ceeddestroy(ceed,err)
178      end
179!-----------------------------------------------------------------------
180