file: Hlib/src/Felisa.c


/*---------------------------------------------------------------------------*
 *                        RCS Information                                    *
 *                                                                           *
 * $Source:$ 
 * $Revision:$ 
 * $Date:$
 * $Author:$
 * $State:$
 *---------------------------------------------------------------------------*/

#include <stdio.h>
#include <hotel.h>
#include <math.h>
#include <veclib.h>
#include <felisa.h>


#define SCALTOL 1e-3
#define TOLV 1E-4

typedef struct fcurve{
  int    id;
  double u;
  struct fcurve *next;
} Fcurve;

typedef struct fsurf{
  int    id;
  double u;
  double v;
  struct fsurf *next;
} Fsurf;

typedef struct fvert{
  int      ncurve;  /* number of curves   that vertex lies in */
  int      nface;   /* number of surfaces that vertex lies in */
  double   x,y,z;   /* global coordinates                     */
  Fcurve  *fcur;    /* curve   information                    */
  Fsurf   *fsur;    /* surface information                    */
} Fvert;

typedef struct splineDat{
  int nu[2];
  double *r;  
} SplineDat;

/* local variables to file */
Fvert      *felvert;
int         nfvert, nfcurve, nfsurf;
SplineDat  *sdat, *cdat;

void Load_Felisa_Surface(char *name){
  register int i,j;
  int      nfface,v[3],err,n;
  double   len,tolg,up,vp,X[3];
  FILE     *infile;
  char     buf[BUFSIZ];
  Fcurve   *fc;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
  Fsurf    *ff;


  if(!(infile = fopen(name,"r"))){
    sprintf(buf,"%s.fro",name);
    if(!(infile = fopen(buf,"r"))){
      fprintf(stderr,"unable to open file %s or %s\n",name,buf);
      exit(-1);
    }
  }

  fgets(buf,BUFSIZ,infile);
  sscanf(buf,"%d%d%*d%*d%d%d",&nfface,&nfvert,&nfcurve,&nfsurf);

  felvert    = (Fvert *)calloc(nfvert,sizeof(Fvert));
  cdat = (SplineDat *)malloc(nfcurve*sizeof(SplineDat));
  sdat = (SplineDat *)malloc(nfsurf*sizeof(SplineDat));

  /* store global vertices as a check */
  for(i = 0; i < nfvert; ++i){
    fgets(buf,BUFSIZ,infile);
    sscanf(buf,"%*d%lf%lf%lf",&felvert[i].x,&felvert[i].y,&felvert[i].z);
  }

  /* can skip this section */
  tolg = 1e8;
  for(i = 0; i < nfface; ++i){
    fgets(buf,BUFSIZ,infile);
    /* calculate minimum distance between vertices */
    sscanf(buf,"%*d%d%d%d",v,v+1,v+2);
    v[0]--; v[1]--;v[2]--;
    for(j = 0; j < 3; ++j){
      len = sqrt(pow(felvert[v[j]].x - felvert[v[(j+1)%3]].x,2) + 
		 pow(felvert[v[j]].y - felvert[v[(j+1)%3]].y,2) + 
		 pow(felvert[v[j]].z - felvert[v[(j+1)%3]].z,2));
      tolg = min(tolg,len);
    }
  }

  settolg(SCALTOL*tolg);        //OVERLOAD CALL: settolg: felisa.h(?), felisa.h(?), felisa.h(?)

  /* assemble curve information into felvertex structure */
  for(i = 0; i < nfcurve; ++i){
    fscanf(infile,"%*d%d",&n);
    for(j = 0; j < n; ++j){
      fscanf(infile,"%d%lf",v,&up);
      v[0]--;

      fc = felvert[v[0]].fcur;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
      felvert[v[0]].fcur = (Fcurve *)malloc(sizeof(Fcurve));
      felvert[v[0]].fcur->id = i;
      felvert[v[0]].fcur->u = up;
      felvert[v[0]].fcur->next = fc;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
      felvert[v[0]].ncurve++;
    }    
  }
  
  /* assemble face information into felvertex structure */
  for(i = 0; i < nfsurf; ++i){
    fscanf(infile,"%*d%d",&n);
    for(j = 0; j < n; ++j){
      fscanf(infile,"%d%lf%lf",v,&up,&vp);
      v[0]--;

      ff = felvert[v[0]].fsur;
      felvert[v[0]].fsur = (Fsurf *)malloc(sizeof(Fsurf));
      felvert[v[0]].fsur->id = i;
      felvert[v[0]].fsur->u  = up;
      felvert[v[0]].fsur->v  = vp;
      felvert[v[0]].fsur->next = ff;
      felvert[v[0]].nface++;
    }    
  }

  /* Load Spline Data for curves  */
  for(i = 0; i < nfcurve; ++i){
    fscanf(infile,"%*d%d",cdat[i].nu);
    cdat[i].r = dvector(0,6*cdat[i].nu[0]-1);
    
    for(j = 0; j < 6*cdat[i].nu[0]; ++j)
      fscanf(infile,"%lf",cdat[i].r+j);
  }
  
  /* Load Spline Data for curves  */
  for(i = 0; i < nfsurf; ++i){
    fscanf(infile,"%*d%d%d",sdat[i].nu,sdat[i].nu+1);
    sdat[i].r = dvector(0,12*sdat[i].nu[0]*sdat[i].nu[1]-1);
    
    for(j = 0; j < 12*sdat[i].nu[0]*sdat[i].nu[1]; ++j)
      fscanf(infile,"%lf",sdat[i].r+j);
  }

  for(i = 0; i < nfvert; ++i){

    for(fc = felvert[i].fcur; fc; fc = fc->next){        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
      fgcurv(0,cdat[fc->id].nu[0],cdat[fc->id].r,fc->u,X,&err);        //OVERLOAD CALL: fgcurv: felisa.h(?), felisa.h(?), felisa.h(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
      
      if((len = sqrt(pow(X[0]-felvert[i].x,2)+pow(X[1]-felvert[i].y,2)
		     +pow(X[2]-felvert[i].z,2))) > TOLV)
	fprintf(stderr,"Felisa vertex %d is represented by a point %lf "        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
		"away in curve %d\n",i+1,len,fc->id+1);        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
    }
    
    for(ff = felvert[i].fsur; ff; ff = ff->next){
      fgsurf(0,sdat[ff->id].nu[0],sdat[ff->id].nu[1],sdat[ff->id].r,        //OVERLOAD CALL: fgsurf: felisa.h(?), felisa.h(?), felisa.h(?)
	     ff->u,ff->v,X,&err);
      
      if((len = sqrt(pow(X[0]-felvert[i].x,2)+pow(X[1]-felvert[i].y,2)
		     +pow(X[2]-felvert[i].z,2))) > TOLV)
	fprintf(stderr,"Felisa vertex %d is represented by a point %lf "        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
		"away in surface %d\n",i+1,len,ff->id+1);
    }
  }
}
       
void Free_Felisa_data(void){
  register int i;
  Fcurve *fc,*fc1;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
  Fsurf  *ff,*ff1;

  for(i = 0; i < nfvert; ++i){
    fc = felvert[i].fcur;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
    while(fc){        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
      fc1 = fc->next;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
      free(fc);        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
      fc = fc1;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
    }
    ff = felvert[i].fsur;
    while(ff){
      ff1 = ff->next;
      free(ff);
      ff = ff1;
    }
  }
  free(felvert);

  for(i = 0; i < nfcurve; ++i)
    free(cdat[i].r);
  free(cdat);

  for(i = 0; i < nfsurf; ++i)
    free(sdat[i].r);
  free(sdat);
}


void genFelFile(Element *E, double *x, double *y, double *z, Curve *curve){
  register int i,j;
  int *fvertid,eflag[Max_Nverts],face,err;
  int qa = E->qa, qb = E->qb,felsu = 0;
  double **xl,**xg,**u,**v,uc,u1,u2,v1,v2,*za,*zb,*w;
  Fcurve  *fc,*fc1;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
  Fsurf   *fs,*fs1;
  int nfv;

  if(E->identify() == Nek_Prism)        //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
    qb = E->qc;
  

  getzw(qa,&za,&w,'a');
  getzw(qb,&zb,&w,'b');
  


  u  = dmatrix(0,qb-1,0,qa-1);
  v  = dmatrix(0,qb-1,0,qa-1);
  xl = dmatrix(0,QGmax-1,0,1);
  xg = dmatrix(0,QGmax-1,0,2);
  
  face  = curve->face;
  fvertid = curve->info.file.vert;
  
  nfv = E->Nfverts(face);        //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(nfv != 3){
    fprintf(stderr, "genFelFile: not setup for quad. faces\n");
    exit(-1);
  }
#if 0
  else if(E->identify() != Nek_Tet){        //OVERLOAD CALL: identify: Hex.h(?), Prism.h(?), Pyr.h(?), Quad.h(?), Tet.h(?), Tri.h(?), nekstruct.h(?)
    fprintf(stderr, "genFelFile: only set up for tetrahedra\n");
    exit(-1);
  }
#endif
  /* check to see if vertices match */
  for(i = 0; i < nfv; ++i)
    if((fabs(felvert[fvertid[i]].x-E->vert[E->vnum(face,i)].x)>TOLV)||        //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
       (fabs(felvert[fvertid[i]].y-E->vert[E->vnum(face,i)].y)>TOLV)||        //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
       (fabs(felvert[fvertid[i]].z-E->vert[E->vnum(face,i)].z)>TOLV)){        //OVERLOAD CALL: vnum: Nums.c(Quad), Nums.c(Tet), Nums.c(Pyr), Nums.c(Prism), Nums.c(Hex), Nums.c(Element), Nums.c(Tri)
      fprintf(stderr,"Vertex %d in element %d does not match vertex %d in"
      " felisa file %s\n",i+1,E->id+1,fvertid[i]+1,curve->info.file.name);
      exit(-1);
    }

  /* first check to see if any edge is on a curve */
  for(i = 0; i < nfv; ++i){
    eflag[i] = 0;
    if(felvert[fvertid[i]].ncurve&&felvert[fvertid[(i+1)%3]].ncurve)
      for(fc = felvert[fvertid[i]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	for(fc1 = felvert[fvertid[(i+1)%3]].fcur; fc1; fc1 = fc1->next)
	  if(fc->id == fc1->id)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    eflag[i] = fc->id;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
  }
  //    if(!eflag[i])
#if 0  
    /* check to see which surface this triangle is in */
     if(!felsu){
	for(fs = felvert[fvertid[i]].fsur; fs; fs = fs->next)
	  for(fs1 = felvert[fvertid[(i+1)%3]].fsur; fs1; fs1 = fs1->next)
	    if(fs->id == fs1->id){
	      felsu = fs->id;
	      break;
	    }
    
	for(err=0,fs = felvert[fvertid[(i+2)%3]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu)
	    err=1;
	
	if(!err){
	  fprintf(stderr,"can not find felisa face for curved "
		  "face %d in element %d\n",face+1,E->id+1);
	  
	  exit(-1);
	}
	
      }	      
#endif
  
    if(!felsu){
      int *iflag = ivector(0, nfsurf-1);
      izero(nfsurf, iflag, 1);
      for(i = 0; i < nfv; ++i){
	
	for(fs = felvert[fvertid[i]].fsur; fs; fs = fs->next)
	  ++iflag[fs->id];
      }
      for(j=0;j<nfsurf;++j)
	if(iflag[j] == 3)
	  felsu = j;
      free(iflag);
    }

  // only set up for tetrahedra
  /* calculate x,y,z, points on any edge along a curve */
  for(i = 0; i < 3; ++i)
    if(eflag[i]){
      /* calculate x y z coordinates based on quadrature distribution in 
	 parametric space */
      switch(i){
      case 0:
	for(fc = felvert[fvertid[0]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	  if(fc->id == eflag[i])        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    u1 = fc->u;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	for(fc = felvert[fvertid[1]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	  if(fc->id == eflag[i])        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    u2 = fc->u;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	
	for(err=0,j = 0; j < qa; ++j){
	  uc = u1*(1-za[j])*0.5 + u2*(1+za[j])*0.5;
	  fgcurv(0,cdat[eflag[i]].nu[0],cdat[eflag[i]].r,uc,xg[j],&err);        //OVERLOAD CALL: fgcurv: felisa.h(?), felisa.h(?), felisa.h(?)
	  
	  if(err == -1){
	    fprintf(stderr,"failed to find coordinate along "
		    "curve in element %d\n",E->id+1);
	    exit(-1);
	  }
	}
        
        locuv(qa,*xg,*xl,sdat[felsu].nu[0],sdat[felsu].nu[1],sdat[felsu].r);        //OVERLOAD CALL: locuv: felisa.h(?), felisa.h(?), felisa.h(?)
	
	for(j = 0; j < qa; ++j){
	  u[0][j] = xl[j][0];
	  v[0][j] = xl[j][1];
	}
	break;
      case 1:
	/* find u value at felvertices */
	for(fc = felvert[fvertid[1]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	  if(fc->id == eflag[i])        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    u1 = fc->u;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	for(fc = felvert[fvertid[2]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	  if(fc->id == eflag[i])        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    u2 = fc->u;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	
	for(err=0,j = 0; j < qb; ++j){
	  uc = u1*(1-zb[j])*0.5 + u2*(1+zb[j])*0.5;
	  fgcurv(0,cdat[eflag[i]].nu[0],cdat[eflag[i]].r,uc,xg[j],&err);        //OVERLOAD CALL: fgcurv: felisa.h(?), felisa.h(?), felisa.h(?)

	  if(err == -1){
	    fprintf(stderr,"failed to find coordinate along "
		    "curve in element %d\n",E->id+1);
	    exit(-1);
	  }
        }
        
        locuv(qb,*xg,*xl,sdat[felsu].nu[0],sdat[felsu].nu[1],sdat[felsu].r);        //OVERLOAD CALL: locuv: felisa.h(?), felisa.h(?), felisa.h(?)
	
	for(j = 0; j < qb; ++j){
	  u[j][qa-1] = xl[j][0];
	  v[j][qa-1] = xl[j][1];
	}
	break;
      case 2:
	/* find u value at felvertices */
	for(fc = felvert[fvertid[0]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	  if(fc->id == eflag[i])        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    u1 = fc->u;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	for(fc = felvert[fvertid[2]].fcur; fc; fc = fc->next)        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?); fc: Coords.c(?), Curvi.c(?)
	  if(fc->id == eflag[i])        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	    u2 = fc->u;        //OVERLOAD CALL: fc: Coords.c(?), Curvi.c(?)
	
	for(err=0,j = 0; j < qb; ++j){
	  uc = u1*(1-zb[j])*0.5 + u2*(1+zb[j])*0.5;
	  fgcurv(0,cdat[eflag[i]].nu[0],cdat[eflag[i]].r,uc,xg[j],&err);        //OVERLOAD CALL: fgcurv: felisa.h(?), felisa.h(?), felisa.h(?)

	  if(err == -1){
	    fprintf(stderr,"failed to find coordinate along "
		    "curve in element %d\n",E->id+1);
	    exit(-1);
	  }
        }
        
        locuv(qb,*xg,*xl,sdat[felsu].nu[0],sdat[felsu].nu[1],sdat[felsu].r);        //OVERLOAD CALL: locuv: felisa.h(?), felisa.h(?), felisa.h(?)
	
	for(j = 0; j < qb; ++j){
	  u[j][0] = xl[j][0];
	  v[j][0] = xl[j][1];
	}
	break;
      }
    }
    else{ 
      /* calcuate u,v as linear fit between vertices */
      switch(i){
      case 0:
	/* find vertex data */
	for(fs = felvert[fvertid[0]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu){
	    u1 = fs->u;
	    v1 = fs->v;
	  }
	for(fs = felvert[fvertid[1]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu){
	    u2 = fs->u;
	    v2 = fs->v;
	  }
	
	for(j = 0; j < qa; ++j){
	  u[0][j] = u1*(1-za[j])*0.5 + u2*(1+za[j])*0.5;
	  v[0][j] = v1*(1-za[j])*0.5 + v2*(1+za[j])*0.5;
	}
	break;
      case 1:
	/* find felvertex data */
	for(fs = felvert[fvertid[1]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu){
	    u1 = fs->u;
	    v1 = fs->v;
	  }
	for(fs = felvert[fvertid[2]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu){
	    u2 = fs->u;
	    v2 = fs->v;
	  }
	
	for(j = 0; j < qb; ++j){
	  u[j][qa-1] = u1*(1-zb[j])*0.5 + u2*(1+zb[j])*0.5;
	  v[j][qa-1] = v1*(1-zb[j])*0.5 + v2*(1+zb[j])*0.5;
	}
	break;
      case 2:
	/* find felvertex data */
	for(fs = felvert[fvertid[0]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu){
	    u1 = fs->u;
	    v1 = fs->v;
	  }
	for(fs = felvert[fvertid[2]].fsur; fs; fs = fs->next)
	  if(fs->id == felsu){
	    u2 = fs->u;
	    v2 = fs->v;
	  }
	for(j = 0; j < qb; ++j){
	  u[j][0] = u1*(1-zb[j])*0.5 + u2*(1+zb[j])*0.5;
	  v[j][0] = v1*(1-zb[j])*0.5 + v2*(1+zb[j])*0.5;
	}
	break;
      default:
	break;
      }
  }

  /* get top vertices for blending */
  for(fs = felvert[fvertid[2]].fsur; fs; fs = fs->next)
    if(fs->id == felsu){
      u2 = fs->u;
      v2 = fs->v;
    }
  
  /* blend u,v in the interior based on edge values */
  for(i = 1; i < qa-1; ++i)
    for(j = 1; j < qb; ++j){
      u[j][i] = u[0][i]*(1-zb[j])*0.5 + u[j][0]*(1-za[i])*0.5 +
	        u[j][qa-1]*(1+za[i])*0.5 - (1-za[i])*(1-zb[j])*0.25*u[0][0] -
                (1+za[i])*(1-zb[j])*0.25*u[0][qa-1];
      v[j][i] = v[0][i]*(1-zb[j])*0.5 + v[j][0]*(1-za[i])*0.5 +
	        v[j][qa-1]*(1+za[i])*0.5 - (1-za[i])*(1-zb[j])*0.25*v[0][0] -
                (1+za[i])*(1-zb[j])*0.25*v[0][qa-1];
    }
  
  /* find x,y,z for u,v, */
  for(i = 0; i < qa; ++i)
    for(j = 0; j < qb; ++j){
      fgsurf(0,sdat[felsu].nu[0],sdat[felsu].nu[1],sdat[felsu].r,        //OVERLOAD CALL: fgsurf: felisa.h(?), felisa.h(?), felisa.h(?)
             u[j][i],v[j][i],*xg,&err);
      if(err == -1){
	fprintf(stderr,"failed to find coordinate along "
		"surface in element %d\n",E->id+1);
	exit(-1);
      }
      x[i + j*qa] = xg[0][0];
      y[i + j*qa] = xg[0][1];
      z[i + j*qa] = xg[0][2];
      
    }
  
  free_dmatrix(u,0,0);  free_dmatrix(v,0,0);
  free_dmatrix(xl,0,0); free_dmatrix(xg,0,0);
}



Back to Source File Index


C++ to HTML Conversion by ctoohtml