/*********************************************************************************
File:            construct.c
Authors:         Eija Achren, Lauri Glad, Pekka Haara, Sarianne Mykra, Tuija Rinne
Date:            11 jun 1992

Description:     Constructs a picture file in UCD file format.

Modifications:

*********************************************************************************/

#include "fem.h"
#include "defin.h"
#include <string.h>
#include <stdio.h>
#include <avs/avs.h>
#include <avs/flow.h>
#include <avs/field.h>
#include <time.h>
#include <sys/file.h>
#include <math.h>
#include <stdlib.h>
#include <errno.h>
#include <avs/ucd_defs.h>


#define WIDGET_COUNT 14	/* number of widgets e.g. datacomponents	*/
#define MAT_TYPE 1	/* material type 				*/
#define MAXNODES 20  	/* max number of nodes / element		*/
#define MAXIM(a,b) ( (a) > (b) ? (a) : (b) )

extern int error_message(char *, char *, char *, int, STime *);
void get_dataset(char *, int *, int *, int *, AllData *);
void find_minmax(int, REAL **, int, int, float *, float *, int *, float *);


/*******************************************************************************
Function:	insufficient_data

Description:	Checks if chosen data components and existing data can
                be put together to construct proper data-array for AVS, 
                if not, the module returns to ask proper data components.

Input:		*changes: contains all the other data

Output:         -

Global variables: bool disp_vec: displacement vector
                  bool disp_x  : displacement x
                  bool disp_y  : displacement y
                  bool disp_z  : displacement z
                  bool stress_vec  : stress vector
                  bool stress_x    : stress x
                  bool stress_y    : stress y
                  bool stress_z    : stress z
                  bool stress_xy   : stress xy
                  bool stress_xz   : stress xz
                  bool stress_yz   : stress yz
                  bool von_Mises   : stress von Mises
                  bool max_shear   : stress max shear
                  bool temperature : temperature

Return:		-
-----------------------------------------------------------------------------*/

void insufficient_data(changes)
 AllData *changes;
{
  char message[130];
  bool separator=FALSE;

  message[0]='\0';
  if(changes->disp.data == NULL && (disp_vec || disp_x || disp_y || disp_z))
  {
    if(separator)
      strcat(message,"and ");
    else
      separator=TRUE;
    strcat(message,"displacement data ");
  }
  if(changes->stress.data == NULL && (stress_vec || stress_x || stress_y || stress_z 
				 || stress_xy || stress_xz || stress_yz))
  {
    if(separator)
      strcat(message,"and ");
    else
      separator=TRUE;
    strcat(message,"stress data ");
  }
  if(changes->von_Mises.data == NULL && von_Mises)
  {
    if(separator)
      strcat(message,"and ");
    else
      separator=TRUE;
    strcat(message,"von Mises stress ");
  }
  if(changes->max_shear.data == NULL && max_shear)
  {
    if(separator)
      strcat(message,"and ");
    else
      separator=TRUE;
    strcat(message,"max shear stress ");
  }
  if(changes->temp.data == NULL && temperature)
  {
    if(separator)
      strcat(message,"and ");
    else
      separator=TRUE;
    strcat(message,"temperatures ");
  }
  if(separator)
    error_message("%s not found.", message, NULL, -1,NULL);
}



/************************************************************************
Function:	  get_dataset

Description:	  Reads widgets (i.e. the components that the user wants) and 
		  constructs lists for datacomponents.

Input:	       	  changes :       data vector from ABAQUS

Output:	          labels :        namelist of data components
                  num_data_comp : number of data components
                  data_comp :     numberlist of data components
                  data_veclen :   lenght of data vector

Global variables: bool disp_vec: displacement vector
                  bool disp_x  : displacement x
                  bool disp_y  : displacement y
                  bool disp_z  : displacement z
                  bool stress_vec  : stress vector
                  bool stress_x    : stress x
                  bool stress_y    : stress y
                  bool stress_z    : stress z
                  bool stress_xy   : stress xy
                  bool stress_xz   : stress xz
                  bool stress_yz   : stress yz
                  bool von_Mises   : stress von Mises
                  bool max_shear   : stress max shear
                  bool temperature : temperature

Return:		  -
-----------------------------------------------------------------------*/

void get_dataset(labels,num_data_comp,data_comp,data_veclen,changes)
char	*labels;      	/* namelist of data components   	*/
int	*num_data_comp; /* number of data components	    	*/
int 	*data_comp;	/* numberlist of data components	*/
int	*data_veclen;	/* lenght of data vector 	    	*/
AllData *changes;       /* vector containing data               */
{

        *num_data_comp=0;
	*data_veclen=0;
	labels[0]='\0';
  	if(disp_vec && (changes->disp.data != NULL))
		{
		 data_comp[*num_data_comp]=3; 	/* a vector -> 3 */
		 strcat(labels,"disp-vec;");
		 *data_veclen+=3;		/* a vector -> 3 */
		 (*num_data_comp)++;
		}
	if(disp_x && (changes->disp.data != NULL))
		{
		 data_comp[*num_data_comp]=1; 	/* a scalar -> 1 */
		 strcat(labels,"disp x;");
		 *data_veclen+=1;		/* a scalar -> 1 */
		 (*num_data_comp)++;
		}
	if(disp_y && (changes->disp.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"disp y;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(disp_z && (changes->disp.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"disp z;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_x && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"stress x;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_y && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"stress y;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_z && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"stress z;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_xy && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"stress xy;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_xz && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"stress xz;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_yz && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"stress yz;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(stress_vec && (changes->stress.data != NULL))
		{
		 data_comp[*num_data_comp]=3;
		 strcat(labels,"stress vec;");
		 *data_veclen+=3;
		 (*num_data_comp)++;
		}
	if(von_Mises && (changes->von_Mises.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"von Mises;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(max_shear && (changes->max_shear.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"max shear;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
		}
	if(temperature && (changes->temp.data != NULL))
		{
		 data_comp[*num_data_comp]=1;
		 strcat(labels,"temperature;");
		 *data_veclen+=1;
		 (*num_data_comp)++;
	        }
	/* labels[strlen(labels)-1]='\0'; 
	   A bug fixed by Sami Saarinen 18-SEP-92 */
	labels[strlen(labels)-2]='\0';
}



/******************************************************************************
Function:       find_minmax

Description:    Finds min and max values of data and constructs a data-array
                from user-defined datacomponents.

Input:          num_nodes:    number of nodes
                data:         data-array from ABAQUS
                data_veclen:  lenght of data vector
                min_index:    where the first value of a particular data
                              component is in it's data-array
                min_array,
                max_array:    list of min and max values of data components
                data_veclen_count:
                              counter for the lenght of data vector

Output:         result:       where the values of the data components are read
             
Global variables:  -

Return:            -
-----------------------------------------------------------------------------*/

static void find_minmax(num_nodes, data, data_veclen, min_index, min_array, max_array,
			data_veclen_count, result)
int num_nodes;
REAL **data;
int data_veclen;
int min_index;
float *min_array;
float *max_array;
int *data_veclen_count;
float *result;

{
  static int index;
  int i,j;

  if(*data_veclen_count == 0) index = 0;
  else index++;

  if(data_veclen == 1)
    min_array[index] = max_array[index] = (float)(data[0][min_index]);


/*---- construct data-array (result) for AVS ------------*/ 
  for (i = 0; i < num_nodes; i++)
    {
      for (j = 0; j < data_veclen; j++)
	{
	  result[num_nodes*(*data_veclen_count)+i*data_veclen+j] = (float)(data[i][min_index+j]);

/*---- find min and max values for each data component ---*/
	  if(data_veclen == 1)
	  {
	    if (result[num_nodes*(*data_veclen_count)+i*data_veclen+j] < min_array[index])
	      min_array[index] = (float)(data[i][min_index+j]);
	    if (result[num_nodes*(*data_veclen_count)+i*data_veclen+j] > max_array[index])
	      max_array[index] = (float)(data[i][min_index+j]);
	  }
	}
    }
*data_veclen_count+=data_veclen;
} /* find_minmax */
  

/**************************************************************************
Function:	  construct_UCD

Description:	  Converts data to UCD-file format.

Input:	       	  original: contains the original coordinate data
		  changes: contains all the other data

Output:           out: contains all the data for one picture

Global variables: bool disp_vec: displacement vector
                  bool disp_x  : displacement x
                  bool disp_y  : displacement y
                  bool disp_z  : displacement z
                  bool stress_vec  : stress vector
                  bool stress_x    : stress x
                  bool stress_y    : stress y
                  bool stress_z    : stress z
                  bool stress_xy   : stress xy
                  bool stress_xz   : stress xz
                  bool stress_yz   : stress yz
                  bool von_Mises   : stress von Mises
                  bool max_shear   : stress max shear
                  bool temperature : temperature

Return:		-
--------------------------------------------------------------------------*/

int construct_UCD(out,original,changes,filename)
UCD_structure **out;
Model *original;
AllData *changes;
char *filename;
{
	/*------------ local variables ---------------------------*/
	char labels[140];		/* namelist of datacomponents 		*/
	int num_data_comp;		/* number of datacomponents		*/
	int data_comp[WIDGET_COUNT];	/* numberlist of datacomponents		*/
	int data_veclen;		/* lenght of datavector			*/

	int num_cells;			/* number of cells			*/
	int num_nodes;			/* number of nodes			*/
	int num_model_data;		/* number of model data (=0)		*/
	int num_cell_data;		/* number of cell data (=0)		*/
	int num_node_data;		/* number of node data (=data_veclen)	*/
	int cell_tsize;			/* size of cell's connectivity list	*/
	int ucd_flags;			/* UCD cell names			*/
	int util_flag;			/* utility flag for general usage	*/
	int node_csize;			/* size of node's connectivity list	*/

	int cell;			/* cell number				*/
	int cell_type;			/* cell to find information		*/
	int name;			/* cell name				*/
	char *element_type;		/* name of the element type		*/
	int material_type;		/* user defined material type		*/
	int mid_edge_flags;		/* does the cell have mid-edge nodes	*/
	int *node_list;			/* array of node numbers		*/
	
	int i,j,k,count;		/* counters				*/
	int data_veclen_count;		/* counter for lenght of datavector	*/
	float *data_array;              /* data in UCD format                   */
	
	float *x,*y,*z;		/* temporary pointers to coordinate vectors	*/
	float *xarr,*yarr,*zarr;	/* scaled coordinate vectors to AVS	*/

	REAL minx,maxx,
	     miny,maxy,minz,maxz; 	/* temporary variables for scaling 	*/ 
   
	static char *cell_str[] =
		     {"pt","line","tri","quad","tet","pyr","prism","hex"};
		     			/* cell names				*/
	static int node_count[2][8]={
	              
       {1,2,3,4, 4, 5, 6, 8},
			     {1,3,6,8,10,13,15,20}};

	float min_array[20],
	      max_array[20];  	 /* arrays of min and max values for every data */ 

	float min_extent[3], max_extent[3]; /* extensions of coordinates */

	/*----------- program code---------------------------------*/
	/*---- general data from the picture ------------------*/
	data_comp[0]=0;
	insufficient_data(changes);
	for(i=0; i<20; i++)
	  min_array[i]=max_array[i]=0.0;
	get_dataset(labels,&num_data_comp,data_comp,&data_veclen,changes);
	num_cells=original->numelems;
	num_nodes=original->numnodes;
	num_model_data=0;
	num_cell_data=0;
	num_node_data=data_veclen;

	cell_tsize = node_csize = num_cells * MAXNODES;

	data_array=(float *)malloc(sizeof(float)*num_nodes*data_veclen);
	if(data_veclen == 0)
	{
	  error_message("No data components selected or selected components do not exist. Select new data components.\n",NULL, NULL,-1,NULL);
	  return ERROR;
	}
	else
	  if(data_array==NULL)
	    error_message("Memory allocation failed (data_array).",NULL,NULL,-3,NULL);

	ucd_flags = UCD_CELL_TYPES | UCD_INT | UCD_CELL_NAMES | UCD_MATERIAL_IDS | UCD_MID_EDGES;
	util_flag=0;

	/*---- allocation of structure out --------------------*/
	if (*out) UCDstructure_free(*out);
	*out=(UCD_structure *)UCDstructure_alloc(filename,num_model_data,
		ucd_flags,num_cells,cell_tsize,num_cell_data,
		num_nodes,node_csize,num_node_data,util_flag);
	if (*out == NULL)
	  error_message("Memory allocation for UCD structure failed.\n",NULL,NULL,-3,NULL);
	   

	/*---- parameters from function get_dataset -----------*/
	UCDstructure_set_node_components(*out,data_comp,num_data_comp);
	UCDstructure_set_node_labels(*out,labels,";");

	/*---- allocate room for coordinates ------------------*/
	xarr=(float *)malloc(sizeof(float)*num_nodes);
	yarr=(float *)malloc(sizeof(float)*num_nodes);
	zarr=(float *)malloc(sizeof(float)*num_nodes);
	if(xarr==NULL || yarr==NULL || zarr==NULL)
	  error_message("Memory allocation failed.",NULL,NULL,-3,NULL);

	x=xarr;y=yarr;z=zarr;

        /*---- find min and  max values for coordinates x, y and z --*/
	/*---- and set extent ----*/
	
	minx=original->node[0].x;
	maxx=original->node[0].x;
	miny=original->node[0].y;
	maxy=original->node[0].y;
	minz=original->node[0].z;
	maxz=original->node[0].z;
      
      
        for (i=1;i<num_nodes;i++)
	{
	       if(original->node[i].x < minx)
		      minx=original->node[i].x;
	       if(original->node[i].x > maxx)
		      maxx=original->node[i].x;
	       if(original->node[i].y < miny)
		      miny=original->node[i].y;
	       if(original->node[i].y > maxy)
		      maxy=original->node[i].y;
	       if(original->node[i].z < minz)
		      minz=original->node[i].z;
	       if(original->node[i].z > maxz)
		      maxz=original->node[i].z;
	}

	min_extent[0] = (float)minx; min_extent[1] = (float)miny; min_extent[2] = (float)minz;
	max_extent[0] = (float)maxx; max_extent[1] = (float)maxy; max_extent[2] = (float)maxz;
	UCDstructure_set_extent( *out, min_extent, max_extent);
       	
        /*---- read node coordinates -----------------------------------*/
	for (i=0;i<num_nodes;i++)
	{
		*x=(float)(original->node[i].x);
		*y=(float)(original->node[i].y);
		*z=(float)(original->node[i].z);
		x++;y++;z++;
	}

	UCDstructure_set_node_positions(*out,xarr,yarr,zarr);

	free(xarr);
	free(yarr);
	free(zarr);

	/*---- read element topology and other general data for elements */
	for(i=0;i<original->numelems;i++)
	{
	        int node_list[MAXNODES];

		cell = i;
		name = original->elem[i].eid; 		 /* element id 	  */
		material_type = MAT_TYPE;			 /* material type */

		if (original->elem[i].geom == GEOM_ILLEGAL)
		  {
		     error_message("The file consists one or more elements of illegal geometry.\n",NULL, NULL,-2,NULL);
		     return ERROR;
		   }
		cell_type = (int)(original->elem[i].geom)-1;
		element_type = cell_str[cell_type];
		if (original->elem[i].order == ORDER_ILLEGAL)
		  {
		    error_message("The order of element shape function is illegal for one or more elements.\n",NULL, NULL,-2,NULL);
		    return ERROR;
		  }
		if (original->elem[i].order == ORDER_CUBIC)
		  {
		    error_message("This module can not process elements with cubic order of element shape function.\n",NULL, NULL,-2,NULL);
		    return ERROR;
		  }
		mid_edge_flags = (int)(original->elem[i].order-1); /* mid-edges ?   */		

		memcpy(node_list, original->elem[i].nodelist, node_count[mid_edge_flags][cell_type]*sizeof(int));   /* copy topology      */

		for(j=0; j<node_count[mid_edge_flags][cell_type]; j++)
		  {
		    for(k=0; original->node[k].nid != node_list[j] && k<original->numnodes; k++)
		      ;
		    node_list[j]=k;
		  }
	        UCDcell_set_information(*out, cell, i, cell_str[cell_type], 1,
					cell_type, mid_edge_flags, node_list);
	}

	/*---- read actual data from structure changes --------------------------*/
	/*---- data is stored in 'data_array' -----------------------------------*/
        /*---- 'data_array' has num_nodes*num_data_comp data elements------------*/
	/*---- displacements are scaled here ------------------------------------*/

 	data_veclen_count=0;
	count=0;

	if (disp_vec && changes->disp.data != NULL)
	    find_minmax(num_nodes,changes->disp.data,3,0,min_array,
			max_array,&data_veclen_count,data_array);

 	if(disp_x && changes->disp.data != NULL)
	    find_minmax(num_nodes,changes->disp.data,1,0,min_array,
			max_array,&data_veclen_count,data_array);

	if(disp_y && changes->disp.data != NULL)
	    find_minmax(num_nodes,changes->disp.data,1,1,min_array,
			max_array,&data_veclen_count,data_array);

	if(disp_z && changes->disp.data != NULL)
	  find_minmax(num_nodes,changes->disp.data,1,2,min_array,
		      max_array,&data_veclen_count,data_array);

	if(stress_x && changes->stress.data != NULL)
	  find_minmax(num_nodes,changes->stress.data,1,0,min_array,
		      max_array,&data_veclen_count,data_array);

	if(stress_y && changes->stress.data != NULL)
	  find_minmax(num_nodes,changes->stress.data,1,1,min_array,
		      max_array,&data_veclen_count,data_array);

	if(stress_z && changes->stress.data != NULL)
	  find_minmax(num_nodes,changes->stress.data,1,2,min_array,
		      max_array,&data_veclen_count,data_array);

	if(stress_xy && changes->stress.data != NULL)
	  find_minmax(num_nodes,changes->stress.data,1,3,min_array,
		      max_array,&data_veclen_count,data_array);

	if(stress_xz && changes->stress.data != NULL)
	  find_minmax(num_nodes,changes->stress.data,1,4,min_array,
		      max_array,&data_veclen_count,data_array);

	if(stress_yz && changes->stress.data != NULL)
	  find_minmax(num_nodes,changes->stress.data,1,5,min_array,
		      max_array,&data_veclen_count,data_array);

 	if(stress_vec && changes->stress.data != NULL)
	    find_minmax(num_nodes,changes->stress.data,3,0,min_array,
			max_array,&data_veclen_count,data_array);

	if(von_Mises && changes->von_Mises.data != NULL)
	  find_minmax(num_nodes,changes->von_Mises.data,1,0,min_array,
		      max_array,&data_veclen_count,data_array);

	if(max_shear && changes->max_shear.data != NULL)
	  find_minmax(num_nodes,changes->max_shear.data,1,0,min_array,
		      max_array,&data_veclen_count,data_array);

	if(temperature && changes->temp.data != NULL)
	  find_minmax(num_nodes,changes->temp.data,1,0,min_array,
		      max_array,&data_veclen_count,data_array);

	  UCDstructure_set_node_data(*out,data_array);
	  UCDstructure_set_node_minmax(*out,min_array,max_array);

      return OK;
}














































