file: Hlib/src/Tri.c


/**************************************************************************/
//                                                                        //
//   Author:    T.Warburton                                               //
//   Design:    T.Warburton && S.Sherwin                                  //
//   Date  :    12/4/96                                                   //
//                                                                        //
//   Copyright notice:  This code shall not be replicated or used without //
//                      the permission of the author.                     //
//                                                                        //
/**************************************************************************/

#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <polylib.h>
#include "veclib.h"
#include "hotel.h"
#include "nekstruct.h"


#define Tri_DIM 2


Tri::Tri(){
  Nverts = 3;
  Nedges = 3;
  Nfaces = 1;
}


Tri::Tri(Element *E){
  if(!Tri_wk)
    Tri_work();
  id      = E->id;
  type    = E->type;
  state   = 'p';
  Nverts = 3;
  Nedges = 3;
  Nfaces = 1;
  
  vert    = (Vert *)calloc(Nverts,sizeof(Vert));
  edge    = (Edge *)calloc(Nedges,sizeof(Edge));
  face    = (Face *)calloc(Nfaces,sizeof(Face));
  lmax    = E->lmax;
  interior_l       = E->interior_l;
  Nmodes  = E->Nmodes;
  Nbmodes = E->Nbmodes;
  qa      = E->qa;
  qb      = E->qb;
  qc      = 0;

  qtot    = E->qtot;

  memcpy(vert,E->vert,Nverts*sizeof(Vert));
  memcpy(edge,E->edge,Nedges*sizeof(Edge));
  memcpy(face,E->face,Nfaces*sizeof(Face));
  
  /* set memory */
  vert[0].hj = (double*)  0;
  face[0].hj = (double**) 0;
  h          = (double**) 0;
  hj_3d = (double***)0;
  h_3d  = (double***)0;

  curve  = E->curve;
  curvX  = E->curvX;
  geom   = E->geom;
  dgL    = E->dgL;

}


Tri::Tri(int i_d, char ty, int L, int Qa, int Qb, int Qc, Coord *X){     //OVERLOAD CALL: Coord: nekstruct.h(?), hotel.h(?)
  int i;

  if(!Tri_wk)
    Tri_work();

  id = i_d;
  type = ty;
  state = 'p';
  Nverts = 3;
  Nedges = 3;
  Nfaces = 1;

  vert = (Vert *)calloc(Nverts,sizeof(Vert));
  edge = (Edge *)calloc(Nedges,sizeof(Edge));
  face = (Face *)calloc(Nfaces,sizeof(Face));
  lmax = L;
  interior_l    = 0;
  Nmodes  = L*(L+1)/2;
  Nbmodes = Nmodes - (L-3)*(L-2)/2;
  qa      = Qa;
  qb      = Qb;
  qc      = Qc;

  qtot    = qa*qb;

  /* set vertex solve mask to 1 by default */
  for(i = 0; i < Nverts; ++i){
    vert[i].id    = i;
    vert[i].eid   = id;
    vert[i].solve = 1;
    vert[i].x     = X->x[i];
    vert[i].y     = X->y[i];
  }
  /* construct edge system */
  for(i = 0; i < Nedges; ++i){
    edge[i].id  = i;
    edge[i].eid = id;
    edge[i].l   = L-2;
  }  

  /* construct face system */
  for(i = 0; i < Nfaces; ++i){
    face[i].id  = i;
    face[i].eid = id;
    face[i].l   = max(0,L-3);
  }

  vert[0].hj = (double*)  0;
  face[0].hj = (double**) 0;
  h          = (double**) 0;
  hj_3d = (double***)0;
  h_3d  = (double***)0;

  curve  = (Curve*)NULL;
  curvX  = (Cmodes*)NULL;
}

Tri::Tri(int , char , int *, int *, Coord *){     //OVERLOAD CALL: Coord: nekstruct.h(?), hotel.h(?)
}


static void CheckVertLoc(Element *U, Element *E, int fac);        //OVERLOAD CALL: CheckVertLoc: Coords.c(?), Tri.c(?)
  
/* make sure that any vertex touching 'face' in element E has the 
   same 'xy' co-ordinates */

static void CheckVertLoc(Element *U, Element *E, int face){
  int     gid1,gid2,i;
  double  x1,x2,y1,y2;
  
  gid1 = E->vert[face].gid;
  x1   = E->vert[face].x;
  y1   = E->vert[face].y;
  gid2 = E->vert[(face+1)%E->Nverts].gid;
  x2   = E->vert[(face+1)%E->Nverts].x;
  y2   = E->vert[(face+1)%E->Nverts].y;

  for(;U; U=U->next)
    for(i = 0; i < U->Nverts; ++i){
      if(U->vert[i].gid == gid1){
	U->vert[i].x = x1;
	U->vert[i].y = y1;
      }
      if(U->vert[i].gid == gid2){
	U->vert[i].x = x2;
	U->vert[i].y = y2;
      }
    }
}
#if 1
static double ***Tri_iplap = (double***)0;
static int Tri_ipmodes = 0;

// Assume fixed order 
int mode_pos(Element *T, int m){
#if 1
  if(m < T->Nverts+T->edge[0].l)
    return m;
  
  if(m < T->Nverts+T->edge[0].l+T->edge[1].l)
    return m-T->edge[0].l+LGmax-2;
  
  if(m < T->Nbmodes)
    return m-T->edge[0].l-T->edge[1].l + 2*(LGmax-2);
  
  int j,n;
  for(j=0,n=T->Nbmodes+T->face[0].l;j<T->face[0].l;++j,n+=T->face[0].l-j)
    if(m < n)
      return 
	m-T->Nbmodes+T->Nverts+(LGmax-2)*T->Nedges+j*(LGmax-3-T->face[0].l);

  return -100;
#else
  return m;
#endif
}

double Tri_iprodhelm(Element *T, int i, int j, double lambda){
  double d;
  double jac = T->geom->jac.d;
  double rx  = T->geom->rx.d, sx = T->geom->sx.d;
  double ry  = T->geom->ry.d, sy = T->geom->sy.d;
  int    pa  = mode_pos(T,i);
  int    pb  = mode_pos(T,j);
  
  d  = (rx*rx+ry*ry)*Tri_iplap[0][pa][pb];
  d += (rx*sx+ry*sy)*Tri_iplap[1][pa][pb];
  d += (sx*sx+sy*sy)*Tri_iplap[2][pa][pb];
  d +=        lambda*Tri_iplap[3][pa][pb];
  d *= jac;

  return d;
}
void Tri_setup_iprodlap();

void  Tri_HelmMat(Element *T, LocMat *helm, double lambda){ 

  if(!Tri_iplap)
    Tri_setup_iprodlap();    
  
  double **a = helm->a,
         **b = helm->b,
         **c = helm->c;
  int    nbl    = T->Nbmodes; 
  int    Nmodes = T->Nmodes; 
  int    i,j;

  /* `A' matrix */
  for(i = 0; i < nbl; ++i)
    for(j = i; j < nbl; ++j)
      a[i][j] = a[j][i] = Tri_iprodhelm(T, i, j, lambda);
  
  /* `B' matrix */
  for(i = 0; i < nbl; ++i)
    for(j = nbl; j < Nmodes; ++j)
      b[i][j-nbl] = Tri_iprodhelm(T, i, j, lambda);
  
  /* 'C' matrix */
  for(i = nbl; i < Nmodes; ++i)
    for(j = nbl; j < Nmodes; ++j)
      c[i-nbl][j-nbl] = c[j-nbl][i-nbl] = Tri_iprodhelm(T, i, j, lambda);

}

void Tri_setup_iprodlap(){

  int qa = QGmax;
#ifndef COMPRESS
  int qb = QGmax-1;
#else
  int qb = QGmax;
#endif
  int qc = 0;
  int L  = LGmax;
  
  // Set up dummy element with maximum quadrature/edge order
  
  Coord X;        //OVERLOAD CALL: Coord: nekstruct.h(?), hotel.h(?)
  X.x = dvector(0,NTri_verts-1);
  X.y = dvector(0,NTri_verts-1);
  
  X.x[0] = 0.0;  X.x[1] = 1.0;  X.x[2] = 0.0;
  X.y[0] = 0.0;  X.y[1] = 0.0;  X.y[2] = 1.0;
  
  Tri *T = (Tri*) new Tri(0,'T', L, qa, qb, qc, &X);        //OVERLOAD CALL: Tri: Tri.c(Tri), Tri.c(Tri), Tri.c(Tri), Tri.c(Tri)

  free(X.x);  free(X.y);

  int i,j,k,n;
  int facs = Tri_DIM*Tri_DIM;        //OVERLOAD CALL: Tri_DIM: hotel.h(?), Tri.c(?); Tri_DIM: hotel.h(?), Tri.c(?)
  Tri_ipmodes = T->Nmodes;
  Tri_iplap = (double***) calloc(facs,sizeof(double**));
  
  for(i = 0; i < facs; ++i){
    Tri_iplap[i] = dmatrix(0, T->Nmodes-1, 0, T->Nmodes-1);
    dzero(T->Nmodes*T->Nmodes, Tri_iplap[i][0], 1);
  }

  // Set up gradient basis
  
  Basis   *B,*DB;
  Mode    *w,*m,*m1,*md,*md1,**gb,**gb1,*fac;
  double  *z;

  B      = T->getbasis();        //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
  DB     = T->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 **) malloc(T->Nmodes*sizeof(Mode *));
  gb[0]  = mvecset(0,Tri_DIM*T->Nmodes,qa, qb, qc);        //OVERLOAD CALL: Tri_DIM: hotel.h(?), Tri.c(?)
  gb1    = (Mode **) malloc(T->Nmodes*sizeof(Mode *));
  gb1[0] = mvecset(0,Tri_DIM*T->Nmodes,qa, qb, qc);        //OVERLOAD CALL: Tri_DIM: hotel.h(?), Tri.c(?)

  for(i = 1; i < T->Nmodes; ++i) gb[i]  = gb[i-1]+Tri_DIM;        //OVERLOAD CALL: Tri_DIM: hotel.h(?), Tri.c(?)
  for(i = 1; i < T->Nmodes; ++i) gb1[i] = gb1[i-1]+Tri_DIM;        //OVERLOAD CALL: Tri_DIM: hotel.h(?), Tri.c(?)

  getzw(qa,&z,&w[0].a,'a');
  getzw(qb,&z,&w[0].b,'b');
  
  /* fill gb with basis info for laplacian calculation */

  // vertex modes
  m  =  B->vert;
  md = DB->vert;
  for(i = 0,n=0; i < T->Nverts; ++i,++n)
    T->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)

  // edge modes
  for(i = 0; i < T->Nedges; ++i){
    m1  = B ->edge[i];
    md1 = DB->edge[i];
    for(j = 0; j < T->edge[i].l; ++j,++n)
      T->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)
  }
  
  // face modes
  for(i = 0; i < T->Nfaces; ++i)
    for(j = 0; j < T->face[0].l; ++j){
      m1  = B ->face[i][j];
      md1 = DB->face[i][j];
      for(k = 0; k < T->face[0].l-j; ++k,++n)
	T->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(i = 0; i < T->Nmodes; ++i){
    Tri_mvmul2d(qa,qb,qc,gb[i]  ,w,gb1[i]);
    Tri_mvmul2d(qa,qb,qc,gb[i]+1,w,gb1[i]+1);
  }
  
  // Calculate Laplacian inner products
  
  double s1, s2, s3, s4, s5, s6, s7, s8;
  double *tmp = dvector(0, QGmax-1); 
  
  fac = B->vert+1;

  for(i = 0; i < T->Nmodes; ++i)
    for(j = 0; j < T->Nmodes; ++j){

      s1  = ddot(qa,gb[i][0].a,1,gb1[j][0].a,1);
      s1 *= ddot(qb,gb[i][0].b,1,gb1[j][0].b,1);
      
      dvmul(qa, fac->a, 1, gb[i][0].a, 1, tmp, 1);

      s2  = ddot(qa,       tmp,1,gb1[j][0].a,1);
      s2 *= ddot(qb,gb[i][0].b,1,gb1[j][0].b,1);

      s3  = ddot(qa,       tmp,1,gb1[j][1].a,1);
      s3 *= ddot(qb,gb[i][0].b,1,gb1[j][1].b,1);

      dvmul(qa, fac->a, 1, tmp, 1, tmp, 1);
      
      s4  = ddot(qa,       tmp,1,gb1[j][0].a,1);
      s4 *= ddot(qb,gb[i][0].b,1,gb1[j][0].b,1);
      
      s5  = ddot(qa,gb[i][0].a,1,gb1[j][1].a,1);
      s5 *= ddot(qb,gb[i][0].b,1,gb1[j][1].b,1);
      
      s6  = ddot(qa,gb[i][1].a,1,gb1[j][0].a,1);
      s6 *= ddot(qb,gb[i][1].b,1,gb1[j][0].b,1);
      
      dvmul(qa, fac->a, 1, gb[i][1].a, 1, tmp, 1);

      s7  = ddot(qa,       tmp,1,gb1[j][0].a,1);
      s7 *= ddot(qb,gb[i][1].b,1,gb1[j][0].b,1);

      s8  = ddot(qa,gb[i][1].a,1,gb1[j][1].a,1);
      s8 *= ddot(qb,gb[i][1].b,1,gb1[j][1].b,1);
      
      Tri_iplap[0][i][j] = s1;
      Tri_iplap[1][i][j] = s5 + 2*s2 + s6;
      Tri_iplap[2][i][j] = s3 + s4 + s7 + s8;
    }

  /* fill gb with basis info for mass matrix calculation */

  // vertex modes
  m  = B->vert;
  for(i = 0,n=0; i < T->Nverts; ++i,++n){
    dcopy(qa, m[i].a, 1, gb[n]->a, 1);
    dcopy(qb, m[i].b, 1, gb[n]->b, 1);
  }

  // edge modes
  for(i = 0; i < T->Nedges; ++i){
    m1 = B ->edge[i];
    for(j = 0; j < T->edge[i].l; ++j,++n){
      dcopy(qa, m1[j].a, 1, gb[n]->a, 1);
      dcopy(qb, m1[j].b, 1, gb[n]->b, 1);
    }
  }
  
  // face modes
  for(i = 0; i < T->Nfaces; ++i)
    for(j = 0; j < T->face[i].l; ++j){
      m1  = B ->face[i][j];
      for(k = 0; k < T->face[i].l-j; ++k,++n){
	dcopy(qa, m1[k].a, 1, gb[n]->a, 1);
	dcopy(qb, m1[k].b, 1, gb[n]->b, 1);
      }
    }
  
  /* multiply by weights */
  for(i = 0; i < T->Nmodes; ++i)
    Tri_mvmul2d(qa,qb,qc,gb[i]  ,w,gb1[i]);
  
  for(i = 0; i < T->Nmodes; ++i)
    for(j = 0; j < T->Nmodes; ++j){
      
      s1  = ddot(qa,gb[i][0].a,1,gb1[j][0].a,1);
      s1 *= ddot(qb,gb[i][0].b,1,gb1[j][0].b,1);
      
      Tri_iplap[3][i][j] = s1;
    }

  free_mvec(gb[0]) ; free((char *) gb);
  free_mvec(gb1[0]); free((char *) gb1);
  free((char*)w);

  free(tmp);
}

#endif







Back to Source File Index


C++ to HTML Conversion by ctoohtml