file: Hlib/src/Matrix.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 <math.h>
#include <veclib.h>
#include "hotel.h"
#include "Tri.h"
#include "Quad.h"
static void Bsystem_mem (Bsystem *Ubsys, Element *U);
static void invert_a (Bsystem *Ubsys); //OVERLOAD CALL: invert_a: Matrix.c(?), Matrix_Stokes.c(?)
void MemPrecon (Element *U, Bsystem *B); //OVERLOAD CALL: MemPrecon: Matrix.c(?), Matrix_Stokes.c(?)
static void FillPrecon(Element *E, Bsystem *B); //OVERLOAD CALL: FillPrecon: Matrix.c(?), Matrix_Stokes.c(?)
static void InvtPrecon(Bsystem *B); //OVERLOAD CALL: InvtPrecon: Matrix.c(?), Matrix_Stokes.c(?)
#if (defined(DUMPMASS) || defined(DUMPHELM))
static void dumpmat(Element *E,LocMat *mat);
#endif
int bandwidth(Element *E, Bsystem *Bsys);
void Tri_HelmMat(Element *T, LocMat *helm, double lambda);
static void Add_robin_matrix(Element *E, LocMat *helm, Bndry *rbc);
void GenMat(Element_List *U, Bndry *Ubc, Bsystem *Ubsys, Metric *lambda,
SolveType Stype){
LocMat *mass = 0,*helm;
const int nel = U->nel;
int na,nb,nc;
int iter = (int) Ubsys->smeth;
Bndry *Bc;
Element *E;
if(!Ubsys->signchange)
setup_signchange(U,Ubsys);
Bsystem_mem(Ubsys,U->fhead);
if (Ubsys->rslv) MemRecur(Ubsys,0,'p');
#if 0
if (Ubsys->rslv) MemRecur(Ubsys,0);
#endif
if (Ubsys->smeth == iterative) MemPrecon(U->fhead,Ubsys); //OVERLOAD CALL: MemPrecon: Matrix.c(?), Matrix_Stokes.c(?)
switch(Stype){
case Mass:
for(E = U->fhead; E;E=E->next){
mass = E->mat_mem (); //OVERLOAD CALL: mat_mem: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
if(!iter||!Ubsys->Gmat->invc[E->geom->id]){
if(!E->curvX)
E->MassMat(mass); //OVERLOAD CALL: MassMat: Matrix.c(Tri), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex)
else
E->MassMatC(mass); //OVERLOAD CALL: MassMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
E->condense(mass,Ubsys); //OVERLOAD CALL: condense: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
}
if(!iter) E->project(mass,Ubsys); //OVERLOAD CALL: project: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
E->mat_free (mass); //OVERLOAD CALL: mat_free: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
}
break;
case Helm:
for(E = U->fhead; E;E=E->next){
helm = E->mat_mem (); //OVERLOAD CALL: mat_mem: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
if(lambda){
mass = E->mat_mem (); //OVERLOAD CALL: mat_mem: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
na = mass->asize*mass->asize;
nb = mass->asize*mass->csize;
nc = mass->csize*mass->csize;
}
if(Ubsys->rslv||
(Ubsys->smeth == direct)||
!Ubsys->Gmat->invc[E->geom->id]){
if(E->curvX || lambda[E->id].p)
E->HelmMatC (helm,lambda+E->id); //OVERLOAD CALL: HelmMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
else{
if(E->identify()==Nek_Tri) //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
Tri_HelmMat(E, helm, lambda[E->id].d);
else
{
E->LapMat(helm); //OVERLOAD CALL: LapMat: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex)
if(lambda[E->id].d){
E->MassMat(mass); //OVERLOAD CALL: MassMat: Matrix.c(Tri), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex)
daxpy(na,lambda[E->id].d,*mass->a,1,*helm->a,1);
if(nc){
daxpy(nb,lambda[E->id].d,*mass->b,1,*helm->b,1);
daxpy(nc,lambda[E->id].d,*mass->c,1,*helm->c,1);
}
}
}
}
for(Bc=Ubc;Bc;Bc=Bc->next)
if(Bc->type == 'R' && Bc->elmt->id == E->id)
Add_robin_matrix(E, helm, Bc);
E->condense (helm,Ubsys); //OVERLOAD CALL: condense: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
}
if(Ubsys->smeth == direct || Ubsys->rslv)
E->project (helm,Ubsys); //OVERLOAD CALL: project: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
E->mat_free (helm); //OVERLOAD CALL: mat_free: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
if(lambda) E->mat_free (mass); //OVERLOAD CALL: mat_free: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
}
#ifdef PARALLEL
DO_PARALLEL {
if((!lambda || (!lambda->p && !lambda->d))
&& Ubsys->pll->nsolve == Ubsys->pll->nglobal)
Ubsys->singular = 1;
}
else
#endif
if((!lambda || (!lambda->p && !lambda->d))
&&(Ubsys->nsolve == Ubsys->nglobal)){
Ubsys->singular = Ubsys->bmap[iparam("IESING")][iparam("IVSING")]+1;
#ifdef DEBUG
// fprintf(stderr, "GenMat: matrix singular, removed node\n");
#endif
/* check to see if singular point is in inner boundary solve */ //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
if(Ubsys->rslv){
int i,k;
int cstart = Ubsys->rslv->rdata[Ubsys->rslv->nrecur-1].cstart;
if(Ubsys->singular > cstart){
for(k = iparam("IESING"); k < nel; ++k)
for(i = 0; i < U->flist[k]->Nverts; ++i)
if(Ubsys->bmap[k][i] < cstart){
Ubsys->singular = Ubsys->bmap[k][i] + 1;
goto finished;
}
for(k = 0; k < iparam("IESING"); ++k)
for(i = 0; i < U->flist[k]->Nverts; ++i)
if(Ubsys->bmap[k][i] < cstart){
Ubsys->singular = Ubsys->bmap[k][i] + 1;
goto finished;
}
finished:;
}
}
}
break;
}
if(Ubsys->rslv){ /* recurssive solver setup */
register int i,j;
double **a;
Recur *rdata;
#if 1
Condense_Recur(Ubsys,0,'p');
#else
Condense_Recur(Ubsys,0);
#endif
for(i = 0; i < Ubsys->rslv->nrecur-1; ++i){
rdata = Ubsys->rslv->rdata+i;
a = Ubsys->rslv->A.a;
#if 1
MemRecur(Ubsys,i+1,'p');
#else
MemRecur(Ubsys,i+1);
#endif
for(j = 0; j < rdata->npatch; ++j)
Project_Recur(j,rdata->patchlen_a[j],a[j],rdata->map[j],i+1,Ubsys);
#if 1
Condense_Recur(Ubsys,i+1,'p');
#else
Condense_Recur(Ubsys,i+1);
#endif
for(j = 0; j < rdata->npatch; ++j) free(a[j]); free(a);
}
#if 1
if(Ubsys->smeth == direct)
Rinvert_a(Ubsys,'p');
else
Rpack_a (Ubsys,'l');
#else
if(Ubsys->smeth == direct)
Rinvert_a(Ubsys);
else
Rpack_a (Ubsys);
#endif
#ifdef DEBUG
{
int memory = 0;
int nrecur = Ubsys->rslv->nrecur;
Recur *r = Ubsys->rslv->rdata;
int nsolve = r[nrecur-1].cstart;
for(i = 0 ; i < nrecur; ++i)
for(j = 0; j < r[i].npatch; ++j)
memory += r[i].patchlen_c[j]*(r[i].patchlen_c[j]+1)/2
+ r[i].patchlen_c[j]*r[i].patchlen_a[j];
if(Ubsys->smeth == direct){
int bwidth = Ubsys->rslv->Ainfo.bwidth_a;
if(2*bwidth < nsolve)
memory += bwidth*nsolve;
else
memory += nsolve*(nsolve+1)/2;
fprintf(stderr,"\nBandwidth a = %d [%d]\n",bwidth,nsolve);
}
else{
int mem_iter = 0;
for(i = 0; i < r[nrecur-1].npatch; ++i)
mem_iter +=r[nrecur-1].patchlen_a[i]*(r[nrecur-1].patchlen_a[i]+1)/2;
memory += mem_iter;
fprintf(stderr,"\n boundary system size: %d\n",mem_iter);
}
fprintf(stderr,"\nMemory = %d\n", memory);
}
#endif
}
else /* standard statically condensed set up */
if(Ubsys->smeth == direct) invert_a (Ubsys); //OVERLOAD CALL: invert_a: Matrix.c(?), Matrix_Stokes.c(?)
else{
FillPrecon(U->fhead, Ubsys); //OVERLOAD CALL: FillPrecon: Matrix.c(?), Matrix_Stokes.c(?)
InvtPrecon(Ubsys); //OVERLOAD CALL: InvtPrecon: Matrix.c(?), Matrix_Stokes.c(?)
}
/* zero U so don't start with funny values */
U->zerofield();
}
static void invert_a(Bsystem *Ubsys){
const int bwidth = Ubsys->Gmat->bwidth_a;
const int nsolve = Ubsys->nsolve;
int info=0;
if(nsolve){
if(2*bwidth < nsolve){ /* banded */
if(Ubsys->singular){
register int k;
int gid = Ubsys->singular-1;
dzero(bwidth,Ubsys->Gmat->inva[gid],1);
for(k = 0; k < min(bwidth,gid+1); ++k)
Ubsys->Gmat->inva[gid-k][k] = 0.0;
Ubsys->Gmat->inva[gid][0] = 1.0;
}
dpbtrf('L', nsolve, bwidth-1, *Ubsys->Gmat->inva, bwidth, info);
}
else{
if(Ubsys->singular){ /* symmetric */
register int k;
int gid = Ubsys->singular-1;
double *s;
s = *Ubsys->Gmat->inva+gid;
for(k = 0; k < gid; ++k, s += nsolve-k )
s[0] = 0.0;
dzero(nsolve-gid,s,1);
s[0] = 1.0;
}
dpptrf('L', nsolve, *Ubsys->Gmat->inva, info);
}
}
if(info) fprintf(stderr,"invertA: info not zero\n");
}
static void Bsystem_mem(Bsystem *Ubsys, Element *U){
const int nsolve = Ubsys->nsolve, nfam = Ubsys->families;
if(!Ubsys->rslv) /* if not a recursive solver */
if(Ubsys->smeth == iterative)
Ubsys->Gmat->a = (double **)calloc(nfam,sizeof(double *));
else{
if(nsolve){ /* banded solver */
Ubsys->Gmat->bwidth_a = bandwidth(U,Ubsys);
if(2*Ubsys->Gmat->bwidth_a < nsolve){
Ubsys->Gmat->inva = dmatrix(0, nsolve-1, 0, Ubsys->Gmat->bwidth_a-1);
dzero(nsolve*Ubsys->Gmat->bwidth_a,*Ubsys->Gmat->inva,1);
}
else{ /* symmetric solver */
Ubsys->Gmat->inva = dmatrix(0,0,0,nsolve*(nsolve+1)/2-1);
dzero(nsolve*(nsolve+1)/2,*Ubsys->Gmat->inva,1);
}
}
else
Ubsys->Gmat->inva = (double **)malloc(sizeof(double *));
}
Ubsys->Gmat->bwidth_c = ivector(0,nfam-1);
Ubsys->Gmat-> invc = (double **)calloc(nfam,sizeof(double *));
Ubsys->Gmat->binvc = (double **)calloc(nfam,sizeof(double *));
return;
}
/* Matrix preconditioning system */
void MemPrecon(Element *U, Bsystem *B){
register int i;
int gid,ll;
const int nvs = B->nv_solve, nes = B->ne_solve;
const int nfs = B->nf_solve;
Element *E;
MatPre *M = B->Pmat = (MatPre *) calloc(1,sizeof(MatPre)); //OVERLOAD CALL: MatPre: nekstruct.h(?), nekstruct.h(?); MatPre: nekstruct.h(?), nekstruct.h(?); MatPre: nekstruct.h(?), nekstruct.h(?)
switch(B->Precon){
case Pre_Diag:
M->info.diag.ndiag = B->nsolve;
M->info.diag.idiag = dvector(0,B->nsolve-1);
dzero(B->nsolve,M->info.diag.idiag,1);
break;
case Pre_Block:
M->info.block.nvert = nvs;
M->info.block.iedge = (double **) calloc((nes)?nes:1,sizeof(double *));
M->info.block.Ledge = ivector(0,(nes)?nes-1:0);
izero(nes, M->info.block.Ledge, 1);
/* set edges */
for(E=U;E;E=E->next)
for(i = 0; i < E->Nedges; ++i){
gid = E->edge[i].gid;
if(gid < nes)
M->info.block.Ledge[gid] = E->edge[i].l;
}
if(U->dim() == 3){ //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
M->info.block.iface = (double **) calloc((nfs)?nfs:1,sizeof(double *));
M->info.block.Lface = ivector(0,(nfs)?nfs-1:0);
izero(nfs, M->info.block.Lface, 1);
/* faces */
for(E=U;E;E=E->next)
for(i = 0; i < E->Nfaces; ++i){
gid = E->face[i].gid;
ll = E->face[i].l;
ll = (E->Nfverts(i) == 3) ? ll*(ll+1)/2 : ll*ll; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
if(gid < nfs)
M->info.block.Lface[gid] = ll;
}
}
/* declare memory */
M->info.block.ivert = dvector(0,nvs-1);
dzero(nvs,M->info.block.ivert,1);
for(i = 0; i < nes; ++i){
ll = M->info.block.Ledge[i];
ll = ll*(ll+1)/2;
M->info.block.iedge[i] = dvector(0,ll-1);
dzero(ll,M->info.block.iedge[i],1);
}
if(U->dim() == 3) //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
for(i = 0; i < nfs; ++i){
ll = M->info.block.Lface[i];
ll = ll*(ll+1)/2;
M->info.block.iface[i] = dvector(0,ll-1);
dzero(ll,M->info.block.iface[i],1);
}
break;
case Pre_Overlap:
M->info.overlap.ndiag = B->nsolve;
M->info.overlap.idiag = dvector(0,B->nsolve-1);
dzero(B->nsolve,M->info.overlap.idiag,1);
break;
case Pre_None:
break;
}
}
/* add A to preconditioner */
static void FillPrecon(Element *U, Bsystem *B){
register int i,j;
const int nvs = B->nv_solve;
int gid,pos,l,nbl;
MatPre *M = B->Pmat; //OVERLOAD CALL: MatPre: nekstruct.h(?), nekstruct.h(?)
double *A;
Element *E;
switch(B->Precon){
case Pre_Diag:
dzero(B->nsolve,M->info.diag.idiag,1); // RECENT FIX
break;
case Pre_Block:
dzero(nvs,M->info.block.ivert,1);
for(i = 0; i < B->ne_solve; ++i){
l = M->info.block.Ledge[i];
l = l*(l+1)/2;
dzero(l,M->info.block.iedge[i],1);
}
if(U->dim() == 3) //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
for(i = 0; i < B->nf_solve; ++i){
l = M->info.block.Lface[i];
l = l*(l+1)/2;
dzero(l,M->info.block.iface[i],1);
}
break;
case Pre_None:
break;
case Pre_Overlap:
dzero(B->nsolve,M->info.overlap.idiag,1);
break;
}
switch(B->Precon){
case Pre_Diag:
for(E=U;E;E=E->next){
A = B->Gmat->a[E->geom->id];
nbl = E->Nbmodes;
#if 0
for(i = 0,pos = 0; i < E->Nverts; ++i,pos+=nbl--){
gid = E->vert[i].gid;
if(gid < nvs)
M->info.diag.idiag[gid] += A[pos];
}
for(i = 0; i < E->Nedges; ++i){
gid = E->edge[i].gid;
l = E->edge[i].l;
if(gid < B->ne_solve){
for(j = 0; j < l; ++j,pos+=nbl--)
M->info.diag.idiag[B->edge[gid]+j] += A[pos];
}
else
for(j = 0; j < l; ++j,pos+=nbl--);
}
if(E->dim() == 3) //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
for(i = 0; i < E->Nfaces; ++i){
gid = E->face[i].gid;
l = E->face[i].l;
l = (E->Nfverts(i) == 3) ? l*(l+1)/2: l*l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
if(gid < B->nf_solve){
for(j = 0; j < l; ++j,pos+=nbl--){
M->info.diag.idiag[B->face[gid]+j] += A[pos];
}
}
else
for(j = 0; j < l; ++j,pos+=nbl--);
}
#endif
pos = 0;
for(i=0;i<nbl;pos += nbl-i, ++i)
if(B->bmap[E->id][i] < B->nsolve)
M->info.diag.idiag[B->bmap[E->id][i]] += A[pos];
}
#ifdef PARALLEL
DO_PARALLEL
parallel_gather(M->info.diag.idiag,B);
#endif
break;
case Pre_Block:
fprintf(stderr,"Block Precondition is not set up in H_Matrix.C\n");
exit(-1);
break;
case Pre_None:
break;
case Pre_Overlap:
#if 0
for(E=U;E;E=E->next){
A = B->Gmat->a[E->geom->id];
nbl = E->Nbmodes;
for(i = 0,pos = 0; i < E->Nverts; ++i,pos+=nbl--){
gid = E->vert[i].gid;
if(gid < nvs)
M->info.overlap.idiag[gid] += A[pos];
}
for(i = 0; i < E->Nedges; ++i){
gid = E->edge[i].gid;
l = E->edge[i].l;
if(gid < B->ne_solve){
for(j = 0; j < l; ++j,pos+=nbl--)
M->info.overlap.idiag[B->edge[gid]+j] += A[pos];
}
else
for(j = 0; j < l; ++j,pos+=nbl--);
}
if(E->dim() == 3) //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
for(i = 0; i < E->Nfaces; ++i){
gid = E->face[i].gid;
l = E->face[i].l;
l = (E->Nfverts(i) == 3) ? l*(l+1)/2: l*l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
if(gid < B->nf_solve){
for(j = 0; j < l; ++j,pos+=nbl--){
M->info.overlap.idiag[B->face[gid]+j] += A[pos];
}
}
else
for(j = 0; j < l; ++j,pos+=nbl--);
}
}
#endif
for(E=U;E;E=E->next){
A = B->Gmat->a[E->geom->id];
nbl = E->Nbmodes;
pos = 0;
for(i=0;i<nbl;pos += nbl-i, ++i)
if(B->bmap[E->id][i] < B->nsolve)
M->info.diag.idiag[B->bmap[E->id][i]] += A[pos];
}
#ifdef PARALLEL
DO_PARALLEL
parallel_gather(M->info.diag.idiag,B);
#endif
break;
}
}
static void InvtPrecon(Bsystem *B){
switch(B->Precon){
case Pre_Diag:
dvrecp(B->Pmat->info.diag.ndiag,B->Pmat->info.diag.idiag,
1,B->Pmat->info.diag.idiag,1);
break;
case Pre_Block:
break;
case Pre_Overlap:
dvrecp(B->Pmat->info.overlap.ndiag,B->Pmat->info.overlap.idiag,
1,B->Pmat->info.overlap.idiag,1);
break;
case Pre_None:
break;
}
}
/* pack 'a' matrix into 'b' */
void PackMatrix(double **a, int n, double *b, int bwidth){
register int i;
if(n>2*bwidth){ /* banded symmetric lower triangular form */
double *s;
for(i = 0,s=b; i < n-bwidth; ++i,s+=bwidth)
dcopy(bwidth,a[i]+i,1,s,1);
for(i = n-bwidth; i < n; ++i,s+=bwidth)
dcopy(n-i,a[i]+i,1,s,1);
}
else{ /* symmetric lower triangular form */
register int j;
for(i = 0,j=0; i < n; j+=n-i++)
dcopy(n-i, a[i]+i, 1, b+j, 1);
}
}
/* factor positive symmertric and banded matrices */
void FacMatrix(double *a, int n, int bwidth){
int info=0;
if(n>2*bwidth)
dpbtrf('L',n,bwidth-1,a,bwidth,info);
else
dpptrf('L',n,a,info);
if(info) {
fprintf(stderr, "FacMatrix: n=%d bw=%d info=%d\n",
n,bwidth,info);
error_msg(FacMatrixP: info not zero); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
}
}
#if (defined(DUMPMASS) || defined(DUMPHELM))
static void dumpmat(Element *E,LocMat *mat){
register int i,j;
int asize;
int csize;
csize = E->Nbmodes;
asize = E->Nmodes - csize;
for(i=0; i< asize; ++i)
for(j = 0; j < asize; ++j)
if(fabs(mat->a[i][j]) > 1e-10)
fprintf(stdout,"%d %d %lg \n",i,j,mat->a[i][j]);
for(i=0; i< csize; ++i)
for(j = 0; j < csize; ++j)
if(fabs(mat->c[i][j]) > 1e-10)
fprintf(stdout,"%d %d %lg \n",i+asize,j+asize,mat->c[i][j]);
for(i = 0; i < asize; ++i)
for(j = 0; j < csize; ++j)
if(fabs(mat->b[i][j])> 1e-10){
fprintf(stdout,"%d %d %lg \n",i,j+asize,mat->b[i][j]);
fprintf(stdout,"%d %d %lg \n",j+asize,i,mat->b[i][j]);
}
}
#endif
int bandwidth(Element *E, Bsystem *Bsys){
register int i;
const int nes = Bsys->ne_solve;
int high,low,bwidth=0;
int l;
const int nfs = Bsys->nf_solve;
for(;E; E = E->next){
low = 1000000; high = 0;
for(i = 0; i < E->Nverts; ++i)
if(E->vert[i].solve){
low = min(low, E->vert[i].gid);
high = max(high,E->vert[i].gid);
}
for(i = 0; i < E->Nedges; ++i)
if(E->edge[i].gid < nes){
low = min(low, Bsys->edge[E->edge[i].gid]);
high = max(high,Bsys->edge[E->edge[i].gid] + E->edge[i].l);
}
// needs to be fixed
if(E->dim() == 3) //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
for(i = 0; i < E->Nfaces; ++i)
if(E->face[i].gid < nfs){
l = E->face[i].l;
l = (E->Nfverts(i) == 3) ? l*(l+1)/2 : l*l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
low = min(low, Bsys->face[E->face[i].gid]);
high = max(high,Bsys->face[E->face[i].gid] + l);
}
bwidth = ((high-low+1) > bwidth)? (high-low+1) : bwidth;
}
return bwidth;
}
void Bsystem_mem_free(Bsystem *Ubsys, Element_List *U){
const int nsolve = Ubsys->nsolve, nfam = Ubsys->families;
int i,j;
if(!Ubsys->rslv) /* if not a recursive solver */
if(Ubsys->smeth == iterative){
for(i=0;i<nfam;++i){
free(Ubsys->Gmat->a[i]);
Ubsys->Gmat->a[i] = (double*)0;
}
free(Ubsys->Gmat->a);
Ubsys->Gmat->a = (double**)0;
}
else{
if(nsolve){ /* banded solver */
Ubsys->Gmat->bwidth_a = 0;
free_dmatrix(Ubsys->Gmat->inva,0,0);
Ubsys->Gmat->inva = (double**)0;
}
else{
free(Ubsys->Gmat->inva); // MIGHT NEED TO FIX
Ubsys->Gmat->inva = (double**)0;
}
}
if(nfam>1){
free(Ubsys->Gmat->bwidth_c);
Ubsys->Gmat->bwidth_c = (int*)0;
}
for(i=0;i<nfam;++i){
free(Ubsys->Gmat->invc[i]);
free(Ubsys->Gmat->binvc[i]);
Ubsys->Gmat->invc[i] = (double*)0;
Ubsys->Gmat->binvc[i] = (double*)0;
}
free(Ubsys->Gmat->invc);
free(Ubsys->Gmat->binvc);
Ubsys->Gmat->invc = (double**)0;
Ubsys->Gmat->binvc = (double**)0;
if (Ubsys->smeth == iterative) {
#if 1
switch(Ubsys->Precon){
case Pre_Diag:
free(Ubsys->Pmat->info.diag.idiag);
free(Ubsys->Pmat);
break;
case Pre_Block:
break;
case Pre_None:
break;
case Pre_Overlap:
free(Ubsys->Pmat->info.overlap.idiag);
free(Ubsys->Pmat);
}
#else
register int i;
int l,gid;
const int nvs = Ubsys->nv_solve, nes = Ubsys->ne_solve;
Element *E;
free(Ubsys->Pmat->ivert);
Ubsys->Pmat->ivert = (double*)0;
free(Ubsys->Pmat->iedge);
Ubsys->Pmat->iedge = (double**)0;
if(U->fhead->dim() == 3){ //OVERLOAD CALL: dim: Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), Hex.h(?), Prism.h(?)
free(Ubsys->Pmat->Lface);
free(Ubsys->Pmat->iface);
Ubsys->Pmat->iface = (double **)0;
Ubsys->Pmat->Lface = (int*)0;
}
free(Ubsys->Pmat->Ledge); Ubsys->Pmat->Ledge = (int*)0;
#endif
}
if (Ubsys->rslv) {
Rsolver *rslv = Ubsys->rslv;
Recur *rdata;
for(j=0;j<rslv->nrecur;++j){
rdata = rslv->rdata + j;
for(i = 0; i < rdata->npatch; ++i){
if(rdata->patchlen_a[i])
free(rdata->binvc[i]);
free(rdata->invc[i]);
}
free(rdata->binvc); rdata->binvc=NULL;
free(rdata->invc); rdata->invc =NULL;
free(rdata->bwidth_c); rdata->bwidth_c=NULL;
}
}
return;
}
static void Add_robin_matrix(Element *E, LocMat *helm, Bndry *rbc){
int edge = rbc->face;
int qa = E->qa;
double *tmp,*z,*w, d;
int i,j;
Basis *b = E->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)
double sigma = dparam("SIGMA");
int start, vn1, vn2, eni, enj;
for(i=0, start=E->Nverts;i<edge;++i)
start += E->edge[i].l;
getzw(qa,&z,&w,'a');
tmp = dvector(0,qa-1);
/* Inner product of vertex modes (0,0)*/
vn1 = E->vnum(edge,0); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
dvmul(qa,b->vert[0].a,1,w,1,tmp,1); // 0->i
if(E->curve)
dvmul(qa,rbc->sjac.p,1,tmp,1,tmp,1);
else
dsmul(qa,rbc->sjac.d,tmp,1,tmp,1);
d = ddot(qa,tmp,1,b->vert[0].a,1);
helm->a[vn1][vn1] += sigma*d;
/* Inner product of vertex modes (1,1)*/
vn2 = E->vnum(edge,1); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
dvmul(qa,b->vert[1].a,1,w,1,tmp,1);
if(E->curve)
dvmul(qa,rbc->sjac.p,1,tmp,1,tmp,1);
else
dsmul(qa,rbc->sjac.d,tmp,1,tmp,1);
d = ddot(qa,tmp,1,b->vert[1].a,1);
helm->a[vn2][vn2] += sigma*d;
/* Inner product of vertex modes (0,1)*/
dvmul(qa,b->vert[0].a,1,w,1,tmp,1);
if(E->curve)
dvmul(qa,rbc->sjac.p,1,tmp,1,tmp,1);
else
dsmul(qa,rbc->sjac.d,tmp,1,tmp,1);
d = ddot(qa,tmp,1,b->vert[1].a,1);
helm->a[vn1][vn2] += sigma*d;
helm->a[vn2][vn1] += sigma*d;
/* Inner product of vertex0 with edge modes */
dvmul(qa,b->vert[0].a,1,w,1,tmp,1);
if(E->curve)
dvmul(qa,rbc->sjac.p,1,tmp,1,tmp,1);
else
dsmul(qa,rbc->sjac.d,tmp,1,tmp,1);
for(j = 0; j < E->edge[edge].l; ++j){
enj = start+j;
d = ddot(qa,tmp,1,b->edge[0][j].a,1);
helm->a[enj][vn1] += sigma*d;
helm->a[vn1][enj] += sigma*d;
}
/* Inner product of vertex1 with edge modes */
dvmul(qa,b->vert[1].a,1,w,1,tmp,1);
if(E->curve)
dvmul(qa,rbc->sjac.p,1,tmp,1,tmp,1);
else
dsmul(qa,rbc->sjac.d,tmp,1,tmp,1);
for(j = 0; j < E->edge[edge].l; ++j){
enj = start+j;
d = ddot(qa,tmp,1,b->edge[0][j].a,1);
helm->a[enj][vn2] += sigma*d;
helm->a[vn2][enj] += sigma*d;
}
/* Inner product of edge modes with edge modes */
for(i = 0; i < E->edge[edge].l; ++i){
dvmul(qa,b->edge[0][i].a,1,w,1,tmp,1);
if(E->curve)
dvmul(qa,rbc->sjac.p,1,tmp,1,tmp,1);
else
dsmul(qa,rbc->sjac.d,tmp,1,tmp,1);
eni = start+i;
for(j = i; j < E->edge[edge].l; ++j){
enj = start+j;
d = ddot(qa,tmp,1,b->edge[0][j].a,1);
helm->a[eni][enj] += sigma*d;
if(eni!=enj)
helm->a[enj][eni] += sigma*d;
}
}
// plotdmatrix(helm->a, 0, helm->asize-1, 0, helm->asize-1);
free(tmp);
}
void CorrectRobinFluxes(Bndry *B){
register int i;
const int face = B->face;
const int qa = B->elmt->qa, l = B->elmt->edge[face].l;
double *wa,*tmp,*fi, fac;
Mode *m;
Element *E = B->elmt;
Basis *Base = E->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)
int gid, vid, vn1, vn2;
Bndry *bc, *Ubc;
fi = dvector(0,qa-1);
tmp = dvector(0,qa-1);
for(Ubc=B;Ubc;Ubc=Ubc->next)
if(Ubc->type=='R'){
E = B->elmt;
vn1 = E->vnum(B->face,0); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
vn2 = E->vnum(B->face,1); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
if(E->vert[vn1].solve && E->vert[vn2].solve) break;
gid = E->vert[vn1].gid;
getzw(qa,&wa,&wa,'a');
/* Add shared, Dirichlet, vertex modes */
if(!E->vert[vn1].solve)
for(bc = Ubc; bc; bc=bc->next){
E = bc->elmt;
for(i=0;i<E->Nfverts(bc->face);++i){ //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
vid = E->vnum(bc->face,i); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
if(gid == E->vert[vid].gid &&
(bc->type == 'V' || bc->type == 'W')){
fac = -bc->bvert[i]*dparam("SIGMA");
dsmul(qa,fac,Base->vert[0].a,1,fi,1);
break;
}
}
}
gid = E->vert[vn2].gid;
if(!E->vert[vn2].solve)
for(bc = Ubc; bc; bc=bc->next){
E = bc->elmt;
for(i=0;i<E->Nfverts(bc->face);++i){ //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
vid = E->vnum(bc->face,i); //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
if(gid == E->vert[vid].gid &&
(bc->type == 'V' || bc->type == 'W')){
fac = -bc->bvert[i]*dparam("SIGMA");
dsvtvp(qa, fac, Base->vert[1].a, 1, fi, 1, fi, 1);
break;
}
}
}
tmp = dvector(0,qa-1);
m = Base->edge[0];
if(B->elmt->curve)
dvmul(qa,B->sjac.p,1,fi,1,fi,1);
else
dsmul(qa,B->sjac.d,fi,1,fi,1);
/* calculate inner product over surface */
dvmul(qa,Base->vert[0].a,1,wa,1,tmp,1);
B->bvert[0] += ddot(qa,tmp,1,fi,1);
dvmul(qa,Base->vert[1].a,1,wa,1,tmp,1);
B->bvert[1] += ddot(qa,tmp,1,fi,1);
for(i = 0; i < l; ++i){
dvmul(qa,m[i].a,1,wa,1,tmp,1);
B->bedge[0][i] += ddot(qa,tmp,1,fi,1);
}
}
free(tmp); free(fi);
}
void Setup_Preconditioner(Element_List *U, Bndry *Ubc, Bsystem *Ubsys, Metric *lambda, SolveType Stype){
switch(Ubsys->Precon){
case Pre_Diag:
if(Ubsys->Pmat->info.diag.idiag)
free(Ubsys->Pmat->info.diag.idiag);
break;
}
if (Ubsys->smeth == iterative) MemPrecon(U->fhead,Ubsys); //OVERLOAD CALL: MemPrecon: Matrix.c(?), Matrix_Stokes.c(?)
FillPrecon(U->fhead, Ubsys); //OVERLOAD CALL: FillPrecon: Matrix.c(?), Matrix_Stokes.c(?)
InvtPrecon(Ubsys); //OVERLOAD CALL: InvtPrecon: Matrix.c(?), Matrix_Stokes.c(?)
}
/* factor symmertric and banded matrices which are not nec. positive definite*/
void FacMatrixSym(double *a, int *pivot, int n, int bwidth){
int info;
if(n>2*bwidth){
error_msg(not implemented banded non-positive definite factorisation); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
dpbtrf('L',n,bwidth-1,a,bwidth,info);
}
else
dsptrf('L',n,a,pivot,info);
if(info) {
fprintf(stderr, "FacMatrixSym: n=%d bw=%d info=%d\n",
n,bwidth,info);
error_msg(FacMatrixP: info not zero); //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
}
}
// Helmholtz only
void UpdateMatList(Element_List *U, Bndry *Ubc, Bsystem *Ubsys,
Metric *lambda, SolveType Stype,
int *list, int init, int id){
LocMat *mass = 0,*helm, *oldhelm;
int na,nb,nc;
Element *E;
static double ****alist;
double *dsave = dvector(0, LGmax*LGmax*LGmax-1);
// If init is non-zero -> setup storage for a-b^t c^-1 b for each elmt
// also setup general matrix storage
if(init && Ubsys->smeth == direct){
if(!alist)
alist = (double****) calloc(MAXFIELDS, sizeof(double***));
if(!alist[id])
alist[id] = (double***) calloc(U->nel, sizeof(double**));
Bsystem_mem(Ubsys,U->fhead);
if (Ubsys->smeth == iterative)
MemPrecon(U->fhead,Ubsys); //OVERLOAD CALL: MemPrecon: Matrix.c(?), Matrix_Stokes.c(?)
}
// setup dummy storage for project routine for unchanging elmts
oldhelm = (LocMat*) calloc(1, sizeof(LocMat));
for(E = U->fhead; E;E=E->next){
if(list[E->id] || init){
// save field values
dcopy(E->Nmodes, E->vert[0].hj, 1, dsave, 1);
// setup local storage for matrices
helm = E->mat_mem (); //OVERLOAD CALL: mat_mem: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
if(lambda){
mass = E->mat_mem (); //OVERLOAD CALL: mat_mem: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
na = mass->asize*mass->asize;
nb = mass->asize*mass->csize;
nc = mass->csize*mass->csize;
}
// generate elmtal matrices
if(E->curvX || lambda[E->id].p)
E->HelmMatC (helm,lambda+E->id); //OVERLOAD CALL: HelmMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
else{
if(E->identify()==Nek_Tri) //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
Tri_HelmMat(E, helm, lambda[E->id].d);
else
{
E->LapMat(helm); //OVERLOAD CALL: LapMat: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex)
if(lambda[E->id].d){
E->MassMat(mass); //OVERLOAD CALL: MassMat: Matrix.c(Tri), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex)
daxpy(na,lambda[E->id].d,*mass->a,1,*helm->a,1);
if(nc){
daxpy(nb,lambda[E->id].d,*mass->b,1,*helm->b,1);
daxpy(nc,lambda[E->id].d,*mass->c,1,*helm->c,1);
}
}
}
}
// calculate a-b^t c^-1 b
E->condense (helm,Ubsys); //OVERLOAD CALL: condense: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
// store result
if(init && Ubsys->smeth == direct){
alist[id][E->id] = dmatrix(0, helm->asize-1, 0, helm->asize-1);
dcopy(helm->asize*helm->asize, helm->a[0], 1,
alist[id][E->id][0], 1);
}
// free mass matrix
if(lambda) E->mat_free (mass); //OVERLOAD CALL: mat_free: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
// project local bndry matrix to global matrix
if(Ubsys->smeth == direct || Ubsys->rslv)
E->project (helm,Ubsys); //OVERLOAD CALL: project: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
// free remaining temporary local matrices
E->mat_free (helm); //OVERLOAD CALL: mat_free: Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet)
// restore field values
dcopy(E->Nmodes, dsave, 1, E->vert[0].hj, 1);
}
else{
// project stored a-b^t c^-1 b to global matrix
if(Ubsys->smeth == direct || Ubsys->rslv){
oldhelm->a = alist[id][E->id];
E->project (oldhelm,Ubsys); //OVERLOAD CALL: project: Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element), Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet)
}
}
}
#ifdef PARALLEL
DO_PARALLEL {
if((!lambda || (!lambda->p && !lambda->d))
&& Ubsys->pll->nsolve == Ubsys->pll->nglobal)
Ubsys->singular = 1;
}
else
#endif
if((!lambda || (!lambda->p && !lambda->d))
&&(Ubsys->nsolve == Ubsys->nglobal)){
Ubsys->singular = Ubsys->bmap[iparam("IESING")][iparam("IVSING")]+1;
}
if(Ubsys->smeth == direct)
invert_a (Ubsys); //OVERLOAD CALL: invert_a: Matrix.c(?), Matrix_Stokes.c(?)
else{
FillPrecon(U->fhead, Ubsys); //OVERLOAD CALL: FillPrecon: Matrix.c(?), Matrix_Stokes.c(?)
InvtPrecon(Ubsys); //OVERLOAD CALL: InvtPrecon: Matrix.c(?), Matrix_Stokes.c(?)
}
free(dsave);
free(oldhelm);
}
/*
Function name: Element::MassMat //OVERLOAD CALL: MassMat: Matrix.c(?), Matrix.c(?)
Function Purpose:
Argument 1: LocMat *mass
Purpose:
Function Notes:
*/
void Tri::MassMat(LocMat *mass){
register int i,i1,j,k,n,n1;
int H;
int ll;
Basis *B;
Mode *w,*mtmp,*m,*m1;
double *z;
double jac = geom->jac.d;
double **a = mass->a,
**b = mass->b,
**c = mass->c;
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)
w = mvector(0,0);
mtmp = mvecset(0,0,qa,qb,qc);
getzw(qa,&z,&w[0].a,'a');
getzw(qb,&z,&w[0].b,'b');
/* sort A and B matrices */
m = B->vert;
for(i = 0,n=0; i < Nverts; ++i,++n){
Tri_mvmul2d(qa,qb,qc,m+i,w,mtmp);
/* vertex with vertex */
for(j = i,n1=n; j < Nverts; ++j,++n1)
a[n1][n] = a[n][n1] = jac*iprod(m+j,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
/* vertex with edge */
for(j = 0; j < Nedges; ++j){
m1 = B->edge[j];
for(k = 0; k < edge[j].l; ++k,++n1)
a[n1][n] = a[n][n1] = jac*iprod(m1+k,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
}
/* vertices with interior */
ll = face->l;
for(k = 0, n1=0; k < ll; ++k){
m1 = B->face[0][k];
for(H = 0; H < ll-k; ++H,++n1)
b[n][n1] = jac*iprod(m1+H,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
}
}
for(i = 0; i < Nedges; ++i){
m = B->edge[i];
for(i1 = 0; i1 < edge[i].l; ++i1,++n){
Tri_mvmul2d(qa,qb,qc,m+i1,w,mtmp);
/* edge with edge */
for(k = i1,n1=n; k < edge[i].l; ++k,++n1)
a[n1][n] = a[n][n1] = jac*iprod(m+k,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
for(j = i+1; j < Nedges; ++j){
m1 = B->edge[j];
for(k = 0; k < edge[j].l; ++k,++n1)
a[n1][n] = a[n][n1] = jac*iprod(m1+k,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
}
ll = face[0].l;
for(k = 0,n1=0; k < ll; ++k){
m1 = B->face[0][k];
for(H = 0; H < ll-k; ++H,++n1)
b[n][n1] = jac*iprod(m1+H,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
}
}
}
/* fill c */
ll = face[0].l;
for(i = 0,n=0; i < ll; ++i){
m = B->face[0][i];
for(i1=0; i1 < ll-i;++i1,++n){
Tri_mvmul2d(qa,qb,qc,m+i1,w,mtmp);
for(H = i1,n1=n; H < ll-i; ++H,++n1)
c[n1][n] = c[n][n1] = jac*iprod(m+H,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
for(k = i+1; k < ll; ++k){
m1 = B->face[0][k];
for(H = 0; H < ll-k; ++H,++n1)
c[n1][n] = c[n][n1] = jac*iprod(m1+H,mtmp); //OVERLOAD CALL: iprod: InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism)
}
}
}
// recent change
free((char*)w);
free_mvec(mtmp);
}
void Quad::MassMat (LocMat *mass){
MassMatC(mass); //OVERLOAD CALL: MassMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
}
void Tet::MassMat(LocMat *mass){
MassMatC(mass); //OVERLOAD CALL: MassMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
}
void Pyr::MassMat(LocMat *mass){
MassMatC(mass); //OVERLOAD CALL: MassMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
}
void Prism::MassMat(LocMat *mass){
MassMatC(mass); //OVERLOAD CALL: MassMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
}
void Hex::MassMat(LocMat *mass){
MassMatC(mass); //OVERLOAD CALL: MassMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
}
void Element::MassMat (LocMat *){ERR;} // return mass-matri(ces) //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
/*
Function name: Element::MassMatC
Function Purpose:
Argument 1: LocMat *mass
Purpose:
Function Notes:
*/
void Tri::MassMatC(LocMat *mass){
register int i,j,n;
const int nbl = Nbmodes, N = Nmodes - nbl, csize = mass->csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
/* fill A */
for(i = 0, n = 0; i < Nverts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
for(i = 0; i < Nedges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
/* fill B and C */
L = face->l;
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L-i; ++j, ++n){
fillElmt(B->face[0][i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N ,vert->hj+nbl,1,c[n],1);
}
}
void Quad::MassMatC(LocMat *mass){
double *save = dvector(0, qtot+Nmodes-1);
dcopy(qtot, h[0], 1, save, 1);
dcopy(Nmodes, vert[0].hj, 1, save+qtot, 1);
register int i,j,n;
const int nbl = Nbmodes, N = Nmodes - nbl, csize = mass->csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
/* fill A */
for(i = 0, n = 0; i < Nverts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
for(i = 0; i < Nedges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
/* fill B and C */
L = face->l;
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L; ++j, ++n){
fillElmt(B->face[0][i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N ,vert->hj+nbl,1,c[n],1);
}
dcopy(qtot, save, 1, h[0], 1);
dcopy(Nmodes, save+qtot, 1, vert[0].hj, 1);
free(save);
// dump_sc(mass->asize, mass->csize, mass->a, mass->b, mass->c);
}
void Tet::MassMatC(LocMat *mass){
register int i,j,k,n;
const int nbl = Nbmodes, N = Nmodes - nbl, csize = mass->csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
/* fill A */
for(i = 0, n = 0; i < Nverts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
for(i = 0; i < Nedges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
for(i = 0; i < Nfaces; ++i)
for(j = 0; j < (L = face[i].l); ++j)
for(k = 0; k < L-j; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,a[n],1);
}
/* fill B and C */
L = interior_l;
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L-i; ++j)
for(k = 0; k < L-i-j; ++k, ++n){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N, vert->hj+nbl,1,c[n],1);
}
// dump_sc(nbl, csize, a, b, c);
}
void Pyr::MassMatC(LocMat *mass){
register int i,j,k,n;
int kk;
int asize, csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
asize = Nbmodes;
csize = Nmodes - Nbmodes;
/* fill A */
for(i = 0, n = 0; i < NPyr_verts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPyr_edges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPyr_faces; ++i)
for(j = 0; j < (L = face[i].l); ++j){
kk = (Nfverts(i) == 3) ? face[i].l-j:face[i].l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(k = 0; k < kk; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
}
/* fill B and C */
L = interior_l;
if(csize)
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L-i; ++j)
for(k = 0; k < L-i-j; ++k){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,*b+n,csize);
dcopy(csize,vert->hj+asize,1,c[n],1);
++n;
}
// dump_sc(mass->asize, mass->csize, mass->a, mass->b, mass->c);
}
void Prism::MassMatC(LocMat *mass){
register int i,j,k,n;
int kk;
int asize, csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
asize = Nbmodes;
csize = Nmodes - Nbmodes;
/* fill A */
for(i = 0, n = 0; i < NPrism_verts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPrism_edges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPrism_faces; ++i)
for(j = 0; j < (L = face[i].l); ++j){
kk = (Nfverts(i) == 3) ? face[i].l-j:face[i].l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(k = 0; k < kk; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
}
/* fill B and C */
L = interior_l;
for(i = 0,n=0; i < L-1; ++i)
for(j = 0; j < L; ++j)
for(k = 0; k < L-1-i; ++k, ++n){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,*b+n,csize);
dcopy(csize,vert->hj+asize,1,c[n],1);
}
}
void Hex::MassMatC(LocMat *mass){
register int i,j,k,n;
int asize, csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
asize = Nbmodes;
csize = Nmodes - Nbmodes;
/* fill A */
for(i = 0, n = 0; i < NHex_verts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NHex_edges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NHex_faces; ++i)
for(j = 0; j < (L = face[i].l); ++j)
for(k = 0; k < L; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,a[n],1);
}
/* fill B and C */
L = interior_l;
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L; ++j)
for(k = 0; k < L; ++k, ++n){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
Iprod(this); //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
dcopy(asize,vert->hj ,1,*b+n,csize);
dcopy(csize,vert->hj+asize,1,c[n],1);
}
}
void Element::MassMatC(LocMat *){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
/*
Function name: Element::HelmMatC
Function Purpose:
Argument 1: LocMat *mass
Purpose:
Argument 2: Metric *lambda
Purpose:
Function Notes:
*/
void Tri::HelmMatC(LocMat *mass, Metric *lambda){
register int i,j,n;
int nbl = Nbmodes, N = Nmodes - nbl, asize = mass->asize,csize = mass->csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c,
**d = mass->d;
/* fill A */
for(i = 0,n=0; i < Nverts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,*a+n,asize);
if(lambda->wave)
dcopy(N,vert->hj+nbl,1,d[n],1);
}
for(i = 0; i < Nedges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,*a+n,asize);
if(lambda->wave)
dcopy(N,vert->hj+nbl,1,d[n],1);
}
/* fill B and C */
L = face->l;
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L-i; ++j, ++n){
fillElmt(B->face[0][i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N ,vert->hj+nbl,1,*c+n,csize);
}
}
void Quad::HelmMatC(LocMat *mass, Metric *lambda){
double *save = dvector(0, qtot+Nmodes-1);
dcopy(qtot, h[0], 1, save, 1);
dcopy(Nmodes, vert[0].hj, 1, save+qtot, 1);
register int i,j,n;
int nbl = Nbmodes, N = Nmodes - nbl, csize = mass->csize,asize = mass->asize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c,
**d = mass->d;
char orig_state = state;
/* fill A */
for(i = 0,n=0; i < Nverts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,*a+n,asize);
if(lambda->wave)
dcopy(N,vert->hj+nbl,1,d[n],1);
}
for(i = 0; i < Nedges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,*a+n,asize);
if(lambda->wave)
dcopy(N,vert->hj+nbl,1,d[n],1);
}
/* fill B and C */
L = face->l;
for(i = 0,n=0; i < L;++i)
for(j = 0; j < L; ++j,++n){
fillElmt(B->face[0][i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N ,vert->hj+nbl,1,*c+n,csize);
}
// dump_sc(mass->asize, mass->csize, mass->a, mass->b, mass->c);
state = orig_state;
dcopy(qtot, save, 1, h[0], 1);
dcopy(Nmodes, save+qtot, 1, vert[0].hj, 1);
free(save);
}
void Tet::HelmMatC(LocMat *mass, Metric *lambda){
register int i,j,k,n;
int nbl = Nbmodes, N = Nmodes - nbl, csize = mass->csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
/* fill A */
for(i = 0,n=0; i < Nverts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,a[n],1);
}
for(i = 0; i < Nedges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,a[n],1);
}
for(i = 0; i < Nfaces; ++i)
for(j = 0; j < (L = face[i].l); ++j)
for(k = 0; k < L-j; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,a[n],1);
}
/* fill B and C */
L = interior_l;
for(i = 0,n = 0; i < L; ++i)
for(j = 0; j < L-i; ++j)
for(k = 0; k < L-i-j; ++k,++n){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N ,vert->hj+nbl,1,c[n],1);
}
// dump_sc(Nbmodes, Nmodes-Nbmodes, a, b, c);
}
void Pyr::HelmMatC(LocMat *mass, Metric *lambda){
register int i,j,k,n;
int kk;
int asize, csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
asize = Nbmodes;
csize = Nmodes - Nbmodes;
/* fill A */
for(i = 0, n = 0; i < NPyr_verts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPyr_edges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPyr_faces; ++i)
for(j = 0; j < (L = face[i].l); ++j){
kk = (Nfverts(i) == 3) ? face[i].l-j:face[i].l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(k = 0; k < kk; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,a[n],1);
}
}
/* fill B and C */
L = interior_l;
if(csize)
for(i = 0,n=0; i < L; ++i)
for(j = 0; j < L-i; ++j)
for(k = 0; k < L-i-j; ++k){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,*b+n,csize);
dcopy(csize,vert->hj+asize,1,c[n],1);
++n;
}
// dump_sc(mass->asize, mass->csize, mass->a, mass->b, mass->c);
}
void Prism::HelmMatC(LocMat *mass, Metric *lambda){
register int i,j,k,n;
int kk;
int asize, csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
asize = Nbmodes;
csize = Nmodes - Nbmodes;
/* fill A */
for(i = 0, n = 0; i < NPrism_verts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPrism_edges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,a[n],1);
}
for(i = 0; i < NPrism_faces; ++i)
for(j = 0; j < (L = face[i].l); ++j){
kk = (Nfverts(i) == 3) ? face[i].l-j:face[i].l; //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(k = 0; k < kk; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,a[n],1);
}
}
/* fill B and C */
L = interior_l;
for(i = 0,n=0; i < L-1; ++i)
for(j = 0; j < L; ++j)
for(k = 0; k < L-1-i; ++k, ++n){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(asize,vert->hj ,1,*b+n,csize);
dcopy(csize,vert->hj+asize,1,c[n],1);
}
}
void Hex::HelmMatC(LocMat *mass, Metric *lambda){
register int i,j,k,n;
int nbl = Nbmodes, N = Nmodes - nbl, csize = mass->csize;
int L;
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)
double **a = mass->a,
**b = mass->b,
**c = mass->c;
/* fill A */
for(i = 0,n=0; i < NHex_verts; ++i,++n){
fillElmt(B->vert+i); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,a[n],1);
}
for(i = 0; i < NHex_edges; ++i)
for(j = 0; j < (L = edge[i].l); ++j, ++n){
fillElmt(B->edge[i]+j); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,a[n],1);
}
for(i = 0; i < NHex_faces; ++i)
for(j = 0; j < (L = face[i].l); ++j)
for(k = 0; k < L; ++k, ++n){
fillElmt(B->face[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj,1,a[n],1);
}
/* fill B and C */
L = interior_l;
for(i = 0,n = 0; i < L; ++i)
for(j = 0; j < L; ++j)
for(k = 0; k < L; ++k,++n){
fillElmt(B->intr[i][j]+k); //OVERLOAD CALL: fillElmt: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
HelmHoltz(lambda); //OVERLOAD CALL: HelmHoltz: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri)
dcopy(nbl,vert->hj ,1,*b+n,csize);
dcopy(N ,vert->hj+nbl,1,c[n],1);
}
// dump_sc(Nbmodes, Nmodes-Nbmodes, a, b, c);
}
void Element::HelmMatC(LocMat *, Metric *){ERR;} // return Helmholtz op. //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
/*
Function name: Element::LapMat
Function Purpose:
Argument 1: LocMat *lap
Purpose:
Function Notes:
*/
void Tri::LapMat(LocMat *lap){
register int i,j,k,n;
int ll,nbl;
Basis *B,*DB;
Mode *w,*m,*m1,*md,*md1,**gb,**gb1,*fac;
double *z,jac = geom->jac.d;
double **a = lap->a,
**b = lap->b,
**c = lap->c;
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)
fac = B->vert;
w = mvector(0,0);
gb = (Mode **) calloc(Nmodes,sizeof(Mode *));
gb[0] = mvecset(0,3*Nmodes,qa,qb,qc);
gb1 = (Mode **) calloc(Nmodes,sizeof(Mode *));
gb1[0] = mvecset(0,3*Nmodes,qa,qb,qc);
for(i = 1; i < Nmodes; ++i) gb[i] = gb[i-1]+3;
for(i = 1; i < Nmodes; ++i) gb1[i] = gb1[i-1]+3;
nbl = Nmodes - face[0].l*(face[0].l+1)/2;
getzw(qa,&z,&w[0].a,'a');
getzw(qb,&z,&w[0].b,'b');
/* fill gb with basis info for laplacian calculation */
m = B ->vert;
md = DB->vert;
for(i = 0,n=0; i < Nverts; ++i,++n)
fill_gradbase(gb[n],m+i,md+i,fac); //OVERLOAD CALL: fill_gradbase: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
for(i = 0; i < Nedges; ++i){
m1 = B ->edge[i];
md1 = DB->edge[i];
for(j = 0; j < edge[i].l; ++j,++n)
fill_gradbase(gb[n],m1+j,md1+j,fac); //OVERLOAD CALL: fill_gradbase: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
}
/* faces */
for(i = 0; i < Nfaces; ++i){
ll = face[i].l;
for(j = 0; j < ll; ++j){
m1 = B ->face[i][j];
md1 = DB->face[i][j];
for(k = 0; k < ll-j; ++k,++n)
fill_gradbase(gb[n],m1+k,md1+k,fac); //OVERLOAD CALL: fill_gradbase: Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex), Gradient.c(Element)
}
}
/* multiply by weights */
for(n = 0; n < Nmodes; ++n){
Tri_mvmul2d(qa,qb,qc,gb[n] ,w,gb1[n]);
Tri_mvmul2d(qa,qb,qc,gb[n]+1,w,gb1[n]+1);
}
/* `A' matrix */
for(i = 0; i < nbl; ++i)
for(j = i; j < nbl; ++j)
a[i][j] = a[j][i] = jac*iprodlap(gb[i],gb1[j],fac+1); //OVERLOAD CALL: iprodlap: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri)
/* `B' matrix */
for(i = 0; i < nbl; ++i)
for(j = nbl; j < Nmodes; ++j)
b[i][j-nbl] = jac*iprodlap(gb[i],gb1[j],fac+1); //OVERLOAD CALL: iprodlap: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri)
/* 'C' matrix */
for(i = nbl; i < Nmodes; ++i)
for(j = nbl; j < Nmodes; ++j)
c[i-nbl][j-nbl] = c[j-nbl][i-nbl] = jac*iprodlap(gb[i],gb1[j],fac+1); //OVERLOAD CALL: iprodlap: InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex), InnerProd.c(Element), InnerProd.c(Tri)
free_mvec(gb[0]) ; free((char *) gb);
free_mvec(gb1[0]); free((char *) gb1);
// new addition 7/18/96
free((char*)w);
}
void Quad::LapMat(LocMat *){
fprintf(stderr,"Quad::LapMat should be calling Quad::HelmMatC\n");
}
void Tet::LapMat(LocMat *lap){
Metric *lambda = (Metric*)calloc(1, sizeof(Metric));
HelmMatC(lap, lambda); //OVERLOAD CALL: HelmMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
free(lambda);
}
void Pyr::LapMat(LocMat *lap){
Metric *lambda = (Metric*)calloc(1, sizeof(Metric));
HelmMatC(lap, lambda); //OVERLOAD CALL: HelmMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
free(lambda);
}
void Prism::LapMat(LocMat *lap){
Metric *lambda = (Metric*)calloc(1, sizeof(Metric));
HelmMatC(lap, lambda); //OVERLOAD CALL: HelmMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
free(lambda);
}
void Hex::LapMat(LocMat *lap){
Metric *lambda = (Metric*)calloc(1, sizeof(Metric));
HelmMatC(lap, lambda); //OVERLOAD CALL: HelmMatC: Matrix.c(Tri), Matrix.c(Quad), Matrix.c(Tet), Matrix.c(Pyr), Matrix.c(Prism), Matrix.c(Hex), Matrix.c(Element)
free(lambda);
}
void Element::LapMat (LocMat *){ERR;} // return Laplacian op. //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
/*
Function name: Element::project
Function Purpose:
Argument 1: LocMat *m
Purpose:
Argument 2: Bsystem *Ubsys
Purpose:
Function Notes:
*/
void Tri::project(LocMat *m, Bsystem *Ubsys){
register int i,j;
int pi,pj,nsolve;
int *list;
const int asize = m->asize;
double **a = m->a;
list = Ubsys->bmap[id];
if(Ubsys->rslv){ /* project a to first recurrsion */
/* close pack the a matrix to elimate any extra zero for variable case */
for(i = 1; i < asize; ++i)
dcopy(asize,a[i],1,*a+i*asize,1);
Project_Recur(id, asize, *a, list, 0, Ubsys);
}
else{ /* normal static condenstion */
double **inva = Ubsys->Gmat->inva;
const int bwidth = Ubsys->Gmat->bwidth_a;
if(nsolve = Ubsys->nsolve){
if(2*bwidth < nsolve){
/* store as packed banded system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj >= pi)
inva[pi][pj-pi] += a[i][j];
}
else{
/* store as symmetric system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj <= pi)
inva[0][(pj*(pj+1))/2 + (nsolve-pj)*pj + (pi-pj)] += a[i][j];
}
}
}
}
void Quad::project(LocMat *m, Bsystem *Ubsys){
register int i,j;
int pi,pj,nsolve;
int *list;
const int asize = m->asize;
double **a = m->a;
list = Ubsys->bmap[id];
if(Ubsys->rslv){ /* project a to first recurrsion */
/* close pack the a matrix to elimate any extra zero for variable case */
for(i = 1; i < asize; ++i)
dcopy(asize,a[i],1,*a+i*asize,1);
Project_Recur(id,asize, *a, list, 0, Ubsys);
}
else{ /* normal static condenstion */
const int bwidth = Ubsys->Gmat->bwidth_a;
double **inva = Ubsys->Gmat->inva;
if(nsolve = Ubsys->nsolve){
if(2*bwidth < nsolve){
/* store as packed banded system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj >= pi)
inva[pi][pj-pi] += a[i][j];
}
else{
/* store as symmetric system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj <= pi)
inva[0][(pj*(pj+1))/2 + (nsolve-pj)*pj + (pi-pj)] += a[i][j];
}
}
}
}
void Tet::project(LocMat *m, Bsystem *Ubsys){
register int i,j;
int pi,pj,nsolve;
int *list;
const int asize = m->asize;
const int bwidth = Ubsys->Gmat->bwidth_a;
double **a = m->a, **inva = Ubsys->Gmat->inva;
list = Ubsys->bmap[id];
if(Ubsys->rslv){ /* project a to first recurrsion */
/* close pack the a matrix to elimate any extra zero for variable case */
for(i = 1; i < asize; ++i)
dcopy(asize,a[i],1,*a+i*asize,1);
Project_Recur(id, asize, *a, list, 0, Ubsys);
}
else{ /* normal static condenstion */
if(nsolve = Ubsys->nsolve){
if(2*bwidth < nsolve){
/* store as packed banded system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj >= pi)
inva[pi][pj-pi] += a[i][j];
}
else{
/* store as symmeTetc system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj <= pi)
inva[0][(pj*(pj+1))/2 + (nsolve-pj)*pj + (pi-pj)] += a[i][j];
}
}
}
}
void Pyr::project(LocMat *m, Bsystem *Ubsys){
register int i,j;
int pi,pj,nsolve;
int *list;
const int asize = m->asize;
const int bwidth = Ubsys->Gmat->bwidth_a;
double **a = m->a, **inva = Ubsys->Gmat->inva;
list = Ubsys->bmap[id];
if(Ubsys->rslv){ /* project a to first recurrsion */
/* close pack the a matrix to elimate any extra zero for variable case */
for(i = 1; i < asize; ++i)
dcopy(asize,a[i],1,*a+i*asize,1);
Project_Recur(id, asize, *a, list, 0, Ubsys);
}
else{ /* normal static condenstion */
if(nsolve = Ubsys->nsolve){
if(2*bwidth < nsolve){
/* store as packed banded system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj >= pi)
inva[pi][pj-pi] += a[i][j];
}
else{
/* store as symmePyrc system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj <= pi)
inva[0][(pj*(pj+1))/2 + (nsolve-pj)*pj + (pi-pj)] += a[i][j];
}
}
}
}
void Prism::project(LocMat *m, Bsystem *Ubsys){
register int i,j;
int pi,pj,nsolve;
int *list;
const int asize = m->asize;
const int bwidth = Ubsys->Gmat->bwidth_a;
double **a = m->a, **inva = Ubsys->Gmat->inva;
list = Ubsys->bmap[id];
if(Ubsys->rslv){ /* project a to first recurrsion */
/* close pack the a matrix to elimate any extra zero for variable case */
for(i = 1; i < asize; ++i)
dcopy(asize,a[i],1,*a+i*asize,1);
Project_Recur(id, asize, *a, list, 0, Ubsys);
}
else{ /* normal static condenstion */
if(nsolve = Ubsys->nsolve){
if(2*bwidth < nsolve){
/* store as packed banded system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj >= pi)
inva[pi][pj-pi] += a[i][j];
}
else{
/* store as symmetric system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj <= pi)
inva[0][(pj*(pj+1))/2 + (nsolve-pj)*pj + (pi-pj)] += a[i][j];
}
}
}
}
void Hex::project(LocMat *m, Bsystem *Ubsys){
register int i,j;
int pi,pj,nsolve;
int *list;
const int asize = m->asize;
const int bwidth = Ubsys->Gmat->bwidth_a;
double **a = m->a, **inva = Ubsys->Gmat->inva;
list = Ubsys->bmap[id];
if(Ubsys->rslv){ /* project a to first recurrsion */
/* close pack the a matrix to elimate any extra zero for variable case */
for(i = 1; i < asize; ++i)
dcopy(asize,a[i],1,*a+i*asize,1);
Project_Recur(id, asize, *a, list, 0, Ubsys);
}
else{ /* normal static condenstion */
if(nsolve = Ubsys->nsolve){
if(2*bwidth < nsolve){
/* store as packed banded system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj >= pi)
inva[pi][pj-pi] += a[i][j];
}
else{
/* store as symmeHexc system */
for(i = 0; i < asize; ++i)
if((pi=list[i]) < nsolve)
for(j = 0; j < asize; ++j)
if((pj=list[j]) < nsolve)
if(pj <= pi)
inva[0][(pj*(pj+1))/2 + (nsolve-pj)*pj + (pi-pj)] += a[i][j];
}
}
}
}
void Element::project(LocMat *, Bsystem *){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
/*
Function name: Element::condense
Function Purpose:
Argument 1: LocMat *m
Purpose:
Argument 2: Bsystem *Ubsys
Purpose:
Function Notes:
*/
void Tri::condense(LocMat *m, Bsystem *Ubsys){
register int i,j,k;
int asize,csize,L,info,n;
int *bw = Ubsys->Gmat->bwidth_c, geom_id = geom->id;
double **a = m->a, **b = m->b, **c = m->c;
double *invc, *binvc;
/* calc local boundary matrix size */
asize = Nverts;
for(i = 0; i < Nedges; ++i)
asize += edge[i].l;
csize = Nfmodes(); //OVERLOAD CALL: Nfmodes: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
m->asize = asize;
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
Ubsys->Gmat->a[geom_id] = dvector(0,asize*(asize+1)/2-1);
if(csize){
if(curvX || (Ubsys->lambda->p))
bw[geom_id] = csize;
else{
bw[geom_id] = 2*face->l+1;
}
if(csize > 2*bw[geom_id])
invc = dvector(0,csize*bw[geom_id]-1);
else
invc = dvector(0,csize*(csize+1)/2-1);
binvc = dvector(0,asize*csize-1);
for(i = 0; i < asize; ++i)
dcopy(csize,b[i],1,binvc + i*csize,1);
PackMatrix(c,csize,invc,bw[geom_id]);
/* calc invc(c) */
FacMatrix(invc,csize,bw[geom_id]);
/* calc inv(c)*b^T */
if(csize > 2*bw[geom_id])
dpbtrs('L',csize,bw[geom_id]-1,asize,invc,bw[geom_id],binvc,csize,info);
else
dpptrs('L',csize,asize,invc,binvc,csize,info);
/* A - b inv(c) b^T */
for(i = 0; i < asize; ++i)
for(j = 0; j < asize; ++j)
a[i][j] -= ddot(csize,b[i],1,binvc+j*csize,1);
}
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
PackMatrix(a,asize,Ubsys->Gmat->a[geom_id],asize);
else{
/* sort connectivity issues */
n = Nverts;
for(i = 0; i < Nedges; ++i){
L = edge[i].l;
if(edge[i].con){
for(j = 1; j < L; j+=2){
dscal(asize,-1.0,a[n+j],1);
for(k=0; k < asize; ++k)
a[k][n+j] *= -1.0;
}
}
n += L;
}
}
if(csize)
if(Ubsys->Gmat->invc[geom_id]){
free(invc);
free(binvc);
}
else{
Ubsys->Gmat-> invc[geom_id] = invc;
Ubsys->Gmat->binvc[geom_id] = binvc;
}
else /* just set it to non-zero value */
Ubsys->Gmat->invc[geom_id] = (double *)malloc(1);
}
void Quad::condense(LocMat *m, Bsystem *Ubsys){
register int i,j,k;
int asize,csize,L,info,n;
int *bw = Ubsys->Gmat->bwidth_c, geom_id = geom->id;
double **a = m->a, **b = m->b, **c = m->c;
double *invc, *binvc;
/* calc local boundary matrix size */
asize = Nverts;
for(i = 0; i < Nedges; ++i)
asize += edge[i].l;
csize = face->l*face->l;
m->asize = asize;
// tcew -- store as a symmetrix (full) matrix ??
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
Ubsys->Gmat->a[geom_id] = dvector(0,asize*(asize+1)/2-1);
if(csize){
bw[geom_id] = csize;
#if 0
if(curve->type != T_Straight || (Ubsys->lambda->p))
bw[geom_id] = csize;
else{
bw[geom_id] = 3*face->l+1;
}
#endif
if(csize > 2*bw[geom_id])
invc = dvector(0,csize*bw[geom_id]-1);
else
invc = dvector(0,csize*(csize+1)/2-1);
binvc = dvector(0,asize*csize-1);
for(i = 0; i < asize; ++i)
dcopy(csize,b[i],1,binvc + i*csize,1);
PackMatrix(c,csize,invc,bw[geom_id]);
/* calc invc(c) */
FacMatrix(invc,csize,bw[geom_id]);
/* calc inv(c)*b^T */
if(csize > 2*bw[geom_id])
dpbtrs('L',csize,bw[geom_id]-1,asize,invc,bw[geom_id],binvc,csize,info);
else
dpptrs('L',csize,asize,invc,binvc,csize,info);
/* A - b inv(c) b^T */
for(i = 0; i < asize; ++i)
for(j = 0; j < asize; ++j)
a[i][j] -= ddot(csize,b[i],1,binvc+j*csize,1);
}
// Possible
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
PackMatrix(a,asize,Ubsys->Gmat->a[geom_id],asize);
else{
/* sort connectivity issues */
n = Nverts;
for(i = 0; i < Nedges; ++i){
L = edge[i].l;
if(edge[i].con){
for(j = 1; j < L; j+=2){
dscal(asize,-1.0,a[n+j],1);
for(k=0; k < asize; ++k)
a[k][n+j] *= -1.0;
}
}
n += L;
}
}
if(csize)
if(Ubsys->Gmat->invc[geom_id]){
free(invc);
free(binvc);
}
else{
Ubsys->Gmat-> invc[geom_id] = invc;
Ubsys->Gmat->binvc[geom_id] = binvc;
}
else /* just set it to non-zero value */
Ubsys->Gmat->invc[geom_id] = (double *)malloc(1);
}
void Tet::condense(LocMat *m, Bsystem *Ubsys){
register int i,j,k,j1,n1;
int asize,csize,L,info,n;
int *bw = Ubsys->Gmat->bwidth_c, geomid = geom->id;
double **a = m->a, **b = m->b, **c = m->c;
double *invc, *binvc;
/* calc local boundary matrix size */
asize = Nbmodes;
csize = Nmodes-Nbmodes;
m->asize = asize;
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
Ubsys->Gmat->a[geomid] = dvector(0,asize*(asize+1)/2-1);
if(csize){
if(curvX)
bw[geomid] = csize;
else
bw[geomid] = interior_l*(interior_l+1);
if(csize > 2*bw[geomid])
invc = dvector(0,csize*bw[geomid]-1);
else
invc = dvector(0,csize*(csize+1)/2-1);
binvc = dvector(0,asize*csize-1);
for(i = 0; i < asize; ++i)
dcopy(csize,b[i],1,binvc + i*csize,1);
PackMatrix(c,csize,invc,bw[geomid]);
/* calc invc(c) */
FacMatrix(invc,csize,bw[geomid]);
/* calc inv(c)*b^T */
if(csize > 2*bw[geomid])
dpbtrs('L',csize,bw[geomid]-1,asize,invc,bw[geomid],binvc,csize,info);
else
dpptrs('L',csize,asize,invc,binvc,csize,info);
/* A - b inv(c) b^T */
for(i = 0; i < asize; ++i)
for(j = 0; j < asize; ++j)
a[i][j] -= ddot(csize,b[i],1,binvc+j*csize,1);
}
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
PackMatrix(a,asize,Ubsys->Gmat->a[geomid],asize);
else{
/* sort connectivity issues */
n = Nverts;
for(i = 0; i < Nedges; ++i){
L = edge[i].l;
if(edge[i].con){
for(j = 1; j < L; j+=2){
dscal(asize,-1.0,a[n+j],1);
for(k=0; k < asize; ++k)
a[k][n+j] *= -1.0;
}
}
n += L;
}
for(i=0;i<Nfaces;++i){
L = face[i].l;
if(face[i].con){
for(j=1,n1=n+L;j<L;j+=2){
for(j1=0;j1<L-j;++j1){
dsmul(asize,-1.0,a[n1+j1],1,a[n1+j1],1);
for(k=0;k<asize;++k)
a[k][n1+j1] *= -1.0;
}
n1+=2*(L-j)-1;
}
}
n += L*(L+1)/2;
}
}
if(csize)
if(Ubsys->Gmat->invc[geomid]){
free(invc);
free(binvc);
}
else{
Ubsys->Gmat-> invc[geomid] = invc;
Ubsys->Gmat->binvc[geomid] = binvc;
}
else /* just set it to non-zero value */
Ubsys->Gmat->invc[geomid] = (double *)malloc(1);
}
void Pyr::condense(LocMat *m, Bsystem *Ubsys){
register int i,j,k;
int asize,csize,L,info,n,n1,j1;
int *bw = Ubsys->Gmat->bwidth_c, geom_id = geom->id;
double **a = m->a, **b = m->b, **c = m->c;
double *invc, *binvc;
/* calc local boundary matrix size */
asize = Nbmodes;
csize = Nmodes-Nbmodes;
m->asize = asize;
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
Ubsys->Gmat->a[geom_id] = dvector(0,asize*(asize+1)/2-1);
if(csize){
bw[geom_id] = csize;
if(csize > bw[geom_id])
invc = dvector(0,csize*bw[geom_id]-1);
else
invc = dvector(0,csize*(csize+1)/2-1);
binvc = dvector(0,asize*csize-1);
for(i = 0; i < asize; ++i)
dcopy(csize,b[i],1,binvc + i*csize,1);
PackMatrix(c,csize,invc,bw[geom_id]);
/* calc invc(c) */
FacMatrix(invc,csize,bw[geom_id]);
/* calc inv(c)*b^T */
if(csize > 2*bw[geom_id])
dpbtrs('L',csize,bw[geom_id]-1,asize,invc,bw[geom_id],binvc,csize,info);
else
dpptrs('L',csize,asize,invc,binvc,csize,info);
/* A - b inv(c) b^T */
for(i = 0; i < asize; ++i)
for(j = 0; j < asize; ++j)
a[i][j] -= ddot(csize,b[i],1,binvc+j*csize,1);
}
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
PackMatrix(a,asize,Ubsys->Gmat->a[geom_id],asize);
else{
/* sort connectivity issues */
n = NPyr_verts;
for(i = 0; i < NPyr_edges; ++i){
L = edge[i].l;
if(edge[i].con){
for(j = 1; j < L; j+=2){
dscal(asize,-1.0,a[n+j],1);
for(k=0; k < asize; ++k)
a[k][n+j] *= -1.0;
}
}
n += L;
}
#if 1
// needs to be fixed
for(i=0;i<NPyr_faces;++i){
L = face[i].l;
// fix for triangle faces
if(Nfverts(i) == 4){ //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
switch(face[i].con ){
case 0:
// nothing to be done
break;
case 1:
// negate in 'a' direction
for(k=0;k<L;++k)
for(j=1;j<L;j+=2){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
case 2:
// negate in 'b' direction
for(k=1;k<L;k+=2)
for(j=0;j<L;++j){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
case 3:
// negate in 'a' direction
for(k=0;k<L;++k)
for(j=1;j<L;j+=2){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
// negate in 'b' direction
for(k=1;k<L;k+=2)
for(j=0;j<L;++j){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
}
n += L*L;
}
else{
if(face[i].con){
for(j=1,n1=n+L;j<L;j+=2){
for(j1=0;j1<L-j;++j1){
dsmul(asize,-1.0,a[n1+j1],1,a[n1+j1],1);
for(k=0;k<asize;++k)
a[k][n1+j1] *= -1.0;
}
n1+=2*(L-j)-1;
}
}
n += L*(L+1)/2;
}
}
#endif
}
if(csize)
if(Ubsys->Gmat->invc[geom_id]){
free(invc);
free(binvc);
}
else{
Ubsys->Gmat-> invc[geom_id] = invc;
Ubsys->Gmat->binvc[geom_id] = binvc;
}
else /* just set it to non-zero value */
Ubsys->Gmat->invc[geom_id] = (double *)malloc(1);
}
void Prism::condense(LocMat *m, Bsystem *Ubsys){
register int i,j,k;
int asize,csize,L,info,n,n1,j1;
int *bw = Ubsys->Gmat->bwidth_c, geom_id = geom->id;
double **a = m->a, **b = m->b, **c = m->c;
double *invc, *binvc;
/* calc local boundary matrix size */
asize = Nbmodes;
csize = Nmodes-Nbmodes;
m->asize = asize;
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
Ubsys->Gmat->a[geom_id] = dvector(0,asize*(asize+1)/2-1);
if(csize){
bw[geom_id] = csize;
if(csize > bw[geom_id])
invc = dvector(0,csize*bw[geom_id]-1);
else
invc = dvector(0,csize*(csize+1)/2-1);
binvc = dvector(0,asize*csize-1);
for(i = 0; i < asize; ++i)
dcopy(csize,b[i],1,binvc + i*csize,1);
PackMatrix(c,csize,invc,bw[geom_id]);
/* calc invc(c) */
FacMatrix(invc,csize,bw[geom_id]);
/* calc inv(c)*b^T */
if(csize > 2*bw[geom_id])
dpbtrs('L',csize,bw[geom_id]-1,asize,invc,bw[geom_id],binvc,csize,info);
else
dpptrs('L',csize,asize,invc,binvc,csize,info);
/* A - b inv(c) b^T */
for(i = 0; i < asize; ++i)
for(j = 0; j < asize; ++j)
a[i][j] -= ddot(csize,b[i],1,binvc+j*csize,1);
}
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
PackMatrix(a,asize,Ubsys->Gmat->a[geom_id],asize);
else{
/* sort connectivity issues */
n = NPrism_verts;
for(i = 0; i < NPrism_edges; ++i){
L = edge[i].l;
if(edge[i].con){
for(j = 1; j < L; j+=2){
dscal(asize,-1.0,a[n+j],1);
for(k=0; k < asize; ++k)
a[k][n+j] *= -1.0;
}
}
n += L;
}
// needs to be fixed
for(i=0;i<NPrism_faces;++i){
L = face[i].l;
// fix for triangle faces
if(Nfverts(i) == 4){ //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
switch(face[i].con ){
case 0:
// nothing to be done
break;
case 1:
// negate in 'a' direction
for(k=0;k<L;++k)
for(j=1;j<L;j+=2){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
case 2:
// negate in 'b' direction
for(k=1;k<L;k+=2)
for(j=0;j<L;++j){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
case 3:
// negate in 'a' direction
for(k=0;k<L;++k)
for(j=1;j<L;j+=2){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
// negate in 'b' direction
for(k=1;k<L;k+=2)
for(j=0;j<L;++j){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
}
n += L*L;
}
else{
if(face[i].con){
for(j=1,n1=n+L;j<L;j+=2){
for(j1=0;j1<L-j;++j1){
dsmul(asize,-1.0,a[n1+j1],1,a[n1+j1],1);
for(k=0;k<asize;++k)
a[k][n1+j1] *= -1.0;
}
n1+=2*(L-j)-1;
}
}
n += L*(L+1)/2;
}
}
}
if(csize)
if(Ubsys->Gmat->invc[geom_id]){
free(invc);
free(binvc);
}
else{
Ubsys->Gmat-> invc[geom_id] = invc;
Ubsys->Gmat->binvc[geom_id] = binvc;
}
else /* just set it to non-zero value */
Ubsys->Gmat->invc[geom_id] = (double *)malloc(1);
}
void Hex::condense(LocMat *m, Bsystem *Ubsys){
register int i,j,k;
int asize,csize,L,info,n;
int *bw = Ubsys->Gmat->bwidth_c, geom_id = geom->id;
double **a = m->a, **b = m->b, **c = m->c;
double *invc, *binvc;
/* calc local boundary matrix size */
asize = Nbmodes;
csize = Nmodes-Nbmodes;
m->asize = asize;
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
Ubsys->Gmat->a[geom_id] = dvector(0,asize*(asize+1)/2-1);
if(csize){
bw[geom_id] = csize;
if(csize > bw[geom_id])
invc = dvector(0,csize*bw[geom_id]-1);
else
invc = dvector(0,csize*(csize+1)/2-1);
binvc = dvector(0,asize*csize-1);
for(i = 0; i < asize; ++i)
dcopy(csize,b[i],1,binvc + i*csize,1);
PackMatrix(c,csize,invc,bw[geom_id]);
/* calc invc(c) */
FacMatrix(invc,csize,bw[geom_id]);
/* calc inv(c)*b^T */
if(csize > 2*bw[geom_id])
dpbtrs('L',csize,bw[geom_id]-1,asize,invc,bw[geom_id],binvc,csize,info);
else
dpptrs('L',csize,asize,invc,binvc,csize,info);
/* A - b inv(c) b^T */
for(i = 0; i < asize; ++i)
for(j = 0; j < asize; ++j)
a[i][j] -= ddot(csize,b[i],1,binvc+j*csize,1);
}
if(Ubsys->smeth == iterative&&!Ubsys->rslv)
PackMatrix(a,asize,Ubsys->Gmat->a[geom_id],asize);
else{
/* sort connectivity issues */
n = NHex_verts;
for(i = 0; i < NHex_edges; ++i){
L = edge[i].l;
if(edge[i].con){
for(j = 1; j < L; j+=2){
dscal(asize,-1.0,a[n+j],1);
for(k=0; k < asize; ++k)
a[k][n+j] *= -1.0;
}
}
n += L;
}
#if 1
// needs to be fixed
for(i=0;i<NHex_faces;++i){
L = face[i].l;
switch(face[i].con ){
case 0:
// nothing to be done
break;
case 1:
// negate in 'a' direction
for(k=0;k<L;++k)
for(j=1;j<L;j+=2){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
case 2:
// negate in 'b' direction
for(k=1;k<L;k+=2)
for(j=0;j<L;++j){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
case 3:
// negate in 'a' direction
for(k=0;k<L;++k)
for(j=1;j<L;j+=2){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
// negate in 'b' direction
for(k=1;k<L;k+=2)
for(j=0;j<L;++j){
dscal(asize,-1.0,a[n+j+L*k],1);
dscal(asize,-1.0,a[0]+n+j+L*k, asize);
}
break;
}
n += L*L;
}
#endif
}
if(csize)
if(Ubsys->Gmat->invc[geom_id]){
free(invc);
free(binvc);
}
else{
Ubsys->Gmat-> invc[geom_id] = invc;
Ubsys->Gmat->binvc[geom_id] = binvc;
}
else /* just set it to non-zero value */
Ubsys->Gmat->invc[geom_id] = (double *)malloc(1);
}
void Element::condense(LocMat *, Bsystem *){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
double Tri_mass_mprod(Element *E, Mode *m, double *wvec);
double Quad_mass_mprod(Element *E, Mode *m, double *wvec);
double Tet_mass_mprod(Element *E, Mode *m, double *wvec);
double Pyr_mass_mprod(Element *E, Mode *m, double *wvec);
double Prism_mass_mprod(Element *E, Mode *m, double *wvec);
double Hex_mass_mprod(Element *E, Mode *m, double *wvec);
/*
Function name: Element::fill_diag_massmat
Function Purpose:
Function Notes:
*/
void Tri::fill_diag_massmat(){
Basis *b;
double *wa, *wb;
int i,j,k;
double *wvec = dvector(0, qtot-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)
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'b');
Mode mw, *m;
mw.a = wa; mw.b = wb;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Tri_mass_mprod(this, m, wvec);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Tri_mass_mprod(this, m, wvec);
}
for(i=0;i<Nfaces;++i)
for(j=0;j<face[i].l;++j)
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Tri_mass_mprod(this, m, wvec);
}
free(wvec);
}
void Quad::fill_diag_massmat(){
Basis *b;
double *wa, *wb;
int i,j,k;
double *wvec = dvector(0, qtot-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)
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
Mode mw, *m;
mw.a = wa; mw.b = wb;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Quad_mass_mprod(this, m, wvec);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Quad_mass_mprod(this, m, wvec);
}
for(i=0;i<Nfaces;++i)
for(j=0;j<face[i].l;++j)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Quad_mass_mprod(this, m, wvec);
}
free(wvec);
}
void Tet::fill_diag_massmat(){
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)
double *wa, *wb, *wc;
int i,j,k;
double *wvec = dvector(0, qtot-1);
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'b');
getzw(qc,&wc,&wc,'c');
Mode mw, *m;
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Tet_mass_mprod(this, m, wvec);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Tet_mass_mprod(this, m, wvec);
}
for(i=0;i<Nfaces;++i){
for(j=0;j<face[i].l;++j){
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Tet_mass_mprod(this, m, wvec);
}
}
}
for(i=0;i<interior_l;++i)
for(j=0;j<interior_l-i;++j)
for(k=0;k<interior_l-i-j;++k){
// fix
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Tet_mass_mprod(this, m, wvec);
}
free(wvec);
}
void Pyr::fill_diag_massmat(){
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)
double *wa, *wb, *wc;
int i,j,k,nfv;
double *wvec = dvector(0, qtot-1);
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
getzw(qc,&wc,&wc,'c');
Mode mw, *m;
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Pyr_mass_mprod(this, m, wvec);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Pyr_mass_mprod(this, m, wvec);
}
for(i=0;i<Nfaces;++i){
nfv = Nfverts(i); //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(j=0;j<face[i].l;++j){
if(nfv == 4)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Pyr_mass_mprod(this, m, wvec);
}
else
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Pyr_mass_mprod(this, m, wvec);
}
}
}
for(i=0;i<interior_l;++i)
for(j=0;j<interior_l-i;++j)
for(k=0;k<interior_l-i-j;++k){
// fix
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Pyr_mass_mprod(this, m, wvec);
}
free(wvec);
}
void Prism::fill_diag_massmat(){
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)
double *wa, *wb, *wc;
int i,j,k,nfv;
double *wvec = dvector(0, qtot-1);
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
getzw(qc,&wc,&wc,'b');
Mode mw, *m;
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Prism_mass_mprod(this, m, wvec);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Prism_mass_mprod(this, m, wvec);
}
for(i=0;i<Nfaces;++i){
nfv = Nfverts(i); //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(j=0;j<face[i].l;++j){
if(nfv == 4)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Prism_mass_mprod(this, m, wvec);
}
else
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Prism_mass_mprod(this, m, wvec);
}
}
}
for(i=0;i<interior_l-1;++i)
for(j=0;j<interior_l;++j)
for(k=0;k<interior_l-i-1;++k){
// fix
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Prism_mass_mprod(this, m, wvec);
}
free(wvec);
}
void Hex::fill_diag_massmat(){
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)
double *wa, *wb, *wc;
int i,j,k;
double *wvec = dvector(0, qtot-1);
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
getzw(qc,&wc,&wc,'a');
Mode mw, *m;
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Hex_mass_mprod(this, m, wvec);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Hex_mass_mprod(this, m, wvec);
}
for(i=0;i<Nfaces;++i)
for(j=0;j<face[i].l;++j)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Hex_mass_mprod(this, m, wvec);
}
for(i=0;i<interior_l;++i)
for(j=0;j<interior_l;++j)
for(k=0;k<interior_l;++k){
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Hex_mass_mprod(this, m, wvec);
}
free(wvec);
}
void Element::fill_diag_massmat(){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
double Tri_mass_mprod(Element *E, Mode *m, double *wvec){
static double *tmp = 0;
if(!tmp)
tmp = dvector(0, QGmax*QGmax-1);
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
double d = ddot(E->qtot, tmp, 1, wvec, 1);
return d;
}
double Quad_mass_mprod(Element *E, Mode *m, double *wvec){
static double *tmp = 0;
if(!tmp)
tmp = dvector(0, QGmax*QGmax-1);
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
double d = ddot(E->qtot, tmp, 1, wvec, 1);
return d;
}
double Tet_mass_mprod(Element *E, Mode *m, double *wvec){
static double *tmp = 0;
if(!tmp)
tmp = dvector(0, QGmax*QGmax*QGmax-1);
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
double d = ddot(E->qtot, tmp, 1, wvec, 1);
return d;
}
double Pyr_mass_mprod(Element *E, Mode *m, double *wvec){
static double *tmp = 0;
if(!tmp)
tmp = dvector(0, QGmax*QGmax*QGmax-1);
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
double d = ddot(E->qtot, tmp, 1, wvec, 1);
return d;
}
double Prism_mass_mprod(Element *E, Mode *m, double *wvec){
static double *tmp = 0;
if(!tmp)
tmp = dvector(0, QGmax*QGmax*QGmax-1);
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
double d = ddot(E->qtot, tmp, 1, wvec, 1);
return d;
}
double Hex_mass_mprod(Element *E, Mode *m, double *wvec){
static double *tmp = 0;
if(!tmp)
tmp = dvector(0, QGmax*QGmax*QGmax-1);
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
double d = ddot(E->qtot, tmp, 1, wvec, 1);
return d;
}
double Tri_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda);
double Quad_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda);
double Tet_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda);
double Pyr_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda);
double Prism_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda);
double Hex_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda);
/*
Function name: Element::fill_diag_helmmat
Function Purpose:
Argument 1: Metric *lambda
Purpose:
Function Notes:
*/
void Tri::fill_diag_helmmat(Metric *lambda){
int i, j, k;
Basis *b;
double *wa, *wb;
double *wvec = dvector(0, qtot-1);
Mode mw, *m;
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)
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'b');
mw.a = wa; mw.b = wb;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Tri_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Tri_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nfaces;++i)
for(j=0;j<face[i].l;++j)
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Tri_helm_mprod(this, m, wvec, lambda);
}
free(wvec);
}
void Quad::fill_diag_helmmat( Metric *lambda){
int i, j, k;
Basis *b;
double *wa, *wb;
double *wvec = dvector(0, qtot-1);
Mode mw, *m;
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)
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
mw.a = wa; mw.b = wb;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Quad_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Quad_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nfaces;++i)
for(j=0;j<face[i].l;++j)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Quad_helm_mprod(this, m, wvec, lambda);
}
free(wvec);
}
void Tet::fill_diag_helmmat(Metric *lambda){
int i, j, k;
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)
double *wa, *wb, *wc;
double *wvec = dvector(0, qtot-1);
Mode mw, *m;
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'b');
getzw(qc,&wc,&wc,'c');
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Tet_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Tet_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nfaces;++i){
for(j=0;j<face[i].l;++j){
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Tet_helm_mprod(this, m, wvec, lambda);
}
}
}
for(i=0;i<interior_l;++i)
for(j=0;j<interior_l-i;++j)
for(k=0;k<interior_l-i-j;++k){
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Tet_helm_mprod(this, m, wvec, lambda);
}
free(wvec);
}
void Pyr::fill_diag_helmmat(Metric *lambda){
int i, j, k, nfv;
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)
double *wa, *wb, *wc;
double *wvec = dvector(0, qtot-1);
Mode mw, *m;
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
getzw(qc,&wc,&wc,'c');
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Pyr_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
// fix
m = b->edge[i]+j;
edge[i].hj[j] = Pyr_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nfaces;++i){
nfv = Nfverts(i); //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(j=0;j<face[i].l;++j){
if(nfv == 4)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Pyr_helm_mprod(this, m, wvec, lambda);
}
else
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Pyr_helm_mprod(this, m, wvec, lambda);
}
}
}
for(i=0;i<interior_l;++i)
for(j=0;j<interior_l-i;++j)
for(k=0;k<interior_l-i-j;++k){
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Pyr_helm_mprod(this, m, wvec, lambda);
}
free(wvec);
}
void Prism::fill_diag_helmmat(Metric *lambda){
int i, j, k, nfv;
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)
double *wa, *wb, *wc;
double *wvec = dvector(0, qtot-1);
Mode mw, *m;
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
getzw(qc,&wc,&wc,'b');
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Prism_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
// fix
m = b->edge[i]+j;
edge[i].hj[j] = Prism_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nfaces;++i){
nfv = Nfverts(i); //OVERLOAD CALL: Nfverts: Nums.c(Tri), Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element)
for(j=0;j<face[i].l;++j){
if(nfv == 4)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Prism_helm_mprod(this, m, wvec, lambda);
}
else
for(k=0;k<face[i].l-j;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Prism_helm_mprod(this, m, wvec, lambda);
}
}
}
for(i=0;i<interior_l-1;++i)
for(j=0;j<interior_l;++j)
for(k=0;k<interior_l-i-1;++k){
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Prism_helm_mprod(this, m, wvec, lambda);
}
free(wvec);
}
void Hex::fill_diag_helmmat(Metric *lambda){
int i, j, k;
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)
double *wa, *wb, *wc;
double *wvec = dvector(0, qtot-1);
Mode mw, *m;
getzw(qa,&wa,&wa,'a');
getzw(qb,&wb,&wb,'a');
getzw(qc,&wc,&wc,'a');
mw.a = wa; mw.b = wb; mw.c = wc;
fillvec(&mw, wvec); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
if(curvX)
dvmul(qtot, wvec, 1, geom->jac.p, 1, wvec, 1);
else
dscal(qtot, geom->jac.d, wvec, 1);
for(i=0;i<Nverts;++i){
m = b->vert+i;
vert[i].hj[0] = Hex_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nedges;++i)
for(j=0;j<edge[i].l;++j){
m = b->edge[i]+j;
edge[i].hj[j] = Hex_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<Nfaces;++i)
for(j=0;j<face[i].l;++j)
for(k=0;k<face[i].l;++k){
m = b->face[i][j]+k;
face[i].hj[j][k] = Hex_helm_mprod(this, m, wvec, lambda);
}
for(i=0;i<interior_l;++i)
for(j=0;j<interior_l;++j)
for(k=0;k<interior_l;++k){
m = b->intr[i][j]+k;
hj_3d[i][j][k] = Hex_helm_mprod(this, m, wvec, lambda);
}
free(wvec);
}
void Element::fill_diag_helmmat(Metric *){ERR;} //OVERLOAD CALL: ERR: Element.c(?), hotel.h(?)
double Tri_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda){
static double *tmp = 0;
static double *dx = 0;
static double *dy = 0;
double *nul = 0;
if(!tmp){
tmp = dvector(0, QGmax*QGmax-1);
dx = dvector(0, QGmax*QGmax-1);
dy = dvector(0, QGmax*QGmax-1);
}
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
E->Grad_h(tmp, dx, dy, nul, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (E->qtot, dx, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dy, 1, dy, 1, dx, 1, dx, 1);
if(lambda && lambda->p)
dvmul(E->qtot, lambda->p, 1, dx, 1, dx, 1);
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
daxpy (E->qtot, lambda->d, tmp, 1, dx, 1);
double d = ddot(E->qtot, dx, 1, wvec, 1);
return d;
}
double Quad_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda){
static double *tmp = 0;
static double *dx = 0;
static double *dy = 0;
double *nul = 0;
if(!tmp){
tmp = dvector(0, QGmax*QGmax-1);
dx = dvector(0, QGmax*QGmax-1);
dy = dvector(0, QGmax*QGmax-1);
}
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
E->Grad_h(tmp, dx, dy, nul, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (E->qtot, dx, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dy, 1, dy, 1, dx, 1, dx, 1);
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
if(lambda && lambda->p)
dvmul(E->qtot, lambda->p, 1, dx, 1, dx, 1);
daxpy (E->qtot, lambda->d, tmp, 1, dx, 1);
double d = ddot(E->qtot, dx, 1, wvec, 1);
return d;
}
double Tet_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda){
static double *tmp = 0;
static double *dx = 0;
static double *dy = 0;
static double *dz = 0;
if(!tmp){
tmp = dvector(0, QGmax*QGmax*QGmax-1);
dx = dvector(0, QGmax*QGmax*QGmax-1);
dy = dvector(0, QGmax*QGmax*QGmax-1);
dz = dvector(0, QGmax*QGmax*QGmax-1);
}
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
E->Grad_h(tmp, dx, dy, dz, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (E->qtot, dx, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dy, 1, dy, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dz, 1, dz, 1, dx, 1, dx, 1);
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
if(lambda->p)
dvvtvp(E->qtot, lambda->p, 1, tmp, 1, dx, 1, dx, 1);
else
daxpy (E->qtot, lambda->d, tmp, 1, dx, 1);
double d = ddot(E->qtot, dx, 1, wvec, 1);
return d;
}
double Pyr_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda){
static double *tmp = 0;
static double *dx = 0;
static double *dy = 0;
static double *dz = 0;
if(!tmp){
tmp = dvector(0, QGmax*QGmax*QGmax-1);
dx = dvector(0, QGmax*QGmax*QGmax-1);
dy = dvector(0, QGmax*QGmax*QGmax-1);
dz = dvector(0, QGmax*QGmax*QGmax-1);
}
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
E->Grad_h(tmp, dx, dy, dz, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (E->qtot, dx, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dy, 1, dy, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dz, 1, dz, 1, dx, 1, dx, 1);
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
if(lambda->p)
dvvtvp (E->qtot, lambda->p, 1, tmp, 1, dx, 1, dx, 1);
else
daxpy (E->qtot, lambda->d, tmp, 1, dx, 1);
double d = ddot(E->qtot, dx, 1, wvec, 1);
return d;
}
double Prism_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda){
static double *tmp = 0;
static double *dx = 0;
static double *dy = 0;
static double *dz = 0;
if(!tmp){
tmp = dvector(0, QGmax*QGmax*QGmax-1);
dx = dvector(0, QGmax*QGmax*QGmax-1);
dy = dvector(0, QGmax*QGmax*QGmax-1);
dz = dvector(0, QGmax*QGmax*QGmax-1);
}
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
E->Grad_h(tmp, dx, dy, dz, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (E->qtot, dx, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dy, 1, dy, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dz, 1, dz, 1, dx, 1, dx, 1);
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
if(lambda->p)
dvvtvp (E->qtot, lambda->p, 1, tmp, 1, dx, 1, dx, 1);
else
daxpy (E->qtot, lambda->d, tmp, 1, dx, 1);
double d = ddot(E->qtot, dx, 1, wvec, 1);
return d;
}
double Hex_helm_mprod(Element *E, Mode *m, double *wvec, Metric *lambda){
static double *tmp = 0;
static double *dx = 0;
static double *dy = 0;
static double *dz = 0;
if(!tmp){
tmp = dvector(0, QGmax*QGmax*QGmax-1);
dx = dvector(0, QGmax*QGmax*QGmax-1);
dy = dvector(0, QGmax*QGmax*QGmax-1);
dz = dvector(0, QGmax*QGmax*QGmax-1);
}
E->fillvec(m , tmp); //OVERLOAD CALL: fillvec: Misc.c(Tri), Misc.c(Quad), Misc.c(Tet), Misc.c(Pyr), Misc.c(Prism), Misc.c(Hex), Misc.c(Element)
E->Grad_h(tmp, dx, dy, dz, 'a'); //OVERLOAD CALL: Grad_h: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
dvmul (E->qtot, dx, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dy, 1, dy, 1, dx, 1, dx, 1);
dvvtvp(E->qtot, dz, 1, dz, 1, dx, 1, dx, 1);
dvmul (E->qtot, tmp, 1, tmp, 1, tmp, 1);
if(lambda->p)
dvvtvp (E->qtot, lambda->p, 1, tmp, 1, dx, 1, dx, 1);
else
daxpy (E->qtot, lambda->d, tmp, 1, dx, 1);
double d = ddot(E->qtot, dx, 1, wvec, 1);
return d;
}
C++ to HTML Conversion by ctoohtml