file: Hlib/src/InnerProd.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 <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <polylib.h>
#include "veclib.h"
#include "hotel.h"
static void Hex_map_to_hj(Element *E, double *d);
/*
Function name: Element::Iprod
Function Purpose:
Argument 1: Element *T
Purpose:
Function Notes:
*/
void Tri::Iprod(Element *T){
Basis *ba = T->getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Basis *bb = ba;
Iprod_d(T, ba,bb); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
void Quad::Iprod(Element *T){
Basis *ba = T->getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Basis *bb = ba;
Iprod_d(T, ba,bb); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
void Tet::Iprod(Element *T){
Basis *ba = T->getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Basis *bb = ba;
Basis *bc = ba;
Tet::Iprod_3d(T, ba,bb,bc); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
void Pyr::Iprod(Element *T){
Basis *ba = T->getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Basis *bb = ba;
Basis *bc = ba;
Pyr::Iprod_3d(T, ba,bb,bc); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
void Prism::Iprod(Element *T){
Basis *ba = T->getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Basis *bb = ba;
Basis *bc = ba;
Prism::Iprod_3d(T, ba,bb,bc); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
void Hex::Iprod(Element *T){
Basis *ba = T->getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Basis *bb = ba;
Basis *bc = ba;
Hex::Iprod_3d(T, ba,bb,bc); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
}
void Element::Iprod(Element *){ERR;} // Inner product to Elmt //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
void Tri::Iprod_d(Element *T, Basis *ba, Basis *bb){
register int i,L=0;
double *wa,*wb,*H, *wk_a;
Edge *e = T->edge;
Face *f = T->face;
size_t d = sizeof(double);
T->state = 't';
getzw(T->qa,&wa,&wa,'a');
getzw(T->qb,&wb,&wb,'b');
/* multiply physical solution by jacobean and weights*/
if(T->curvX)
dvmul(qa*qb,T->geom->jac.p,1,*h,1,*T->h,1);
else
dsmul(qa*qb,T->geom->jac.d,*h,1,*T->h,1);
for(i = 0; i < qb; ++i) dvmul(qa,wa,1, T->h[i],1, T->h[i],1);
for(i = 0; i < qb; ++i) dscal(qa,wb[i],T->h[i],1);
dgemm('T','N',qb,lmax,qa,1.0,T->h[0],qa,ba->vert[1].a,qa,0.0,Tri_wk,qb);
#if 0
H = *T->h;
dgemv('T',qb,lmax,1.0,bb->vert[2].b,qb,Tri_wk ,1,0.0,H ,1);
dgemv('T',qb,lmax,1.0,bb->vert[2].b,qb,Tri_wk+qb,1,0.0,H+lmax,1);
H += 2*lmax;
Tri_wk_a = Tri_wk+2*qb;
for(i = 0; i < lmax-2; Tri_wk_a+=qb,H+=lmax-2-i, ++i)
dgemv('T',qb,lmax-2-i,1.0,bb->edge[0][i].b,qb,1,Tri_wk_a,1,0.0,H,1);
#else
H = *T->h;
mxva(bb->vert[2].b,qb,1,Tri_wk ,1,H ,1,lmax,qb);
mxva(bb->vert[2].b,qb,1,Tri_wk+qb,1,H+lmax,1,lmax,qb);
H += 2*lmax;
wk_a = Tri_wk+2*qb;
for(i = 0; i < lmax-2; wk_a+=qb,H+=lmax-2-i, ++i)
mxva(bb->edge[0][i].b,qb,1,wk_a,1,H,1,lmax-2-i,qb);
#endif
/* unpack h */
/* copy field into h for moment */
H = *T->h;
T->vert[2].hj[0] = *H; H++;
T->vert[1].hj[0] = *H; H++;
memcpy(e[1].hj,H,e[1].l*d); H += lmax-2;
T->vert[2].hj[0] += *H; H++;
T->vert[0].hj[0] = *H; H++ ;
memcpy(e[2].hj,H,e[2].l*d); H += lmax-2;
for(i = 0, L = 0; i < e[0].l; L += lmax-2-i, ++i)
e[0].hj[i] = H[L];
++H;
for(i = 0; i < f[0].l; H += lmax-2-i, ++i)
memcpy(f[0].hj[i],H,(f[0].l-i)*d);
}
void Quad::Iprod_d(Element *T, Basis *ba, Basis *bb){
register int i;
double *wa,*wb,*H, *wk;
Edge *e = T->edge;
Face *f = T->face;
size_t d = sizeof(double);
T->state = 't';
getzw(T->qa,&wa,&wa,'a');
getzw(T->qb,&wb,&wb,'a');
/* multiply physical solution by jacobean and weights */
H = Quad_Iprod_wk;
wk = Quad_wk;
dvmul(qa*qb,T->geom->jac.p,1,h[0],1,H,1);
for(i = 0; i < qb; ++i) dvmul(qa,wa,1, H+i*qa,1, H+i*qa,1);
for(i = 0; i < qa; ++i) dvmul(qb,wb,1, H+i, qa, H+i, qa);
dgemm('T','N',lmax, qb,qa,1.0, ba->vert[1].a, qa, H, qa, 0.0, wk, lmax);
dgemm('N','N',lmax,lmax,qb,1.0, wk, lmax, bb->vert[2].b, qb, 0.0, H, lmax);
/* unpack h */
/* copy field into T->h for moment */
T->vert[2].hj[0] = *H; ++H;
T->vert[3].hj[0] = *H; ++H;
memcpy(e[2].hj,H,e[2].l*d); H += lmax-2;
T->vert[1].hj[0] = *H; ++H;
T->vert[0].hj[0] = *H; ++H;
memcpy(e[0].hj,H,e[0].l*d); H += lmax-2;
dcopy(e[1].l, H, lmax, e[1].hj, 1);++H;
dcopy(e[3].l, H, lmax, e[3].hj, 1);++H;
for(i = 0; i < f[0].l; ++H, ++i)
dcopy(f[0].l, H, lmax, f[0].hj[i], 1);
}
void Tet::Iprod_d(Element *, Basis *, Basis *){
fprintf(stderr,"Tet::Iprod_d not implemented\n");
}
void Pyr::Iprod_d(Element *, Basis *, Basis *){
fprintf(stderr,"Pyr::Iprod_d not implemented\n");
}
void Prism::Iprod_d(Element *, Basis *, Basis *){
fprintf(stderr,"Prism::Iprod_d not implemented\n");
}
void Hex::Iprod_d(Element *, Basis *, Basis *){
fprintf(stderr,"Hex::Iprod_d not implemented\n");
}
void Tet::Iprod_3d(Element *E, Basis *ba, Basis *bb, Basis *bc){
register int i,j,k;
const int qbc = qb*qc;
int ll;
double *wa,*wb,*wc,*H,*fh,*fha,*fa, *wk = Tet_Iprod_wk; //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?)
Edge *e = E->edge;
Face *f = E->face;
Mode *m,**m1;
E->state = 't';
/* multiply physical solution by jacobean */
if(E->curvX)
dvmul(qtot,E->geom->jac.p,1,h_3d[0][0],1,E->h_3d[0][0],1);
else
dsmul(qtot,E->geom->jac.d,h_3d[0][0],1,E->h_3d[0][0],1);
/* assuming at present that quadrature order if fixed */
getzw(E->qa,&wa,&wa,'a');
getzw(E->qb,&wb,&wb,'b');
getzw(E->qc,&wc,&wc,'c');
for(i = 0; i < qb*qc; ++i)
dvmul(qa, wa, 1, E->h_3d[0][i], 1, E->h_3d[0][i], 1);
for(i = 0; i < qc; ++i)
for(j = 0; j < qb; ++j)
dsmul(qa, wb[j], E->h_3d[i][j], 1, E->h_3d[i][j], 1);
for(i = 0; i < qc; ++i)
dsmul(qa*qb, wc[i], E->h_3d[i][0], 1, E->h_3d[i][0], 1);
fh = wk;
fha = fh + (LGmax+1)*QGmax*QGmax;
/* calculate fh[l][j][k] */
/* edge 6, vertex C and D */
for(k = 0, H = E->h_3d[0][0]; k < qc; ++k)
for(j = 0; j < qb; ++j, H+=qa)
fh[j*qc + k] = ddot(qa,ba->vert[2].a,1,H,1);
/* face 3, edges 2 and 5, vertex B */
fh += qbc;
for(k = 0, H = E->h_3d[0][0]; k < qc; ++k)
for(j = 0; j < qb; ++j, H+=qa)
fh[j*qc+k] = ddot(qa,ba->vert[1].a,1,H,1);
/* face 4, edges 3 and 4, vertex A */
fh += qbc;
for(k = 0, H = E->h_3d[0][0]; k < qc; ++k)
for(j = 0; j < qb; ++j, H += qa)
fh[j*qc+k] = ddot(qa,ba->vert[0].a,1,H,1);
/* interior, faces 1 and 2, edge 1 */
ll= max(E->interior_l,max(f[0].l,max(f[1].l,e[0].l)));
fh += qbc;
m = ba->edge[0];
for(i = 0; i < ll; ++i,fh+=qbc)
for(k = 0,H = E->h_3d[0][0]; k < qc; ++k)
for(j = 0; j < qb; ++j,H+=qa)
fh[j*qc+k] = ddot(qa,m[i].a,1,H,1);
/* calculate fhat[l][m][k] needs to be done mulitply for each *
* different type of g^2(b) since l may vary */
/* interior and face 1 */
ll = max(f[0].l,E->interior_l);
m1 = bb->face[0];
fh = wk+3*qbc;
for(i = 0,fa = fha; i < ll; fa += (ll-i)*qc, ++i) //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?); fa: Curvi.c(?), Coords.c(?)
for(j = 0; j < ll-i; ++j)
for(k = 0; k < qc; ++k)
fa[j*qc+k] = ddot(qb,m1[i][j].b,1,fh+i*qbc+k,qc); //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?)
/* inner product of interior */
for(i = 0,fa = fha; i < E->interior_l; fa += (ll-i)*qc, ++i) //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?); fa: Curvi.c(?), Coords.c(?)
for(j = 0; j < E->interior_l-i; ++j)
for(k = 0; k < E->interior_l-i-j; ++k)
E->hj_3d[i][j][k] = ddot(qc,fa+j*qc,1,bc->intr[i][j][k].c,1); //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?)
/* inner product of face 1 */
for(i = 0,fa = fha; i < f[0].l; fa += (ll-i)*qc, ++i) //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?); fa: Curvi.c(?), Coords.c(?)
for(j = 0; j < f[0].l-i; ++j)
f[0].hj[i][j] = ddot(qc,fa+j*qc,1,bc->face[0][i][j].c,1); //OVERLOAD CALL: fa: Curvi.c(?), Coords.c(?)
/* face 2 and edge 1*/
ll = max(f[1].l,e[0].l);
m = bb->edge[0];
for(i = 0; i < ll; ++i,fh+=qbc)
for(k = 0; k < qc; ++k)
fha[i*qc+k] = ddot(qb,m[i].b,1,fh+k,qc);
/* inner product of face 2 */
for(i = 0; i < f[1].l; ++i)
for(k = 0; k < f[1].l-i; ++k)
f[1].hj[i][k] = ddot(qc,fha+i*qc,1,bc->face[1][i][k].c,1);
/* inner product of edge 1 */
for(i = 0; i < e[0].l; ++i)
e[0].hj[i] = ddot(qc,fha+i*qc,1,bc->edge[0][i].c,1);
/* faces 3, edge 2*/
ll = max(f[2].l,e[1].l);
m = bb->edge[1];
fh = wk + qbc;
for(i = 0; i < ll; ++i)
for(k = 0; k < qc; ++k)
fha[i*qc+k] = ddot(qb,m[i].b,1,fh+k,qc);
/* inner product of face 3 */
for(i = 0; i < f[2].l; ++i)
for(k = 0; k < f[2].l-i; ++k)
f[2].hj[i][k] = ddot(qc,fha+i*qc,1,bc->face[2][i][k].c,1);
/* inner product of edge 2 */
for(i = 0; i < e[1].l; ++i)
e[1].hj[i] = ddot(qc,fha+i*qc,1,bc->edge[1][i].c,1);
/* f1 for faces 4, edges 3*/
ll= max(f[3].l,e[2].l);
m = bb->edge[2];
fh += qbc;
for(i = 0; i < ll; ++i)
for(k = 0; k < qc; ++k)
fha[i*qc+k] = ddot(qb,m[i].b,1,fh+k,qc);
/* inner product of face 4 */
for(i = 0; i < f[3].l; ++i)
for(k = 0; k < f[3].l-i; ++k)
f[3].hj[i][k] = ddot(qc,fha+i*qc,1,bc->face[3][i][k].c,1);
/* inner product of edge 3 */
for(i = 0; i < e[2].l; ++i)
e[2].hj[i] = ddot(qc,fha+i*qc,1,bc->edge[2][i].c,1);
/* f1 for edge 4, vertex A*/
for(k = 0; k < qc; ++k)
fha[k] = ddot(qb,bb->vert[0].b,1,fh+k,qc);
/* inner product of edge 4 */
for(i = 0; i < e[3].l; ++i)
e[3].hj[i] = ddot(qc,fha,1,bc->edge[3][i].c,1);
/* inner product of vertex A */
E->vert[0].hj[0] = ddot(qc,fha,1,bc->vert[0].c,1);
/* f1 for edge 5, vertex B*/
fh = wk + qbc;
for(k = 0; k < qc; ++k)
fha[k] = ddot(qb,bb->vert[0].b,1,fh+k,qc);
/* inner product of edge 5 */
for(i = 0; i < e[4].l; ++i)
e[4].hj[i] = ddot(qc,fha,1,bc->edge[4][i].c,1);
/* inner product of vertex B */
E->vert[1].hj[0] = ddot(qc,fha,1,bc->vert[1].c,1);
/* f1 for edge 6 and vertex C */
for(k = 0; k < qc; ++k)
fha[k] = ddot(qb,bb->vert[2].b,1,wk+k,qc);
/* calculate interior product of edge 6 */
for(i = 0; i < e[5].l; ++i)
e[5].hj[i] = ddot(qc,fha,1,bc->edge[5][i].c,1);
/* inner product of vertex C */
E->vert[2].hj[0] = ddot(qc,fha,1,bc->vert[2].c,1);
/* f1 for vertex D */
for(k = 0; k < qc; ++k)
fha[k] = ddot(qb,bb->vert[3].b,1,wk+k,qc);
/* inner product of vertex D */
E->vert[3].hj[0] = ddot(qc,fha,1,bc->vert[3].c,1);
}
void Pyr::Iprod_3d(Element *E, Basis *ba, Basis *bb, Basis *bc){
register int i,j,k;
double *wa,*wb,*wc;
E->state = 't';
dvmul(qtot,geom->jac.p,1,h_3d[0][0],1,E->h_3d[0][0],1);
getzw(E->qa,&wa,&wa,'a');
getzw(E->qb,&wb,&wb,'a');
getzw(E->qc,&wc,&wc,'c'); // changed
for(i = 0; i < qb*qc; ++i)
dvmul(qa, wa, 1, E->h_3d[0][i], 1, E->h_3d[0][i], 1);
for(i = 0; i < qc; ++i)
for(j = 0; j < qb; ++j)
dsmul(qa, wb[j], E->h_3d[i][j], 1, E->h_3d[i][j], 1);
for(i = 0; i < qc; ++i)
dsmul(qa*qb, wc[i], E->h_3d[i][0], 1, E->h_3d[i][0], 1);
Mode *bv = bc->vert;
Mode **be = bc->edge;
Mode ***bf = bc->face;
Mode ***bi = bc->intr;
Vert *v = E->vert;
Edge *e = E->edge;
Face *f = E->face;
double *H = Pyr_Iprod_wk;
dgemm('T','N',lmax,qb*qc,qa,1.0,ba->vert[0].a,qa,**E->h_3d,qa,0.0,H,lmax);
for(i=0;i<qc;++i)
dgemm('N','N',lmax,lmax,qb,1.,H+i*lmax*qb,lmax, bb->vert[0].b, qb,0.0, (**E->h_3d)+i*lmax*lmax,lmax);
H = (**E->h_3d);
v[0].hj[0] = ddot(qc, bv[0].c, 1, H, lmax*lmax);
v[4].hj[0] = ddot(qc, bv[4].c, 1, H, lmax*lmax);
for(i=0;i<e[4].l;++i)
e[4].hj[i] = ddot(qc, be[4][i].c, 1, H, lmax*lmax);
H += 1;
for(i=0;i<e[0].l;++i)
e[0].hj[i] = ddot(qc, be[0][i].c, 1, H+i, lmax*lmax);
for(i=0;i<f[1].l;++i)
for(j=0;j<f[1].l-i;++j)
f[1].hj[i][j] = ddot(qc, bf[1][i][j].c, 1, H+i, lmax*lmax);
H += lmax-2;
v[1].hj[0] = ddot(qc, bv[1].c, 1, H, lmax*lmax);
v[4].hj[0] += ddot(qc, bv[4].c, 1, H, lmax*lmax);
for(i=0;i<e[5].l;++i)
e[5].hj[i] = ddot(qc, be[5][i].c, 1, H, lmax*lmax);
H = (**E->h_3d) + lmax;
for(i=0;i<e[3].l;++i)
e[3].hj[i] = ddot(qc, be[3][i].c, 1, H+i*lmax, lmax*lmax);
for(i=0;i<f[4].l;++i)
for(j=0;j<f[4].l-i;++j)
f[4].hj[i][j] = ddot(qc, bf[4][i][j].c, 1, H+i*lmax, lmax*lmax);
H += 1;
int il = interior_l;
for(i=0;i<il;++i)
for(j=0;j<il-i;++j)
for(k=0;k<il-i-j;++k)
hj_3d[i][j][k] = ddot(qc, bi[i][j][k].c, 1, H+i+j*lmax, lmax*lmax);
for(i=0;i<f[0].l;++i)
for(j=0;j<f[0].l;++j)
f[0].hj[i][j] = ddot(qc, bf[0][i][j].c, 1, H+j+i*lmax, lmax*lmax);
H += (lmax-2);
for(i=0;i<e[1].l;++i)
e[1].hj[i] = ddot(qc, be[1][i].c, 1, H+i*lmax, lmax*lmax);
for(i=0;i<f[2].l;++i)
for(j=0;j<f[2].l-i;++j)
f[2].hj[i][j] = ddot(qc, bf[2][i][j].c, 1, H+i*lmax, lmax*lmax);
H = (**E->h_3d) + lmax*(lmax-1); // fill in top layer
v[3].hj[0] = ddot(qc, bv[3].c, 1, H, lmax*lmax);
v[4].hj[0] += ddot(qc, bv[4].c, 1, H, lmax*lmax);
for(i=0;i<e[7].l;++i)
e[7].hj[i] = ddot(qc, be[7][i].c, 1, H, lmax*lmax);
H += 1;
for(i=0;i<e[2].l;++i)
e[2].hj[i] = ddot(qc, be[2][i].c, 1, H+i, lmax*lmax);
for(i=0;i<f[3].l;++i)
for(j=0;j<f[3].l-i;++j)
f[3].hj[i][j] = ddot(qc, bf[3][i][j].c, 1, H+i, lmax*lmax);
H += (lmax-2);
v[2].hj[0] = ddot(qc, bv[2].c, 1, H, lmax*lmax);
v[4].hj[0] += ddot(qc, bv[4].c, 1, H, lmax*lmax);
for(i=0;i<e[6].l;++i)
e[6].hj[i] = ddot(qc, be[6][i].c, 1, H, lmax*lmax);
}
void Prism::Iprod_3d(Element *E, Basis *ba, Basis *bb, Basis *bc){
register int i,j,k;
double *wa,*wb,*wc;
E->state = 't';
dvmul(qtot,geom->jac.p,1,h_3d[0][0],1,E->h_3d[0][0],1);
getzw(E->qa,&wa,&wa,'a');
getzw(E->qb,&wb,&wb,'a');
getzw(E->qc,&wc,&wc,'b'); // changed
for(i = 0; i < qb*qc; ++i)
dvmul(qa, wa, 1, E->h_3d[0][i], 1, E->h_3d[0][i], 1);
for(i = 0; i < qc; ++i)
for(j = 0; j < qb; ++j)
dsmul(qa, wb[j], E->h_3d[i][j], 1, E->h_3d[i][j], 1);
for(i = 0; i < qc; ++i)
dsmul(qa*qb, wc[i], E->h_3d[i][0], 1, E->h_3d[i][0], 1);
Mode *bv = bc->vert;
Mode **be = bc->edge;
Mode ***bf = bc->face;
Mode ***bi = bc->intr;
Vert *v = E->vert;
Edge *e = E->edge;
Face *f = E->face;
double *H = Prism_Iprod_wk;
dgemm('T','N',lmax,qb*qc,qa,1.0,ba->vert[0].a,qa,**E->h_3d,qa,0.0,H,lmax);
for(i=0;i<qc;++i)
dgemm('N','N',lmax,lmax,qb,1.,H+i*lmax*qb,lmax, bb->vert[0].b, qb,0.0, (**E->h_3d)+i*lmax*lmax,lmax);
H = (**E->h_3d);
v[0].hj[0] = ddot(qc, bv[0].c, 1, H, lmax*lmax);
v[4].hj[0] = ddot(qc, bv[4].c, 1, H, lmax*lmax);
for(i=0;i<e[4].l;++i)
e[4].hj[i] = ddot(qc, be[4][i].c, 1, H, lmax*lmax);
H += 1;
for(i=0;i<e[0].l;++i)
e[0].hj[i] = ddot(qc, be[0][i].c, 1, H+i, lmax*lmax);
for(i=0;i<f[1].l;++i)
for(j=0;j<f[1].l-i;++j)
f[1].hj[i][j] = ddot(qc, bf[1][i][j].c, 1, H+i, lmax*lmax);
H += lmax-2;
v[1].hj[0] = ddot(qc, bv[1].c, 1, H, lmax*lmax);
v[4].hj[0] += ddot(qc, bv[4].c, 1, H, lmax*lmax);
for(i=0;i<e[5].l;++i)
e[5].hj[i] = ddot(qc, be[5][i].c, 1, H, lmax*lmax);
H = (**E->h_3d) + lmax;
for(i=0;i<e[3].l;++i)
e[3].hj[i] = ddot(qc, be[3][i].c, 1, H+i*lmax, lmax*lmax);
for(i=0;i<e[8].l;++i)
e[8].hj[i] = ddot(qc, be[8][i].c, 1, H+i*lmax, lmax*lmax);
for(i=0;i<f[4].l;++i)
for(j=0;j<f[4].l;++j)
f[4].hj[i][j] = ddot(qc, bf[4][i][j].c, 1, H+j*lmax, lmax*lmax);
H += 1;
int il = interior_l;
for(i=0;i<il-1;++i)
for(j=0;j<il;++j)
for(k=0;k<il-1-i;++k)
hj_3d[i][j][k] = ddot(qc, bi[i][j][k].c, 1, H+i+j*lmax, lmax*lmax);
for(i=0;i<f[0].l;++i)
for(j=0;j<f[0].l;++j)
f[0].hj[i][j] = ddot(qc, bf[0][i][j].c, 1, H+j+i*lmax, lmax*lmax);
H += (lmax-2);
for(i=0;i<e[1].l;++i)
e[1].hj[i] = ddot(qc, be[1][i].c, 1, H+i*lmax, lmax*lmax);
for(i=0;i<e[8].l;++i)
e[8].hj[i] += ddot(qc, be[8][i].c, 1, H+i*lmax, lmax*lmax);
for(i=0;i<f[2].l;++i)
for(j=0;j<f[2].l;++j)
f[2].hj[i][j] = ddot(qc, bf[2][i][j].c, 1, H+j*lmax, lmax*lmax);
H = (**E->h_3d) + lmax*(lmax-1); // fill in top layer
v[3].hj[0] = ddot(qc, bv[3].c, 1, H, lmax*lmax);
v[5].hj[0] = ddot(qc, bv[5].c, 1, H, lmax*lmax);
for(i=0;i<e[7].l;++i)
e[7].hj[i] = ddot(qc, be[7][i].c, 1, H, lmax*lmax);
H += 1;
for(i=0;i<e[2].l;++i)
e[2].hj[i] = ddot(qc, be[2][i].c, 1, H+i, lmax*lmax);
for(i=0;i<f[3].l;++i)
for(j=0;j<f[3].l-i;++j)
f[3].hj[i][j] = ddot(qc, bf[3][i][j].c, 1, H+i, lmax*lmax);
H += (lmax-2);
v[2].hj[0] = ddot(qc, bv[2].c, 1, H, lmax*lmax);
v[5].hj[0] += ddot(qc, bv[5].c, 1, H, lmax*lmax);
for(i=0;i<e[6].l;++i)
e[6].hj[i] = ddot(qc, be[6][i].c, 1, H, lmax*lmax);
}
void Hex::Iprod_3d(Element *E, Basis *ba, Basis *bb, Basis *bc){
register int i,j;
double *wa,*wb,*wc;
E->state = 't';
dvmul(qtot,geom->jac.p,1,h_3d[0][0],1,E->h_3d[0][0],1);
getzw(E->qa,&wa,&wa,'a');
getzw(E->qb,&wb,&wb,'a');
getzw(E->qb,&wc,&wc,'a');
for(i = 0; i < qb*qc; ++i)
dvmul(qa, wa, 1, E->h_3d[0][i], 1, E->h_3d[0][i], 1);
for(i = 0; i < qc; ++i)
for(j = 0; j < qb; ++j)
dsmul(qa, wb[j], E->h_3d[i][j], 1, E->h_3d[i][j], 1);
for(i = 0; i < qc; ++i)
dsmul(qa*qb, wc[i], E->h_3d[i][0], 1, E->h_3d[i][0], 1);
double *H = Hex_Iprod_wk;
dgemm('T','N',lmax,qb*qc,qa,1.0,ba->vert[0].a,qa,**E->h_3d,qa,0.0,H,lmax);
for(i=0;i<qc;++i)
dgemm('N','N',lmax,lmax,qb,1.,H+i*lmax*qb,lmax, bb->vert[0].b, qb,0.0, (**E->h_3d)+i*lmax*lmax,lmax);
dgemm('N','N',lmax*lmax,lmax,qc,1.0, **E->h_3d, lmax*lmax, bc->vert[0].c, qc, 0.0, H, lmax*lmax);
Hex_map_to_hj(this, H);
}
static void Hex_map_to_hj(Element *E, double *d){
int lm = E->lmax;
int i,j,k;
E->vert[0].hj[0] = d[0];
E->vert[1].hj[0] = d[lm-1];
E->vert[2].hj[0] = d[lm*lm-1];
E->vert[3].hj[0] = d[lm*(lm-1)];
E->vert[4].hj[0] = d[lm*lm*(lm-1)+0];
E->vert[5].hj[0] = d[lm*lm*(lm-1)+lm-1];
E->vert[6].hj[0] = d[lm*lm*(lm-1)+lm*lm-1];
E->vert[7].hj[0] = d[lm*lm*(lm-1)+lm*(lm-1)];
dcopy(E->edge[0].l, d+1, 1, E->edge[0].hj, 1);
dcopy(E->edge[1].l, d+lm+lm-1, lm, E->edge[1].hj, 1);
dcopy(E->edge[2].l, d+lm*(lm-1)+1, 1, E->edge[2].hj, 1);
dcopy(E->edge[3].l, d+lm, lm, E->edge[3].hj, 1);
dcopy(E->edge[4].l, d+lm*lm, lm*lm, E->edge[4].hj, 1);
dcopy(E->edge[5].l, d+lm*lm+lm-1, lm*lm, E->edge[5].hj, 1);
dcopy(E->edge[6].l, d+lm*lm+lm*lm-1, lm*lm, E->edge[6].hj, 1);
dcopy(E->edge[7].l, d+lm*lm+lm*(lm-1), lm*lm, E->edge[7].hj, 1);
dcopy(E->edge[8].l, d+lm*lm*(lm-1)+1, 1, E->edge[8].hj, 1);
dcopy(E->edge[9].l, d+lm*lm*(lm-1)+lm+lm-1, lm, E->edge[9].hj, 1);
dcopy(E->edge[10].l, d+lm*lm*(lm-1)+lm*(lm-1)+1, 1, E->edge[10].hj, 1);
dcopy(E->edge[11].l, d+lm*lm*(lm-1)+lm, lm, E->edge[11].hj, 1);
// faces
for(i=0;i<E->face[0].l;++i)
dcopy(E->face[0].l, d+lm+1+i*lm, 1, E->face[0].hj[i], 1);
for(i=0;i<E->face[1].l;++i)
dcopy(E->face[1].l, d+lm*lm+1+i*lm*lm, 1, E->face[1].hj[i], 1);
for(i=0;i<E->face[2].l;++i)
dcopy(E->face[2].l, d+lm*lm+lm+lm-1+i*lm*lm, lm,E->face[2].hj[i], 1 );
for(i=0;i<E->face[3].l;++i)
dcopy(E->face[3].l, d+lm*lm+lm*(lm-1)+1+i*lm*lm, 1, E->face[3].hj[i], 1);
for(i=0;i<E->face[4].l;++i)
dcopy(E->face[4].l, d+lm*lm+lm+i*lm*lm, lm, E->face[4].hj[i], 1);
for(i=0;i<E->face[5].l;++i)
dcopy(E->face[5].l, d+lm*lm*(lm-1)+lm+1+i*lm, 1, E->face[5].hj[i], 1);
// interior
for(i=0;i<E->interior_l;++i)
for(j=0;j<E->interior_l;++j)
for(k=0;k<E->interior_l;++k)
E->hj_3d[i][j][k] = d[lm*lm+lm+1+k+j*lm+i*lm*lm];
}
/*
Function name: Element::iprod
Function Purpose:
Argument 1: Mode *x
Purpose:
Argument 2: Mode *y
Purpose:
Function Notes:
*/
double Tri::iprod(Mode *x, Mode *y){
double d;
d = ddot(qa,x->a,1,y->a,1);
d *= ddot(qb,x->b,1,y->b,1);
return d;
}
double Quad::iprod(Mode *x, Mode *y){
double d = geom->jac.d;
d *= ddot(qa,x->a,1,y->a,1);
d *= ddot(qb,x->b,1,y->b,1);
return d;
}
double Tet::iprod(Mode *x, Mode *y){
double d;
d = ddot(qa,x->a,1,y->a,1);
d *= ddot(qb,x->b,1,y->b,1);
d *= ddot(qc,x->c,1,y->c,1);
return d;
}
double Pyr::iprod(Mode *x, Mode *y){
double d;
d = ddot(qa,x->a,1,y->a,1);
d *= ddot(qb,x->b,1,y->b,1);
d *= ddot(qc,x->c,1,y->c,1);
return d;
}
double Prism::iprod(Mode *x, Mode *y){
double d;
d = ddot(qa,x->a,1,y->a,1);
d *= ddot(qb,x->b,1,y->b,1);
d *= ddot(qc,x->c,1,y->c,1);
return d;
}
double Hex::iprod(Mode *x, Mode *y){
double d;
d = ddot(qa,x->a,1,y->a,1);
d *= ddot(qb,x->b,1,y->b,1);
d *= ddot(qc,x->c,1,y->c,1);
return d;
}
double Element::iprod(Mode *, Mode *){return 0.0;}
/*
Function name: Element::iprodlap
Function Purpose:
Argument 1: Mode *x
Purpose:
Argument 2: Mode *y
Purpose:
Argument 3: Mode *fac
Purpose:
Function Notes:
*/
double Tri::iprodlap(Mode *x, Mode *y, Mode *fac){
double prod,s1,s2,s3,s4,s5;
double g1 = geom->rx.d*geom->rx.d + geom->ry.d*geom->ry.d;
double g2 = geom->sx.d*geom->sx.d + geom->sy.d*geom->sy.d;
double g4 = geom->rx.d*geom->sx.d + geom->ry.d*geom->sy.d;
double g3 = 0,g5 = 0,g6=0;
double *tmp = dvector(0,QGmax-1);
prod = g1*ddot(qa,x[0].a,1,y[0].a,1);
dvmul(qa,fac[0].a,1,x[0].a,1,tmp,1);
prod += 2*(g4+g5)*ddot(qa,tmp,1,y[0].a,1);
dvmul(qa,fac[0].a,1,tmp,1,tmp,1);
prod += (g2+2*g6+g3)*ddot(qa,tmp,1,y[0].a,1);
prod *= ddot(qb,x[0].b,1,y[0].b,1);
s1 = ddot(qa,x[0].a,1,y[1].a,1);
dvmul(qa,fac[0].a,1,x[0].a,1,tmp,1);
s2 = ddot(qa,tmp,1,y[1].a,1);
s3 = ddot(qa,y[0].a,1,x[1].a,1);
dvmul(qa,fac[0].a,1,y[0].a,1,tmp,1);
s4 = ddot(qa,tmp,1,x[1].a,1);
s5 = ddot(qa,x[1].a,1,y[1].a,1);
prod += (g4*s1 + (g2+g6)*s2)*ddot(qb,x[0].b,1,y[1].b,1);
prod += (g4*s3 + (g2+g6)*s4)*ddot(qb,y[0].b,1,x[1].b,1);
prod += g2*ddot(qb,x[1].b,1,y[1].b,1)*s5;
free(tmp);
return prod;
}
double Quad::iprodlap(Mode *x, Mode *y, Mode *fac){
double prod,s1,s2,s3,s4,s5;
double g1 = geom->rx.d*geom->rx.d + geom->ry.d*geom->ry.d;
double g2 = geom->sx.d*geom->sx.d + geom->sy.d*geom->sy.d;
double g4 = geom->rx.d*geom->sx.d + geom->ry.d*geom->sy.d;
double g3 = 0,g5 = 0,g6=0;
double *tmp =dvector(0,QGmax-1);
prod = g1*ddot(qa,x[0].a,1,y[0].a,1);
dvmul(qa,fac[0].a,1,x[0].a,1,tmp,1);
prod += 2*(g4+g5)*ddot(qa,tmp,1,y[0].a,1);
dvmul(qa,fac[0].a,1,tmp,1,tmp,1);
prod += (g2+2*g6+g3)*ddot(qa,tmp,1,y[0].a,1);
prod *= ddot(qb,x[0].b,1,y[0].b,1);
s1 = ddot(qa,x[0].a,1,y[1].a,1);
dvmul(qa,fac[0].a,1,x[0].a,1,tmp,1);
s2 = ddot(qa,tmp,1,y[1].a,1);
s3 = ddot(qa,y[0].a,1,x[1].a,1);
dvmul(qa,fac[0].a,1,y[0].a,1,tmp,1);
s4 = ddot(qa,tmp,1,x[1].a,1);
s5 = ddot(qa,x[1].a,1,y[1].a,1);
prod += (g4*s1 + (g2+g6)*s2)*ddot(qb,x[0].b,1,y[1].b,1);
prod += (g4*s3 + (g2+g6)*s4)*ddot(qb,y[0].b,1,x[1].b,1);
prod += g2*ddot(qb,x[1].b,1,y[1].b,1)*s5;
free(tmp);
return prod;
}
double Tet::iprodlap(Mode *, Mode *, Mode *){
fprintf(stderr,"Tet::iprodlap not implemented\n");
return 0.0;
}
double Pyr::iprodlap(Mode *, Mode *, Mode *){
fprintf(stderr,"Pyr::iprodlap not implemented\n");
return 0.0;
}
double Prism::iprodlap(Mode *, Mode *, Mode *){
fprintf(stderr,"Prism::iprodlap not implemented\n");
return 0.0;
}
double Hex::iprodlap(Mode *, Mode *, Mode *){
fprintf(stderr,"Hex::iprodlap not implemented\n");
return 0.0;
}
double Element::iprodlap(Mode *, Mode *, Mode *){return 0.0;}
/*
Function name: Element::HelmHoltz
Function Purpose:
Argument 1: Metric *lambda
Purpose:
Function Notes:
*/
void Tri::HelmHoltz(Metric *lambda){
int Nm,qt;
Basis *b,*db;
double *store,*u1,*u2;
store = Tri_wk+QGmax*QGmax;
u1 = Tri_wk+2*QGmax*QGmax;
u2 = Tri_wk+3*QGmax*QGmax;
Nm = Nmodes;
b = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
db = derbasis(); //OVERLOAD CALL: derbasis: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
qt = qa*qb;
Grad_d(u1, u2, 0, 'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
if(lambda){
if(lambda->wave){
dsmul(qt,lambda->d, *h,1,*h,1);
dvvtvp(qt,lambda->wave[0],1,u1,1,*h,1,*h,1);
dvvtvp(qt,lambda->wave[1],1,u2,1,*h,1,*h,1);
Iprod_d(this,b,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(Nm,vert->hj,1,store,1);
}
else{
Iprod_d(this,b,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dsmul(Nm, lambda->d, vert->hj, 1, store, 1);
}
}
else if(lambda->wave){ /* zero lambda case */
dvmul(qt,lambda->wave[0],1,u1,1,*h,1);
dvvtvp(qt,lambda->wave[1],1,u2,1,*h,1,*h,1);
Iprod_d(this,b,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(Nm,vert->hj,1,store,1);
}
else
dzero(Nm,store,1);
if(option("AXISYM")){
// assume lambda->p holds r
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
dvmul(qtot, lambda->p, 1, u2, 1, u2, 1);
}
else if(lambda && lambda->p)
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
form_diprod(u1,u2,NULL,b->vert+1); //OVERLOAD CALL: form_diprod: InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet)
dcopy(qt, u1, 1, *h, 1);
Iprod_d(this,db,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u2,1,*h,1);
Iprod_d((Element*)this,b,db); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,vert->hj,1);
state = 't';
// free(store); free(u1);
}
void Quad::HelmHoltz(Metric *lambda){
int Nm,qt;
Basis *b,*db;
double *store,*u1,*u2;
store = dvector(0,Nmodes-1);
u1 = dvector(0,2*QGmax*QGmax-1);
u2 = u1 + QGmax*QGmax;
Nm = Nmodes;
b = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
db = derbasis(); //OVERLOAD CALL: derbasis: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
qt = qa*qb;
Grad_d(u1, u2, 0, 'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
if(lambda){
if(lambda->wave){
dsmul(qt,lambda->d, *h,1,*h,1);
dvvtvp(qt,lambda->wave[0],1,u1,1,*h,1,*h,1);
dvvtvp(qt,lambda->wave[1],1,u2,1,*h,1,*h,1);
Iprod_d(this,b,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(Nm,vert->hj,1,store,1);
}
else{
Iprod_d(this,b,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dsmul(Nm, lambda->d, vert->hj, 1, store, 1);
}
}
else if(lambda->wave){ /* zero lambda case */
dvmul(qt,lambda->wave[0],1,u1,1,*h,1);
dvvtvp(qt,lambda->wave[1],1,u2,1,*h,1,*h,1);
Iprod_d(this,b,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(Nm,vert->hj,1,store,1);
}
else
dzero(Nm,store,1);
if(option("AXISYM")){
// assume lambda->p holds r
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
dvmul(qtot, lambda->p, 1, u2, 1, u2, 1);
}
else if(lambda&&lambda->p) // add sponge
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
form_diprod(u1,u2,NULL,(Mode*)NULL); //OVERLOAD CALL: form_diprod: InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet)
dcopy(qt,u1,1,*h,1);
Iprod_d(this,db,b); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u2,1,*h,1);
Iprod_d(this,b,db); //OVERLOAD CALL: Iprod_d: InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,vert->hj,1);
state = 't';
free(store); free(u1);
}
void Tet::HelmHoltz(Metric *lambda){
int Nm,qt;
Basis *b,*db;
double *store,*u1,*u2,*u3;
store = Tet_HelmHoltz_wk;
u1 = store+QGmax*QGmax*QGmax;
u2 = store+2*QGmax*QGmax*QGmax;
u3 = store+3*QGmax*QGmax*QGmax;
dzero(3*QGmax*QGmax*QGmax, u1, 1);
Nm = Nmodes;
b = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
db = derbasis(); //OVERLOAD CALL: derbasis: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
qt = qtot;
Grad_d(u1, u2, u3, 'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
// add sponge
if(lambda->p)
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
form_diprod(u1,u2,u3,b->vert+1); //OVERLOAD CALL: form_diprod: InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet)
if(lambda){
Tet::Iprod_3d(this,b,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dsmul(Nm, lambda->d, vert->hj, 1, store, 1);
}
else
dzero(Nm,store,1);
dcopy(qt,u1,1,*h_3d[0],1);
Tet::Iprod_3d(this,db,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u2,1,*h_3d[0],1);
Tet::Iprod_3d(this,b,db,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u3,1,*h_3d[0],1);
Tet::Iprod_3d(this,b,b,db); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,vert->hj,1);
state = 't';
}
void Pyr::HelmHoltz(Metric *lambda){
int Nm,qt;
Basis *b,*db;
double *store,*u1,*u2,*u3;
store = Pyr_HelmHoltz_wk;
u1 = Pyr_HelmHoltz_wk+QGmax*QGmax*QGmax;
u2 = Pyr_HelmHoltz_wk+2*QGmax*QGmax*QGmax;
u3 = Pyr_HelmHoltz_wk+3*QGmax*QGmax*QGmax;
dzero(3*QGmax*QGmax*QGmax, u1, 1);
Nm = Nmodes;
b = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
db = derbasis(); //OVERLOAD CALL: derbasis: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
qt = qtot;
Grad_d(u1, u2, u3, 'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
// add sponge
if(lambda->p)
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
form_diprod(u1,u2,u3,NULL); //OVERLOAD CALL: form_diprod: InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet)
if(lambda){
Pyr::Iprod_3d(this,b,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dsmul(Nm, lambda->d, vert->hj, 1, store, 1);
}
else
dzero(Nm,store,1);
dcopy(qt,u1,1,*h_3d[0],1);
Pyr::Iprod_3d(this,db,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u2,1,*h_3d[0],1);
Pyr::Iprod_3d(this,b,db,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u3,1,*h_3d[0],1);
Pyr::Iprod_3d(this,b,b,db); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,vert->hj,1);
state = 't';
}
void Prism::HelmHoltz(Metric *lambda){
Basis *b,*db;
double *store,*u1,*u2,*u3;
store = Prism_HelmHoltz_wk;
u1 = Prism_HelmHoltz_wk+QGmax*QGmax*QGmax;
u2 = Prism_HelmHoltz_wk+2*QGmax*QGmax*QGmax;
u3 = Prism_HelmHoltz_wk+3*QGmax*QGmax*QGmax;
dzero(3*QGmax*QGmax*QGmax, u1, 1);
b = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
db = derbasis(); //OVERLOAD CALL: derbasis: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
Grad_d(u1, u2, u3, 'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
// add sponge
if(lambda->p)
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
form_diprod(u1,u2,u3,NULL); //OVERLOAD CALL: form_diprod: InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet)
if(lambda){
Prism::Iprod_3d(this,b,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dsmul(Nmodes, lambda->d, vert->hj, 1, store, 1);
}
else
dzero(Nmodes,store,1);
dcopy(qtot,u1,1,*h_3d[0],1);
Prism::Iprod_3d(this,db,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nmodes,vert->hj,1,store,1,store,1);
dcopy(qtot,u2,1,*h_3d[0],1);
Prism::Iprod_3d(this,b,db,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nmodes,vert->hj,1,store,1,store,1);
dcopy(qtot,u3,1,*h_3d[0],1);
Prism::Iprod_3d(this,b,b,db); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nmodes,vert->hj,1,store,1,vert->hj,1);
state = 't';
}
void Hex::HelmHoltz(Metric *lambda){
int Nm,qt;
Basis *b,*db;
double *store,*u1,*u2,*u3;
store = Hex_HelmHoltz_wk;
u1 = store+QGmax*QGmax*QGmax;
u2 = store+2*QGmax*QGmax*QGmax;
u3 = store+3*QGmax*QGmax*QGmax;
dzero(3*QGmax*QGmax*QGmax, u1, 1);
Nm = Nmodes;
b = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
db = derbasis(); //OVERLOAD CALL: derbasis: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
qt = qtot;
Grad_d(u1, u2, u3, 'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
// add sponge
if(lambda->p)
dvmul(qtot, lambda->p, 1, u1, 1, u1, 1);
form_diprod(u1,u2,u3,NULL); //OVERLOAD CALL: form_diprod: InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet)
if(lambda){
Hex::Iprod_3d(this,b,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dsmul(Nm, lambda->d, vert->hj, 1, store, 1);
}
else
dzero(Nm,store,1);
dcopy(qt,u1,1,*h_3d[0],1);
Hex::Iprod_3d(this,db,b,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u2,1,*h_3d[0],1);
Hex::Iprod_3d(this,b,db,b); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,store,1);
dcopy(qt,u3,1,*h_3d[0],1);
Hex::Iprod_3d(this,b,b,db); //OVERLOAD CALL: Iprod_3d: InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dvadd(Nm,vert->hj,1,store,1,vert->hj,1);
state = 't';
}
void Element::HelmHoltz(Metric *){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
/*
Function name: Element::form_diprod
Function Purpose:
Argument 1: double *u1
Purpose:
Argument 2: double *u2
Purpose:
Argument 3: double *
Purpose:
Argument 4: Mode *m
Purpose:
Function Notes:
*/
void Tri::form_diprod(double *u1, double *u2, double *, Mode *m)
{
register int i;
const int qt = qa*qb;
Geom *G = geom;
double *tmp = Tri_wk+4*QGmax*QGmax; // sort out later
double *Ur = tmp;
// Grad_d(u1,u2,(double*)NULL,'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
if(curvX){
dvmul (qt,G->rx.p,1,u1,1,Ur,1);
dvvtvp(qt,G->ry.p,1,u2,1,Ur,1,Ur,1);
dvmul (qt,G->sy.p,1,u2,1,u2,1);
dvvtvp(qt,G->sx.p,1,u1,1,u2,1,u2,1);
}
else{
if(G->rx.d) dsmul(qt,G->rx.d,u1,1,Ur,1);
else dzero(qt,Ur,1);
if(G->ry.d) daxpy(qt,G->ry.d,u2,1,Ur,1);
if(G->sy.d) dscal(qt,G->sy.d,u2,1);
else dzero(qt,u2,1);
if(G->sx.d) daxpy(qt,G->sx.d,u1,1,u2,1);
}
for(i = 0; i < qa; ++i)
dsmul(qb,m[0].a[i],u2+i,qa,u1+i,qa);
dvadd(qt,Ur,1,u1,1,u1,1);
for(i = 0; i < qb; ++i)
dscal(qa,1.0/m[0].b[i],u1+i*qa,1);
// free(tmp);
}
void Quad::form_diprod(double *u1, double *u2, double *, Mode *)
{
const int qt = qa*qb;
Geom *G = geom;
double *Ur = dvector(0,qt-1);
// Grad_d(u1,u2,(double*)NULL,'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (qt,G->rx.p,1,u1,1,Ur,1);
dvvtvp(qt,G->ry.p,1,u2,1,Ur,1,Ur,1);
dvmul (qt,G->sy.p,1,u2,1,u2,1);
dvvtvp(qt,G->sx.p,1,u1,1,u2,1,u2,1);
dcopy (qt,Ur,1,u1,1);
free(Ur);
}
void Tet::form_diprod(double *u1, double *u2, double *u3, Mode *m)
{
register int i,j;
Geom *G = geom;
double *Ur = Tet_form_diprod_wk;
double *Us = Ur+QGmax*QGmax*QGmax;
double *us;
// Grad_d(u1,u2,u3,'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
if(curvX){
dvmul (qtot,G->rx.p,1,u1,1,Ur,1);
dvvtvp(qtot,G->ry.p,1,u2,1,Ur,1,Ur,1);
dvvtvp(qtot,G->rz.p,1,u3,1,Ur,1,Ur,1);
dvmul (qtot,G->sx.p,1,u1,1,Us,1);
dvvtvp(qtot,G->sy.p,1,u2,1,Us,1,Us,1);
dvvtvp(qtot,G->sz.p,1,u3,1,Us,1,Us,1);
dvmul (qtot,G->tz.p,1,u3,1,u3,1);
dvvtvp(qtot,G->ty.p,1,u2,1,u3,1,u3,1);
dvvtvp(qtot,G->tx.p,1,u1,1,u3,1,u3,1);
}
else{
if(G->rx.d) dsmul(qtot,G->rx.d,u1,1,Ur,1);
else dzero(qtot,Ur,1);
if(G->ry.d) daxpy(qtot,G->ry.d,u2,1,Ur,1);
if(G->rz.d) daxpy(qtot,G->rz.d,u3,1,Ur,1);
if(G->sx.d) dsmul(qtot,G->sx.d,u1,1,Us,1);
else dzero(qtot,Us,1);
if(G->sy.d) daxpy(qtot,G->sy.d,u2,1,Us,1);
if(G->sz.d) daxpy(qtot,G->sz.d,u3,1,Us,1);
if(G->tz.d) dsmul(qtot,G->tz.d,u3,1,u3,1);
else dzero(qtot,u3,1);
if(G->ty.d) daxpy(qtot,G->ty.d,u2,1,u3,1);
if(G->tx.d) daxpy(qtot,G->tx.d,u1,1,u3,1);
}
dvadd(qtot,Us,1,u3,1,u1,1);
for(i = 0; i < qa; ++i)
dsmul(qa*qb,m[0].a[i],u1+i,qa,u1+i,qa);
dvadd(qtot,Ur,1,u1,1,u1,1);
for(i = 0; i < qc; ++i)
dsmul(qa*qb,1.0/m[0].c[i],u1+i*qa*qb,1,u1+i*qa*qb,1);
dvrecp(qb,m[0].b,1,Ur,1);
for(i = 0,us=Us; i < qc; ++i)
for(j = 0; j < qb; ++j,u1+=qa,u3+=qa,us+=qa){
dsmul(qa,Ur[j],u1,1,u1,1);
daxpy(qa,m[1].b[j],u3,1,us,1);
}
for(i = 0; i < qc; ++i)
dsmul(qa*qb,1.0/m[0].c[i],Us+i*qa*qb,1,u2+i*qa*qb,1);
}
void Pyr::form_diprod(double *u1, double *u2, double *u3, Mode *)
{
#if 1
register int i,j;
Geom *G = geom;
double *Ur = Pyr_form_diprod_wk;
double *Us = Pyr_form_diprod_wk+QGmax*QGmax*QGmax;
double *Ut = Pyr_form_diprod_wk+2*QGmax*QGmax*QGmax;
Basis *B = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
// Grad_d(u1,u2,u3,'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (qtot,G->rx.p,1,u1,1,Ur,1);
dvvtvp(qtot,G->ry.p,1,u2,1,Ur,1,Ur,1);
dvvtvp(qtot,G->rz.p,1,u3,1,Ur,1,Ur,1);
dvmul (qtot,G->sx.p,1,u1,1,Us,1);
dvvtvp(qtot,G->sy.p,1,u2,1,Us,1,Us,1);
dvvtvp(qtot,G->sz.p,1,u3,1,Us,1,Us,1);
dvmul (qtot,G->tx.p,1,u1,1,Ut,1);
dvvtvp(qtot,G->ty.p,1,u2,1,Ut,1,Ut,1);
dvvtvp(qtot,G->tz.p,1,u3,1,Ut,1,Ut,1);
// u1
for(i=0;i<qb*qc;++i)
dvvtvp(qa, B->vert[1].a, 1, Ut+i*qa, 1, Ur+i*qa, 1, Ur+i*qa, 1);
for(i=0;i<qc;++i)
dscal(qa*qb, 1./B->vert[0].c[i], Ur+i*qa*qb, 1);
// u2
for(i=0;i<qc;++i)
for(j=0;j<qb;++j)
daxpy(qa, B->vert[1].a[j], Ut+j*qa+i*qa*qb, 1, Us+j*qa+i*qa*qb, 1);
for(i=0;i<qc;++i)
dscal(qa*qb, 1./B->vert[0].c[i], Us+i*qa*qb, 1);
dcopy(qtot, Ur, 1, u1, 1);
dcopy(qtot, Us, 1, u2, 1);
dcopy(qtot, Ut, 1, u3, 1);
#endif
}
void Prism::form_diprod(double *u1, double *u2, double *u3, Mode *)
{
#if 1
register int i;
Geom *G = geom;
double *Ur = Prism_form_diprod_wk;
double *Us = Prism_form_diprod_wk+QGmax*QGmax*QGmax;
double *Ut = Prism_form_diprod_wk+2*QGmax*QGmax*QGmax;
Basis *B = getbasis(); //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
// Grad_d(u1,u2,u3,'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (qtot,G->rx.p,1,u1,1,Ur,1);
dvvtvp(qtot,G->ry.p,1,u2,1,Ur,1,Ur,1);
dvvtvp(qtot,G->rz.p,1,u3,1,Ur,1,Ur,1);
dvmul (qtot,G->sx.p,1,u1,1,Us,1);
dvvtvp(qtot,G->sy.p,1,u2,1,Us,1,Us,1);
dvvtvp(qtot,G->sz.p,1,u3,1,Us,1,Us,1);
dvmul (qtot,G->tx.p,1,u1,1,Ut,1);
dvvtvp(qtot,G->ty.p,1,u2,1,Ut,1,Ut,1);
dvvtvp(qtot,G->tz.p,1,u3,1,Ut,1,Ut,1);
for(i=0;i<qb*qc;++i)
dvvtvp(qa, B->vert[1].a, 1, Ut+i*qa, 1, Ur+i*qa, 1, Ur+i*qa, 1);
for(i=0;i<qc;++i)
dscal(qa*qb, 1./B->vert[0].c[i], Ur+i*qa*qb, 1);
dcopy(qtot, Ur, 1, u1, 1);
dcopy(qtot, Us, 1, u2, 1);
dcopy(qtot, Ut, 1, u3, 1);
#endif
}
void Hex::form_diprod(double *u1, double *u2, double *u3, Mode *)
{
Geom *G = geom;
double *Ur = Hex_form_diprod_wk;
double *Us = Hex_form_diprod_wk+QGmax*QGmax*QGmax;
double *Ut = Hex_form_diprod_wk+2*QGmax*QGmax*QGmax;
// Grad_d(u1,u2,u3,'a'); //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (qtot,G->rx.p,1,u1,1,Ur,1);
dvvtvp(qtot,G->ry.p,1,u2,1,Ur,1,Ur,1);
dvvtvp(qtot,G->rz.p,1,u3,1,Ur,1,Ur,1);
dvmul (qtot,G->sx.p,1,u1,1,Us,1);
dvvtvp(qtot,G->sy.p,1,u2,1,Us,1,Us,1);
dvvtvp(qtot,G->sz.p,1,u3,1,Us,1,Us,1);
dvmul (qtot,G->tx.p,1,u1,1,Ut,1);
dvvtvp(qtot,G->ty.p,1,u2,1,Ut,1,Ut,1);
dvvtvp(qtot,G->tz.p,1,u3,1,Ut,1,Ut,1);
dcopy(qtot, Ur, 1, u1, 1);
dcopy(qtot, Us, 1, u2, 1);
dcopy(qtot, Ut, 1, u3, 1);
}
void Element::form_diprod(double *, double *, double *, Mode *){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
C++ to HTML Conversion by ctoohtml