xref: /libCEED/tests/t314-basis-f.f90 (revision a7bd39da094468e8222391a6dcb47d808ad9b491) 
1*a7bd39daSjeremylt!-----------------------------------------------------------------------
2*a7bd39daSjeremylt      subroutine eval(dimn,x,rslt)
3*a7bd39daSjeremylt      integer dimn
4*a7bd39daSjeremylt      real*8 x(3)
5*a7bd39daSjeremylt      real*8 rslt
6*a7bd39daSjeremylt
7*a7bd39daSjeremylt      rslt=tanh(x(1)+0.1)
8*a7bd39daSjeremylt      if (dimn>1) then
9*a7bd39daSjeremylt        rslt=rslt+atan(x(2)+0.2)
10*a7bd39daSjeremylt      endif
11*a7bd39daSjeremylt      if (dimn>2) then
12*a7bd39daSjeremylt        rslt=rslt+exp(-(x(3)+0.3)*(x(3)+0.3))
13*a7bd39daSjeremylt      endif
14*a7bd39daSjeremylt
15*a7bd39daSjeremylt      end
16*a7bd39daSjeremylt!-----------------------------------------------------------------------
17*a7bd39daSjeremylt      program test
18*a7bd39daSjeremylt
19*a7bd39daSjeremylt      include 'ceedf.h'
20*a7bd39daSjeremylt
21*a7bd39daSjeremylt      integer ceed,err
22*a7bd39daSjeremylt      integer x,xq,u,uq,ones,gtposeones
23*a7bd39daSjeremylt      integer bxl,bug
24*a7bd39daSjeremylt      integer dimn,d
25*a7bd39daSjeremylt      integer i
26*a7bd39daSjeremylt      integer p
27*a7bd39daSjeremylt      integer q
28*a7bd39daSjeremylt      parameter(p=8)
29*a7bd39daSjeremylt      parameter(q=7)
30*a7bd39daSjeremylt      integer maxdim
31*a7bd39daSjeremylt      parameter(maxdim=3)
32*a7bd39daSjeremylt      integer qdimnmax
33*a7bd39daSjeremylt      parameter(qdimnmax=q**maxdim)
34*a7bd39daSjeremylt      integer pdimnmax
35*a7bd39daSjeremylt      parameter(pdimnmax=p**maxdim)
36*a7bd39daSjeremylt      integer xdimmax
37*a7bd39daSjeremylt      parameter(xdimmax=2**maxdim)
38*a7bd39daSjeremylt      integer pdimn,qdimn,xdim
39*a7bd39daSjeremylt
40*a7bd39daSjeremylt      real*8 xx(xdimmax*maxdim)
41*a7bd39daSjeremylt      real*8 xxx(maxdim)
42*a7bd39daSjeremylt      real*8 xxq(pdimnmax*maxdim)
43*a7bd39daSjeremylt      real*8 uuq(qdimnmax*maxdim)
44*a7bd39daSjeremylt      real*8 uu(pdimnmax)
45*a7bd39daSjeremylt      real*8 ggtposeones(pdimnmax)
46*a7bd39daSjeremylt      real*8 sum1
47*a7bd39daSjeremylt      real*8 sum2
48*a7bd39daSjeremylt      integer dimxqdimn
49*a7bd39daSjeremylt      integer*8 uoffset,xoffset,offset1,offset2,offset3
50*a7bd39daSjeremylt
51*a7bd39daSjeremylt      character arg*32
52*a7bd39daSjeremylt
53*a7bd39daSjeremylt      call getarg(1,arg)
54*a7bd39daSjeremylt      call ceedinit(trim(arg)//char(0),ceed,err)
55*a7bd39daSjeremylt
56*a7bd39daSjeremylt      do dimn=1,maxdim
57*a7bd39daSjeremylt        qdimn=q**dimn
58*a7bd39daSjeremylt        pdimn=p**dimn
59*a7bd39daSjeremylt        xdim=2**dimn
60*a7bd39daSjeremylt        dimxqdimn=dimn*qdimn
61*a7bd39daSjeremylt        sum1=0
62*a7bd39daSjeremylt        sum2=0
63*a7bd39daSjeremylt
64*a7bd39daSjeremylt        do d=0,dimn-1
65*a7bd39daSjeremylt          do i=1,xdim
66*a7bd39daSjeremylt            if ((mod(i-1,2**(dimn-d))/(2**(dimn-d-1))).ne.0) then
67*a7bd39daSjeremylt              xx(d*xdim+i)=1
68*a7bd39daSjeremylt            else
69*a7bd39daSjeremylt              xx(d*xdim+i)=-1
70*a7bd39daSjeremylt            endif
71*a7bd39daSjeremylt          enddo
72*a7bd39daSjeremylt        enddo
73*a7bd39daSjeremylt
74*a7bd39daSjeremylt        call ceedvectorcreate(ceed,xdim*dimn,x,err)
75*a7bd39daSjeremylt        xoffset=0
76*a7bd39daSjeremylt        call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,xx,xoffset,err)
77*a7bd39daSjeremylt        call ceedvectorcreate(ceed,pdimn*dimn,xq,err)
78*a7bd39daSjeremylt        call ceedvectorsetvalue(xq,0.d0,err)
79*a7bd39daSjeremylt        call ceedvectorcreate(ceed,pdimn,u,err)
80*a7bd39daSjeremylt        call ceedvectorcreate(ceed,qdimn*dimn,uq,err)
81*a7bd39daSjeremylt        call ceedvectorsetvalue(uq,0.d0,err)
82*a7bd39daSjeremylt        call ceedvectorcreate(ceed,qdimn*dimn,ones,err)
83*a7bd39daSjeremylt        call ceedvectorsetvalue(ones,1.d0,err)
84*a7bd39daSjeremylt        call ceedvectorcreate(ceed,pdimn,gtposeones,err)
85*a7bd39daSjeremylt        call ceedvectorsetvalue(gtposeones,0.d0,err)
86*a7bd39daSjeremylt
87*a7bd39daSjeremylt        call ceedbasiscreatetensorh1lagrange(ceed,dimn,dimn,2,p,&
88*a7bd39daSjeremylt     &   ceed_gauss_lobatto,bxl,err)
89*a7bd39daSjeremylt        call ceedbasisapply(bxl,1,ceed_notranspose,ceed_eval_interp,x,xq,err)
90*a7bd39daSjeremylt
91*a7bd39daSjeremylt        call ceedvectorgetarrayread(xq,ceed_mem_host,xxq,offset1,err)
92*a7bd39daSjeremylt        do i=1,pdimn
93*a7bd39daSjeremylt          do d=0,dimn-1
94*a7bd39daSjeremylt            xxx(d+1)=xxq(d*pdimn+i+offset1)
95*a7bd39daSjeremylt          enddo
96*a7bd39daSjeremylt          call eval(dimn,xxx,uu(i))
97*a7bd39daSjeremylt        enddo
98*a7bd39daSjeremylt        call ceedvectorrestorearrayread(xq,xxq,offset1,err)
99*a7bd39daSjeremylt        uoffset=0
100*a7bd39daSjeremylt        call ceedvectorsetarray(u,ceed_mem_host,ceed_use_pointer,uu,uoffset,err)
101*a7bd39daSjeremylt
102*a7bd39daSjeremylt        call ceedbasiscreatetensorh1lagrange(ceed,dimn,1,p,q,ceed_gauss,bug,err)
103*a7bd39daSjeremylt
104*a7bd39daSjeremylt        call ceedbasisapply(bug,1,ceed_notranspose,ceed_eval_grad,u,uq,err)
105*a7bd39daSjeremylt        call ceedbasisapply(bug,1,ceed_transpose,ceed_eval_grad,ones,&
106*a7bd39daSjeremylt     &   gtposeones,err)
107*a7bd39daSjeremylt
108*a7bd39daSjeremylt        call ceedvectorgetarrayread(gtposeones,ceed_mem_host,ggtposeones,&
109*a7bd39daSjeremylt     &   offset1,err)
110*a7bd39daSjeremylt        call ceedvectorgetarrayread(u,ceed_mem_host,uu,offset2,err)
111*a7bd39daSjeremylt        call ceedvectorgetarrayread(uq,ceed_mem_host,uuq,offset3,err)
112*a7bd39daSjeremylt        do i=1,pdimn
113*a7bd39daSjeremylt          sum1=sum1+ggtposeones(i+offset1)*uu(i+offset2)
114*a7bd39daSjeremylt        enddo
115*a7bd39daSjeremylt        do i=1,dimxqdimn
116*a7bd39daSjeremylt          sum2=sum2+uuq(i+offset3)
117*a7bd39daSjeremylt        enddo
118*a7bd39daSjeremylt        call ceedvectorrestorearrayread(gtposeones,ggtposeones,offset1,err)
119*a7bd39daSjeremylt        call ceedvectorrestorearrayread(u,uu,offset2,err)
120*a7bd39daSjeremylt        call ceedvectorrestorearrayread(uq,uuq,offset3,err)
121*a7bd39daSjeremylt        if(abs(sum1-sum2) > 1.0D-10) then
122*a7bd39daSjeremylt          write(*,'(A,I1,A,F12.6,A,F12.6)')'[',dimn,'] Error: ',sum1,  ' != ',&
123*a7bd39daSjeremylt     &     sum2
124*a7bd39daSjeremylt        endif
125*a7bd39daSjeremylt
126*a7bd39daSjeremylt        call ceedvectordestroy(x,err)
127*a7bd39daSjeremylt        call ceedvectordestroy(xq,err)
128*a7bd39daSjeremylt        call ceedvectordestroy(u,err)
129*a7bd39daSjeremylt        call ceedvectordestroy(uq,err)
130*a7bd39daSjeremylt        call ceedvectordestroy(ones,err)
131*a7bd39daSjeremylt        call ceedvectordestroy(gtposeones,err)
132*a7bd39daSjeremylt        call ceedbasisdestroy(bxl,err)
133*a7bd39daSjeremylt        call ceedbasisdestroy(bug,err)
134*a7bd39daSjeremylt      enddo
135*a7bd39daSjeremylt
136*a7bd39daSjeremylt      call ceeddestroy(ceed,err)
137*a7bd39daSjeremylt      end
138*a7bd39daSjeremylt!-----------------------------------------------------------------------
139