file: Hlib/src/Basis.c


/**************************************************************************/
// // 
// Author: T.Warburton // // Design: T.Warburton && S.Sherwin     // 
// Date : 12/4/96                                                         
// // // // Copyright notice: This code shall not
be replicated or used without // // the permission of the author. // //
// /**************************************************************************/
#include <stdio.h>
#include <math.h>
#include <veclib.h>
#include <polylib.h> #include "hotel.h"
/* zeros weights info */ typedef struct 
zwinfo
{ int q;
 /* number of quadrature points */ char dir;
 /* direction of
zeros/weights either 'a', 'b' or 'c' */ double *z;
 /* zeros */ double *w;

/* weigths */ struct zwinfo *next;
 } 
ZWinfo;

/*-------------------------------------------------------------------*
* link list for zeros and weights for jacobi polynomial in [0,1] * * at
the gauss points. Note: alpha > -1 , beta > -1 * *-------------------------------------------------------------------*/
static ZWinfo *zwinf,*zwbase;
 static ZWinfo
*addzw(int, char);
 void 
getzw(int
q, double **z, double **w, char dir){ /* check link list */ for(zwinf
= zwbase;
 zwinf;
 zwinf = zwinf->next) if(zwinf->q
== q) if(zwinf->dir == dir){ *z = zwinf->z;

*w = zwinf->w;
 return;
 } /* else add new zeros and weights */ zwinf = zwbase;

zwbase = addzw(q,dir);
 zwbase->next = zwinf;
 *z =
zwbase->z;
 *w = zwbase->w;
 return;
 } /*--------------------------------------------------------------------------*
* This function gets the zeros and weights in [-1,1]. In the 'a' direction
* * alpha = 0, beta = 0. In the 'b' alpha = 1, beta = 0 and in the 'c'
* * direction alpha = 2, beta = 0. * * Note: that the 'b','c' direction
weights are scaled by (1/2)^alpha * * to be consistent with the jacobean
d(r,s)/d(t,s) * *--------------------------------------------------------------------------*/
static ZWinfo *addzw(int n,
char dir){ ZWinfo *zw = (ZWinfo
*)calloc(1,sizeof(ZWinfo));
 zw->q = n;
 zw->dir =
dir;
 zw->z = dvector(0,n-1);
 zw->w = dvector(0,n-1);
 #ifndef COMPRESS switch(dir){
case 'a': zwgll(zw->z,zw->w,n);
 break;

case 'b': zwgrj(zw->z,zw->w,n,1.0,0.0);

dscal(n,0.5,zw->w,1);
 break;
 case 'c': zwgrj(zw->z,zw->w,n,2.0,0.0);

dscal(n,0.25,zw->w,1);
 break;
 case 'g': /* just Gauss points */ zwgl
(zw->z,zw->w,n);
 break;
 case 'h': /* Gauss points with 1,0 weight */
zwgj (zw->z,zw->w,n,1.0,0.0);
 dscal(n,0.5,zw->w,1);

break;
 case 'i': zwgj (zw->z,zw->w,n,2.0,0.0);

dscal(n,0.25,zw->w,1);
 break;
 default: fprintf(stderr,"addzw - incorect
direction specified");
 break;
 } #else int i;
 switch(dir){ case 'a': zwgll(zw->z,zw->w,n);

break;
 case 'b': zwgrj(zw->z,zw->w,n,0.0,0.0);

for(i = 0;
 i <
 n;
 ++i) zw->w[i] *= (1-zw->z[i])*0.5;
 break;
 case 'c':
zwgrj(zw->z,zw->w,n,0.0,0.0);
 for(i
= 0;
 i <
 n;
 ++i){ zw->w[i] *= (1-zw->z[i])*0.5;
 zw->w[i] *= (1-zw->z[i])*0.5;

} break;
 case 'g': /* just Gauss points */ zwgl
(zw->z,zw->w,n);
 break;
 case 'h': /* Gauss points with 1,0 weight */
zwgj (zw->z,zw->w,n,0.0,0.0);
 for(i =
0;
 i <
 n;
 ++i) zw->w[i] *= (1-zw->z[i])*0.5;
 break;
 case 'i': zwgj
(zw->z,zw->w,n,0.0,0.0);
 for(i = 0;
 i <
 n;
 ++i){ zw->w[i] *= (1-zw->z[i])*0.5;

zw->w[i] *= (1-zw->z[i])*0.5;
 } break;
 case 'j': zwgrj(zw->z,zw->w,n,0.0,0.0);

break;
 default: fprintf(stderr,"addzw - incorect direction specified");

break;
 } #endif return zw;
 } typedef struct 
iminfo
{ int q1;
 /* number of quadrature points */ int q2;
 InterpDir
dir;
 /* direction of zeros/weights either 'a', 'b' or 'c' */ double **im;

/* zeros */ struct iminfo *next;
 } 
IMinfo;

/*---------------------------------------------------------------------*
* link list for interpolation matrices: * * InterpDir
: a2b go from gll points to grj_(1,0) points * * b2a go from grj_(1,0)
points to gll points * * b2c go from grj_(1,0) points to grj_(2,0) points
* * c2b go from grj_(2,0) points to grj_(2,0) points * *---------------------------------------------------------------------*/
/*---------------------------------------------------------------------*
* Matrix storage for interpolation matrices: * * InterpDir
: a2a go from gll points to gll points * * a2b go from gll points to grj_(1,0)
points * * a2g go from gll points to gl points * * b2a go from grj_(1,0)
points to gll points * * b2g go from grj_(1,0) points to glj_(2,0) points
* * b2c go from grj_(1,0) points to grj_(2,0) points * * b2g go from grj_(1,0)
points to gl points * * c2b go from grj_(2,0) points to grj_(1,0) points
* * c2c go from grj_(2,0) points to grj_(2,0) points * * c2g go from grj_(2,0)
points to gl points * *---------------------------------------------------------------------*/
static IMinfo *iminf,*imbase;
 static IMinfo
*addim(int,int,InterpDir);

//OVERLOAD CALL: addim: Basis.c(?), Basis.c(?)
void 
getim(int q1, int q2, double ***im, InterpDir
dir){ /* check link list */ for(iminf = imbase;
 iminf;
 iminf = iminf->next)
if(iminf->q1 == q1) if(iminf->q2
== q2) if(iminf->dir == dir){ *im = iminf->im;

return;
 } /* else add new interpolation matrix */ iminf = imbase;
 imbase
= addim(q1,q2,dir);
 //OVERLOAD CALL: addim: Basis.c(?),
Basis.c(?) imbase->next = iminf;
 *im =
imbase->im;
 return;
 } #ifndef COMPRESS static IMinfo
*addim(int n1, int n2, InterpDir
dir){ IMinfo *im = (IMinfo
*)malloc(sizeof(IMinfo));
 double *z1,*z2,*w;
 im->q1
= n1;
 im->q2 = n2;
 im->dir = dir;
 im->im = dmatrix(0,n2-1,0,n1-1);
 switch(dir){
case a2a: /* interpolate a to a */ { getzw(n1,&
z1,&
w,'a');

getzw(n2,&
z2,&
w,'a');
 igllm(im->im,z1,z2,n1,n2);

break;
 } case a2b: /* interpolate a to b */ { getzw(n1,&
z1,&
w,'a');

getzw(n2,&
z2,&
w,'b');
 igllm(im->im,z1,z2,n1,n2);

break;
 } case a2g: /* interpolate a to g (gauss) */ { getzw(n1,&
z1,&
w,'a');

getzw(n2,&
z2,&
w,'g');
 igllm(im->im,z1,z2,n1,n2);

break;
 } case b2a: /* interpolate b to a */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'a');
 igrjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } case b2b: /* interpolate b to b */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'b');
 igrjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } case b2c: /* interpolate b to c */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'c');
 igrjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } case b2g: /* interpolate b to g (gauss) */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'g');
 igrjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } case b2h: /* interpolate b (Radau) to h (gauss) */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'h');
 igrjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } case c2b: /* interpolate c to b */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'b');
 igrjm(im->im,z1,z2,n1,n2,2.0,0.0);

break;
 } case c2c: /* interpolate c to c */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'c');
 igrjm(im->im,z1,z2,n1,n2,2.0,0.0);

break;
 } case c2g: /* interpolate c to g (gauss) */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'g');
 igrjm(im->im,z1,z2,n1,n2,2.0,0.0);

break;
 } case c2h: /* interpolate c (Radau1) to h (gauss) */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'h');
 igrjm(im->im,z1,z2,n1,n2,2.0,0.0);

break;
 } case g2g: /* interpolate g (gauss) to g (gauss) */ { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'g');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case g2a: /* interpolate g (gauss) to g (gauss) */ { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'a');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case g2b: { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'b');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case h2b: { getzw(n1,&
z1,&
w,'h');

getzw(n2,&
z2,&
w,'b');
 igjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } case h2c: { getzw(n1,&
z1,&
w,'h');

getzw(n2,&
z2,&
w,'c');
 igjm(im->im,z1,z2,n1,n2,1.0,0.0);

break;
 } default: fprintf(stderr,"addim - incorect direction specified");

break;
 } return im;
 } #else static IMinfo *addim(int
n1, int n2, InterpDir dir){ IMinfo
*im = (IMinfo *)malloc(sizeof(IMinfo));

double *z1,*z2,*w;
 im->q1 = n1;
 im->q2 = n2;
 im->dir = dir;
 im->im = dmatrix(0,n2-1,0,n1-1);

switch(dir){ case a2a: /* interpolate a to a */ { getzw(n1,&
z1,&
w,'a');

getzw(n2,&
z2,&
w,'a');
 igllm(im->im,z1,z2,n1,n2);

break;
 } case a2b: /* interpolate a to b */ { getzw(n1,&
z1,&
w,'a');

getzw(n2,&
z2,&
w,'b');
 igllm(im->im,z1,z2,n1,n2);

break;
 } case a2g: /* interpolate a to g (gauss) */ { getzw(n1,&
z1,&
w,'a');

getzw(n2,&
z2,&
w,'g');
 igllm(im->im,z1,z2,n1,n2);

break;
 } case b2a: /* interpolate b to a */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'a');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case b2b: /* interpolate b to b */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'b');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case b2c: /* interpolate b to c */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'c');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case b2g: /* interpolate b to g (gauss) */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'g');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case b2h: /* interpolate b (Radau) to h (gauss) */ { getzw(n1,&
z1,&
w,'b');

getzw(n2,&
z2,&
w,'h');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case c2b: /* interpolate c to b */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'b');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case c2c: /* interpolate c to c */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'c');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case c2g: /* interpolate c to g (gauss) */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'g');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case c2h: /* interpolate c (Radau1) to h (gauss) */ { getzw(n1,&
z1,&
w,'c');

getzw(n2,&
z2,&
w,'h');
 igrjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case g2g: /* interpolate g (gauss) to g (gauss) */ { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'g');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case g2a: /* interpolate g (gauss) to g (gauss) */ { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'a');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case g2b: { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'b');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case g2c: { getzw(n1,&
z1,&
w,'g');

getzw(n2,&
z2,&
w,'c');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case h2b: { getzw(n1,&
z1,&
w,'h');

getzw(n2,&
z2,&
w,'b');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } case h2c: { getzw(n1,&
z1,&
w,'h');

getzw(n2,&
z2,&
w,'c');
 igjm(im->im,z1,z2,n1,n2,0.0,0.0);

break;
 } default: fprintf(stderr,"addim - incorect direction specified");

break;
 } return im;
 } #endif // L_2 normalise a void 
normalise_mode(int
q, double *v, char dir){ if(option("TWOTWO")){
double *z, *w, d, hd;
 double *tmp = dvector(0, q-1);
 double *tmp1 = dvector(0,
q-1);
 getzw(q, &
z, &
w, dir);
 dvmul(q, v, 1,
v, 1, tmp, 1);
 d = ddot(q, tmp, 1, w, 1);
 hd = sqrt(d);
 dscal(q, 1./hd,
v, 1);
 free(tmp);
 free(tmp1);
 } return;
 } void 
normalised_iprod_1d(int
q, double *h, int L, double *m, double *res, double *invnorm, char
dir){ double *wk = dvector(0, q-1);
 double *z, *w;
 getzw(q,
&
z, &
w, dir);
 dvmul(q, w, 1, h, 1, wk, 1);
 dgemv('T', q, L, 1.,
m, q, wk, 1, 0.0, res, 1);
 dvmul(L, invnorm, 1, res, 1, res, 1);
 free(wk);

} // assumes 'a'x'a' // interior modes.. void 
quad_normalised_iprod_2d(int
q, double *h, int L, double *m, double *res, double *invnorm) { int
i;
 double *wk = dvector(0, q*q-1);
 double *wk1 = dvector(0, q*q-1);
 double
*z, *w;
 getzw(q, &
z, &
w, 'a');
 // multiply
field by weights for(i = 0;
 i <
 q;
 ++i) dvmul(q, w, 1, h+i*q,1, wk+i*q,1);

for(i = 0;
 i <
 q;
 ++i) dvmul(q, w, 1, wk+i,q, wk+i,q);
 dgemm('T','N',
L, q, q, 1.0, m, q, wk, q, 0.0, wk1, L);
 // 'a' direction dgemm('N','N',
L, L, q, 1.0, wk, L, m, q, 0.0, res, L);
 // 'b' direction dvmul(L*L, invnorm,
1, res, 1, res, 1);
 free(wk);
 free(wk1);
 } // assumes 'a'x'b' // interior
modes.. void 
tri_normalised_iprod_2d(int
qa, char dira, char dirb, int qb, double *h, int L, double *ma, double
*mb, double *res, double *invnorm) { int i;
 double *wk = dvector(0,
max(qa, qb)*max(qa,qb)
-1);
 double *wk1 = dvector(0, max(qa, qb)*max(qa,qb)
-1);
 double *z, *w, *wk1_a, *res_a;
 // multiply field by weights getzw(qa,
&
z, &
w, dira);
 for(i = 0;
 i <
 qb;
 ++i) dvmul(qa, w,1, h+i*qa,
1, wk+i*qa, 1);
 getzw(qb, &
z, &
w, dirb);
 for(i
= 0;
 i <
 qb;
 ++i) dscal(qa, w[i], wk+i*qa, 1);
 // 'a' direction integrals
dgemm('T','N',qb, L, qa, 1.0, wk, qa, ma, qa, 0.0, wk1, qb);
 // 'b' direction
integrals wk1_a = wk1;
 res_a = res;
 // +1 -- ignore first edge mode
for(i = 0;
 i <
 L;
 mb += qb*(L-i+1), wk1_a += qb, res_a+=L-i, ++i) dgemv('T',
qb, L-i, 1.0, mb, qb, wk1_a, 1, 0.0, res_a, 1);
 dvmul(L*(L+1)/2, invnorm,
1, res, 1, res, 1);
 free(wk);
 free(wk1);
 } static double 
norm(int
np, double *d, char dir){ double fac =0.;
 double *tmp = dvector(0,
np-1);
 double *z, *w;
 getzw(np, &
z, &
w, dir);

dvmul(np, d, 1, d, 1, tmp, 1);
 fac = sqrt(ddot(np, tmp, 1, w, 1));
 free(tmp);

return fac;
 } /* interpolation matrix */ static double ***BasisA_G, ***BasisA_GL,
***BasisA_GR;
 static double ****BasisB_G, ****BasisB_GR;
 static double
****BasisC_GR, ****BasisC_G;
 /* This function initialises the orthogonal
basis statically defined matrices based
on QGmax. These arrays may be accessed form
[1-QGmax]*/ void 
init_ortho_basis(void){
BasisA_G = (double ***) calloc(QGmax, sizeof(double **));
 BasisA_G -= 1;

BasisA_GL = (double ***) calloc(QGmax, sizeof(double **));
 BasisA_GL-=
1;
 BasisA_GR = (double ***) calloc(QGmax, sizeof(double **));
 BasisA_GR-=
1;
 BasisB_G = (double ****)calloc(QGmax, sizeof(double ***));
 BasisB_G
-= 1;
 BasisB_GR = (double ****)calloc(QGmax, sizeof(double ***));
 BasisB_GR-=
1;
 BasisC_GR = (double ****)calloc(QGmax, sizeof(double ***));
 BasisC_GR-=
1;
 BasisC_G = (double ****)calloc(QGmax, sizeof(double ***));
 BasisC_G-=
1;
 } /* Gets the values of "mode" (0..LGmax)
degree Leg. poly on "size" Gauss points, stores
the result in BasisA_G[size][mode] array
*/ void 
get_moda_G (int size, double ***mat){
if (!BasisA_G[size]){ double *w, *z;
 int
i;
 BasisA_G[size] = dmatrix(0,LGmax-1,0,size-1);
 getzw(size,&
z,&
w,'g');

for (i=0;
 i <
 LGmax;
 ++i){ jacobf(size,
z, BasisA_G[size][i], i, 0.0, 0.0);
 /* normalise so that (g,g)=1*/ #ifndef
NONORM dscal(size,sqrt(0.5*(2.0*i+1.0)),BasisA_G[size][i],1);
 #endif }
} *mat = BasisA_G[size];
 return;
 } /* computes the values of a_mode[mode_a][j],
mode_a=0..LGmax, j=0..size-1, at size Gauss-Lob. points, stores the result
in BasisA_GL[size][mode_a] array */ void 
get_moda_GL
(int size, double ***mat){ if (!BasisA_GL[size]){
double *w, *z;
 int i;
 double *tmp = dvector(0, size-1);
 BasisA_GL[size]
= dmatrix(0,LGmax-1,0,size-1);
 getzw(size,&
z,&
w,'a');

for (i=0;
 i <
 LGmax;
 ++i){ jacobf(size,
z, BasisA_GL[size][i], i, 0.0, 0.0);
 #ifndef NONORM /* normalise so that
(g,g)=1*/ dscal(size,sqrt(0.5*(2.0*i+1.0)),BasisA_GL[size][i],1);
 #endif
} free(tmp);
 } *mat = BasisA_GL[size];
 return;
 } /* computes the values
of a_mode[mode_a][j], mode_a=0..LGmax, j=0..size-1, at size Gauss-Lob.
points, stores the result in BasisA_GL[size][mode_a] array */ void 
get_moda_GR
(int size, double ***mat){ if (!BasisA_GR[size]){
double *w, *z;
 int i;
 double *tmp = dvector(0, size-1);
 BasisA_GR[size]
= dmatrix(0,LGmax-1,0,size-1);
 getzw(size,&
z,&
w,'b');

for (i=0;
 i <
 LGmax;
 ++i){ jacobf(size,
z, BasisA_GR[size][i], i, 0.0, 0.0);
 #ifndef NONORM /* normalise so that
(g,g)=1*/ dscal(size,sqrt(0.5*(2.0*i+1.0)),BasisA_GR[size][i],1);
 #endif
} free(tmp);
 } *mat = BasisA_GR[size];
 return;
 } /* computes the values
of b dependent modes at either the gauss or gauss radau points depending
on the function and stores the result in BasisB_G
array */ static double ***get_Bmodes(int size,
double *z);
 static double ***get_Cmodes(int size,
double *z);
 void 
get_modb_G (int size,
double ****mat){ if (!BasisB_G[size]){
double *w,*z;
 getzw (size,&
z,&
w,'h');
 BasisB_G[size]
= get_Bmodes(size,z);
 } *mat = BasisB_G[size];

return;
 } void 
get_modb_GR (int size,
double ****mat){ if (!BasisB_GR[size]){
double *w,*z;
 getzw (size,&
z,&
w,'b');
 BasisB_GR[size]
= get_Bmodes(size,z);
 } *mat = BasisB_GR[size];

return;
 } static double ***get_Bmodes(int size,
double *z){ double ***mat;
 double fac;
 int i,j,cnt;
 /* declare matrix
*/ mat = (double ***)malloc(LGmax*sizeof(double **));
 mat[0] = (double
**) malloc(LGmax*(LGmax+1)/2*sizeof(double *));
 mat[0][0] = (double *)
malloc(LGmax*(LGmax+1)/2*size*sizeof(double));
 for(i = 0,cnt=0;
 i <

LGmax;
 cnt+=LGmax-i,++i) mat[i] = mat[0] + cnt;
 for(i = 0,cnt=0;
 i <

LGmax;
 ++i) for(j = 0;
 j <
 LGmax-i;
 ++j,++cnt) mat[i][j] = mat[0][0]
+ size*cnt;
 dfill(size,1.0,mat[0][0],1);
 for(j = 1;
 j <
 LGmax;
 ++j){
jacobf(size,z,mat[0][j],j,1.0,0.0);

dvmul(size,mat[0][0],1,mat[0][j],1,mat[0][j],1);
 } if(LGmax
> 1){ for(i = 0;
 i <
 size;
 ++i) mat[1][0][i] = 1.0 - z[i];
 for(j = 1;

j <
 LGmax-1;
 ++j){ jacobf(size,z,mat[1][j],j,3.0,0.0);

dvmul(size,mat[1][0],1,mat[1][j],1,mat[1][j],1);
 } for(i = 2;
 i <
 LGmax;

++i){ dvmul(size,mat[1][0],1,mat[i-1][0],1,mat[i][0],1);
 for(j = 1;
 j <

LGmax-i;
 ++j){ jacobf(size,z,mat[i][j],j,2.0*i+1.0,0.0);

dvmul(size,mat[i][0],1,mat[i][j],1,mat[i][j],1);
 } } } #ifndef NONORM /*
normalise (recalling factor of
2 for weight (1-b)/2 */ for(i =0, fac = 1.0;
 i <
 LGmax;
 ++i, fac *=
0.25) for(j = 0;
 j <
 LGmax-i;
 ++j) dscal(size,sqrt((i+j+1.0)*fac),mat[i][j],1);

#endif return mat;
 } void 
get_modc_GR
(int size, double ****mat){ if (!BasisC_GR[size]){
double *w,*z;
 getzw (size,&
z,&
w,'c');
 BasisC_GR[size]
= get_Cmodes(size,z);
 } *mat = BasisC_GR[size];

return;
 } void 
get_modc_G (int size, double
****mat){ if (!BasisC_G[size]){ double
*w,*z;
 getzw (size,&
z,&
w,'h');
 BasisC_G[size]
= get_Cmodes(size,z);
 } *mat = BasisC_G[size];

return;
 } static double ***get_Cmodes(int size,
double *z){ double ***mat;
 double fac;
 int i,j,cnt;
 /* declare matrix
*/ mat = (double ***)malloc(LGmax*sizeof(double **));
 mat[0] = (double
**) malloc(LGmax*(LGmax+1)/2*sizeof(double *));
 mat[0][0] = (double *)
malloc(LGmax*(LGmax+1)/2*size*sizeof(double));
 for(i = 0,cnt=0;
 i <

LGmax;
 cnt+=LGmax-i,++i) mat[i] = mat[0] + cnt;
 for(i = 0,cnt=0;
 i <

LGmax;
 ++i) for(j = 0;
 j <
 LGmax-i;
 ++j,++cnt) mat[i][j] = mat[0][0]
+ size*cnt;
 dfill(size,1.0,mat[0][0],1);
 for(j = 1;
 j <
 LGmax;
 ++j){
jacobf(size,z,mat[0][j],j,2.0,0.0);

dvmul(size,mat[0][0],1,mat[0][j],1,mat[0][j],1);
 } if(LGmax
> 1){ for(i = 0;
 i <
 size;
 ++i) mat[1][0][i] = 1.0 - z[i];
 for(j = 1;

j <
 LGmax-1;
 ++j){ jacobf(size,z,mat[1][j],j,4.0,0.0);

dvmul(size,mat[1][0],1,mat[1][j],1,mat[1][j],1);
 } for(i = 2;
 i <
 LGmax;

++i){ dvmul(size,mat[1][0],1,mat[i-1][0],1,mat[i][0],1);
 for(j = 1;
 j <

LGmax-i;
 ++j){ jacobf(size,z,mat[i][j],j,2.0*i+2.0,0.0);

dvmul(size,mat[i][0],1,mat[i][j],1,mat[i][j],1);
 } } } #ifndef NONORM /*
normalise (recalling factor of
2 for weight (1-b)/2 */ for(i =0, fac = 1.;
 i <
 LGmax;
 ++i, fac *= 0.25)
for(j = 0;
 j <
 LGmax-i;
 ++j) dscal(size,sqrt((i+j+1.5)*fac),mat[i][j],1);

#endif return mat;
 } void Tri_reset_basis(void);

void Quad_reset_basis(void);
 void Tet_reset_basis(void);

void Pyr_reset_basis(void);
 void Prism_reset_basis(void);

void Hex_reset_basis(void);
 void 
reset_bases(){
Tri_reset_basis();
 Quad_reset_basis();

Tet_reset_basis();
 Pyr_reset_basis();

Prism_reset_basis();
 Hex_reset_basis();

} static Basis ****Tri_bbase=NULL;
 static
Basis ****Quad_bbase=NULL;
 static Basis
*Tet_bbinf=NULL,*Tet_bbase=NULL;
 static Basis
*Pyr_bbinf=NULL,*Pyr_bbase=NULL;
 static Basis
*Prism_bbinf=NULL,*Prism_bbase=NULL;
 static Basis
*Hex_bbinf=NULL,*Hex_bbase=NULL;
 static Basis
*Tri_addbase(int L, int qa, int qb, int qc);

static Basis *Quad_addbase(int
L, int qa, int qb, int qc);
 static Basis
*Tet_addbase(int L, int qa, int qb, int qc);

static Basis *Pyr_addbase(int
L, int qa, int qb, int qc);
 static Basis
*Prism_addbase(int L, int qa, int qb, int
qc);
 static Basis *Hex_addbase(int
L, int qa, int qb, int qc);
 static void Tri_normalise(Basis
*b, int qa, int qb, int qc);
 static void Quad_normalise(Basis
*b, char type );
 static void Tet_normalise(Basis
*b, int , int , int );
 static void Pyr_normalise(Basis
*, int , int , int );
 static void Prism_normalise(Basis
*, int , int , int );
 static void Hex_normalise(Basis
*, int , int , int );
 /* Function name: Element::getbasis
Function Purpose: Function Notes: */ Basis
*Tri::getbasis(){ int
i,j;
 if(!Tri_bbase){ Tri_bbase = (Basis****)
calloc(LGmax+1, sizeof(Basis***));
 for(i=0;
i<
LGmax+1;
++i){
Tri_bbase[i] = (Basis***) calloc(QGmax+1,
sizeof(Basis**));
 for(j=0;
j<
QGmax+1;
++j)
Tri_bbase[i][j] = (Basis**) calloc(QGmax+1,
sizeof(Basis*));
 } } if(!Tri_bbase[lmax][qa][qb])
Tri_bbase[lmax][qa][qb] = Tri_addbase(lmax,qa,qb,qc);

return Tri_bbase[lmax][qa][qb];
 } Basis
*Quad::getbasis(){
int i,j;
 if(!Quad_bbase){ Quad_bbase
= (Basis****) calloc(LGmax+1, sizeof(Basis***));

for(i=0;
i<
LGmax+1;
++i){ Quad_bbase[i] = (Basis***)
calloc(QGmax+1, sizeof(Basis**));
 for(j=0;
j<
QGmax+1;
++j)
Quad_bbase[i][j] = (Basis**) calloc(QGmax+1,
sizeof(Basis*));
 } } if(!Quad_bbase[lmax][qa][qb])
Quad_bbase[lmax][qa][qb] = Quad_addbase(lmax,qa,qb,qc);

return Quad_bbase[lmax][qa][qb];
 } Basis
*Tet::getbasis(){ /*
check link list */ for(Tet_bbinf = Tet_bbase;
 Tet_bbinf;
 Tet_bbinf = Tet_bbinf->next)
if(Tet_bbinf->id == lmax){ return Tet_bbinf;

} /* else add new zeros and weights */ Tet_bbinf = Tet_bbase;
 Tet_bbase
= Tet_addbase(lmax,qa,qb,qc);
 Tet_bbase->next
= Tet_bbinf;
 return Tet_bbase;
 } Basis
*Pyr::getbasis(){ /*
check link list */ for(Pyr_bbinf = Pyr_bbase;
 Pyr_bbinf;
 Pyr_bbinf = Pyr_bbinf->next)
if(Pyr_bbinf->id == lmax){ return Pyr_bbinf;

} /* else add new zeros and weights */ Pyr_bbinf = Pyr_bbase;
 Pyr_bbase
= Pyr_addbase(lmax,qa,qb,qc);
 Pyr_bbase->next
= Pyr_bbinf;
 return Pyr_bbase;
 } Basis
*Prism::getbasis(){
/* check link list */ for(Prism_bbinf = Prism_bbase;
 Prism_bbinf;
 Prism_bbinf
= Prism_bbinf->next) if(Prism_bbinf->id
== lmax){ return Prism_bbinf;
 } /* else add new zeros and weights */ Prism_bbinf
= Prism_bbase;
 Prism_bbase = Prism_addbase(lmax,qa,qb,qc);

Prism_bbase->next = Prism_bbinf;
 return Prism_bbase;
 } Basis
*Hex::getbasis(){ /*
check link list */ for(Hex_bbinf = Hex_bbase;
 Hex_bbinf;
 Hex_bbinf = Hex_bbinf->next)
if(Hex_bbinf->id == lmax){ return Hex_bbinf;

} /* else add new zeros and weights */ Hex_bbinf = Hex_bbase;
 Hex_bbase
= Hex_addbase(lmax,qa,qb,qc);
 Hex_bbase->next
= Hex_bbinf;
 return Hex_bbase;
 } Basis
*Element::getbasis(){return
(Basis*)NULL;
} static Basis
*Tri_addbase(int L, int qa, int qb, int qc){
Basis *b = Tri_mem_base(L,qa,qb,qc);

/* set up storage */ Tri_mem_modes(b);
 /*
set up vertices */ Tri_set_vertices(b->vert,qa,'a');

Tri_set_vertices(b->vert,qb,'b');
 if(L>2){
Tri_set_edges(b,qa,'a');
 Tri_set_edges(b,qb,'b');

} if(L>3) Tri_set_faces(b,qb,'b');

if(option("Normalise")) Tri_normalise(b,qa,qb,qc);

/* normalise basis so that L2 of each
mode = 1 */ return b;
 } static Basis
*Quad_addbase(int L, int qa, int qb, int qc){
Basis *b = Quad_mem_base(L,qa,qb,qc);

/* set up storage */ Quad_mem_modes(b);
 /*
set up vertices */ Quad_set_vertices(b,qa,'a');

Quad_set_vertices(b,qb,'b');
 if(L>2){
Quad_set_edges(b,qa,'a');
 Quad_set_edges(b,qb,'b');

} if(L>2) Quad_set_faces(b,qb,'b');

/* normalise basis so that L2 of each
mode = 1 */ if(option("Normalise"))
Quad_normalise(b,'B');
 return b;
 } static
Basis *Tet_addbase(int
L, int qa, int qb, int qc){ Basis
*b = Tet_mem_base(L,qa,qb,qc);
 /* set up storage
*/ Tet_mem_modes(b);
 /* set up vertices */
Tet_set_vertices(b->vert,qa,'a');
 Tet_set_vertices(b->vert,qb,'b');

Tet_set_vertices(b->vert,qc,'c');
 if(L>2){
Tet_set_edges(b,qa,'a');
 Tet_set_edges(b,qb,'b');

Tet_set_edges(b,qc,'c');
 } if(L>3){
Tet_set_faces(b,qa,'a');
 Tet_set_faces(b,qb,'b');

Tet_set_faces(b,qc,'c');
 } if(option("Normalise"))
Tet_normalise(b,qa,qb,qc);
 /* normalise basis
so that L2 of each mode = 1 */ return
b;
 } static Basis *Pyr_addbase(int
L, int qa, int qb, int qc){ Basis
*b = Pyr_mem_base(L,qa,qb,qc);
 /* set up storage
*/ Pyr_mem_modes(b);
 /* set up vertices */
Pyr_set_vertices(b->vert,qa,'a');
 Pyr_set_vertices(b->vert,qb,'b');

Pyr_set_vertices(b->vert,qc,'c');
 if(L>2){
Pyr_set_edges(b,qa,'a');
 Pyr_set_edges(b,qa,'b');

Pyr_set_edges(b,qc,'c');
 } if(L>3){
Pyr_set_faces(b,qa,'a');
 Pyr_set_faces(b,qa,'b');

Pyr_set_faces(b,qc,'c');
 } if(option("Normalise"))
Pyr_normalise(b,qa,qb,qc);
 /* normalise basis
so that L2 of each mode = 1 */ return
b;
 } static Basis *Prism_addbase(int
L, int qa, int qb, int qc){ Basis
*b = Prism_mem_base(L,qa,qb,qc);
 /* set up
storage */ Prism_mem_modes(b);
 /* set up
vertices */ Prism_set_vertices(b->vert,qa,'a');

Prism_set_vertices(b->vert,qb,'b');
 Prism_set_vertices(b->vert,qc,'c');

if(L>2){ Prism_set_edges(b,qa,'a');

Prism_set_edges(b,qa,'b');
 Prism_set_edges(b,qc,'c');

} if(L>3){ Prism_set_faces(b,qa,'a');

Prism_set_faces(b,qa,'b');
 Prism_set_faces(b,qc,'c');

} if(option("Normalise")) Prism_normalise(b,qa,qb,qc);

/* normalise basis so that L2 of each
mode = 1 */ return b;
 } static Basis
*Hex_addbase(int L, int qa, int qb, int qc){
Basis *b = Hex_mem_base(L,qa,qb,qc);

/* set up storage */ Hex_mem_modes(b);
 /*
set up vertices */ Hex_set_vertices(b->vert,qa,'a');

if(L>2) Hex_set_edges(b,qa,'a');

if(L>2) Hex_set_faces(b,qa,'a');

if(option("Normalise")) Hex_normalise(b,qa,qb,qc);

/* normalise basis so that L2 of each
mode = 1 */ return b;
 } void 
Tri_reset_basis(void){
int i,j,k;
 if(Tri_bbase){ for(i=0;
i<
LGmax+1;
++i)
for(j=0;
j<
QGmax+1;
++j) for(k=0;
k<
QGmax+1;
++k) if(Tri_bbase[i][j][k]){
free(Tri_bbase[i][j][k]->vert[2].a);
 free(Tri_bbase[i][j][k]->vert[2].b);

free(Tri_bbase[i][j][k]->vert);
 free(Tri_bbase[i][j][k]->edge);
 free(Tri_bbase[i][j][k]->face);

free(Tri_bbase[i][j][k]);
 } for(i=0;
i<
LGmax+1;
++i){ for(j=0;
j<
QGmax+1;
++j)
free(Tri_bbase[i][j]);
 free(Tri_bbase[i]);
 } free(Tri_bbase);
 Tri_bbase
= NULL;
 } } void 
Quad_reset_basis(void){
int i,j,k;
 if(Quad_bbase){ for(i=0;
i<
LGmax+1;
++i)
for(j=0;
j<
QGmax+1;
++j) for(k=0;
k<
QGmax+1;
++k) if(Quad_bbase[i][j][k]){
free(Quad_bbase[i][j][k]->vert[1].a);
 free(Quad_bbase[i][j][k]->vert[2].b);

free(Quad_bbase[i][j][k]->vert);
 free(Quad_bbase[i][j][k]->edge);
 free(Quad_bbase[i][j][k]->face);

free(Quad_bbase[i][j][k]);
 } for(i=0;
i<
LGmax+1;
++i){ for(j=0;
j<
QGmax+1;
++j)
free(Quad_bbase[i][j]);
 free(Quad_bbase[i]);
 } free(Quad_bbase);
 Quad_bbase
= NULL;
 } } void 
Tet_reset_basis(void){
for(Tet_bbinf = Tet_bbase;
 Tet_bbinf;
 Tet_bbinf = Tet_bbinf->next){ free(Tet_bbinf->vert[2].a);

free(Tet_bbinf->vert[3].b);
 free(Tet_bbinf->vert[3].c);
 free(Tet_bbinf->vert);

free(Tet_bbinf->edge);
 free(Tet_bbinf->face);
 free(Tet_bbinf->intr);
 }
if(Tet_bbase) while(Tet_bbase){ Tet_bbinf
= Tet_bbase->next;
 free(Tet_bbase);
 Tet_bbase = Tet_bbinf;
 } Tet_bbase
= (Basis *)NULL;
 } void 
Pyr_reset_basis(void){
for(Pyr_bbinf = Pyr_bbase;
 Pyr_bbinf;
 Pyr_bbinf = Pyr_bbinf->next){ free(Pyr_bbinf->vert[4].a);

// small leak here free(Pyr_bbinf->vert);
 free(Pyr_bbinf->edge);
 free(Pyr_bbinf->face);

free(Pyr_bbinf->intr);
 } if(Pyr_bbase)
while(Pyr_bbase){ Pyr_bbinf = Pyr_bbase->next;
 free(Pyr_bbase);
 Pyr_bbase
= Pyr_bbinf;
 } Pyr_bbase = (Basis *)NULL;

} void 
Prism_reset_basis(void){
for(Prism_bbinf = Prism_bbase;
 Prism_bbinf;
 Prism_bbinf = Prism_bbinf->next){
free(Prism_bbinf->vert[4].a);
 // small leak here free(Prism_bbinf->vert);

free(Prism_bbinf->edge);
 free(Prism_bbinf->face);
 free(Prism_bbinf->intr);

} if(Prism_bbase) while(Prism_bbase){ Prism_bbinf
= Prism_bbase->next;
 free(Prism_bbase);
 Prism_bbase = Prism_bbinf;
 } Prism_bbase
= (Basis *)NULL;
 } void 
Hex_reset_basis(void){
for(Hex_bbinf = Hex_bbase;
 Hex_bbinf;
 Hex_bbinf = Hex_bbinf->next){ free(Hex_bbinf->vert[0].a);

free(Hex_bbinf->vert);
 free(Hex_bbinf->edge);
 free(Hex_bbinf->face);
 }
if(Hex_bbase) while(Hex_bbase){ Hex_bbinf
= Hex_bbase->next;
 free(Hex_bbase);
 Hex_bbase = Hex_bbinf;
 } Hex_bbase
= (Basis *)NULL;
 } /* set up pointer
structure for basis */ Basis
*Tri_mem_base(int L, int qa, int qb, int qc)
{ register int i,j;
 Basis *b = (Basis
*)calloc(1,sizeof(Basis));
 Mode
*m;
 b->id = L;
 b->qa = qa;
 b->qb = qb;
 b->qc = qc;
 /* set up basis
so that all mode structure * * are consequative
from Basis->v[0] */ b->vert = m = mvector(0,NTri_verts
+ NTri_edges*(L-2) + (L-3)*(L-2)/2);

m += NTri_verts;
 /* set up edges
*/ b->edge = (Mode **)calloc(NTri_edges,sizeof(Mode
*));
 if(L>2) for(i = 0 ;
 i <
 NTri_edges;

++i,m+=(L-2)) b->edge[i] = m;
 /* set up faces */ b->face = (Mode
***)calloc(NTri_faces,sizeof(Mode
**));
 if(L>3){ for(i = 0;
 i <
 NTri_faces;

++i){ b->face[i] = (Mode **)calloc((L-3),sizeof(Mode
*));
 for(j = 0;
 j <
 L-3;
 m+= (L-3-j),++j) b->face[i][j] = m;
 } } return
b;
 } // vertex == 1 mode // edge == L-2
modes // face == (L-2)*(L-2) modes /* set up pointer structure for basis
*/ Basis *Quad_mem_base(int
L, int qa, int qb, int qc) { register int i,j;
 Basis
*b = (Basis *)malloc(sizeof(Basis));

Mode *m;
 b->id = L;
 b->qa = qa;
 b->qb
= qb;
 b->qc = qc;
 /* set up basis so
that all mode structure * * are consequative
from Basis->v[0] */ b->vert = m = mvector(0,NQuad_verts
+ NQuad_edges*(L-2) + (L-2)*(L-2));

m += NQuad_verts;
 /* set up edges
*/ b->edge = (Mode **)malloc(NQuad_edges*sizeof(Mode
*));
 if(L>2) for(i = 0 ;
 i <
 NQuad_edges;

++i,m+=(L-2)) b->edge[i] = m;
 /* set up faces */ b->face = (Mode
***)malloc(NQuad_faces*sizeof(Mode
**));
 if(L>2){ for(i = 0;
 i <
 NQuad_faces;

++i){ b->face[i] = (Mode **)malloc((L-2)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-2;
 m+= (L-2),++j) b->face[i][j] = m;
 } } return
b;
 } /* set up pointer structure for basis
*/ Basis *Tet_mem_base(int
L, int qa, int qb, int qc) { register int i,j;
 Basis
*b = (Basis *)malloc(sizeof(Basis));

Mode *m;
 b->id = L;
 b->qa = qa;
 b->qb
= qb;
 b->qc = qc;
 /* set up basis so
that all mode structure * * are consequative
from Basis->v[0] */ b->vert = m = mvector(0,NTet_verts
+ NTet_edges*(L-2) + NTet_faces*(L-3)*(L-2)/2
+ (L-4)*(L-3)*(L-2)/6-1);
 m += NTet_verts;

/* set up edges */ b->edge = (Mode **)malloc(NTet_edges*sizeof(Mode
*));
 if(L>2) for(i = 0 ;
 i <
 NTet_edges;

++i,m+=(L-2)) b->edge[i] = m;
 /* set up faces */ b->face = (Mode
***)malloc(NTet_faces*sizeof(Mode
**));
 if(L>3){ for(i = 0;
 i <
 NTet_faces;

++i){ b->face[i] = (Mode **)malloc((L-3)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-3;
 m+= (L-3-j),++j) b->face[i][j] = m;
 } } /*
set up interior (3D) */ if(L>4){ b->intr
= (Mode ***)malloc((L-4)*sizeof(Mode
**));
 for(i = 0;
 i <
 L-4;
 ++i){ b->intr[i] = (Mode
**)malloc((L-4)*sizeof(Mode *));
 for(j
= 0;
 j <
 L-4-i;
 m+= (L-4-i-j),++j) b->intr[i][j] = m;
 } } return b;

} /* set up pointer structure for basis
*/ Basis *Pyr_mem_base(int
L, int qa, int qb, int qc) { register int i,j;
 Basis
*b = (Basis *)malloc(sizeof(Basis));

Mode *m;
 b->id = L;
 b->qa = qa;
 b->qb
= qb;
 b->qc = qc;
 /* set up basis so
that all mode structure * * are consequative
from Basis->v[0] */ b->vert = m = mvector(0,NPyr_verts
+ NPyr_edges*(L-2) + (L-2)*(L-2)
+ 4*(L-3)*(L-2)/2 + (L-3)*(L-2)*(L-1)/6);
 // interior m += NPyr_verts;

/* set up edges */ b->edge = (Mode **)malloc(NPyr_edges*sizeof(Mode
*));
 if(L>2) for(i = 0 ;
 i <
 NPyr_edges;

++i,m+=(L-2)) b->edge[i] = m;
 /* set up faces */ b->face = (Mode
***)malloc(NPyr_faces*sizeof(Mode
**));
 if(L>2){ b->face[0] = (Mode
**)malloc((L-2)*sizeof(Mode *));
 for(j
= 0;
 j <
 L-2;
 m+= (L-2),++j) b->face[0][j] = m;
 for(i=1;
i<
NPyr_faces;
++i){
b->face[i] = (Mode **)malloc((L-3)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-3;
 m+= (L-3-j),++j) b->face[i][j] = m;
 } } /*
set up interior (3D) */ // here if(L>3){
b->intr = (Mode ***)malloc((L-3)*sizeof(Mode
**));
 for(i = 0;
 i <
 L-3;
 ++i){ b->intr[i] = (Mode
**)malloc((L-3)*sizeof(Mode *));
 for(j
= 0;
 j <
 L-3-i;
 m += L-3-i-j, ++j){ b->intr[i][j] = m;
 } } } return
b;
 } /* set up pointer structure for basis
*/ Basis *Prism_mem_base(int
L, int qa, int qb, int qc) { register int i,j;
 Basis
*b = (Basis *)calloc(1,sizeof(Basis));

Mode *m;
 b->id = L;
 b->qa = qa;
 b->qb
= qb;
 b->qc = qc;
 /* set up basis so
that all mode structure * * are consequative
from Basis->v[0] */ b->vert = m = mvector(0,NPrism_verts
+ NPrism_edges*(L-2) + 3*(L-2)*(L-2)
+ 2*(L-3)*(L-2)/2 + (L-2)*(L-3)*(L-2)/2-1);
 m += NPrism_verts;

/* set up edges */ b->edge = (Mode **)calloc(NPrism_edges,sizeof(Mode
*));
 if(L>2) for(i = 0 ;
 i <
 NPrism_edges;

++i,m+=(L-2)) b->edge[i] = m;
 /* set up faces */ b->face = (Mode
***)calloc(NPrism_faces,sizeof(Mode
**));
 if(L>2){ b->face[0] = (Mode
**)malloc((L-2)*sizeof(Mode *));
 for(j
= 0;
 j <
 L-2;
 m+= (L-2),++j) b->face[0][j] = m;
 if(L>3){
b->face[1] = (Mode **)malloc((L-3)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-3;
 m+= (L-3-j),++j) b->face[1][j] = m;
 } b->face[2]
= (Mode **)malloc((L-2)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-2;
 m+= (L-2),++j) b->face[2][j] = m;
 if(L>3){
b->face[3] = (Mode **)malloc((L-3)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-3;
 m+= (L-3-j),++j) b->face[3][j] = m;
 } b->face[4]
= (Mode **)malloc((L-2)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-2;
 m+= (L-2),++j) b->face[4][j] = m;
 } /* set
up interior (3D) */ if(L>3){ b->intr =
(Mode ***)malloc((L-3)*sizeof(Mode
**));
 for(i = 0;
 i <
 L-3;
 ++i){ b->intr[i] = (Mode
**)malloc((L-2)*sizeof(Mode *));
 for(j
= 0;
 j <
 L-2;
 m+= L-3-i,++j) b->intr[i][j] = m;
 } } return b;
 } /* set
up pointer structure for basis */ Basis
*Hex_mem_base(int L, int qa, int qb, int qc)
{ register int i,j;
 Basis *b = (Basis
*)malloc(sizeof(Basis));
 Mode
*m;
 b->id = L;
 b->qa = qa;
 b->qb = qb;
 b->qc = qc;
 /* set up basis
so that all mode structure * * are consequative
from Basis->v[0] */ b->vert = m = mvector(0,NHex_verts
+ NHex_edges*(L-2) + NHex_faces*(L-2)*(L-2)
+ (L-2)*(L-2)*(L-2)-1);
 m += NHex_verts;

/* set up edges */ b->edge = (Mode **)malloc(NHex_edges*sizeof(Mode
*));
 if(L>2) for(i = 0 ;
 i <
 NHex_edges;

++i,m+=(L-2)) b->edge[i] = m;
 /* set up faces */ b->face = (Mode
***)malloc(NHex_faces*sizeof(Mode
**));
 if(L>2){ for(i = 0;
 i <
 NHex_faces;

++i){ b->face[i] = (Mode **)malloc((L-2)*sizeof(Mode
*));
 for(j = 0;
 j <
 L-2;
 m+= (L-2),++j) b->face[i][j] = m;
 } } /* set
up interior (3D) */ if(L>2){ b->intr =
(Mode ***)malloc((L-2)*sizeof(Mode
**));
 for(i = 0;
 i <
 L-2;
 ++i){ b->intr[i] = (Mode
**)malloc((L-2)*sizeof(Mode *));
 for(j
= 0;
 j <
 L-2;
 m+= L-2,++j) b->intr[i][j] = m;
 } } return b;
 } /*----------------------------------------------------------------------------*
* Set up memory for modal storage.
Not all modes need to have independent * * storage since the modes are
repeated. The independent terms in 'a' 'b' * * * * 'a' 'b' * * 1.0 1.0
(3d only ) * * (1+a)/2 (1+b)/2 * * (1-a)/2 (1-b)/2 * * (1-a)/2 (1+a)/2
(1-b)/2 (1+b)/2 P_{m-1}^{1,1}(b) * * (1-a)/2 (1+a)/2 P_1^{1,1}(a) {(1-b)/2}^l
* * (1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a) {(1-b)/2}^2 (1+b)/2 P_{m-1}^{3,1}(b)
* * {(1-b)/2}^l (1+b)/2 P_{m-1}^{2l-1,1}(b) * * * * Each polynomial in
a,b,c are stored contiguously in this order and then * * the appropriate
pointers are assigned to the correct position in this * * vector.
This allows matrix operations in sumfactorisation routines * //OVERLOAD
CALL: vector: Coords.c(?), Curvi.c(?)
*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*
* this routine sets up the storage in continuous vectors as described above
* * it also sets up the pointer system so that the modes can be accessed
* * individually * *---------------------------------------------------------------------------*/
void 
Tri_mem_modes(Basis
*b){ register int i,j;
 const int L = b->id, qa = b->qa, qb = b->qb;

double *s;
 Mode *v = b->vert, **e = b->edge,
***f = b->face;
 /* allocate memory
for 'a' polynomials */ s = dvector(0,qa*(L+1)-1);
 v[2].a = s;
 s += qa;

/* set up pointers */ v[1].a = s;
 s += qa;
 v[0].a = s;
 s += qa;
 for(i =
0;
 i <
 L-2;
 ++i, s+=qa){ e[0][i].a = s;
 e[1][i].a = v[1].a;
 e[2][i].a
= v[0].a;
 } for(i = 0;
 i <
 L-3;
 ++i) for(j = 0;
 j <
 L-3-i;
 ++j){
f[0][i][j].a = e[0][i].a;
 } /* allocation for b */ s = dvector(0,qb*(L+(L-2)*(L-1)/2)-1);

v[2].b = s;
 s += qb;
 v[0].b = s;
 v[1].b = s;
 s += qb;
 for(i = 0;
 i <

L-2;
 ++i, s+=qb){ e[1][i].b = s;
 e[2][i].b = e[1][i].b;
 } for(i = 0;
 i
<
 L-2;
 ++i){ e[0][i].b = s;
 s+=qb;
 for(j = 0;
 j <
 L-3-i;
 ++j,s+=qb){
f[0][i][j].b = s;
 } } } /*--------------------------------------------------------------------------*
* Set up memory for modal storage.
Not all modes need to have independent * * storage since the modes are
repeated. The independent terms in 'a' 'b' * * 'a' 'b' * * (1+a)/2 (1+b)/2
* * (1-a)/2 (1-b)/2 * * (1-a)/2 (1+a)/2 (1-b)/2 (1+b)/2 * * (1-a)/2 (1+a)/2
P_1^{1,1}(a) (1-b)/2 (1+b)/2 P_1^{1,1}(b) * * (1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a)
(1-b)/2 (1+b)/2 P_{l-2}^{1,1}(b) * * * * * * Each polynomial in a,b,c are
stored contiguously in this order and then * * the appropriate pointers
are assigned to the correct position in this * * vector. This allows
matrix operations in sumfactorisation routines * //OVERLOAD CALL: vector:
Coords.c(?), Curvi.c(?)
*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*
* this routine sets up the storage in continuous vectors as described above*
* it also sets up the pointer system so that the modes can be accessed
* * individually * *--------------------------------------------------------------------------*/
// DONE -- checked void 
Quad_mem_modes(Basis
*b){ register int i,j;
 const int L = b->id, qa = b->qa, qb = b->qb;

double *s;
 Mode *v = b->vert, **e = b->edge,
***f = b->face;
 // count up number of different modes /* allocate memory
for 'a' polynomials */ s = dvector(0,qa*L-1);
 /* set up pointers */ v[1].a
= v[2].a = s;
 s += qa;
 v[0].a = v[3].a = s;
 s += qa;
 for(i = 0;
 i <

L-2;
 ++i, s+=qa){ e[0][i].a = s;
 e[1][i].a = v[1].a;
 e[2][i].a = s;
 e[3][i].a
= v[0].a;
 } for(i = 0;
 i <
 L-2;
 ++i) for(j = 0;
 j <
 L-2;
 ++j) f[0][i][j].a
= e[0][i].a;
 /* allocation for b */ s = dvector(0,qb*L-1);
 v[3].b = v[2].b
= s;
 s += qb;
 v[1].b = v[0].b = s;
 s += qb;
 for(i = 0;
 i <
 L-2;
 ++i,
s+=qb){ e[0][i].b = v[1].b;
 e[1][i].b = s;
 e[2][i].b = v[2].b;
 e[3][i].b
= s;
 } for(i = 0;
 i <
 L-2;
 ++i) for(j = 0;
 j <
 L-2;
 ++j) f[0][i][j].b
= e[1][j].b;
 } /*----------------------------------------------------------------------------*
* Set up memory for modal storage.
Not all modes need to have independent * * storage since the modes are
repeated. The independent terms in 'a' 'b' * * and 'c' are: * * * * 'a'
'b' * * 1.0 1.0 (3d only ) * * (1+a)/2 (1+b)/2 * * (1-a)/2 (1-b)/2 * *
(1-a)/2 (1+a)/2 (1-b)/2 (1+b)/2 P_{m-1}^{1,1}(b) * * (1-a)/2 (1+a)/2 P_1^{1,1}(a)
{(1-b)/2}^l * * (1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a) {(1-b)/2}^2 (1+b)/2 P_{m-1}^{3,1}(b)
* * {(1-b)/2}^l (1+b)/2 P_{m-1}^{2l-1,1}(b) * * * * 'c' * * (1+c)/2 * *
(1-c)/2 * * {(1-c)/2}^2 * * {(1-c)/2}^{l+m} * * (1-c)/2 (1+c)/2 P_{n-1}^{1,1}(c)
* * {(1-c)/2}^2 (1+c)/2 P_{n-1}^{3,1}(c) * * {(1-c)/2}^{l+m} (1+c)/2 P_{n-1}^{2L+2m-1,1}(c)
* * * * Each polynomial in a,b,c are stored contiguously in this order
and then * * the appropriate pointers are assigned to the correct position
in this * * vector. This allows matrix operations in sumfactorisation
routines * //OVERLOAD CALL: vector: Coords.c(?),
Curvi.c(?) *---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*
* this routine sets up the storage in continuous vectors as described above
* * it also sets up the pointer system so that the modes can be accessed
* * individually * *---------------------------------------------------------------------------*/
void 
Tet_mem_modes(Basis
*b){ register int i,j;
 const int L = b->id, qa = b->qa, qb = b->qb;

double *s;
 Mode *v = b->vert, **e = b->edge,
***f = b->face;
 register int k;
 int qc = b->qc;
 Mode
***in = b->intr;
 /* allocate memory
for 'a' polynomials */ s = dvector(0,qa*(L+1)-1);
 v[2].a = s;
 s += qa;

/* set up pointers */ v[1].a = s;
 s += qa;
 v[0].a = s;
 s += qa;
 v[3].a
= v[2].a;
 for(i = 0;
 i <
 L-2;
 ++i, s+=qa){ e[0][i].a = s;
 e[1][i].a
= v[1].a;
 e[2][i].a = v[0].a;
 e[3][i].a = v[0].a;
 e[4][i].a = v[1].a;
 e[5][i].a
= v[2].a;
 } for(i = 0;
 i <
 L-3;
 ++i) for(j = 0;
 j <
 L-3-i;
 ++j){
f[0][i][j].a = e[0][i].a;
 f[1][i][j].a = e[0][i].a;
 f[2][i][j].a = e[1][i].a;

f[3][i][j].a = e[2][i].a;
 } for(i = 0;
 i <
 L-4;
 ++i) for(j = 0;
 j <

L-4-i;
 ++j) for(k = 0;
 k <
 L-4-i-j;
 ++k) in[i][j][k].a = f[0][i][j].a;

/* allocation for b */ s = dvector(0,qb*(L+(L-2)*(L-1)/2+1)-1);
 /* set
up modal pointers */ v[3].b = s;
 s += qb;
 v[2].b = s;
 s += qb;
 v[0].b =
v[1].b = s;
 s+=qb;
 for(i = 0;
 i <
 L-2;
 ++i, s+=qb){ e[1][i].b = s;
 e[2][i].b
= e[1][i].b;
 e[3][i].b = v[0].b;
 e[4][i].b = v[1].b;
 e[5][i].b = v[2].b;

} for(i = 0;
 i <
 L-2;
 ++i){ e[0][i].b = s;
 s+=qb;
 for(j = 0;
 j <

L-3-i;
 ++j,s+=qb){ f[0][i][j].b = s;
 f[1][i][j].b = e[0][i].b;
 f[2][i][j].b
= e[1][i].b;
 f[3][i][j].b = e[2][i].b;
 } } for(i = 0;
 i <
 L-4;
 ++i)
for(j = 0;
 j <
 L-4-i;
 ++j) for(k = 0;
 k <
 L-4-i-j;
 ++k) in[i][j][k].b
= f[0][i][j].b;
 /* allocation for c */ s = dvector(0,qc*(L+(L-2)*(L-1)/2)-1);

v[3].c = s;
 s += qc;
 v[0].c = v[1].c = v[2].c = s;
 s += qc;
 for(i = 0;

i <
 L-2;
 ++i, s+=qc){ e[0][i].c = s;
 e[1][i].c = e[0][i].c;
 e[2][i].c
= e[0][i].c;
 } for(i = 0;
 i <
 L-2;
 ++i, s+=qc){ e[3][i].c = s;
 e[4][i].c
= e[3][i].c;
 e[5][i].c = e[3][i].c;
 } for(i = 0;
 i <
 L-3;
 ++i) for(j
= 0;
 j <
 L-3-i;
 ++j,s+=qc){ f[1][i][j].c = s;
 f[2][i][j].c = f[1][i][j].c;

f[3][i][j].c = f[2][i][j].c;
 } for(i = 0;
 i <
 L-3;
 ++i) for(j = 0;
 j
<
 (L-3-i);
 ++j) f[0][i][j].c = e[0][j+i+1].c;
 for(i = 0;
 i <
 L-4;

++i) for(j = 0;
 j <
 L-4-i;
 ++j) for(k = 0;
 k <
 L-4-i-j;
 ++k) in[i][j][k].c
= f[1][j+i+1][k].c;
 } /*----------------------------------------------------------------------------*
* Set up memory for modal storage.
Not all modes need to have independent * * storage since the modes are
repeated. The independent terms in 'a' 'b' * * and 'c' are: * * * * 'a'
'b' * * (1-a)/2 * (1+a)/2 (1-b)/2 * * (1-a)/2 (1+a)/2 (1-b)/2 (1+b)/2 P_{m-1}^{1,1}(b)
* * (1-a)/2 (1+a)/2 P_1^{1,1}(a) {(1-b)/2}^l * * (1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a)
{(1-b)/2}^2 (1+b)/2 P_{m-1}^{3,1}(b) * * {(1-b)/2}^l (1+b)/2 P_{m-1}^{2l-1,1}(b)
* * * * 'c' * * (1+c)/2 * * (1-c)/2 * * {(1-c)/2}^2 * * {(1-c)/2}^{l+m}
* * (1-c)/2 (1+c)/2 P_{n-1}^{1,1}(c) * * {(1-c)/2}^2 (1+c)/2 P_{n-1}^{3,1}(c)
* * {(1-c)/2}^{l+m} (1+c)/2 P_{n-1}^{2L+2m-1,1}(c) * * * * Each polynomial
in a,b,c are stored contiguously in this order and then * * the appropriate
pointers are assigned to the correct position in this * * vector.
This allows matrix operations in sumfactorisation routines * //OVERLOAD
CALL: vector: Coords.c(?), Curvi.c(?)
*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*
* this routine sets up the storage in continuous vectors as described above
* * it also sets up the pointer system so that the modes can be accessed
* * individually * *---------------------------------------------------------------------------*/
void 
Pyr_mem_modes(Basis
*b){ register int i,j,k;
 const int L = b->id, qa = b->qa, qc = b->qc;

double *sa, *onea, *v1a, *v2a, *eia;
 double *sc, *onec, *v1c, *v2c, *eic;

Mode *v = b->vert, **e = b->edge, ***f
= b->face;
 Mode ***intr = b->intr;
 /*
allocate memory for 'a' polynomials
*/ sa = dvector(0,qa*(L+1)-1);
 sc = dvector(0,qc*(1+L+L-2+ (L-2)*(L-2)+(L-3)*(L-2)/2+
(L-3)*(L-2)*(L-1)/6)-1);
 // interior onea = sa;
 sa += qa;
 onec = sc;
 sc+=qc;

v1a = sa;
 sa += qa;
 v1c = sc;
 sc+=qc;
 eia = sa;
 sa += (L-2)*qa;
 v2a = sa;

sa += qa;
 v2c = sc;
 sc+=qc;
 // done v[0].a = v1a;
 v[0].b = v1a;
 v[0].c
= v1c;
 v[1].a = v2a;
 v[1].b = v1a;
 v[1].c = v1c;
 v[2].a = v2a;
 v[2].b =
v2a;
 v[2].c = v1c;
 v[3].a = v1a;
 v[3].b = v2a;
 v[3].c = v1c;
 v[4].a = onea;

v[4].b = onea;
 v[4].c = v2c;
 sa = eia;
 for(i = 0;
 i <
 L-2;
 ++i, sa+=qa,
sc += qc){ eia = sa;
 eic = sc;
 e[0][i].a = eia;
 e[0][i].b = v1a;
 e[1][i].a
= v2a;
 e[1][i].b = eia;
 e[2][i].a = eia;
 e[2][i].b = v2a;
 e[3][i].a = v1a;

e[3][i].b = eia;
 e[4][i].a = v1a;
 e[4][i].b = v1a;
 e[4][i].c = eic;
 e[5][i].a
= v2a;
 e[5][i].b = v1a;
 e[5][i].c = eic;
 e[6][i].a = v2a;
 e[6][i].b = v2a;

e[6][i].c = eic;
 e[7][i].a = v1a;
 e[7][i].b = v2a;
 e[7][i].c = eic;
 } for(i
= 0;
 i <
 L-2;
 ++i, sc += qc) e[3][i].c = e[2][i].c = e[1][i].c = e[0][i].c
= sc;
 for(i = 0;
 i <
 L-2;
 ++i) for(j = 0;
 j <
 L-2;
 ++j,sc+=qc){ f[0][i][j].a
= e[0][j].a;
 f[0][i][j].b = e[0][i].a;
 f[0][i][j].c = sc;
 } for(i = 0;

i <
 L-3;
 ++i) for(j = 0;
 j <
 L-3-i;
 ++j, sc += qc){ f[1][i][j].a
= e[0][i].a;
 f[1][i][j].b = v1a;
 f[1][i][j].c = sc;
 f[3][i][j].a = e[0][i].a;

f[3][i][j].b = v2a;
 f[3][i][j].c = sc;
 f[2][i][j].a = v2a;
 f[2][i][j].b
= e[0][i].a;
 f[2][i][j].c = sc;
 f[4][i][j].a = v1a;
 f[4][i][j].b = e[0][i].a;

f[4][i][j].c = sc;
 } // interior // here for(i = 0;
 i <
 L-3;
 ++i) for(j
= 0;
 j <
 L-3-i;
 ++j) for(k = 0;
 k <
 L-3-i-j;
 ++k,sc+=qc){ intr[i][j][k].a
= e[0][i].a;
 intr[i][j][k].b = e[0][j].a;
 intr[i][j][k].c = sc;
 } } /*----------------------------------------------------------------------------*
* Set up memory for modal storage.
Not all modes need to have independent * * storage since the modes are
repeated. The independent terms in 'a' 'b' * * and 'c' are: * * * * 'a'
'b' * * (1-a)/2 * (1+a)/2 (1-b)/2 * * (1-a)/2 (1+a)/2 (1-b)/2 (1+b)/2 P_{m-1}^{1,1}(b)
* * (1-a)/2 (1+a)/2 P_1^{1,1}(a) {(1-b)/2}^l * * (1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a)
{(1-b)/2}^2 (1+b)/2 P_{m-1}^{3,1}(b) * * {(1-b)/2}^l (1+b)/2 P_{m-1}^{2l-1,1}(b)
* * * * 'c' * * (1+c)/2 * * (1-c)/2 * * {(1-c)/2}^2 * * {(1-c)/2}^{l+m}
* * (1-c)/2 (1+c)/2 P_{n-1}^{1,1}(c) * * {(1-c)/2}^2 (1+c)/2 P_{n-1}^{3,1}(c)
* * {(1-c)/2}^{l+m} (1+c)/2 P_{n-1}^{2L+2m-1,1}(c) * * * * Each polynomial
in a,b,c are stored contiguously in this order and then * * the appropriate
pointers are assigned to the correct position in this * * vector.
This allows matrix operations in sumfactorisation routines * //OVERLOAD
CALL: vector: Coords.c(?), Curvi.c(?)
*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*
* this routine sets up the storage in continuous vectors as described above
* * it also sets up the pointer system so that the modes can be accessed
* * individually * *---------------------------------------------------------------------------*/
void 
Prism_mem_modes(Basis
*b){ register int i,j,k;
 const int L = b->id, qa = b->qa, qc = b->qc;

double *sa, *onea, *v1a, *v2a, *eia;
 double *sc, *onec, *v1c, *v2c, *eic;

Mode *v = b->vert, **e = b->edge, ***f
= b->face;
 Mode ***intr = b->intr;
 /*
allocate memory for 'a' polynomials
*/ sa = dvector(0,qa*(L+1)-1);
 sc = dvector(0,qc*(1+L+L-2+(L-3)*(L-2)/2)-1);

onea = sa;
 sa += qa;
 onec = sc;
 sc+=qc;
 v1a = sa;
 sa += qa;
 v1c = sc;
 sc+=qc;

eia = sa;
 sa += (L-2)*qa;
 v2a = sa;
 sa += qa;
 v2c = sc;
 sc+=qc;
 // done
v[0].a = v1a;
 v[0].b = v1a;
 v[0].c = v1c;
 v[1].a = v2a;
 v[1].b = v1a;
 v[1].c
= v1c;
 v[2].a = v2a;
 v[2].b = v2a;
 v[2].c = v1c;
 v[3].a = v1a;
 v[3].b =
v2a;
 v[3].c = v1c;
 v[4].a = onea;
 v[4].b = v1a;
 v[4].c = v2c;
 v[5].a =
onea;
 v[5].b = v2a;
 v[5].c = v2c;
 sa = eia;
 for(i = 0;
 i <
 L-2;
 ++i,
sa+=qa, sc += qc){ eia = sa;
 eic = sc;
 e[0][i].a = eia;
 e[0][i].b = v1a;

e[1][i].a = v2a;
 e[1][i].b = eia;
 e[1][i].c = v1c;
 e[2][i].a = eia;
 e[2][i].b
= v2a;
 e[3][i].a = v1a;
 e[3][i].b = eia;
 e[3][i].c = v1c;
 e[4][i].a = v1a;

e[4][i].b = v1a;
 e[4][i].c = eic;
 e[5][i].a = v2a;
 e[5][i].b = v1a;
 e[5][i].c
= eic;
 e[6][i].a = v2a;
 e[6][i].b = v2a;
 e[6][i].c = eic;
 e[7][i].a = v1a;

e[7][i].b = v2a;
 e[7][i].c = eic;
 e[8][i].a = onea;
 e[8][i].b = eia;
 e[8][i].c
= v2c;
 } for(i = 0;
 i <
 L-2;
 ++i, sc += qc) e[2][i].c = e[0][i].c =
sc;
 for(i = 0;
 i <
 L-2;
 ++i) for(j = 0;
 j <
 L-2;
 ++j){ eia = e[0][j].a;

eic = e[4][i].c;
 f[0][i][j].a = eia;
 f[0][j][i].b = eia;
 f[0][i][j].c =
e[0][j].c;
 f[2][i][j].a = v2a;
 f[2][i][j].b = eia;
 f[2][i][j].c = eic;

f[4][i][j].b = eia;
 f[4][i][j].a = v1a;
 f[4][i][j].c = eic;
 } for(i = 0;

i <
 L-3;
 ++i) for(j = 0;
 j <
 L-3-i;
 ++j, sc += qc){ f[1][i][j].a
= e[0][i].a;
 f[1][i][j].b = v1a;
 f[1][i][j].c = sc;
 f[3][i][j].a = e[0][i].a;

f[3][i][j].b = v2a;
 f[3][i][j].c = sc;
 } for(i = 0;
 i <
 L-3;
 ++i) for(j
= 0;
 j <
 L-2;
 ++j) for(k = 0;
 k <
 L-3-i;
 ++k){ intr[i][j][k].a =
e[0][i].a;
 intr[i][j][k].b = e[0][j].a;
 intr[i][j][k].c = f[1][i][k].c;

} } /*----------------------------------------------------------------------------*
* Set up memory for modal storage.
Not all modes need to have independent * * storage since the modes are
repeated. The independent terms in 'a' 'b' * * and 'c' are: * * * * 'a'
'b' * * (1-a)/2 * (1+a)/2 (1-b)/2 * * (1-a)/2 (1+a)/2 (1-b)/2 (1+b)/2 P_{m-1}^{1,1}(b)
* * (1-a)/2 (1+a)/2 P_1^{1,1}(a) {(1-b)/2}^l * * (1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a)
{(1-b)/2}^2 (1+b)/2 P_{m-1}^{3,1}(b) * * {(1-b)/2}^l (1+b)/2 P_{m-1}^{2l-1,1}(b)
* * * * 'c' * * (1+c)/2 * * (1-c)/2 * * {(1-c)/2}^2 * * {(1-c)/2}^{l+m}
* * (1-c)/2 (1+c)/2 P_{n-1}^{1,1}(c) * * {(1-c)/2}^2 (1+c)/2 P_{n-1}^{3,1}(c)
* * {(1-c)/2}^{l+m} (1+c)/2 P_{n-1}^{2L+2m-1,1}(c) * * * * Each polynomial
in a,b,c are stored contiguously in this order and then * * the appropriate
pointers are assigned to the correct position in this * * vector.
This allows matrix operations in sumfactorisation routines * //OVERLOAD
CALL: vector: Coords.c(?), Curvi.c(?)
*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*
* this routine sets up the storage in continuous vectors as described above
* * it also sets up the pointer system so that the modes can be accessed
* * individually * *---------------------------------------------------------------------------*/
void 
Hex_mem_modes(Basis
*b){ register int i,j,k;
 const int L = b->id, qa = b->qa;
 double *s,
*v1, *v2, *ei, *eia, *eja;
 Mode *v =
b->vert, **e = b->edge, ***f = b->face;
 Mode
***intr = b->intr;
 /* allocate memory
for 'a' polynomials */ s = dvector(0,qa*L-1);
 v1 = s;
 s+= qa;
 ei = s;
 s+=
qa*(L-2);
 v2 = s;
 s+= qa;
 v[0].a = v1;
 v[0].b = v1;
 v[0].c = v1;
 v[1].a
= v2;
 v[1].b = v1;
 v[1].c = v1;
 v[2].a = v2;
 v[2].b = v2;
 v[2].c = v1;

v[3].a = v1;
 v[3].b = v2;
 v[3].c = v1;
 v[4].a = v1;
 v[4].b = v1;
 v[4].c
= v2;
 v[5].a = v2;
 v[5].b = v1;
 v[5].c = v2;
 v[6].a = v2;
 v[6].b = v2;

v[6].c = v2;
 v[7].a = v1;
 v[7].b = v2;
 v[7].c = v2;
 for(i = 0;
 i <
 L-2;

++i, ei+=qa){ e[0][i].a = ei;
 e[0][i].b = v1;
 e[0][i].c = v1;
 e[1][i].b
= ei;
 e[1][i].a = v2;
 e[1][i].c = v1;
 e[2][i].a = ei;
 e[2][i].b = v2;
 e[2][i].c
= v1;
 e[3][i].b = ei;
 e[3][i].a = v1;
 e[3][i].c = v1;
 e[4][i].c = ei;
 e[4][i].a
= v1;
 e[4][i].b = v1;
 e[5][i].c = ei;
 e[5][i].a = v2;
 e[5][i].b = v1;
 e[6][i].c
= ei;
 e[6][i].a = v2;
 e[6][i].b = v2;
 e[7][i].c = ei;
 e[7][i].a = v1;
 e[7][i].b
= v2;
 e[8][i].a = ei;
 e[8][i].b = v1;
 e[8][i].c = v2;
 e[9][i].b = ei;
 e[9][i].a
= v2;
 e[9][i].c = v2;
 e[10][i].a = ei;
 e[10][i].b = v2;
 e[10][i].c = v2;

e[11][i].b = ei;
 e[11][i].a = v1;
 e[11][i].c = v2;
 } for(i = 0;
 i <

L-2;
 ++i) for(j = 0;
 j <
 L-2;
 ++j){ eia = e[0][j].a;
 eja = e[0][i].a;

f[0][i][j].a = eia;
 f[0][i][j].b = eja;
 f[0][i][j].c = v1;
 f[1][i][j].a
= eia;
 f[1][i][j].b = v1;
 f[1][i][j].c = eja;
 f[2][i][j].a = v2;
 f[2][i][j].b
= eia;
 f[2][i][j].c = eja;
 f[3][i][j].a = eia;
 f[3][i][j].b = v2;
 f[3][i][j].c
= eja;
 f[4][i][j].a = v1;
 f[4][i][j].b = eia;
 f[4][i][j].c = eja;
 f[5][i][j].a
= eia;
 f[5][i][j].b = eja;
 f[5][i][j].c = v2;
 } for(i = 0;
 i <
 L-2;

++i) for(j = 0;
 j <
 L-2;
 ++j) for(k = 0;
 k <
 L-2;
 ++k){ intr[i][j][k].a
= e[0][k].a;
 intr[i][j][k].b = e[0][j].a;
 intr[i][j][k].c = e[0][i].a;

} } void 
Tri_set_vertices(Mode
*v, int q, char dir){ double one = 1.0,*z,*w;
 getzw(q,&
z,&
w,dir);

switch (dir){ case 'a': dvsub(q,&
one,0,z,1,v[0].a,1);
 dscal(q,0.5,v[0].a,1);

dsadd(q,one,z,1,v[1].a,1);
 dscal(q,0.5,v[1].a,1);
 dfill(q,one,v[2].a,1);

break;
 case 'b': dvsub(q,&
one,0,z,1,v[0].b,1);
 dscal(q,0.5,v[0].b,1);

dsadd(q,one,z,1,v[2].b,1);
 dscal(q,0.5,v[2].b,1);
 break;
 } } void 
Quad_set_vertices(Basis
*b, int q, char dir){ double one = 1.0,*z,*w;
 Mode
*v = b->vert;
 getzw(q,&
z,&
w,'a');
 switch (dir){
case 'a': dvsub(q,&
one,0,z,1,v[0].a,1);
 dscal(q,0.5,v[0].a,1);
 dsadd(q,one,z,1,v[1].a,1);

dscal(q,0.5,v[1].a,1);
 break;
 case 'b': dvsub(q,&
one,0,z,1,v[1].b,1);

dscal(q,0.5,v[1].b,1);
 dsadd(q,one,z,1,v[3].b,1);
 dscal(q,0.5,v[3].b,1);

break;
 } } void 
Tet_set_vertices(Mode
*v, int q, char dir){ double one = 1.0,*z,*w;
 getzw(q,&
z,&
w,dir);

switch (dir){ case 'a': dvsub(q,&
one,0,z,1,v[0].a,1);
 dsmul(q,0.5,v[0].a,1,v[0].a,1);

dsadd(q,one,z,1,v[1].a,1);
 dsmul(q,0.5,v[1].a,1,v[1].a,1);
 dfill(q,one,v[2].a,1);

break;
 case 'b': dvsub(q,&
one,0,z,1,v[0].b,1);
 dsmul(q,0.5,v[0].b,1,v[0].b,1);

dsadd(q,one,z,1,v[2].b,1);
 dsmul(q,0.5,v[2].b,1,v[2].b,1);
 dfill(q,one,v[3].b,1);

break;
 case 'c': dvsub(q,&
one,0 ,z,1,v[0].c,1);
 dsmul(q,0.5,v[0].c,1,v[0].c,1);

dsadd(q,one,z,1,v[3].c,1);
 dsmul(q,0.5,v[3].c,1,v[3].c,1);
 break;
 } } void 
Pyr_set_vertices(Mode
*v, int q, char dir){ double one = 1.0,*z,*w;
 switch (dir){ case 'a':
getzw(q,&
z,&
w,'a');
 dvsub(q,&
one,0,z,1,v[0].a,1);

dscal(q,0.5,v[0].a,1);
 // (1-a)/2 dsadd(q,one,z,1,v[1].a,1);
 dscal(q,0.5,v[1].a,1);

// (1+a)/2 dfill(q,one,v[4].a,1);
 // 1 break;
 case 'b': // already set
up in 'a' break;
 case 'c': getzw(q,&
z,&
w,'c');

dvsub(q,&
one,0,z,1,v[0].c,1);
 dscal(q,0.5,v[0].c,1);
 // (1-c)/2 dsadd(q,one,z,1,v[4].c,1);

dscal(q,0.5,v[4].c,1);
 // (1+c)/2 break;
 } } void 
Prism_set_vertices(Mode
*v, int q, char dir){ double one = 1.0,*z,*w;
 switch (dir){ case 'a':
getzw(q,&
z,&
w,'a');
 dvsub(q,&
one,0,z,1,v[0].a,1);

dscal(q,0.5,v[0].a,1);
 // (1-a)/2 dsadd(q,one,z,1,v[1].a,1);
 dscal(q,0.5,v[1].a,1);

// (1+a)/2 dfill(q,one,v[4].a,1);
 // 1 break;
 case 'b': // already set
up in 'a' break;
 case 'c': getzw(q,&
z,&
w,'b');

dvsub(q,&
one,0,z,1,v[0].c,1);
 dscal(q,0.5,v[0].c,1);
 // (1-c)/2 dsadd(q,one,z,1,v[4].c,1);

dscal(q,0.5,v[4].c,1);
 // (1+c)/2 break;
 } } void 
Hex_set_vertices(Mode
*v, int q, char dir){ double one = 1.0,*z,*w;
 getzw(q,&
z,&
w,'a');

switch (dir){ case 'a': dvsub(q,&
one,0,z,1,v[0].a,1);
 dsmul(q,0.5,v[0].a,1,v[0].a,1);

dsadd(q,one,z,1,v[1].a,1);
 dsmul(q,0.5,v[1].a,1,v[1].a,1);
 break;
 case
'b': break;
 case 'c': break;
 } } void 
Tri_set_edges(Basis
*b, int q, char dir){ register int i;
 int L = b->id-2;
 double *z,*w;

Mode *v = b->vert, **e = b->edge;
 getzw(q,&
z,&
w,dir);

switch (dir){ case 'a': dvmul(q,v[0].a,1,v[1].a,1,e[0][0].a,1);
 for(i =
1;
 i <
 L;
 ++i){ jacobf(q,z,e[0][i].a,i,1.0,1.0);

dvmul (q,e[0][0].a,1,e[0][i].a,1,e[0][i].a,1);
 } break;
 case 'b': dvmul(q,v[0].b,1,v[0].b,1,e[0][0].b,1);

dvmul(q,v[1].b,1,v[2].b,1,e[1][0].b,1);
 for(i = 1;
 i <
 L;
 ++i){ dvmul
(q,e[0][i-1].b,1,v[0].b,1,e[0][i].b,1);
 jacobf(q,z,e[1][i].b,i,1.0,1.0);

dvmul (q,e[1][0].b,1,e[1][i].b,1,e[1][i].b,1);
 } break;
 } } void 
Quad_set_edges(Basis
*b, int q, char dir){ register int i;
 int L = b->id-2;
 double *z,*w;

Mode *v = b->vert, **e = b->edge;
 getzw(q,&
z,&
w,'a');

switch (dir){ case 'a': dvmul(q,v[0].a,1,v[1].a,1,e[0][0].a,1);
 for(i =
1;
 i <
 L;
 ++i){ jacobf(q,z,e[0][i].a,i,1.0,1.0);

dvmul (q,e[0][0].a,1,e[0][i].a,1,e[0][i].a,1);
 } break;
 case 'b': dvmul(q,v[0].b,1,v[2].b,1,e[1][0].b,1);

for(i = 1;
 i <
 L;
 ++i){ jacobf(q,z,e[1][i].b,i,1.0,1.0);

dvmul (q,e[1][0].b,1,e[1][i].b,1,e[1][i].b,1);
 } break;
 } } void 
Tet_set_edges(Basis
*b, int q, char dir){ register int i;
 int L = b->id-2;
 double *z,*w;

Mode *v = b->vert, **e = b->edge;
 double
alpha = 1.;
 double beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } getzw(q,&
z,&
w,dir);

switch (dir){ case 'a': dvmul(q,v[0].a,1,v[1].a,1,e[0][0].a,1);
 for(i =
1;
 i <
 L;
 ++i){ jacobf(q,z,e[0][i].a,i,alpha,beta);

dvmul (q,e[0][0].a,1,e[0][i].a,1,e[0][i].a,1);
 } break;
 case 'b': dvmul(q,v[0].b,1,v[0].b,1,e[0][0].b,1);

dvmul(q,v[1].b,1,v[2].b,1,e[1][0].b,1);
 for(i = 1;
 i <
 L;
 ++i){ dvmul
(q,e[0][i-1].b,1,v[0].b,1,e[0][i].b,1);
 jacobf(q,z,e[1][i].b,i,alpha,beta);

dvmul (q,e[1][0].b,1,e[1][i].b,1,e[1][i].b,1);
 } break;
 case 'c': dvmul(q,v[0].c,1,v[0].c,1,e[0][0].c,1);

dvmul(q,v[2].c,1,v[3].c,1,e[3][0].c,1);
 for(i = 1;
 i <
 L;
 ++i){ dvmul
(q,e[0][i-1].c,1,v[0].c,1,e[0][i].c,1);
 jacobf(q,z,e[3][i].c,i,alpha,beta);

dvmul (q,e[3][0].c,1,e[3][i].c,1,e[3][i].c,1);
 } break;
 } } void 
Pyr_set_edges(Basis
*b, int q, char dir){ register int i;
 int L = b->id-2;
 double *z,*w;

Mode *v = b->vert, **e = b->edge;
 double
alpha = 1., beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } switch (dir){ case 'a': getzw(q,&
z,&
w,'a');

dvmul(q,v[0].a,1,v[1].a,1,e[0][0].a,1);
 for(i = 1;
 i <
 L;
 ++i){ jacobf(q,z,e[0][i].a,i,alpha,beta);

dvmul (q,e[0][0].a,1,e[0][i].a,1,e[0][i].a,1);
 } break;
 case 'b': break;

case 'c': getzw(q,&
z,&
w,'c');
 // dfill(q,one,e[1][0].c,
1);
 dvmul(q, v[0].c, 1, v[0].c, 1, e[0][0].c, 1);
 for(i=1;
i<
L;
++i) dvmul(q,
v[0].c, 1, e[0][i-1].c, 1, e[0][i].c, 1);
 dvmul(q,v[0].c,1,v[4].c,1,e[4][0].c,1);

for(i = 1;
 i <
 L;
 ++i){ jacobf(q,z,e[4][i].c,i,alpha,beta);

dvmul (q,e[4][0].c,1,e[4][i].c,1,e[4][i].c,1);
 } break;
 } } void 
Prism_set_edges(Basis
*b, int q, char dir){ register int i;
 int L = b->id-2;
 double *z,*w;

Mode *v = b->vert, **e = b->edge;
 double
alpha = 1.;
 double beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } switch (dir){ case 'a': getzw(q,&
z,&
w,'a');

dvmul(q,v[0].a,1,v[1].a,1,e[0][0].a,1);
 for(i = 1;
 i <
 L;
 ++i){ jacobf(q,z,e[0][i].a,i,alpha,beta);

dvmul (q,e[0][0].a,1,e[0][i].a,1,e[0][i].a,1);
 } break;
 case 'b': break;

case 'c': getzw(q,&
z,&
w,'b');
 // dfill(q,one,e[1][0].c,
1);
 dvmul(q, v[0].c, 1, v[0].c, 1, e[0][0].c, 1);
 for(i=1;
i<
L;
++i) dvmul(q,
v[0].c, 1, e[0][i-1].c, 1, e[0][i].c, 1);
 dvmul(q,v[0].c,1,v[4].c,1,e[4][0].c,1);

for(i = 1;
 i <
 L;
 ++i){ jacobf(q,z,e[4][i].c,i,alpha,beta);

dvmul (q,e[4][0].c,1,e[4][i].c,1,e[4][i].c,1);
 } break;
 } } void 
Hex_set_edges(Basis
*b, int q, char dir){ register int i;
 int L = b->id-2;
 double *z,*w;

Mode *v = b->vert, **e = b->edge;
 double
alpha = 1.;
 double beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } getzw(q,&
z,&
w,'a');

switch (dir){ case 'a': dvmul(q,v[0].a,1,v[1].a,1,e[0][0].a,1);
 for(i =
1;
 i <
 L;
 ++i){ jacobf(q,z,e[0][i].a,i,alpha,beta);

dvmul (q,e[0][0].a,1,e[0][i].a,1,e[0][i].a,1);
 } break;
 case 'b': break;

case 'c': break;
 } } void 
Tri_set_faces(Basis
*b, int q, char dir){ register int i,j;
 int L = b->id;
 double *z,*w;

Mode *v = b->vert, **e = b->edge, ***f
= b->face;
 switch (dir){ case 'a': break;
 case 'b': getzw(q,&
z,&
w,dir);

for(i = 0;
 i <
 L-3;
 ++i) dvmul(q,v[2].b,1,e[0][i].b,1,f[0][i][0].b,1);

for(i = 0;
 i <
 L-3;
 ++i) for(j = 1;
 j <
 L-3-i;
 ++j){ jacobf(q,z,f[0][i][j].b,j,2.0*i+3.0,1.0);

dvmul (q,f[0][i][0].b,1,f[0][i][j].b,1,f[0][i][j].b,1);
 } break;
 } } void 
Quad_set_faces(Basis
*, int , char ){ // already set by memory
pointer to edge modes } void 
Tet_set_faces(Basis
*b, int q, char dir){ register int i,j;
 int L = b->id;
 double *z,*w;

Mode *v = b->vert, **e = b->edge, ***f
= b->face;
 double alpha = 1.;
 double beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } switch (dir){ case 'a': break;
 case 'b': getzw(q,&
z,&
w,dir);

for(i = 0;
 i <
 L-3;
 ++i){ dvmul(q,v[2].b,1,e[0][i].b,1,f[0][i][0].b,1);

} for(i = 0;
 i <
 L-3;
 ++i) for(j = 1;
 j <
 L-3-i;
 ++j){ alpha = 2.0*i+3.;

beta = 1.;
 if(option("TWOTWO")){ alpha
+= 1.;
 beta += 1.;
 } jacobf(q,z,f[0][i][j].b,j,alpha,beta);

dvmul (q,f[0][i][0].b,1,f[0][i][j].b,1,f[0][i][j].b,1);
 } break;
 case 'c':
getzw(q,&
z,&
w,dir);
 for(i = 0;
 i <
 L-3;

++i){ dvmul(q,v[3].c,1,e[0][i].c,1,f[1][i][0].c,1);
 } for(i = 0;
 i <

L-3;
 ++i) for(j = 1;
 j <
 L-3-i;
 ++j){ alpha = 2.0*i+3.;
 beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } jacobf(q,z,f[1][i][j].c,j,alpha,beta);

dvmul (q,f[1][i][0].c,1,f[1][i][j].c,1,f[1][i][j].c,1);
 } break;
 } } void 
Pyr_set_faces(Basis
*b, int q, char dir){ int L = b->id, i,j,k;
 Mode
*v = b->vert, **e = b->edge, ***f = b->face, ***in = b->intr;
 double *w,
*z;
 int cnt = 0;
 double alpha, beta;
 switch (dir){ case 'a': break;
 case
'b': break;
 case 'c': getzw(q,&
z,&
w,'c');

for(i = 0;
 i <
 L-2;
 ++i) for(j = 0;
 j <
 L-2;
 ++j){ k = min(L,
i+j);
 dcopy(q, e[0][k].c, 1, f[0][i][j].c, 1);
 } // dvmul(q, e[0][j].c,
1, e[0][i].c, 1,f[0][i][j].c,1);
 for(i = 0;
 i <
 L-3;
 ++i){ dvmul(q,v[4].c,1,e[0][i].c,1,f[1][i][0].c,1);

} for(i = 0;
 i <
 L-3;
 ++i) for(j = 1;
 j <
 L-3-i;
 ++j){ alpha = 2.*i+3.;

beta = 1.;
 if(option("TWOTWO")){ alpha
+= 1.;
 beta += 1.;
 } jacobf(q,z,f[1][i][j].c,j,alpha,beta);

dvmul (q,f[1][i][0].c,1,f[1][i][j].c,1,f[1][i][j].c,1);
 } // interior for(i
= 0;
 i <
 L-3;
 ++i) for(j = 0;
 j <
 L-3-i;
 ++j) for(k=0;
k <
 L-3-i-j;
++k){
alpha = 2.*(i+j)+3.;
 beta = 1.;
 if(option("TWOTWO")){
alpha += 1.;
 beta += 1.;
 } // jacobf(q,z,in[i][j][k].c,k,2.0*(i+j)+3.0,1.0);

jacobf(q,z,in[i][j][k].c,k,alpha,beta);

dvmul (q,f[0][i][j].c,1,in[i][j][k].c,1,in[i][j][k].c,1);
 dvmul (q, v[4].c,1,in[i][j][k].c,1,in[i][j][k].c,1);

++cnt;
 } break;
 } } void 
Prism_set_faces(Basis
*b, int q, char dir){ int L = b->id, i,j;
 Mode
*v = b->vert, **e = b->edge, ***f = b->face;
 double *w, *z;
 double alpha,
beta;
 switch (dir){ case 'a': break;
 case 'b': break;
 case 'c': getzw(q,&
z,&
w,'b');

for(i = 0;
 i <
 L-3;
 ++i){ dvmul(q,v[4].c,1,e[0][i].c,1,f[1][i][0].c,1);

} for(i = 0;
 i <
 L-3;
 ++i) for(j = 1;
 j <
 L-3-i;
 ++j){ alpha = 2.0*i+3.;

beta = 1.;
 if(option("TWOTWO")){ alpha
+= 1.;
 beta += 1.;
 } jacobf(q,z,f[1][i][j].c,j,alpha,beta);

dvmul (q,f[1][i][0].c,1,f[1][i][j].c,1,f[1][i][j].c,1);
 } break;
 } } void 
Hex_set_faces(Basis
*b, int q, char dir){ switch (dir){ case 'a': break;
 case 'b': break;

case 'c': break;
 } } /* normalised basis
in an L2 fashion. * * don't want to normalise vertices since have a collocation
* * property that would be destroyed by it. Want normalisation * * of edges
and faces to match so that can match expansion * * coefficients for C^0
continuity */ static void 
Tri_normalise(Basis
*b, int qa, int qb, int qc){ register int i;
 double *wa,*wb,*s,integral;

double *tmp;
 const int L = b->id;
 qc=qc;
 // compiler fix getzw(qa,&
wa,&
wa,'a');

getzw(qb,&
wb,&
wb,'b');
 tmp = dvector(0,QGmax-1);

/* sort out edge normalisation using side 1 `a' component */ if(L>2){
s = b->edge[0][0].a;
 for(i = 0;
 i <
 L-2;
 ++i, s+=qa){ dvmul(qa,s,1,wa,1,tmp,1);

integral = 1.0/sqrt(ddot(qa,s,1,tmp,1));
 dscal(qa,integral,s,1);
 dscal(qb,integral,b->edge[1][i].b,1);

} } /* normalise face modes using `b' component */ if(L>3){
s = b->face[0][0][0].b;
 for(i = 0;
 i <
 (L-3)*(L-2)/2;
 ++i,s+=qb){ dvmul(qb,s,1,wb,1,tmp,1);

integral = 1.0/sqrt(ddot(qb,s,1,tmp,1));
 dscal(qb,integral,s,1);
 } } free(tmp);

} static void 
Quad_normalise(Basis
*b, char type ){ register int i;
 int L = b->id;
 int qa = b->qa;
 int
qb = b->qb;
 switch(type){ case 'B': /* babuska and Szabo's normalisation
*/ if(L > 2){ for(i = 0;
 i <
 L-2;
++i){
dscal(qa,sqrt((2*i+3)*2)/(double)(i+1),b->edge[0][i].a,1);
 dscal(qa,sqrt((2*i+3)*2)/(double)(i+1),b->edge[2][i].a,1);

dscal(qb,sqrt((2*i+3)*2)/(double)(i+1),b->edge[1][i].b,1);
 dscal(qb,sqrt((2*i+3)*2)/(double)(i+1),b->edge[3][i].b,1);

} } break;
 default: error_msg("type not known");
 //OVERLOAD CALL:
error_msg: gen_utils.h(?), hotel.h(?)
} } /* normalised basis in an L2 fashion.
* * don't want to normalise vertices since have a collocation * * property
that would be destroyed by it. Want normalisation * * of edges and faces
to match so that can match expansion * * coefficients for C^0 continuity
*/ static void 
Tet_normalise(Basis
*b, int , int , int ){ } /* normalised basis
in an L2 fashion. * * don't want to normalise vertices since have a collocation
* * property that would be destroyed by it. Want normalisation * * of edges
and faces to match so that can match expansion * * coefficients for C^0
continuity */ static void 
Pyr_normalise(Basis
*, int , int , int ){ } /* normalised basis
in an L2 fashion. * * don't want to normalise vertices since have a collocation
* * property that would be destroyed by it. Want normalisation * * of edges
and faces to match so that can match expansion * * coefficients for C^0
continuity */ static void 
Prism_normalise(Basis
*, int , int , int ){ } /* normalised basis
in an L2 fashion. * * don't want to normalise vertices since have a collocation
* * property that would be destroyed by it. Want normalisation * * of edges
and faces to match so that can match expansion * * coefficients for C^0
continuity */ static void 
Hex_normalise(Basis
*, int , int , int ){ } /* Function name: Element::fillElmt
Function Purpose: Argument 1: Mode *v
Purpose: Function Notes: */ void Tri::fillElmt(Mode
*v){ register int i;
 for(i = 0;
 i <
 qb;
 ++i) dcopy(qa,v->a,1,h[i],1);

for(i = 0;
 i <
 qb;
 ++i) dscal(qa,v->b[i],h[i],1);
 } void Quad::fillElmt(Mode
*v){ register int i;
 for(i = 0;
 i <
 qb;
 ++i) dcopy(qa,v->a,1,h[i],1);

for(i = 0;
 i <
 qb;
 ++i) dscal(qa,v->b[i],h[i],1);
 } void Tet::fillElmt(Mode
*v){ register int i, j;
 for(i = 0;
 i <
 qb*qc;
 ++i) dcopy(qa,v->a,1,(*h_3d)[i],1);

for(i = 0;
 i <
 qc;
 ++i) for(j = 0;
 j <
 qb;
 ++j) dsmul(qa,v->b[j],h_3d[i][j],1,h_3d[i][j],1);

for(i = 0;
 i <
 qc;
 ++i) dsmul(qa*qb,v->c[i],*(h_3d[i]),1,*(h_3d[i]),1);

} void Pyr::fillElmt(Mode
*v){ register int i, j;
 for(i = 0;
 i <
 qb*qc;
 ++i) dcopy(qa,v->a,1,(*h_3d)[i],1);

for(i = 0;
 i <
 qc;
 ++i) for(j = 0;
 j <
 qb;
 ++j) dsmul(qa,v->b[j],h_3d[i][j],1,h_3d[i][j],1);

for(i = 0;
 i <
 qc;
 ++i) dsmul(qa*qb,v->c[i],*(h_3d[i]),1,*(h_3d[i]),1);

} void Prism::fillElmt(Mode
*v){ register int i, j;
 for(i = 0;
 i <
 qb*qc;
 ++i) dcopy(qa,v->a,1,(*h_3d)[i],1);

for(i = 0;
 i <
 qc;
 ++i) for(j = 0;
 j <
 qb;
 ++j) dsmul(qa,v->b[j],h_3d[i][j],1,h_3d[i][j],1);

for(i = 0;
 i <
 qc;
 ++i) dsmul(qa*qb,v->c[i],*(h_3d[i]),1,*(h_3d[i]),1);

} void Hex::fillElmt(Mode
*v){ register int i, j;
 for(i = 0;
 i <
 qb*qc;
 ++i) dcopy(qa,v->a,1,(*h_3d)[i],1);

for(i = 0;
 i <
 qc;
 ++i) for(j = 0;
 j <
 qb;
 ++j) dsmul(qa,v->b[j],h_3d[i][j],1,h_3d[i][j],1);

for(i = 0;
 i <
 qc;
 ++i) dsmul(qa*qb,v->c[i],*(h_3d[i]),1,*(h_3d[i]),1);

} void Element::fillElmt(Mode
*){ERR;
} //OVERLOAD CALL: ERR: Element.c(?),
hotel.h(?)




Back to Source File Index



C++ to HTML Conversion by ctoohtml