Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
lob_bnd_test.c
Go to the documentation of this file.
1 #include <stddef.h>
2 #include <stdlib.h>
3 #include <stdio.h>
4 #include <math.h>
5 #include <string.h>
6 #include "c99.h"
7 #include "types.h"
8 #include "name.h"
9 #include "fail.h"
10 #include "mem.h"
11 #include "tensor.h"
12 #include "poly.h"
13 #include "lob_bnd.h"
14 
15 
16 #define RESFAC 4
17 #define N 12
18 #define NY 9
19 #define NZ 4
20 #define REPEAT 1000000
21 
22 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923
23 
24 
25 int main()
26 {
27  int failure=0;
28  uint i,r;
29  double p[NZ*NY*N];
30  double lb[2*(RESFAC*NZ)*(RESFAC*NY)*(RESFAC*N)];
31  double work[2*(RESFAC*N)*(RESFAC*NY)*(NZ+1)];
32 
33  double *ld_N = tmalloc(double,N+gll_lag_size(N));
34  lagrange_fun *lag_N = gll_lag_setup(ld_N+N,N);
35 
36  double *ld_NY= tmalloc(double,NY+gll_lag_size(NY));
37  lagrange_fun *lag_NY = gll_lag_setup(ld_NY+NY,NY);
38 
39  double *ld_NZ= tmalloc(double,NZ+gll_lag_size(NZ));
40  lagrange_fun *lag_NZ = gll_lag_setup(ld_NZ+NZ,NZ);
41 
42  double *lb_N = tmalloc(double,lob_bnd_size(N ,RESFAC*N ));
43  double *lb_NY = tmalloc(double,lob_bnd_size(NY,RESFAC*NY));
44  double *lb_NZ = tmalloc(double,lob_bnd_size(NZ,RESFAC*NZ));
45  lob_bnd_setup(lb_N , N ,RESFAC*N );
46  lob_bnd_setup(lb_NY, NY,RESFAC*NY );
47  lob_bnd_setup(lb_NZ, NZ,RESFAC*NZ );
48  /*for(i=0;i<NY*N;++i) p[i]=rand()/(double)RAND_MAX;
49 
50  lob_bnd_lin_1(lb, lb_N,N,RESFAC*N, p,NY); */
51 
52  /* 1D */
53  for(r=0;r<REPEAT;++r) {
54  int m = RESFAC*N;
55  double x = (rand()/(double)RAND_MAX)*2-1;
56  /* x = cos((m-1-j)*PI/(m-1)) */
57  int j = -1 + m - 1 - (int) (acos(x) * (m-1) / PI);
58  double f = (x - cos((m-1-j)*PI/(m-1))) /
59  (cos((m-1-(j+1))*PI/(m-1)) - cos((m-1-j)*PI/(m-1)));
60 
61  if(r%256==0) {
62  for(i=0;i<NY*N;++i) p[i]=rand()/(double)RAND_MAX;
63 
64  lob_bnd_lin_1(lb, lb_N,N,RESFAC*N, p,NY);
65  }
66 
67  if(r<3)
68  printf("%g <= %g <= %g, f = %g\n",
69  cos((m-1-j)*PI/(m-1)), x, cos((m-1-(j+1))*PI/(m-1)), f);
70  lag_N(ld_N,ld_N+N,N,0,x);
71  for(i=0;i<NY;++i) {
72  double lo = (1-f)*lb[(i*m+j)*2 ] + f*lb[(i*m+j+1)*2 ],
73  up = (1-f)*lb[(i*m+j)*2+1] + f*lb[(i*m+j+1)*2+1],
74  px = tensor_dot(ld_N,p+i*N,N);
75  if(r<3 || px < lo || up < px)
76  printf("p_%02d(%g) = %g in [%g,%g]\n",i,x,px,lo,up);
77  if(px<lo || up<px) {failure=1; break;}
78  }
79  if(i!=NY) break;
80  }
81 
82  /* x = cos((m-1-j)*PI/(m-1)) */
83  #define GET_JF(x) \
84  int j##x = -1 + m##x - 1 - (int) (acos(x) * (m##x-1) / PI); \
85  double f##x = (x - cos((m##x-1-j##x)*PI/(m##x-1))) / \
86  (cos((m##x-1-(j##x+1))*PI/(m##x-1)) \
87  - cos((m##x-1-j##x)*PI/(m##x-1)))
88 
89  /* 2D */
90  for(r=0;r<REPEAT;++r) {
91  int mx = RESFAC*N, my = RESFAC*NY;
92  double x = (rand()/(double)RAND_MAX)*2-1,
93  y = (rand()/(double)RAND_MAX)*2-1;
94  GET_JF(x); GET_JF(y);
95 
96  if(r%256==0) {
97  for(i=0;i<NZ*NY*N;++i) p[i]=rand()/(double)RAND_MAX;
98 
99  lob_bnd_lin_2(lb, lb_N,N,mx, lb_NY,NY,my, p,NZ, work);
100  }
101 
102  if(r<3)
103  printf("x: %g <= %g <= %g, f = %g\n",
104  cos((mx-1-jx)*PI/(mx-1)), x, cos((mx-1-(jx+1))*PI/(mx-1)), fx),
105  printf("y: %g <= %g <= %g, f = %g\n",
106  cos((my-1-jy)*PI/(my-1)), y, cos((my-1-(jy+1))*PI/(my-1)), fy);
107  lag_N (ld_N ,ld_N +N ,N ,0,x);
108  lag_NY(ld_NY,ld_NY+NY,NY,0,y);
109 
110  for(i=0;i<NZ;++i) {
111  double lo = (1-fx)*(1-fy)*lb[((i*mx+jx )*my+jy )*2 ]
112  + fx *(1-fy)*lb[((i*mx+jx+1)*my+jy )*2 ]
113  + (1-fx)* fy *lb[((i*mx+jx )*my+jy+1)*2 ]
114  + fx * fy *lb[((i*mx+jx+1)*my+jy+1)*2 ],
115  up = (1-fx)*(1-fy)*lb[((i*mx+jx )*my+jy )*2+1]
116  + fx *(1-fy)*lb[((i*mx+jx+1)*my+jy )*2+1]
117  + (1-fx)* fy *lb[((i*mx+jx )*my+jy+1)*2+1]
118  + fx * fy *lb[((i*mx+jx+1)*my+jy+1)*2+1],
119  pxy = tensor_i2(ld_N,N, ld_NY,NY, p+i*N*NY, work);
120  if(r<3 || pxy < lo || up < pxy)
121  printf("p_%02d(%g,%g) = %g in [%g,%g]\n",i,x,y,pxy,lo,up);
122  if(pxy<lo || up<pxy) {failure=1; break;}
123  }
124  if(i!=NZ) break;
125 
126  }
127 
128  /* 3D */
129  for(r=0;r<REPEAT;++r) {
130  int mx = RESFAC*N, my = RESFAC*NY, mz = RESFAC*NZ;
131  double x = (rand()/(double)RAND_MAX)*2-1,
132  y = (rand()/(double)RAND_MAX)*2-1,
133  z = (rand()/(double)RAND_MAX)*2-1;
134  GET_JF(x); GET_JF(y); GET_JF(z);
135  if(r%256==0) {
136  for(i=0;i<NZ*NY*N;++i) p[i]=rand()/(double)RAND_MAX;
137 
138  lob_bnd_lin_3(lb, lb_N,N,mx, lb_NY,NY,my, lb_NZ,NZ,mz, p,1, work);
139  }
140 
141  if(r<3)
142  printf("x: %g <= %g <= %g, f = %g\n",
143  cos((mx-1-jx)*PI/(mx-1)), x, cos((mx-1-(jx+1))*PI/(mx-1)), fx),
144  printf("y: %g <= %g <= %g, f = %g\n",
145  cos((my-1-jy)*PI/(my-1)), y, cos((my-1-(jy+1))*PI/(my-1)), fy),
146  printf("z: %g <= %g <= %g, f = %g\n",
147  cos((mz-1-jz)*PI/(mz-1)), z, cos((mz-1-(jz+1))*PI/(mz-1)), fz);
148  lag_N (ld_N ,ld_N +N ,N ,0,x);
149  lag_NY(ld_NY,ld_NY+NY,NY,0,y);
150  lag_NZ(ld_NZ,ld_NZ+NZ,NZ,0,z);
151 
152  {
153  double lo =
154  + (1-fx)*(1-fy)*(1-fz)*lb[(((jx )*my+jy )*mz+jz )*2 ]
155  + fx *(1-fy)*(1-fz)*lb[(((jx+1)*my+jy )*mz+jz )*2 ]
156  + (1-fx)* fy *(1-fz)*lb[(((jx )*my+jy+1)*mz+jz )*2 ]
157  + fx * fy *(1-fz)*lb[(((jx+1)*my+jy+1)*mz+jz )*2 ]
158  + (1-fx)*(1-fy)* fz *lb[(((jx )*my+jy )*mz+jz+1)*2 ]
159  + fx *(1-fy)* fz *lb[(((jx+1)*my+jy )*mz+jz+1)*2 ]
160  + (1-fx)* fy * fz *lb[(((jx )*my+jy+1)*mz+jz+1)*2 ]
161  + fx * fy * fz *lb[(((jx+1)*my+jy+1)*mz+jz+1)*2 ],
162  up =
163  + (1-fx)*(1-fy)*(1-fz)*lb[(((jx )*my+jy )*mz+jz )*2+1]
164  + fx *(1-fy)*(1-fz)*lb[(((jx+1)*my+jy )*mz+jz )*2+1]
165  + (1-fx)* fy *(1-fz)*lb[(((jx )*my+jy+1)*mz+jz )*2+1]
166  + fx * fy *(1-fz)*lb[(((jx+1)*my+jy+1)*mz+jz )*2+1]
167  + (1-fx)*(1-fy)* fz *lb[(((jx )*my+jy )*mz+jz+1)*2+1]
168  + fx *(1-fy)* fz *lb[(((jx+1)*my+jy )*mz+jz+1)*2+1]
169  + (1-fx)* fy * fz *lb[(((jx )*my+jy+1)*mz+jz+1)*2+1]
170  + fx * fy * fz *lb[(((jx+1)*my+jy+1)*mz+jz+1)*2+1],
171  pxyz = tensor_i3(ld_N,N, ld_NY,NY, ld_NZ,NZ, p, work);
172  if(r<3 || pxyz < lo || up < pxyz)
173  printf("p(%g,%g,%g) = %g in [%g,%g]\n",x,y,z,pxyz,lo,up);
174  if(pxyz<lo || up<pxyz) failure=1;
175  }
176  if(failure) break;
177 
178  }
179 
180  free(lb_NZ), free(lb_NY), free(lb_N), free(ld_NZ), free(ld_NY), free(ld_N);
181 
182  printf("Tests %s\n", failure?"failed":"successful");
183 
184  return failure;
185 }
main
int main()
Definition: lob_bnd_test.c:25
uint
#define uint
Definition: types.h:70
tensor_i2
static double tensor_i2(const double *Jr, uint nr, const double *Js, uint ns, const double *u, double *work)
Definition: tensor.h:87
tmalloc
#define tmalloc(type, count)
Definition: mem.h:91
lagrange_fun
void lagrange_fun(double *restrict p, double *restrict data, unsigned n, int d, double x)
Definition: poly.c:20
RESFAC
#define RESFAC
Definition: lob_bnd_test.c:16
lob_bnd_lin_2
#define lob_bnd_lin_2
Definition: lob_bnd.c:15
x
#define x
GET_JF
#define GET_JF(x)
NZ
#define NZ
Definition: lob_bnd_test.c:19
tensor_dot
#define tensor_dot
Definition: tensor.c:7
gll_lag_setup
#define gll_lag_setup
Definition: poly.c:18
p
p
Definition: xxt_test2.m:1
gll_lag_size
#define gll_lag_size
Definition: poly.c:17
lob_bnd_lin_3
#define lob_bnd_lin_3
Definition: lob_bnd.c:16
PI
#define PI
Definition: lob_bnd_test.c:22
lob_bnd.h
i
for i
Definition: xxt_test.m:74
lob_bnd_size
static unsigned lob_bnd_size(unsigned n, unsigned m)
Definition: lob_bnd.h:64
NY
#define NY
Definition: lob_bnd_test.c:18
tensor_i3
static double tensor_i3(const double *Jr, uint nr, const double *Js, uint ns, const double *Jt, uint nt, const double *u, double *work)
Definition: tensor.h:96
lob_bnd_setup
#define lob_bnd_setup
Definition: lob_bnd.c:13
lob_bnd_lin_1
#define lob_bnd_lin_1
Definition: lob_bnd.c:14
REPEAT
#define REPEAT
Definition: lob_bnd_test.c:20
work
static double work[TNR *NS]
Definition: findpts_el_2_test2.c:50
name.h
establishes some macros to establish naming conventions
y
static double y[NR *NS *NT *N]
Definition: obbox_test.c:31
N
#define N
Definition: lob_bnd_test.c:17
z
static double z[NR *NS *NT *N]
Definition: obbox_test.c:31
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8