/* This module will connect points together in a mesh */
#include "express.h"
#include "math.h"
#include <stdio.h>
#include <stdlib.h>
#define DEBUG 0
// the tricky part is to define the subset.
// if we have 48 planes (really 48+1), and we have 481 radial slices, then
// subset ={0,48,0,481} if we read in all of the data....

int connect_zlin_mesh_quad(OMobj_id connect_zlin_mesh_quad_id, OMevent_mask event_mask, int seq_num){
  /***********************/
  /*  Declare variables  */
  /***********************/
  int  nplane;
  int  num_radial_slices;
  int  n_slice_size = 0;
  int *n_slice = NULL; 
  int  subset_size =0;
  int *subset = NULL;
  int  x1_size = 0;
  float *x1 = NULL; 
  int  x2_size = 0;
  float *x2 = NULL; 
  int  x3_size = 0;
  float *x3 = NULL; 
  int *tmp_connections_quad = NULL;
  int nt;
  int *first_point;
  int num_cells=0;
  int ishift; 
  int *point_connect;/*problem with the tri's......*/
  int status;
  int mem_quad;//, mem_tri, mem_pyr;
  float *mesh_coord;

  OMobj_id field_id,tri_id,quad_id,pyr_id,cell_id,mesh_cell_set;
  
  int i,j,k,l,m,ii,mp1,lp1;
  int unique_connections=0;
  float rdiff, r, rmin;
  int mem=0; 
  int points_used=0;
  int al_quad;//al_tri,al_pyr;
  
  if (OMget_name_int_val(connect_zlin_mesh_quad_id, OMstr_to_name("nplane"),
			 &nplane) != 1) 
    nplane = 0;    
  if (OMget_name_int_val(connect_zlin_mesh_quad_id, OMstr_to_name("num_radial_slices"),
			 &num_radial_slices) != 1) 
    num_radial_slices = 0;
  
  n_slice = (int *)OMret_name_array_ptr(connect_zlin_mesh_quad_id,
					OMstr_to_name("n_slice"), OM_GET_ARRAY_RD,
					&n_slice_size, NULL);

  subset = (int *)OMret_name_array_ptr(connect_zlin_mesh_quad_id,
					OMstr_to_name("subset"), OM_GET_ARRAY_RD,
				       &subset_size, NULL); 

  //we only care about subset[2],subset[3]... min,max

  x1 = (float *)OMret_name_array_ptr(connect_zlin_mesh_quad_id, OMstr_to_name("x1"),
				     OM_GET_ARRAY_RD, &x1_size, NULL);
   
  x2 = (float *)OMret_name_array_ptr(connect_zlin_mesh_quad_id, OMstr_to_name("x2"),
				     OM_GET_ARRAY_RD, &x2_size, NULL);
  
  x3 = (float *)OMret_name_array_ptr(connect_zlin_mesh_quad_id, OMstr_to_name("x3"),
				     OM_GET_ARRAY_RD, &x3_size, NULL);

  nt=0;
  for (i=subset[2];i<subset[3];i++) {
    nt+=n_slice[i];
  }
  point_connect = (int *)ARRalloc(NULL, DTYPE_INT,x1_size , NULL);
  first_point = (int *) ARRalloc(NULL, DTYPE_INT,num_radial_slices , NULL);

  mesh_coord = (float *) ARRalloc(NULL, DTYPE_FLOAT, 3*x1_size, NULL);
  l = 0;
  for (j=0;j<x1_size;j++) {
    mesh_coord[l] = x1[j];
    l++;
    mesh_coord[l] = x2[j];
    l++;
    mesh_coord[l] = x3[j];
    l++;
  }
  first_point[subset[2]] = 0;
  for (i=subset[2]+1;i<subset[3];i++) {
    first_point[i] = first_point[i-1] + n_slice[i-1];
  }
  ii=0;
  for (k=subset[3]-1;k>subset[2];k--) { /* which radial slice am I on? */
    for (j=0;j<n_slice[k];j++) {
      l = ii*nt + first_point[k] + j;
      // l is where I am at.
      m = ii*nt + first_point[k-1] + 0;
      rmin = (x1[l]- x1[m])*(x1[l]-x1[m])+(x2[l]-x2[m])*(x2[l]-x2[m])+(x3[l]-x3[m])*(x3[l]-x3[m]);
      point_connect[l] = m;
      for (i=0;i<n_slice[k-1];i++) { /* search on the k-1 slice */
	m = ii*nt + first_point[k-1] + i;
	rdiff = (x1[l]- x1[m])*(x1[l]-x1[m])+(x2[l]-x2[m])*(x2[l]-x2[m])+(x3[l]-x3[m])*(x3[l]-x3[m]);
	if (rdiff  <=  rmin) {
	  point_connect[l] = m;
	  rmin = rdiff;
	}
      }//i loop
    }//j loop
  } //k loop
  // }
  /* now I must find the uniques, and output this */
  al_quad = 4 * nt;
  tmp_connections_quad = (int *)ARRalloc(NULL, DTYPE_INT, al_quad, NULL); 

  mem=0;
  mem_quad=0;
  points_used = nt;
  /* form connections, get sides, etc.  */
  i=0;//we just use 1 plane.... any plane
  // for (k=subset[3]-1;k>subset[2];k--) { /* which radial slice am I on? */
  for (k=subset[2]+1;k<subset[3];k++) { /* which radial slice am I on? */
    for (j=0;j<n_slice[k]-1;j++) {
      l = i*nt + first_point[k] + j;
      lp1 = l+1;
      //m = nt+l;
      //mp1 = m + 1;
      if (i == 0 ) {
	ishift = 0;// Why do we need the shift for the 2d????
	l = l + ishift;
	lp1 = l+1;
	if (j +ishift   >= n_slice[k]) l-= n_slice[k];
	if (j +1 + ishift >= n_slice[k]) lp1-= n_slice[k];
      }
      
      
      //tmp_connections_quad[mem_quad] = l;
      //tmp_connections_quad[mem_quad+1] = lp1;
      //tmp_connections_quad[mem_quad+2] = point_connect[lp1];
      //tmp_connections_quad[mem_quad+3] = point_connect[l];
      tmp_connections_quad[mem_quad+2] = l;
      tmp_connections_quad[mem_quad+3] = lp1;
      tmp_connections_quad[mem_quad+0] = point_connect[lp1];
      tmp_connections_quad[mem_quad+1] = point_connect[l];
      mem+=4;
      mem_quad+=4;
      num_cells++;
    }
  }// k loop  
  //}

  field_id =  OMfind_subobj(connect_zlin_mesh_quad_id, OMstr_to_name("field"), OM_OBJ_RD);
  OMuser_destroy_obj(field_id);
  
  field_id =  OMcreate_obj_from_path("FLD.Field_Types.Mesh","field",connect_zlin_mesh_quad_id);
  
  OMset_int_val(OMfind_str_subobj(field_id,"nnodes",OM_OBJ_RW),points_used);
  OMset_int_val(OMfind_str_subobj(field_id,"nspace",OM_OBJ_RW),3);

  FLDadd_cell_set(field_id,"Quad");

  FLDget_cell_set(field_id,0,&mesh_cell_set); 
  /* get the first (0) cell set */
  FLDget_cell_set(field_id,0,&mesh_cell_set); 
  FLDset_ncells(mesh_cell_set,mem_quad/4); 
  FLDset_node_connect(mesh_cell_set,tmp_connections_quad,mem_quad,OM_SET_ARRAY_FREE);
  
  FLDset_coord(field_id,mesh_coord,3*points_used,OM_SET_ARRAY_FREE);
  OMlink_objs(OMfind_subobj(connect_zlin_mesh_quad_id,OMstr_to_name("link"),OM_OBJ_RW),field_id);
  
  if (n_slice != NULL) ARRfree(n_slice);
  if (x1 != NULL) ARRfree(x1);
  if (x2 != NULL) ARRfree(x2);   
  if (x3 != NULL) ARRfree(x3);
  if (first_point !=NULL) ARRfree(first_point);
  if (point_connect !=NULL) ARRfree(point_connect);
  if (subset !=NULL) ARRfree(subset);

  return(1);
}