path: root/src/AmrUcdCompute.cc

/* *****************************************/
/* Includes and Global Defs                */
/* *****************************************/
#include <stdio.h>
#include <math.h>
#include <avs/avs_math.h>

#include <avs/avs.h>
#include <avs/port.h>
#include <avs/field.h>
#include <avs/ucd_defs.h>
#include <stdlib.h>
#include <string.h>
#include <ucd_defs.h>

// #include "config.h"
#include "FlexArrayTmpl.H"
#include <IEEEIO.hh>
#include "AmrNode.hh"
// #include "AmrConvert.hh"
#include "AmrUcdFileReader.hh"
void DumpUCD(UCD_structure *u);

extern "C" 
{
  int amrucd_compute(UCD_structure **Output,char *Filename,int Hold,int Threshold,
		     int GenerateLevelData, int Step,float *Min,float *Max,
		     char *VisibleLevels);
  void amrucd_destroy();
}

void amrucd_destroy(){
  puts("Destroyed AmrUcd"); // erase UCD datastructures?
}

void GenerateFakeUCD(int Step,FlexArray<AmrNode *> &nodes,FlexArray<int> &cells){
  float data[128];
  int i,index;

  //********* Generate some fake data (nsds quads in the z direction)
  for(i=0;i<128;i++) data[i]=(float)i;
  for(index=index=i=0;i<Step+2;i++){
    int z = i;
    for(int y=0;y<2;y++){
      for(int x=0;x<2;x++,index++){
	AmrNode *node=new AmrNode;
	node->location.cartesian.x=(float)x;
	node->location.cartesian.y=(float)y;
	node->location.cartesian.z=-(float)z;
	node->data = data+index;
	node->index=index;
	nodes.append(node);
      }
    }
  }
  for(index=i=0;i<Step+1;i++){
    int connectmap[4]={0,1,3,2};
    for(int y=0;y<4;y++){
      cells.append(connectmap[y]+index);
    }
    index+=4;
    for(y=0;y<4;y++){
      cells.append(connectmap[y]+index);
    } 
  }
  printf("Ncellentries = %u /8 = %u:  nnodes=%u  /4 = %u\n",
	 cells.getSize(),cells.getSize()/8,nodes.getSize(), nodes.getSize()/4);
  /*
   for(i=0;i<nodes.getSize();i++){
     printf("node[%3u]: %f:%f:%f\n",i,
	    (nodes[i])->location.cartesian.x,
	    (nodes[i])->location.cartesian.y,
	    (nodes[i])->location.cartesian.z);
   }
   for(i=0;i<cells.getSize();i++){
     if(!(i%8)) printf("\ncell[%3u]:",i);
     printf("%2u,",cells[i]);
     //if(!((i+1)%8)) puts("");
  }
  */
  puts("");
}

int amrucd_compute (UCD_structure **Output,char *Filename,int Hold,int Threshold,
		    int GenerateLevelData, int Step,float *Min,float *Max,
		    char *VisibleLevels){
  static IObase *datafile; // Need Global Delete
  static AmrUcdFileReader *convertor; // Need Global Delete
  static FlexArray<char> VisibleLevelString;
  // must be pointers to nodes because some are possibly shared...
  FlexArray<AmrNode *> nodes; // the nodes with their data values
  FlexArray<int> cells;       // connectivity between the nodes
  //AVSfield *GenerateFakeData();
  //AVSfield *ReadHDF();
  UCD_structure *ConvertToUCD(FlexArray<AmrNode*> &,FlexArray<int> &,int);
  if(AVSparameter_changed("Filename")){
    int mintime,maxtime;
    IObase *newfile = new IEEEIO(Filename,IObase::Read);
    if(!newfile->isValid()){
      printf("AmrUcd: Error opening datafile %s\n",Filename);
      delete newfile;
      return 0;
    }
    if(datafile) delete datafile;
    datafile=newfile; // make current file
    if(convertor) delete convertor;
    convertor = new AmrUcdFileReader(*datafile);
    // convertor->printGridInfo();
    convertor->getTimeRange(mintime,maxtime);
    if((maxtime-mintime)<=1){
      AVSmodify_parameter("Step",AVS_VALUE|AVS_MINVAL|AVS_MAXVAL,mintime,mintime,maxtime+1);
      AVSparameter_visible("Step",0);
    }
    else {
      AVSmodify_parameter("Step",AVS_MINVAL|AVS_MAXVAL|AVS_VALUE,mintime,mintime,maxtime);
      AVSparameter_visible("Step",1);
    }
    Step=mintime;

    // Now setup the VisibleLevelString
    if(VisibleLevelString.getSize()>0)
      VisibleLevelString.purge();
    for(int i=0;i<convertor->nLevels();i++){
      char buffer[12];
      sprintf(buffer,"*Level %2u:",i);
      
      // Only display first 3 levels
      if(i>=3) { 
	*buffer='.';
	convertor->hideLevel(i);
     }
     
      buffer[10]='\0';
      VisibleLevelString.append(buffer,10);
    }
    VisibleLevelString[convertor->nLevels()*10 - 1]='\0'; // cap the string off
    // Now assign the string for the choice widget
    AVSmodify_parameter("VisibleLevels",AVS_MINVAL|AVS_MAXVAL|AVS_VALUE,
			"*Level  0",VisibleLevelString.getData(),":");
  }
  if(!datafile || !convertor){
    puts("No valid datafile or convertor exists.  No action taken.");
    return 0;
  }
  
  if(AVSparameter_changed("VisibleLevels")){
    int choicenum = AVSchoice_number("VisibleLevels",VisibleLevels) - 1;
    if(choicenum>=0){
      if(VisibleLevelString[choicenum*10]=='*'){ // toggle off
	VisibleLevelString[choicenum*10]='.';
	convertor->hideLevel(choicenum);
      }
      else {
	VisibleLevelString[choicenum*10]='*';
	convertor->showLevel(choicenum);
      }
      AVSmodify_parameter("VisibleLevels",AVS_VALUE|AVS_MINVAL|AVS_MAXVAL,
			  VisibleLevels,VisibleLevelString.getData(),":");
    }
  }
  if(Hold) return 0; // no change as of yet
  //  GenerateFakeUCD(Step,nodes,cells);
  if(*Output){
    puts("FREE UCD");
    DumpUCD(*Output);
    UCDstructure_free(*Output);
  }
  // *Output=NULL;
  convertor->setTime(Step);
  convertor->getUcd(nodes,cells);
  *Output=ConvertToUCD(nodes,cells,GenerateLevelData); // final gen of output
  if(!*Output) return 0;
  // Now if thresholding is there 
  if(Threshold){
    int idx;
    if(GenerateLevelData) idx=1;
    else idx=0;
    (*Output)->min_node_data[idx]=*Min;
    (*Output)->max_node_data[idx]=*Max;
    for(int n=0;n<(*Output)->node_veclen;n++){
      printf("Component[%u] min/max data=%f:%f\n",
	     n,(*Output)->min_node_data[n],(*Output)->max_node_data[n]);
    }
  }
  //for(int i=0;i<nodes.getSize();i++) {delete nodes[i]; nodes[i]=0;}
  return 1;
}

UCD_structure *ConvertToUCD(FlexArray<AmrNode*> &nodes,FlexArray<int> &cells,
			    int GenerateLevelData){
  int i,j,k;
  int util_flag,ucd_flags;/* define which components will be availible in structure */
  int num_cells,num_nodes;
  char model_name[65]; 
  int cell_tsize, node_csize;
  int num_model_data, num_cell_data, num_node_data;
  int node_comp_list[10], num_node_comp,nodelist[8],active_node_list[10];
  float min_ext[3],max_ext[3],node_data_max[10],node_data_min[10],*xc,*yc,*zc;
  int cell,node;
  UCD_structure *output;

  printf("Convert to UCD\n");
  strcpy(model_name,"AMRUCD model");
  /*********************************/
  /* establish the UCD structure   */
  /*********************************/
  /* Options: UCD_INT | UCD_MATERIAL_IDS | UCD_NODE_NAMES 
     | UCD_CELL_NAMES | UCD_CELL_TYPES */
  ucd_flags = UCD_INT;
  num_cells = cells.getSize()/8;
  num_nodes = nodes.getSize();  
  printf("So numcells=%u and numnodes=%u\n",num_cells,num_nodes);
  /* node connectivity, and only scalar nodal data */
  num_model_data = 0; 
  num_cell_data = 0; /* scalar cell data */
  util_flag = 0; /* no cell-edge connectivity */

  /* or should it be total size */
  cell_tsize =  num_cells*8; // This is tested to be correct (will warn) 
  node_csize = 0; // no node connectivity
  if(GenerateLevelData){
    num_node_data = 2; /* how much ??? */
    node_comp_list[0]=1;
    node_comp_list[1]=1;
    num_node_comp = 2;  /* no nodal data */
  }
  else{
    num_node_data = 1; /* how much ??? */
    node_comp_list[0]=1;
    node_comp_list[1]=0;
    num_node_comp = 1; 
  }

  printf("Allocate UCD structure\n");
  printf("Params:\n\tmodel_name: %s\n\t\
num_model_data: %u\n\t\
ucd_flags: %u\n\t\
num_cells: %u\n\t\
cell_tsize: %u\n\t\
num_cell_data: %u\n\t\
num_nodes: %u\n\t\
node_csize: %u\n\t\
num_node_data: %u\n\t\
util_flag: %u\n",model_name,
	 num_model_data, 
	 ucd_flags,
	 num_cells, 
	 cell_tsize,
	 num_cell_data, 
	 num_nodes, 
	 node_csize, 
	 // num_node_data, 
	 num_node_comp,
	 util_flag);

  output = (UCD_structure *)UCDstructure_alloc (model_name,
						num_model_data, 
						ucd_flags,
						num_cells, 
						cell_tsize,
						num_cell_data, 
						num_nodes, 
						node_csize, 
						num_node_data,
						util_flag);
  DumpUCD(output);
  if(!output){
    puts("AmrUcdCompute: UCD structure Allocation Failed!!!!");
    return 0;
  }
  // printf("Address of struct is %u\n",(unsigned)output);
  // should copy directly
  xc = new float[num_nodes];
  yc = new float[num_nodes];
  zc = new float[num_nodes];
  puts("Set Node connectivity and coords");
  for(i=0;i<3;i++)
    min_ext[i]=max_ext[i]=(nodes[0])->location.array[i];

  for(node=0;node < nodes.getSize();node++){
    float x,y,z;
    int clist[8]={0,0,0,0,0,0,0,0};
    AmrNode *n=nodes[node];
    xc[node]=x=n->location.cartesian.x;
    if(x<min_ext[0]) min_ext[0]=x;
    if(x>max_ext[0]) max_ext[0]=x;
    yc[node]=y=n->location.cartesian.y;
    if(y<min_ext[1]) min_ext[1]=y;
    if(y>max_ext[1]) max_ext[1]=y;
    zc[node]=z=n->location.cartesian.z;
    if(z<min_ext[2]) min_ext[2]=z;
    if(z>max_ext[2]) max_ext[2]=z;
    // was (char *) node+1
    if(!UCDnode_set_information(output,node,(char *)(node+1),0,clist)){
      puts("Error setting node information!!!!");
    }
  }
  printf("need to store %u positions\n",node);
  printf("set positions in 3-space minext=%f:%f:%f\tmaxext=%f:%f:%f\n",
	 min_ext[0],min_ext[1],min_ext[2],
	 max_ext[0],max_ext[1],max_ext[2]);
  //for(i=0;i<3;i++) min_ext[i]=max_ext[i]=0.0;
  if(!UCDstructure_set_node_positions(output,xc,yc,zc)){
    puts("Error setting node positions!!!!");
  }
  //printf("Now free the pointers");
  //for(i=0;i<nodes.getSize();i++)
  //  xc[i]=yc[i]=zc[i]=0;
  delete xc;
  delete yc;
  delete zc;
  puts("Set Extents");
  printf("Set Extents-> %u\n",
	 UCDstructure_set_extent(output,min_ext,max_ext));
  /* now cell topology */
  puts("Set cell connectivity");
  int cellidx;
  for(cell=cellidx=0;cell<cells.getSize();cell+=8,cellidx++){
    int nodelist[8];
    for(i=0;i<8;i++) nodelist[i]=cells[cell+i]; // superfluous copy
    UCDcell_set_information(output,cellidx,
			    // (char *)(cellidx+1),"hex", // was (char *)(cellidx+1)
			    (char*)cellidx,"hex",
			    1/*material ID*/,
			    UCD_HEXAHEDRON,
			    0/*no mid-edged nodes*/,
			    nodelist);
    if(cellidx>=num_cells || cell>=cells.getSize()) puts("Error!!!! Too many Cells");
  }
  puts("now set the data");
  /* Order of operations 
     for each node position
     UCDnode_set_information(*output,nodenumber,nodenumber+1,*/
  /* Now Set the node/vertex data */
  /* nodedata is in the field */
  printf("set node Data numcomp=%u compsize=%u:%u\n",
	 num_node_comp,node_comp_list[0],node_comp_list[1]);
  UCDstructure_set_node_components(output,node_comp_list,num_node_comp);

  // Just manually set stuff in the data structure
  printf("util_flag=%u\n",output->util_flag);
  printf("max_ext[2]=%u\n",output->max_extent[2]);

  active_node_list[0]=1; active_node_list[1]=0;
  UCDstructure_set_node_active(output,active_node_list);
  int dataindex;
  if(GenerateLevelData){
    UCDstructure_set_node_labels(output,"level.data",".");
    node_data_min[0]=node_data_max[0]=(float)((nodes[0])->level);
    node_data_min[1]=node_data_max[1]=(float)((nodes[0])->data[0]);
    dataindex=1;
  }
  else {
    UCDstructure_set_node_labels(output,"data",".");
    //node_data_min[0]=node_data_max[0]=(float)((nodes[0])->level);
    node_data_min[0]=node_data_max[0]=(float)((nodes[0])->data[0]);
    dataindex=0;
  }
  // KLUDGE!!! Possibly redundant malloc
  
  float *node_data = new float[nodes.getSize()*num_node_data];
  // Crash!!
  printf("Num Node Data= %u floats\n",nodes.getSize());
  int index=0;
  // Colorize by GridID??? wah?
  if(GenerateLevelData){
    for(i=0,index=0;i<nodes.getSize();i++,index++){
      float val = (float)((nodes[i])->level); // need level later!!!
      // printf("index[%u] valI=%u val=%f\n",i,nodes[i]->gridID,val);
      if(node_data_max[0]<val) node_data_max[0]=val;
      if(node_data_min[0]>val) node_data_min[0]=val;
      if(val>6) printf("index[%u] valI=%u val=%f\n",i,nodes[i]->gridID,val);
      node_data[index]=val;
    } // Crash here
    printf("node_data_max[0]=%f node_data_min[0]=%f\n",node_data_max[0],node_data_min[0]);
    puts("done setting levelinfo");
  }
  
  for(i=0;i<nodes.getSize();i++,index++){
    float *valptr = (nodes[i])->data;
    //printf("DataPtr[%u]=%lu val=%f\r",i,valptr,*valptr);
    float val = *valptr;
    if(node_data_max[dataindex]<val) node_data_max[dataindex]=val;
    if(node_data_min[dataindex]>val) node_data_min[dataindex]=val;
    node_data[index]=val;
  }
  puts("\nset node data");
  printf("UCDstructure_set_node_data->%u\n",
	 UCDstructure_set_node_data(output,node_data));

  printf("After set node data util_flag=%u\n",output->util_flag);
  printf("max_ext[2]=%u\n",output->max_extent[2]);
  
  printf("UCDstructure_set_node_minmax->%u\n",
	 UCDstructure_set_node_minmax(output,node_data_min,node_data_max));
  printf("Set Data min:max %f:%f %f:%f\n",
	 node_data_min[0],node_data_max[0],
	 node_data_min[1],node_data_max[1]);
  // for(i=0;i<nodes.getSize();i++) node_data[i]=0.0; // see if clobbering data
  delete node_data; 
  
  // Just manually set stuff in the data structure
  printf("util_flag=%u\n",output->util_flag);
  printf("max_ext[2]=%u\n",output->max_extent[0]);
  //output->util_flag=0;
  DumpUCD(output);
  return output;
}


void DumpUCD(UCD_structure *u){
  if(u->name) printf("UCD_Structure: name=%s\n",u->name);
  else printf("UCD_Structure: <no name>\n");
  printf("\tnameflag=%u util_flag=%u ncells=%u nnodes=%u\n",
	 u->name_flag,u->util_flag,
	 u->ncells,u->nnodes);
  printf("\tmin:max extent %f,%f,%f:%f,%f,%f\n",
	 u->min_extent[0], u->min_extent[1], u->min_extent[2],
	 u->max_extent[0], u->max_extent[1], u->max_extent[2]);
  printf("\tStructure data=%lu\n",(unsigned long)(u->data));
  if(u->cell_name) printf("Has Cell Names\n");
  else printf("<no cell names>\n");
  printf("\tCell: veclen=%u, node_conn_size=%u last_cell=%u\n",
	 u->cell_veclen,u->node_conn_size,u->ucd_last_cell);
  if(u->node_conn_size){
    if(u->node_list) {
      printf("\thas cell list %lu\n",u->node_list);
    }
    if(u->node_list_ptr){
      printf("\thas node list pointer %lu\n",u->node_list_ptr);
      printf("\tFirstNode %u:%u ++ %u:%u  8th node = %u:%u\n",
	     u->node_list[0],u->node_list_ptr[0],
	     u->node_list[1],u->node_list_ptr[1],
	     u->node_list[8],u->node_list_ptr[8]);
    }
  }
  if(u->node_name) printf("\tHas Node Names\n");
  else printf("\t<no node names>\n");
  printf("\tNode: veclen=%u, cell__conn_size=%u last_node=%u\n",
	 u->node_veclen,u->cell_conn_size,u->ucd_last_node);
  if(u->cell_conn_size){
    if(u->cell_list) {
      printf("\thas cell list %lu\n",u->cell_list);
    }
    if(u->cell_list_ptr){
      printf("\thas cell list pointer %lu\n",u->cell_list_ptr);
      printf("\tFirstCell %u:%u ++ %u:%u  8th cell = %u:%u\n",
	     u->cell_list[0],u->cell_list_ptr[0],
	     u->cell_list[1],u->cell_list_ptr[1],
	     u->cell_list[8],u->cell_list_ptr[8]);
    }
  }
  printf("\trefcount=%u shm_base=%lu\n",u->refcnt,(unsigned long)(u->shm_base));
}