Posted on

LBM Code

main.c

#ifndef _MAIN_C_
#define _MAIN_C_

#include <time.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#include "lbm_model.h"
#include "initialization.h"
#include "initialization_gpu.h"
#include "streaming.h"
#include "collision.h"
#include "boundary.h"
#include "lbm_solver_gpu.h"
#include "visualization.h"
#include "utils.h"

int main(int argc, char *argv[]) {
	float *collide_field=NULL, *stream_field=NULL, *collide_field_d=NULL, *stream_field_d=NULL,
			*swap=NULL, tau, wall_velocity[D_LBM], num_cells, mlups_sum;
	int *flag_field=NULL, *flag_field_d=NULL, xlength, t, timesteps, timesteps_per_plotting,
			gpu_enabled;
	clock_t mlups_time;
	size_t field_size;

	/* process parameters */
	ReadParameters(&xlength, &tau, wall_velocity, &timesteps, &timesteps_per_plotting, argc, argv,
			&gpu_enabled);

	/* check if provided parameters are legitimate */
	ValidateModel(wall_velocity, xlength, tau);

	/* initializing fields */
	num_cells = pow(xlength + 2, D_LBM);
	field_size = Q_LBM*num_cells*sizeof(float);
	collide_field = (float*) malloc(field_size);
	stream_field = (float*) malloc(field_size);
	flag_field = (int*) malloc(num_cells*sizeof(int));
	InitialiseFields(collide_field, stream_field, flag_field, xlength, gpu_enabled);
	InitialiseDeviceFields(collide_field, stream_field, flag_field, xlength, &collide_field_d, &stream_field_d, &flag_field_d);

	for (t = 0; t < timesteps; t++) {
		printf("Time step: #%d\n", t);
		if (gpu_enabled){
			DoIteration(collide_field, stream_field, flag_field, tau, wall_velocity, xlength,
					&collide_field_d, &stream_field_d, &flag_field_d, &mlups_sum);
			/* Copy data from devcice in memory only when we need VTK output */
			if (!(t%timesteps_per_plotting))
				CopyFieldsFromDevice(collide_field, stream_field, xlength,
						&collide_field_d, &stream_field_d);
		} else {
			mlups_time = clock();
			/* Copy pdfs from neighbouring cells into collide field */
			DoStreaming(collide_field, stream_field, flag_field, xlength);
			/* Perform the swapping of collide and stream fields */
			swap = collide_field; collide_field = stream_field; stream_field = swap;
			/* Compute post collision distributions */
			DoCollision(collide_field, flag_field, tau, xlength);
			/* Treat boundaries */
			TreatBoundary(collide_field, flag_field, wall_velocity, xlength);
			mlups_time = clock()-mlups_time;
			/* Print out the MLUPS value */
			mlups_sum += num_cells/(MLUPS_EXPONENT*(float)mlups_time/CLOCKS_PER_SEC);
			if(VERBOSE)
				printf("MLUPS: %f\n", num_cells/(MLUPS_EXPONENT*(float)mlups_time/CLOCKS_PER_SEC));
		}
		/* Print out vtk output if needed */
		if (!(t%timesteps_per_plotting))
			WriteVtkOutput(collide_field, flag_field, "img/lbm-img", t, xlength);
	}

	printf("Average MLUPS: %f\n", mlups_sum/(t+1));

	if (VERBOSE) {
		if(gpu_enabled)
	  	      CopyFieldsFromDevice(collide_field, stream_field, xlength, &collide_field_d, &stream_field_d);
		WriteField(collide_field, "img/collide-field", 0, xlength, gpu_enabled);
		writeFlagField(flag_field, "img/flag-field", xlength, gpu_enabled);
	}

	/* Free memory */
	free(collide_field);
	free(stream_field);
	free(flag_field);

	FreeDeviceFields(&collide_field_d, &stream_field_d, &flag_field_d);

	printf("Simulation complete.\n");
	return 0;
}

#endif

lbm_model.h

#ifndef _LBM_MODEL_H_
#define _LBM_MODEL_H_

#define VERBOSE 0
#define MLUPS_MIN_CELLS 300
#define MLUPS_EXPONENT 1000000
/* simulation parameters */
#define D_LBM 3
#define Q_LBM 19
/* arbitrary cell identifiers */
#define FLUID 0
#define NO_SLIP 1
#define MOVING_WALL 2
/* computation enhancing values */
#define EPS 0.05
#define C_S_POW2_INV 3.0
#define C_S_POW4_INV 9.0
/* reference values, not used in actual computation */
#define SQRT3 1.73205080756887729
static const float C_S = 1.0/SQRT3;
/* predefined simulation parameters */
static const int LATTICE_VELOCITIES[19][3] = {
    {0,-1,-1},{-1,0,-1},{0,0,-1},{1,0,-1},{0,1,-1},{-1,-1,0},{0,-1,0},{1,-1,0},
    {-1,0,0}, {0,0,0},  {1,0,0}, {-1,1,0},{0,1,0}, {1,1,0},  {0,-1,1},{-1,0,1},
    {0,0,1},  {1,0,1},  {0,1,1}
};
static const float LATTICE_WEIGHTS[19] = {
    1.0/36.0, 1.0/36.0, 2.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 2.0/36.0, 1.0/36.0, 
    2.0/36.0, 12.0/36.0,2.0/36.0, 1.0/36.0, 2.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 
    2.0/36.0, 1.0/36.0, 1.0/36.0
};

/* additional cell definitions */
#define LEFT_BOUNDARY 5100
#define RIGHT_BOUNDARY 5900
#define BOTTOM_BOUNDARY 5010
#define TOP_BOUNDARY 5090
#define BACK_BOUNDARY 5001
#define FRONT_BOUNDARY 5009

#define LEFT_BOTTOM_EDGE 5110
#define LEFT_UPPER_EDGE 5190
#define RIGHT_BOTTOM_EDGE 5910
#define RIGHT_UPPER_EDGE 5990

#define BACK_BOTTOM_EDGE 5011
#define BACK_UPPER_EDGE 5091
#define FRONT_BOTTOM_EDGE 5019
#define FRONT_UPPER_EDGE 5099

#define LEFT_BACK_EDGE 5101
#define LEFT_FRONT_EDGE 5109
#define RIGHT_BACK_EDGE 5901
#define RIGHT_FRONT_EDGE 5909

#define LEFT_BOTTOM_BACK_CORNER 5111
#define LEFT_BOTTOM_FRONT_CORNER 5119
#define LEFT_UPPER_BACK_CORNER 5191
#define LEFT_UPPER_FRONT_CORNER 5199
#define RIGHT_BOTTOM_BACK_CORNER 5911
#define RIGHT_BOTTOM_FRONT_CORNER 5919
#define RIGHT_UPPER_BACK_CORNER 5991
#define RIGHT_UPPER_FRONT_CORNER 5999

/**
 * Validates the configured physical model by calculating characteristic numbers
 */
void ValidateModel(float wall_velocity[D_LBM], int domain_size, float tau);

#endif

lbm_model.c

#include <math.h>
#include <stdio.h>

#include "lbm_model.h"
#include "utils.h"


void ValidateModel(float wall_velocity[D_LBM], int domain_size, float tau){
	float u_wall_length, mach_number, reynolds_number;
	/* Compute Mach number and Reynolds number */
	u_wall_length = sqrt(wall_velocity[0]*wall_velocity[0]+
			wall_velocity[1]*wall_velocity[1]+
			wall_velocity[2]*wall_velocity[2]);
	mach_number = u_wall_length*SQRT3;
	reynolds_number = u_wall_length*domain_size*C_S_POW2_INV/(tau-0.5);
	printf("Computed Mach number: %f\n", mach_number);
	printf("Computed Reynolds number: %f\n", reynolds_number);

	/* Check if characteristic numbers are correct */
	if (mach_number >= 1)
		ERROR("Computed Mach number is too large.");
	if (reynolds_number > 500)
		ERROR("Computed Reynolds number is too large for simulation to be run on a laptop/pc.");
}

 

GPU:

initialization_gpu.h

#ifndef _INITIALIZATION_GPU_H_
#define _INITIALIZATION_GPU_H_

/**
 * Copy everything from fields on DRAM in GPU
 */
void InitialiseDeviceFields(float *collide_field, float *stream_field,int *flag_field, int xlength,
		float **collide_field_d, float **stream_field_d,int **flag_field_d);

/**
 * Free GPU memory
 */
void FreeDeviceFields(float **collide_field_d, float **stream_field_d,int **flag_field_d);

/**
 * Copy data from device in global memory.
 */
void CopyFieldsFromDevice(float *collide_field, float *stream_field, int xlength,
		float **collide_field_dd, float **stream_field_dd);

#endif

initialization_gpu.cu

#include "initialization_gpu.h"
#include "lbm_model.h"
#include "utils_gpu.h"


void InitialiseDeviceFields(float *collide_field, float *stream_field,int *flag_field, int xlength,
		float **collide_field_d, float **stream_field_d,int **flag_field_d){
	int num_cells = pow(xlength+2, D_LBM);
	size_t computational_field_size = Q_LBM*num_cells*sizeof(float);
	size_t flag_field_size = num_cells*sizeof(int);

	cudaErrorCheck(cudaMalloc(collide_field_d, computational_field_size));
	cudaErrorCheck(cudaMemcpy(*collide_field_d, collide_field, computational_field_size, cudaMemcpyHostToDevice));
	cudaErrorCheck(cudaMalloc(stream_field_d, computational_field_size));
	cudaErrorCheck(cudaMemcpy(*stream_field_d, stream_field, computational_field_size, cudaMemcpyHostToDevice));
	cudaErrorCheck(cudaMalloc(flag_field_d, flag_field_size));
	cudaErrorCheck(cudaMemcpy(*flag_field_d, flag_field, flag_field_size, cudaMemcpyHostToDevice));
}


void FreeDeviceFields(float **collide_field_d, float **stream_field_d,int **flag_field_d){
	cudaErrorCheck(cudaFree(*collide_field_d));
	cudaErrorCheck(cudaFree(*stream_field_d));
	cudaErrorCheck(cudaFree(*flag_field_d));
}


void CopyFieldsFromDevice(float *collide_field, float *stream_field, int xlength,
		float **collide_field_dd, float **stream_field_dd) {
	int num_cells = pow(xlength+2, D_LBM);
	size_t computational_field_size = Q_LBM*num_cells*sizeof(float);
	cudaErrorCheck(cudaMemcpy(collide_field, *collide_field_dd, computational_field_size, cudaMemcpyDeviceToHost));
	cudaErrorCheck(cudaMemcpy(stream_field, *stream_field_dd, computational_field_size, cudaMemcpyDeviceToHost));
}

lbm_solver_gpu.h

#ifndef _COLLISION_GPU_H_
#define _COLLISION_GPU_H_

/**
 * Performs streaming, collision and boundary treatment.
 */
void DoIteration(float *collide_field, float *stream_field, int *flag_field, float tau,
		float *wall_velocity, int xlength, float **collide_field_d, float **stream_field_d,
		int **flag_field_d, float *mlups_sum);

#endif

lbm_solver_gpu.cu

#include <math.h>
#include <stdio.h>

#include "lbm_solver_gpu.h"
#include "lbm_model.h"
#include "lbm_model_gpu.cuh"
#include "utils.h"
#include "utils_gpu.h"
#include "cell_computation_gpu.cuh"


__constant__ float tau_d, wall_velocity_d[D_LBM];
__constant__ int xlength_d, num_cells_d;
__device__ float *stream_field_d, *collide_field_d;


/**
 * Computes the post-collision distribution functions according to the BGK update rule and
 * stores the results again at the same position.
 */
__device__ void ComputePostCollisionDistributionsGpu(float *current_cell, float *feq){
	current_cell[0]=current_cell[0]-(current_cell[0]-feq[0])/tau_d;
	current_cell[1]=current_cell[1]-(current_cell[1]-feq[1])/tau_d;
	current_cell[2]=current_cell[2]-(current_cell[2]-feq[2])/tau_d;
	current_cell[3]=current_cell[3]-(current_cell[3]-feq[3])/tau_d;
	current_cell[4]=current_cell[4]-(current_cell[4]-feq[4])/tau_d;
	current_cell[5]=current_cell[5]-(current_cell[5]-feq[5])/tau_d;
	current_cell[6]=current_cell[6]-(current_cell[6]-feq[6])/tau_d;
	current_cell[7]=current_cell[7]-(current_cell[7]-feq[7])/tau_d;
	current_cell[8]=current_cell[8]-(current_cell[8]-feq[8])/tau_d;
	current_cell[9]=current_cell[9]-(current_cell[9]-feq[9])/tau_d;
	current_cell[10]=current_cell[10]-(current_cell[10]-feq[10])/tau_d;
	current_cell[11]=current_cell[11]-(current_cell[11]-feq[11])/tau_d;
	current_cell[12]=current_cell[12]-(current_cell[12]-feq[12])/tau_d;
	current_cell[13]=current_cell[13]-(current_cell[13]-feq[13])/tau_d;
	current_cell[14]=current_cell[14]-(current_cell[14]-feq[14])/tau_d;
	current_cell[15]=current_cell[15]-(current_cell[15]-feq[15])/tau_d;
	current_cell[16]=current_cell[16]-(current_cell[16]-feq[16])/tau_d;
	current_cell[17]=current_cell[17]-(current_cell[17]-feq[17])/tau_d;
	current_cell[18]=current_cell[18]-(current_cell[18]-feq[18])/tau_d;
}

/*
 * Performs streaming on cells.
 */
__global__ void DoStreaming(){
	int x = threadIdx.x+blockIdx.x*blockDim.x;
	int y = threadIdx.y+blockIdx.y*blockDim.y;
	int z = threadIdx.z+blockIdx.z*blockDim.z;
	int step=xlength_d+2, idx=x+y*step+z*step*step, nx, ny, nz;

	/* check that indices are within the bounds since there could be more threads than needed */
	if (0<x && x<(step-1) && 0<y && y<(step-1) && 0<z && z<(step-1)){
		nx=x-LATTICE_VELOCITIES_D[0][0];
		ny=y-LATTICE_VELOCITIES_D[0][1];
		nz=z-LATTICE_VELOCITIES_D[0][2];
		stream_field_d[Q_LBM*idx]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)];
		nx=x-LATTICE_VELOCITIES_D[1][0];
		ny=y-LATTICE_VELOCITIES_D[1][1];
		nz=z-LATTICE_VELOCITIES_D[1][2];
		stream_field_d[Q_LBM*idx+1]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+1];
		nx=x-LATTICE_VELOCITIES_D[2][0];
		ny=y-LATTICE_VELOCITIES_D[2][1];
		nz=z-LATTICE_VELOCITIES_D[2][2];
		stream_field_d[Q_LBM*idx+2]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+2];
		nx=x-LATTICE_VELOCITIES_D[3][0];
		ny=y-LATTICE_VELOCITIES_D[3][1];
		nz=z-LATTICE_VELOCITIES_D[3][2];
		stream_field_d[Q_LBM*idx+3]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+3];
		nx=x-LATTICE_VELOCITIES_D[4][0];
		ny=y-LATTICE_VELOCITIES_D[4][1];
		nz=z-LATTICE_VELOCITIES_D[4][2];
		stream_field_d[Q_LBM*idx+4]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+4];
		nx=x-LATTICE_VELOCITIES_D[5][0];
		ny=y-LATTICE_VELOCITIES_D[5][1];
		nz=z-LATTICE_VELOCITIES_D[5][2];
		stream_field_d[Q_LBM*idx+5]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+5];
		nx=x-LATTICE_VELOCITIES_D[6][0];
		ny=y-LATTICE_VELOCITIES_D[6][1];
		nz=z-LATTICE_VELOCITIES_D[6][2];
		stream_field_d[Q_LBM*idx+6]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+6];
		nx=x-LATTICE_VELOCITIES_D[7][0];
		ny=y-LATTICE_VELOCITIES_D[7][1];
		nz=z-LATTICE_VELOCITIES_D[7][2];
		stream_field_d[Q_LBM*idx+7]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+7];
		nx=x-LATTICE_VELOCITIES_D[8][0];
		ny=y-LATTICE_VELOCITIES_D[8][1];
		nz=z-LATTICE_VELOCITIES_D[8][2];
		stream_field_d[Q_LBM*idx+8]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+8];
		nx=x-LATTICE_VELOCITIES_D[9][0];
		ny=y-LATTICE_VELOCITIES_D[9][1];
		nz=z-LATTICE_VELOCITIES_D[9][2];
		stream_field_d[Q_LBM*idx+9]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+9];
		nx=x-LATTICE_VELOCITIES_D[10][0];
		ny=y-LATTICE_VELOCITIES_D[10][1];
		nz=z-LATTICE_VELOCITIES_D[10][2];
		stream_field_d[Q_LBM*idx+10]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+10];
		nx=x-LATTICE_VELOCITIES_D[11][0];
		ny=y-LATTICE_VELOCITIES_D[11][1];
		nz=z-LATTICE_VELOCITIES_D[11][2];
		stream_field_d[Q_LBM*idx+11]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+11];
		nx=x-LATTICE_VELOCITIES_D[12][0];
		ny=y-LATTICE_VELOCITIES_D[12][1];
		nz=z-LATTICE_VELOCITIES_D[12][2];
		stream_field_d[Q_LBM*idx+12]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+12];
		nx=x-LATTICE_VELOCITIES_D[13][0];
		ny=y-LATTICE_VELOCITIES_D[13][1];
		nz=z-LATTICE_VELOCITIES_D[13][2];
		stream_field_d[Q_LBM*idx+13]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+13];
		nx=x-LATTICE_VELOCITIES_D[14][0];
		ny=y-LATTICE_VELOCITIES_D[14][1];
		nz=z-LATTICE_VELOCITIES_D[14][2];
		stream_field_d[Q_LBM*idx+14]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+14];
		nx=x-LATTICE_VELOCITIES_D[15][0];
		ny=y-LATTICE_VELOCITIES_D[15][1];
		nz=z-LATTICE_VELOCITIES_D[15][2];
		stream_field_d[Q_LBM*idx+15]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+15];
		nx=x-LATTICE_VELOCITIES_D[16][0];
		ny=y-LATTICE_VELOCITIES_D[16][1];
		nz=z-LATTICE_VELOCITIES_D[16][2];
		stream_field_d[Q_LBM*idx+16]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+16];
		nx=x-LATTICE_VELOCITIES_D[17][0];
		ny=y-LATTICE_VELOCITIES_D[17][1];
		nz=z-LATTICE_VELOCITIES_D[17][2];
		stream_field_d[Q_LBM*idx+17]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+17];
		nx=x-LATTICE_VELOCITIES_D[18][0];
		ny=y-LATTICE_VELOCITIES_D[18][1];
		nz=z-LATTICE_VELOCITIES_D[18][2];
		stream_field_d[Q_LBM*idx+18]=collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+18];
	}
}


/*
 * Performs collision computation.
 */
__global__ void DoCollision(){
	int x = threadIdx.x+blockIdx.x*blockDim.x;
	int y = threadIdx.y+blockIdx.y*blockDim.y;
	int z = threadIdx.z+blockIdx.z*blockDim.z;
	int step=xlength_d+2, idx=x+y*step+z*step*step;
	float density, velocity[D_LBM], feq[Q_LBM], *current_cell_s;

	/* check that indices are within the bounds since there could be more threads than needed */
	if (0<x && x<(step-1) && 0<y && y<(step-1) && 0<z && z<(step-1)){
		current_cell_s=&collide_field_d[Q_LBM*idx];
		ComputeDensityGpu(current_cell_s,&density);
		ComputeVelocityGpu(current_cell_s,&density,velocity);
		ComputeFeqGpu(&density,velocity,feq);
		ComputePostCollisionDistributionsGpu(current_cell_s,feq);
	}
}


/*
 * Computes proper boundary values.
 */
__global__ void TreatBoundary(int *flag_field_d){
	int x = threadIdx.x+blockIdx.x*blockDim.x;
	int y = threadIdx.y+blockIdx.y*blockDim.y;
	int z = threadIdx.z+blockIdx.z*blockDim.z;

	int step=xlength_d+2, idx=x+y*step+z*step*step, nx, ny, nz, i, boundary_side=0, boundary_idx=100500;
	float density, dot_prod;

	if(idx<num_cells_d) {
		if(flag_field_d[idx] == BOTTOM_BOUNDARY) {
			boundary_side = BOTTOM_BOUNDARY;
			boundary_idx = BOTTOM_BOUNDARY_IDX;
		} else if (flag_field_d[idx] == LEFT_BOUNDARY) {
			boundary_side = LEFT_BOUNDARY;
			boundary_idx = LEFT_BOUNDARY_IDX;
		} else if (flag_field_d[idx] == RIGHT_BOUNDARY) {
			boundary_side = RIGHT_BOUNDARY;
			boundary_idx = RIGHT_BOUNDARY_IDX;
		} else if (flag_field_d[idx] == BACK_BOUNDARY) {
			boundary_side = BACK_BOUNDARY;
			boundary_idx = BACK_BOUNDARY_IDX;
		} else if (flag_field_d[idx] == FRONT_BOUNDARY) {
			boundary_side = FRONT_BOUNDARY;
			boundary_idx = FRONT_BOUNDARY_IDX;
		} else if (flag_field_d[idx] == LEFT_BOTTOM_EDGE) {
			boundary_side = LEFT_BOTTOM_EDGE;
			boundary_idx = 13;
		} else if (flag_field_d[idx] == RIGHT_BOTTOM_EDGE) {
			boundary_side = RIGHT_BOTTOM_EDGE;
			boundary_idx = 11;
		} else if (flag_field_d[idx] == BACK_BOTTOM_EDGE) {
			boundary_side = BACK_BOTTOM_EDGE;
			boundary_idx = 18;
		} else if (flag_field_d[idx] == FRONT_BOTTOM_EDGE) {
			boundary_side = FRONT_BOTTOM_EDGE;
			boundary_idx = 4;
		} else if (flag_field_d[idx] == LEFT_BACK_EDGE) {
			boundary_side = LEFT_BACK_EDGE;
			boundary_idx = 17;
		} else if (flag_field_d[idx] == LEFT_FRONT_EDGE) {
			boundary_side = LEFT_FRONT_EDGE;
			boundary_idx = 3;
		} else if (flag_field_d[idx] == RIGHT_BACK_EDGE) {
			boundary_side = RIGHT_BACK_EDGE;
			boundary_idx = 15;
		} else if (flag_field_d[idx] == RIGHT_FRONT_EDGE) {
			boundary_side = RIGHT_FRONT_EDGE;
			boundary_idx = 1;
		} else if (flag_field_d[idx] == LEFT_UPPER_EDGE) {
			boundary_side = LEFT_UPPER_EDGE;
			boundary_idx = 7;
		} else if (flag_field_d[idx] == RIGHT_UPPER_EDGE) {
			boundary_side = RIGHT_UPPER_EDGE;
			boundary_idx = 5;
		} else if (flag_field_d[idx] == BACK_UPPER_EDGE) {
			boundary_side = BACK_UPPER_EDGE;
			boundary_idx = 14;
		} else if (flag_field_d[idx] == FRONT_UPPER_EDGE) {
			boundary_side = FRONT_UPPER_EDGE;
			boundary_idx = 0;
		} else if (flag_field_d[idx] == TOP_BOUNDARY) {
			boundary_side = TOP_BOUNDARY;
			boundary_idx = TOP_BOUNDARY_IDX;
		}

		if( boundary_side==LEFT_BOUNDARY || boundary_side==RIGHT_BOUNDARY ||
				boundary_side==BOTTOM_BOUNDARY ||
				boundary_side==BACK_BOUNDARY || boundary_side==FRONT_BOUNDARY) {
			i = TREAT_BOUNDARY_INDECES[boundary_idx][0];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			collide_field_d[Q_LBM*(x+y*step+z*step*step)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)];
			i = TREAT_BOUNDARY_INDECES[boundary_idx][1];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			collide_field_d[Q_LBM*(x+y*step+z*step*step)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)];
			i = TREAT_BOUNDARY_INDECES[boundary_idx][2];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			collide_field_d[Q_LBM*(x+y*step+z*step*step)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)];
			i = TREAT_BOUNDARY_INDECES[boundary_idx][3];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			collide_field_d[Q_LBM*(x+y*step+z*step*step)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)];
			i = TREAT_BOUNDARY_INDECES[boundary_idx][4];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			collide_field_d[Q_LBM*(x+y*step+z*step*step)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)];
		} else if (boundary_side == LEFT_BOTTOM_EDGE || boundary_side == RIGHT_BOTTOM_EDGE ||
				boundary_side == BACK_BOTTOM_EDGE || boundary_side == FRONT_BOTTOM_EDGE ||
				boundary_side == LEFT_BACK_EDGE || boundary_side == LEFT_FRONT_EDGE ||
				boundary_side == RIGHT_BACK_EDGE || boundary_side == RIGHT_FRONT_EDGE) {
			nx=x+LATTICE_VELOCITIES_D[boundary_idx][0];
			ny=y+LATTICE_VELOCITIES_D[boundary_idx][1];
			nz=z+LATTICE_VELOCITIES_D[boundary_idx][2];
			collide_field_d[Q_LBM*(x+y*step+z*step*step)+boundary_idx]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(boundary_idx)];
		} else if(boundary_side == LEFT_UPPER_EDGE || boundary_side == RIGHT_UPPER_EDGE ||
				boundary_side == BACK_UPPER_EDGE || boundary_side == FRONT_UPPER_EDGE) {
			i = boundary_idx;
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			/* Compute density in the neighbour cell */
			ComputeDensityGpu(&collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)],&density);
			/* Compute dot product */
			dot_prod=LATTICE_VELOCITIES_D[i][0]*wall_velocity_d[0]+
					LATTICE_VELOCITIES_D[i][1]*wall_velocity_d[1]+
					LATTICE_VELOCITIES_D[i][2]*wall_velocity_d[2];
			/* Assign the boudary cell value */
			collide_field_d[Q_LBM*(idx)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)]+
					2*LATTICE_WEIGHTS_D[i]*density*C_S_POW2_INV*dot_prod;
		} else if(boundary_side == TOP_BOUNDARY) {
			i = TREAT_BOUNDARY_INDECES[boundary_idx][0];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			ComputeDensityGpu(&collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)],&density);
			dot_prod=LATTICE_VELOCITIES_D[i][0]*wall_velocity_d[0]+
					LATTICE_VELOCITIES_D[i][1]*wall_velocity_d[1]+
					LATTICE_VELOCITIES_D[i][2]*wall_velocity_d[2];
			collide_field_d[Q_LBM*(idx)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)]+
					2*LATTICE_WEIGHTS_D[i]*density*C_S_POW2_INV*dot_prod;
			i = TREAT_BOUNDARY_INDECES[boundary_idx][1];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			ComputeDensityGpu(&collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)],&density);
			dot_prod=LATTICE_VELOCITIES_D[i][0]*wall_velocity_d[0]+
					LATTICE_VELOCITIES_D[i][1]*wall_velocity_d[1]+
					LATTICE_VELOCITIES_D[i][2]*wall_velocity_d[2];
			collide_field_d[Q_LBM*(idx)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)]+
					2*LATTICE_WEIGHTS_D[i]*density*C_S_POW2_INV*dot_prod;
			i = TREAT_BOUNDARY_INDECES[boundary_idx][2];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			ComputeDensityGpu(&collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)],&density);
			dot_prod=LATTICE_VELOCITIES_D[i][0]*wall_velocity_d[0]+
					LATTICE_VELOCITIES_D[i][1]*wall_velocity_d[1]+
					LATTICE_VELOCITIES_D[i][2]*wall_velocity_d[2];
			collide_field_d[Q_LBM*(idx)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)]+
					2*LATTICE_WEIGHTS_D[i]*density*C_S_POW2_INV*dot_prod;
			i = TREAT_BOUNDARY_INDECES[boundary_idx][3];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			ComputeDensityGpu(&collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)],&density);
			dot_prod=LATTICE_VELOCITIES_D[i][0]*wall_velocity_d[0]+
					LATTICE_VELOCITIES_D[i][1]*wall_velocity_d[1]+
					LATTICE_VELOCITIES_D[i][2]*wall_velocity_d[2];
			collide_field_d[Q_LBM*(idx)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)]+
					2*LATTICE_WEIGHTS_D[i]*density*C_S_POW2_INV*dot_prod;
			i = TREAT_BOUNDARY_INDECES[boundary_idx][4];
			nx=x+LATTICE_VELOCITIES_D[i][0];
			ny=y+LATTICE_VELOCITIES_D[i][1];
			nz=z+LATTICE_VELOCITIES_D[i][2];
			ComputeDensityGpu(&collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)],&density);
			dot_prod=LATTICE_VELOCITIES_D[i][0]*wall_velocity_d[0]+
					LATTICE_VELOCITIES_D[i][1]*wall_velocity_d[1]+
					LATTICE_VELOCITIES_D[i][2]*wall_velocity_d[2];
			collide_field_d[Q_LBM*(idx)+i]=
					collide_field_d[Q_LBM*(nx+ny*step+nz*step*step)+InvGpu(i)]+
					2*LATTICE_WEIGHTS_D[i]*density*C_S_POW2_INV*dot_prod;
		}
	}
}


/*
 * Performs pointer swap on GPU.
 */
__global__ void DoSwap(){
	float *swap=collide_field_d; collide_field_d=stream_field_d; stream_field_d=swap;
}


void DoIteration(float *collide_field, float *stream_field, int *flag_field, float tau,
		float *wall_velocity, int xlength, float **collide_field_dd, float **stream_field_dd,
		int **flag_field_d, float *mlups_sum){
	int num_cells = pow(xlength+2, D_LBM);
	clock_t mlups_time;

	/* initialize constant data */
	cudaErrorCheck(cudaMemcpyToSymbol(xlength_d, &xlength, sizeof(int), 0, cudaMemcpyHostToDevice));
	cudaErrorCheck(cudaMemcpyToSymbol(num_cells_d, &num_cells, sizeof(int), 0, cudaMemcpyHostToDevice));
	cudaErrorCheck(cudaMemcpyToSymbol(tau_d, &tau, sizeof(float), 0, cudaMemcpyHostToDevice));
	cudaErrorCheck(cudaMemcpyToSymbol(wall_velocity_d, wall_velocity, D_LBM*sizeof(float), 0, cudaMemcpyHostToDevice));

	cudaErrorCheck(cudaMemcpyToSymbol(collide_field_d, collide_field_dd, sizeof(*collide_field_dd), 0, cudaMemcpyHostToDevice));
	cudaErrorCheck(cudaMemcpyToSymbol(stream_field_d, stream_field_dd, sizeof(*stream_field_dd), 0, cudaMemcpyHostToDevice));

	/* define grid structure */
	dim3 block(BLOCK_SIZE, BLOCK_SIZE, BLOCK_SIZE);
	dim3 grid((xlength+2+block.x-1)/block.x, (xlength+2+block.y-1)/block.y, (xlength+2+block.z-1)/block.z);

	mlups_time = clock();

	/* perform streaming */
	DoStreaming<<<grid,block>>>();
	cudaErrorCheck(cudaThreadSynchronize());
	cudaErrorCheck(cudaPeekAtLastError());

	/* Perform the swapping of collide and stream fields */
	DoSwap<<<1,1>>>();
	cudaErrorCheck(cudaThreadSynchronize());
	cudaErrorCheck(cudaPeekAtLastError());

	/* perform collision */
	DoCollision<<<grid,block>>>();
	cudaErrorCheck(cudaThreadSynchronize());
	cudaErrorCheck(cudaPeekAtLastError());

	/* perform boundary treatment */
	TreatBoundary<<<grid,block>>>(*flag_field_d);
	cudaErrorCheck(cudaPeekAtLastError());

	mlups_time = clock()-mlups_time;

	*mlups_sum += num_cells/(MLUPS_EXPONENT*(float)mlups_time/CLOCKS_PER_SEC);
	if(VERBOSE)
		printf("MLUPS: %f\n", num_cells/(MLUPS_EXPONENT*(float)mlups_time/CLOCKS_PER_SEC));

	/* copy data back to host */
	cudaErrorCheck(cudaMemcpyFromSymbol(collide_field_dd, collide_field_d, sizeof(*collide_field_dd), 0, cudaMemcpyDeviceToHost));
	cudaErrorCheck(cudaMemcpyFromSymbol(stream_field_dd, stream_field_d, sizeof(*stream_field_dd), 0, cudaMemcpyDeviceToHost));
}

 

cell_computation_gpu.h

#ifndef _CELL_COMPUTATION_GPU_CUH_
#define _CELL_COMPUTATION_GPU_CUH_


/**
 * Computes the density from the particle distribution functions stored at currentCell.
 * currentCell thus denotes the address of the first particle distribution function of the
 * respective cell. The result is stored in density.
 */
__device__ void ComputeDensityGpu(float *current_cell, float *density){
    *density=0;
    *density+=current_cell[0];
    *density+=current_cell[1];
    *density+=current_cell[2];
    *density+=current_cell[3];
    *density+=current_cell[4];
    *density+=current_cell[5];
    *density+=current_cell[6];
    *density+=current_cell[7];
    *density+=current_cell[8];
    *density+=current_cell[9];
    *density+=current_cell[10];
    *density+=current_cell[11];
    *density+=current_cell[12];
    *density+=current_cell[13];
    *density+=current_cell[14];
    *density+=current_cell[15];
    *density+=current_cell[16];
    *density+=current_cell[17];
    *density+=current_cell[18];
}


/**
 * Computes the velocity within currentCell and stores the result in velocity
 */
__device__ void ComputeVelocityGpu(float *current_cell, float *density, float *velocity){
    velocity[0]=0;
	velocity[0]+=current_cell[0]*LATTICE_VELOCITIES_D[0][0];
	velocity[0]+=current_cell[1]*LATTICE_VELOCITIES_D[1][0];
	velocity[0]+=current_cell[2]*LATTICE_VELOCITIES_D[2][0];
	velocity[0]+=current_cell[3]*LATTICE_VELOCITIES_D[3][0];
	velocity[0]+=current_cell[4]*LATTICE_VELOCITIES_D[4][0];
	velocity[0]+=current_cell[5]*LATTICE_VELOCITIES_D[5][0];
	velocity[0]+=current_cell[6]*LATTICE_VELOCITIES_D[6][0];
	velocity[0]+=current_cell[7]*LATTICE_VELOCITIES_D[7][0];
	velocity[0]+=current_cell[8]*LATTICE_VELOCITIES_D[8][0];
	velocity[0]+=current_cell[9]*LATTICE_VELOCITIES_D[9][0];
	velocity[0]+=current_cell[10]*LATTICE_VELOCITIES_D[10][0];
	velocity[0]+=current_cell[11]*LATTICE_VELOCITIES_D[11][0];
	velocity[0]+=current_cell[12]*LATTICE_VELOCITIES_D[12][0];
	velocity[0]+=current_cell[13]*LATTICE_VELOCITIES_D[13][0];
	velocity[0]+=current_cell[14]*LATTICE_VELOCITIES_D[14][0];
	velocity[0]+=current_cell[15]*LATTICE_VELOCITIES_D[15][0];
	velocity[0]+=current_cell[16]*LATTICE_VELOCITIES_D[16][0];
	velocity[0]+=current_cell[17]*LATTICE_VELOCITIES_D[17][0];
	velocity[0]+=current_cell[18]*LATTICE_VELOCITIES_D[18][0];
	velocity[0]/=*density;

	velocity[1]=0;
	velocity[1]+=current_cell[0]*LATTICE_VELOCITIES_D[0][1];
	velocity[1]+=current_cell[1]*LATTICE_VELOCITIES_D[1][1];
	velocity[1]+=current_cell[2]*LATTICE_VELOCITIES_D[2][1];
	velocity[1]+=current_cell[3]*LATTICE_VELOCITIES_D[3][1];
	velocity[1]+=current_cell[4]*LATTICE_VELOCITIES_D[4][1];
	velocity[1]+=current_cell[5]*LATTICE_VELOCITIES_D[5][1];
	velocity[1]+=current_cell[6]*LATTICE_VELOCITIES_D[6][1];
	velocity[1]+=current_cell[7]*LATTICE_VELOCITIES_D[7][1];
	velocity[1]+=current_cell[8]*LATTICE_VELOCITIES_D[8][1];
	velocity[1]+=current_cell[9]*LATTICE_VELOCITIES_D[9][1];
	velocity[1]+=current_cell[10]*LATTICE_VELOCITIES_D[10][1];
	velocity[1]+=current_cell[11]*LATTICE_VELOCITIES_D[11][1];
	velocity[1]+=current_cell[12]*LATTICE_VELOCITIES_D[12][1];
	velocity[1]+=current_cell[13]*LATTICE_VELOCITIES_D[13][1];
	velocity[1]+=current_cell[14]*LATTICE_VELOCITIES_D[14][1];
	velocity[1]+=current_cell[15]*LATTICE_VELOCITIES_D[15][1];
	velocity[1]+=current_cell[16]*LATTICE_VELOCITIES_D[16][1];
	velocity[1]+=current_cell[17]*LATTICE_VELOCITIES_D[17][1];
	velocity[1]+=current_cell[18]*LATTICE_VELOCITIES_D[18][1];
	velocity[1]/=*density;

	velocity[2]=0;
	velocity[2]+=current_cell[0]*LATTICE_VELOCITIES_D[0][2];
	velocity[2]+=current_cell[1]*LATTICE_VELOCITIES_D[1][2];
	velocity[2]+=current_cell[2]*LATTICE_VELOCITIES_D[2][2];
	velocity[2]+=current_cell[3]*LATTICE_VELOCITIES_D[3][2];
	velocity[2]+=current_cell[4]*LATTICE_VELOCITIES_D[4][2];
	velocity[2]+=current_cell[5]*LATTICE_VELOCITIES_D[5][2];
	velocity[2]+=current_cell[6]*LATTICE_VELOCITIES_D[6][2];
	velocity[2]+=current_cell[7]*LATTICE_VELOCITIES_D[7][2];
	velocity[2]+=current_cell[8]*LATTICE_VELOCITIES_D[8][2];
	velocity[2]+=current_cell[9]*LATTICE_VELOCITIES_D[9][2];
	velocity[2]+=current_cell[10]*LATTICE_VELOCITIES_D[10][2];
	velocity[2]+=current_cell[11]*LATTICE_VELOCITIES_D[11][2];
	velocity[2]+=current_cell[12]*LATTICE_VELOCITIES_D[12][2];
	velocity[2]+=current_cell[13]*LATTICE_VELOCITIES_D[13][2];
	velocity[2]+=current_cell[14]*LATTICE_VELOCITIES_D[14][2];
	velocity[2]+=current_cell[15]*LATTICE_VELOCITIES_D[15][2];
	velocity[2]+=current_cell[16]*LATTICE_VELOCITIES_D[16][2];
	velocity[2]+=current_cell[17]*LATTICE_VELOCITIES_D[17][2];
	velocity[2]+=current_cell[18]*LATTICE_VELOCITIES_D[18][2];
    velocity[2]/=*density;
}


/**
 * Computes the equilibrium distributions for all particle distribution functions of one
 * cell from density and velocity and stores the results in feq.
 */
__device__ void ComputeFeqGpu(float *density, float *velocity, float *feq){
    float s1, s2, s3=velocity[0]*velocity[0]+velocity[1]*velocity[1]+velocity[2]*velocity[2];
	s1 = LATTICE_VELOCITIES_D[0][0]*velocity[0]+LATTICE_VELOCITIES_D[0][1]*velocity[1]+LATTICE_VELOCITIES_D[0][2]*velocity[2];
	s2 = s1*s1;
	feq[0]=LATTICE_WEIGHTS_D[0]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[1][0]*velocity[0]+LATTICE_VELOCITIES_D[1][1]*velocity[1]+LATTICE_VELOCITIES_D[1][2]*velocity[2];
	s2 = s1*s1;
	feq[1]=LATTICE_WEIGHTS_D[1]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[2][0]*velocity[0]+LATTICE_VELOCITIES_D[2][1]*velocity[1]+LATTICE_VELOCITIES_D[2][2]*velocity[2];
	s2 = s1*s1;
	feq[2]=LATTICE_WEIGHTS_D[2]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[3][0]*velocity[0]+LATTICE_VELOCITIES_D[3][1]*velocity[1]+LATTICE_VELOCITIES_D[3][2]*velocity[2];
	s2 = s1*s1;
	feq[3]=LATTICE_WEIGHTS_D[3]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[4][0]*velocity[0]+LATTICE_VELOCITIES_D[4][1]*velocity[1]+LATTICE_VELOCITIES_D[4][2]*velocity[2];
	s2 = s1*s1;
	feq[4]=LATTICE_WEIGHTS_D[4]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[5][0]*velocity[0]+LATTICE_VELOCITIES_D[5][1]*velocity[1]+LATTICE_VELOCITIES_D[5][2]*velocity[2];
	s2 = s1*s1;
	feq[5]=LATTICE_WEIGHTS_D[5]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[6][0]*velocity[0]+LATTICE_VELOCITIES_D[6][1]*velocity[1]+LATTICE_VELOCITIES_D[6][2]*velocity[2];
	s2 = s1*s1;
	feq[6]=LATTICE_WEIGHTS_D[6]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[7][0]*velocity[0]+LATTICE_VELOCITIES_D[7][1]*velocity[1]+LATTICE_VELOCITIES_D[7][2]*velocity[2];
	s2 = s1*s1;
	feq[7]=LATTICE_WEIGHTS_D[7]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[8][0]*velocity[0]+LATTICE_VELOCITIES_D[8][1]*velocity[1]+LATTICE_VELOCITIES_D[8][2]*velocity[2];
	s2 = s1*s1;
	feq[8]=LATTICE_WEIGHTS_D[8]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[9][0]*velocity[0]+LATTICE_VELOCITIES_D[9][1]*velocity[1]+LATTICE_VELOCITIES_D[9][2]*velocity[2];
	s2 = s1*s1;
	feq[9]=LATTICE_WEIGHTS_D[9]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[10][0]*velocity[0]+LATTICE_VELOCITIES_D[10][1]*velocity[1]+LATTICE_VELOCITIES_D[10][2]*velocity[2];
	s2 = s1*s1;
	feq[10]=LATTICE_WEIGHTS_D[10]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[11][0]*velocity[0]+LATTICE_VELOCITIES_D[11][1]*velocity[1]+LATTICE_VELOCITIES_D[11][2]*velocity[2];
	s2 = s1*s1;
	feq[11]=LATTICE_WEIGHTS_D[11]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[12][0]*velocity[0]+LATTICE_VELOCITIES_D[12][1]*velocity[1]+LATTICE_VELOCITIES_D[12][2]*velocity[2];
	s2 = s1*s1;
	feq[12]=LATTICE_WEIGHTS_D[12]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[13][0]*velocity[0]+LATTICE_VELOCITIES_D[13][1]*velocity[1]+LATTICE_VELOCITIES_D[13][2]*velocity[2];
	s2 = s1*s1;
	feq[13]=LATTICE_WEIGHTS_D[13]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[14][0]*velocity[0]+LATTICE_VELOCITIES_D[14][1]*velocity[1]+LATTICE_VELOCITIES_D[14][2]*velocity[2];
	s2 = s1*s1;
	feq[14]=LATTICE_WEIGHTS_D[14]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[15][0]*velocity[0]+LATTICE_VELOCITIES_D[15][1]*velocity[1]+LATTICE_VELOCITIES_D[15][2]*velocity[2];
	s2 = s1*s1;
	feq[15]=LATTICE_WEIGHTS_D[15]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[16][0]*velocity[0]+LATTICE_VELOCITIES_D[16][1]*velocity[1]+LATTICE_VELOCITIES_D[16][2]*velocity[2];
	s2 = s1*s1;
	feq[16]=LATTICE_WEIGHTS_D[16]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[17][0]*velocity[0]+LATTICE_VELOCITIES_D[17][1]*velocity[1]+LATTICE_VELOCITIES_D[17][2]*velocity[2];
	s2 = s1*s1;
	feq[17]=LATTICE_WEIGHTS_D[17]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[18][0]*velocity[0]+LATTICE_VELOCITIES_D[18][1]*velocity[1]+LATTICE_VELOCITIES_D[18][2]*velocity[2];
	s2 = s1*s1;
	feq[18]=LATTICE_WEIGHTS_D[18]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);
}


/**
 * Finds an inverse probability distribution for the provided lattice index.
 */
__device__ int InvGpu(int i){
    return (Q_LBM-1)-i;
}

#endif

cell_computation_gpu.cu

#ifndef _CELL_COMPUTATION_GPU_CUH_
#define _CELL_COMPUTATION_GPU_CUH_


/**
 * Computes the density from the particle distribution functions stored at currentCell.
 * currentCell thus denotes the address of the first particle distribution function of the
 * respective cell. The result is stored in density.
 */
__device__ void ComputeDensityGpu(float *current_cell, float *density){
    *density=0;
    *density+=current_cell[0];
    *density+=current_cell[1];
    *density+=current_cell[2];
    *density+=current_cell[3];
    *density+=current_cell[4];
    *density+=current_cell[5];
    *density+=current_cell[6];
    *density+=current_cell[7];
    *density+=current_cell[8];
    *density+=current_cell[9];
    *density+=current_cell[10];
    *density+=current_cell[11];
    *density+=current_cell[12];
    *density+=current_cell[13];
    *density+=current_cell[14];
    *density+=current_cell[15];
    *density+=current_cell[16];
    *density+=current_cell[17];
    *density+=current_cell[18];
}


/**
 * Computes the velocity within currentCell and stores the result in velocity
 */
__device__ void ComputeVelocityGpu(float *current_cell, float *density, float *velocity){
    velocity[0]=0;
	velocity[0]+=current_cell[0]*LATTICE_VELOCITIES_D[0][0];
	velocity[0]+=current_cell[1]*LATTICE_VELOCITIES_D[1][0];
	velocity[0]+=current_cell[2]*LATTICE_VELOCITIES_D[2][0];
	velocity[0]+=current_cell[3]*LATTICE_VELOCITIES_D[3][0];
	velocity[0]+=current_cell[4]*LATTICE_VELOCITIES_D[4][0];
	velocity[0]+=current_cell[5]*LATTICE_VELOCITIES_D[5][0];
	velocity[0]+=current_cell[6]*LATTICE_VELOCITIES_D[6][0];
	velocity[0]+=current_cell[7]*LATTICE_VELOCITIES_D[7][0];
	velocity[0]+=current_cell[8]*LATTICE_VELOCITIES_D[8][0];
	velocity[0]+=current_cell[9]*LATTICE_VELOCITIES_D[9][0];
	velocity[0]+=current_cell[10]*LATTICE_VELOCITIES_D[10][0];
	velocity[0]+=current_cell[11]*LATTICE_VELOCITIES_D[11][0];
	velocity[0]+=current_cell[12]*LATTICE_VELOCITIES_D[12][0];
	velocity[0]+=current_cell[13]*LATTICE_VELOCITIES_D[13][0];
	velocity[0]+=current_cell[14]*LATTICE_VELOCITIES_D[14][0];
	velocity[0]+=current_cell[15]*LATTICE_VELOCITIES_D[15][0];
	velocity[0]+=current_cell[16]*LATTICE_VELOCITIES_D[16][0];
	velocity[0]+=current_cell[17]*LATTICE_VELOCITIES_D[17][0];
	velocity[0]+=current_cell[18]*LATTICE_VELOCITIES_D[18][0];
	velocity[0]/=*density;

	velocity[1]=0;
	velocity[1]+=current_cell[0]*LATTICE_VELOCITIES_D[0][1];
	velocity[1]+=current_cell[1]*LATTICE_VELOCITIES_D[1][1];
	velocity[1]+=current_cell[2]*LATTICE_VELOCITIES_D[2][1];
	velocity[1]+=current_cell[3]*LATTICE_VELOCITIES_D[3][1];
	velocity[1]+=current_cell[4]*LATTICE_VELOCITIES_D[4][1];
	velocity[1]+=current_cell[5]*LATTICE_VELOCITIES_D[5][1];
	velocity[1]+=current_cell[6]*LATTICE_VELOCITIES_D[6][1];
	velocity[1]+=current_cell[7]*LATTICE_VELOCITIES_D[7][1];
	velocity[1]+=current_cell[8]*LATTICE_VELOCITIES_D[8][1];
	velocity[1]+=current_cell[9]*LATTICE_VELOCITIES_D[9][1];
	velocity[1]+=current_cell[10]*LATTICE_VELOCITIES_D[10][1];
	velocity[1]+=current_cell[11]*LATTICE_VELOCITIES_D[11][1];
	velocity[1]+=current_cell[12]*LATTICE_VELOCITIES_D[12][1];
	velocity[1]+=current_cell[13]*LATTICE_VELOCITIES_D[13][1];
	velocity[1]+=current_cell[14]*LATTICE_VELOCITIES_D[14][1];
	velocity[1]+=current_cell[15]*LATTICE_VELOCITIES_D[15][1];
	velocity[1]+=current_cell[16]*LATTICE_VELOCITIES_D[16][1];
	velocity[1]+=current_cell[17]*LATTICE_VELOCITIES_D[17][1];
	velocity[1]+=current_cell[18]*LATTICE_VELOCITIES_D[18][1];
	velocity[1]/=*density;

	velocity[2]=0;
	velocity[2]+=current_cell[0]*LATTICE_VELOCITIES_D[0][2];
	velocity[2]+=current_cell[1]*LATTICE_VELOCITIES_D[1][2];
	velocity[2]+=current_cell[2]*LATTICE_VELOCITIES_D[2][2];
	velocity[2]+=current_cell[3]*LATTICE_VELOCITIES_D[3][2];
	velocity[2]+=current_cell[4]*LATTICE_VELOCITIES_D[4][2];
	velocity[2]+=current_cell[5]*LATTICE_VELOCITIES_D[5][2];
	velocity[2]+=current_cell[6]*LATTICE_VELOCITIES_D[6][2];
	velocity[2]+=current_cell[7]*LATTICE_VELOCITIES_D[7][2];
	velocity[2]+=current_cell[8]*LATTICE_VELOCITIES_D[8][2];
	velocity[2]+=current_cell[9]*LATTICE_VELOCITIES_D[9][2];
	velocity[2]+=current_cell[10]*LATTICE_VELOCITIES_D[10][2];
	velocity[2]+=current_cell[11]*LATTICE_VELOCITIES_D[11][2];
	velocity[2]+=current_cell[12]*LATTICE_VELOCITIES_D[12][2];
	velocity[2]+=current_cell[13]*LATTICE_VELOCITIES_D[13][2];
	velocity[2]+=current_cell[14]*LATTICE_VELOCITIES_D[14][2];
	velocity[2]+=current_cell[15]*LATTICE_VELOCITIES_D[15][2];
	velocity[2]+=current_cell[16]*LATTICE_VELOCITIES_D[16][2];
	velocity[2]+=current_cell[17]*LATTICE_VELOCITIES_D[17][2];
	velocity[2]+=current_cell[18]*LATTICE_VELOCITIES_D[18][2];
    velocity[2]/=*density;
}


/**
 * Computes the equilibrium distributions for all particle distribution functions of one
 * cell from density and velocity and stores the results in feq.
 */
__device__ void ComputeFeqGpu(float *density, float *velocity, float *feq){
    float s1, s2, s3=velocity[0]*velocity[0]+velocity[1]*velocity[1]+velocity[2]*velocity[2];
	s1 = LATTICE_VELOCITIES_D[0][0]*velocity[0]+LATTICE_VELOCITIES_D[0][1]*velocity[1]+LATTICE_VELOCITIES_D[0][2]*velocity[2];
	s2 = s1*s1;
	feq[0]=LATTICE_WEIGHTS_D[0]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[1][0]*velocity[0]+LATTICE_VELOCITIES_D[1][1]*velocity[1]+LATTICE_VELOCITIES_D[1][2]*velocity[2];
	s2 = s1*s1;
	feq[1]=LATTICE_WEIGHTS_D[1]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[2][0]*velocity[0]+LATTICE_VELOCITIES_D[2][1]*velocity[1]+LATTICE_VELOCITIES_D[2][2]*velocity[2];
	s2 = s1*s1;
	feq[2]=LATTICE_WEIGHTS_D[2]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[3][0]*velocity[0]+LATTICE_VELOCITIES_D[3][1]*velocity[1]+LATTICE_VELOCITIES_D[3][2]*velocity[2];
	s2 = s1*s1;
	feq[3]=LATTICE_WEIGHTS_D[3]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[4][0]*velocity[0]+LATTICE_VELOCITIES_D[4][1]*velocity[1]+LATTICE_VELOCITIES_D[4][2]*velocity[2];
	s2 = s1*s1;
	feq[4]=LATTICE_WEIGHTS_D[4]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[5][0]*velocity[0]+LATTICE_VELOCITIES_D[5][1]*velocity[1]+LATTICE_VELOCITIES_D[5][2]*velocity[2];
	s2 = s1*s1;
	feq[5]=LATTICE_WEIGHTS_D[5]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[6][0]*velocity[0]+LATTICE_VELOCITIES_D[6][1]*velocity[1]+LATTICE_VELOCITIES_D[6][2]*velocity[2];
	s2 = s1*s1;
	feq[6]=LATTICE_WEIGHTS_D[6]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[7][0]*velocity[0]+LATTICE_VELOCITIES_D[7][1]*velocity[1]+LATTICE_VELOCITIES_D[7][2]*velocity[2];
	s2 = s1*s1;
	feq[7]=LATTICE_WEIGHTS_D[7]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[8][0]*velocity[0]+LATTICE_VELOCITIES_D[8][1]*velocity[1]+LATTICE_VELOCITIES_D[8][2]*velocity[2];
	s2 = s1*s1;
	feq[8]=LATTICE_WEIGHTS_D[8]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[9][0]*velocity[0]+LATTICE_VELOCITIES_D[9][1]*velocity[1]+LATTICE_VELOCITIES_D[9][2]*velocity[2];
	s2 = s1*s1;
	feq[9]=LATTICE_WEIGHTS_D[9]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[10][0]*velocity[0]+LATTICE_VELOCITIES_D[10][1]*velocity[1]+LATTICE_VELOCITIES_D[10][2]*velocity[2];
	s2 = s1*s1;
	feq[10]=LATTICE_WEIGHTS_D[10]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[11][0]*velocity[0]+LATTICE_VELOCITIES_D[11][1]*velocity[1]+LATTICE_VELOCITIES_D[11][2]*velocity[2];
	s2 = s1*s1;
	feq[11]=LATTICE_WEIGHTS_D[11]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[12][0]*velocity[0]+LATTICE_VELOCITIES_D[12][1]*velocity[1]+LATTICE_VELOCITIES_D[12][2]*velocity[2];
	s2 = s1*s1;
	feq[12]=LATTICE_WEIGHTS_D[12]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[13][0]*velocity[0]+LATTICE_VELOCITIES_D[13][1]*velocity[1]+LATTICE_VELOCITIES_D[13][2]*velocity[2];
	s2 = s1*s1;
	feq[13]=LATTICE_WEIGHTS_D[13]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[14][0]*velocity[0]+LATTICE_VELOCITIES_D[14][1]*velocity[1]+LATTICE_VELOCITIES_D[14][2]*velocity[2];
	s2 = s1*s1;
	feq[14]=LATTICE_WEIGHTS_D[14]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[15][0]*velocity[0]+LATTICE_VELOCITIES_D[15][1]*velocity[1]+LATTICE_VELOCITIES_D[15][2]*velocity[2];
	s2 = s1*s1;
	feq[15]=LATTICE_WEIGHTS_D[15]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[16][0]*velocity[0]+LATTICE_VELOCITIES_D[16][1]*velocity[1]+LATTICE_VELOCITIES_D[16][2]*velocity[2];
	s2 = s1*s1;
	feq[16]=LATTICE_WEIGHTS_D[16]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[17][0]*velocity[0]+LATTICE_VELOCITIES_D[17][1]*velocity[1]+LATTICE_VELOCITIES_D[17][2]*velocity[2];
	s2 = s1*s1;
	feq[17]=LATTICE_WEIGHTS_D[17]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);

	s1 = LATTICE_VELOCITIES_D[18][0]*velocity[0]+LATTICE_VELOCITIES_D[18][1]*velocity[1]+LATTICE_VELOCITIES_D[18][2]*velocity[2];
	s2 = s1*s1;
	feq[18]=LATTICE_WEIGHTS_D[18]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);
}


/**
 * Finds an inverse probability distribution for the provided lattice index.
 */
__device__ int InvGpu(int i){
    return (Q_LBM-1)-i;
}

#endif


10 AM Meeting

lbm_model_gpu.cuh

#ifndef _LBM_MODEL_GPU_CUH_
#define _LBM_MODEL_GPU_CUH_

/* Restricts 3D blocks to have 512 threads (limits: 512 CC<2.x; 1024 CC>2.x) */
#define BLOCK_SIZE 8

/* Arbitrary boundary identifiers */
#define LEFT_BOUNDARY_IDX 0
#define RIGHT_BOUNDARY_IDX 1
#define BOTTOM_BOUNDARY_IDX 2
#define TOP_BOUNDARY_IDX 3
#define BACK_BOUNDARY_IDX 4
#define FRONT_BOUNDARY_IDX 5


__constant__ static const int LATTICE_VELOCITIES_D[19][3] = {
    {0,-1,-1},{-1,0,-1},{0,0,-1},{1,0,-1},{0,1,-1},{-1,-1,0},{0,-1,0},{1,-1,0},
    {-1,0,0}, {0,0,0},  {1,0,0}, {-1,1,0},{0,1,0}, {1,1,0},  {0,-1,1},{-1,0,1},
    {0,0,1},  {1,0,1},  {0,1,1}
};

__constant__ static const float LATTICE_WEIGHTS_D[19] = {
    1.0/36.0, 1.0/36.0, 2.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 2.0/36.0, 1.0/36.0,
    2.0/36.0, 12.0/36.0,2.0/36.0, 1.0/36.0, 2.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0,
    2.0/36.0, 1.0/36.0, 1.0/36.0
};

/**
 * Stores number of probability distributions which we need to copy in boundary treatment step.
 */
__constant__ static const float TREAT_BOUNDARY_INDECES[6][5] = {
	{3,7,10,13,17},
	{1,5,8,11,15},
	{4,11,12,13,18},
	{0,5,6,7,14},
	{14,15,16,17,18},
	{0,1,2,3,4}
};

#endif

CPU:

initialization.h

#ifndef _INITIALIZATION_H_
#define _INITIALIZATION_H_


/**
 * Reads the parameters for the lid driven cavity scenario from a configuration file
 */
void ReadParameters(int *xlength, float *tau, float *velocity_wall, int *timesteps,
		int *timesteps_per_plotting, int argc, char *argv[], int *gpu_enabled);


/**
 * Initializes the particle distribution functions and the flag field
 */
void InitialiseFields(float *collide_field, float *stream_field,int *flag_field, int xlength,
		int gpu_enabled);

#endif

 

initialization.c

#include <unistd.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

#include "initialization.h"
#include "lbm_model.h"
#include "utils.h"


/**
 * Prints help message that includes program usage instructions and control flags.
 */
void PrintHelpMessage(){
	printf("List of control flags:\n");
	printf("\t -gpu             all computations are to be performed on gpu\n");
	printf("\t -cpu             all computations are to be performed on cpu\n");
	printf("\t -help            prints this help message\n");
	printf("NOTE: Control flags are mutually exclusive and only one flag at a time is allowed\n");
	printf("Example program usage:\n");
	printf("\t./lbm-sim ./data/lbm.dat -gpu\n");
	exit(1);
}


void ReadParameters(int *xlength, float *tau, float *velocity_wall, int *timesteps,
		int *timesteps_per_plotting, int argc, char *argv[], int *gpu_enabled){
    float *velocity_wall_1, *velocity_wall_2, *velocity_wall_3;
    /* printing out help message */
    if(argc<3) PrintHelpMessage();
    if(!strcmp(argv[1], "-help") || !strcmp(argv[2], "-help")) PrintHelpMessage();

    /* checking parameters */
    if(access(argv[1], R_OK) != 0)
        ERROR("Provided configuration file path either doesn't exist or can not be read.");
    if(!strcmp(argv[2], "-gpu")) *gpu_enabled=1; else *gpu_enabled=0;

    /* reading parameters */
    READ_FLOAT(argv[1], *tau);

    velocity_wall_1=&velocity_wall[0];
    velocity_wall_2=&velocity_wall[1];
    velocity_wall_3=&velocity_wall[2];

    READ_FLOAT(argv[1], *velocity_wall_1);
    READ_FLOAT(argv[1], *velocity_wall_2);
    READ_FLOAT(argv[1], *velocity_wall_3);
    
    READ_INT(argv[1], *xlength);
    READ_INT(argv[1], *timesteps);
    READ_INT(argv[1], *timesteps_per_plotting);
}


void InitialiseFields(float *collide_field, float *stream_field, int *flag_field, int xlength,
		int gpu_enabled){
    int x,y,z,i,step=xlength+2;
    
    /* NOTE: We use z=xlength+1 as the moving wall */
    for(x=0;x<step;x++){
        for(y=0;y<step;y++){
            for(z=0;z<step;z++){
                /* Initializing flags */
            	if(gpu_enabled) {
					if(x==0 && y==0 && z==0)
						flag_field[x+y*step+z*step*step] = LEFT_BOTTOM_BACK_CORNER;
					else if(x==0 && y==0 && z==step-1)
						flag_field[x+y*step+z*step*step] = LEFT_BOTTOM_FRONT_CORNER;
					else if(x==0 && y==step-1 && z==0)
						flag_field[x+y*step+z*step*step] = LEFT_UPPER_BACK_CORNER;
					else if(x==0 && y==step-1 && z==step-1)
						flag_field[x+y*step+z*step*step] = LEFT_UPPER_FRONT_CORNER;
					else if(x==step-1 && y==0 && z==0)
						flag_field[x+y*step+z*step*step] = RIGHT_BOTTOM_BACK_CORNER;
					else if(x==step-1 && y==0 && z==step-1)
						flag_field[x+y*step+z*step*step] = RIGHT_BOTTOM_FRONT_CORNER;
					else if(x==step-1 && y==step-1 && z==0)
						flag_field[x+y*step+z*step*step] = RIGHT_UPPER_BACK_CORNER;
					else if(x==step-1 && y==step-1 && z==step-1)
						flag_field[x+y*step+z*step*step] = RIGHT_UPPER_FRONT_CORNER;
					else if(x==0 && y==0)
						flag_field[x+y*step+z*step*step]=LEFT_BOTTOM_EDGE;
					else if(x==xlength+1 && y==0)
						flag_field[x+y*step+z*step*step]=RIGHT_BOTTOM_EDGE;
					else if(y==0 && z==0)
						flag_field[x+y*step+z*step*step]=BACK_BOTTOM_EDGE;
					else if(y==0 && z==xlength+1)
						flag_field[x+y*step+z*step*step]=FRONT_BOTTOM_EDGE;

					else if(x==0 && z==0)
						flag_field[x+y*step+z*step*step]=LEFT_BACK_EDGE;
					else if(x==0 && z==xlength+1)
						flag_field[x+y*step+z*step*step]=LEFT_FRONT_EDGE;
					else if(x==xlength+1 && z==0)
						flag_field[x+y*step+z*step*step]=RIGHT_BACK_EDGE;
					else if(x==xlength+1 && z==xlength+1)
						flag_field[x+y*step+z*step*step]=RIGHT_FRONT_EDGE;
					else if(x==0 && y==xlength+1)
						flag_field[x+y*step+z*step*step]=LEFT_UPPER_EDGE;
					else if(x==xlength+1 && y==xlength+1)
						flag_field[x+y*step+z*step*step]=RIGHT_UPPER_EDGE;
					else if(y==xlength+1 && z==0)
						flag_field[x+y*step+z*step*step]=BACK_UPPER_EDGE;
					else if(y==xlength+1 && z==xlength+1)
						flag_field[x+y*step+z*step*step]=FRONT_UPPER_EDGE;


					else if(y==step-1 && x!=0 && x!=step-1 && z!=0 && z!=step-1)
						flag_field[x+y*step+z*step*step] = TOP_BOUNDARY;
					else if(y==0 && x!=0 && x!=step-1 && z!=0 && z!=step-1)
						flag_field[x+y*step+z*step*step] = BOTTOM_BOUNDARY;
					else if (x==0 && y!=0 && y!=step-1 && z!=0 && z!=step-1)
						flag_field[x+y*step+z*step*step] = LEFT_BOUNDARY;
					else if (x==step-1 && y!=0 && y!=step-1 && z!=0 && z!=step-1)
						flag_field[x+y*step+z*step*step] = RIGHT_BOUNDARY;
					else if (z==0 && y!=0 && y!=step-1 && x!=0 && x!=step-1)
						flag_field[x+y*step+z*step*step] = BACK_BOUNDARY;
					else if (z==step-1 && y!=0 && y!=step-1 && x!=0 && x!=step-1)
						flag_field[x+y*step+z*step*step] = FRONT_BOUNDARY;
					else
						flag_field[x+y*step+z*step*step]=FLUID;
            	} else {
            		if(y == xlength+1)
						flag_field[x+y*step+z*step*step]=MOVING_WALL;
					else if(x == 0 || x == xlength+1 || y == 0 || z == xlength+1 || z == 0)
						flag_field[x+y*step+z*step*step]=NO_SLIP;
					else
						flag_field[x+y*step+z*step*step]=FLUID;
            	}
                
                /* Initializing distributions for stream and collide fields */
                for(i=0;i<Q_LBM;i++){
                    /*
                     * NOTE:Stream field is initilized to 0s because that helps to
                     * track down mistakes and has no impact whatsoever to on the 
                     * computation further on.
                     */
                    stream_field[Q_LBM*(x+y*step+z*step*step)+i]=0;
                    collide_field[Q_LBM*(x+y*step+z*step*step)+i]=LATTICE_WEIGHTS[i];
                }
            }
        }
    }
}

 

streaming.h

#ifndef _STREAMING_H_
#define _STREAMING_H_

/**
 * Carries out the streaming step and writes the respective distribution functions from
 * collideField to streamField.
 */
void DoStreaming(float *collide_field, float *stream_field, int *flag_field, int xlength);

#endif

streaming.c

#include "streaming.h"
#include "lbm_model.h"


void DoStreaming(float *collide_field, float *stream_field, int *flag_field, int xlength){
    int x,nx,y,ny,z,nz,i,step=xlength+2;

    for(x=0;x<step;x++){
        for(y=0;y<step;y++){
            for(z=0;z<step;z++){
                if(flag_field[x+y*step+z*step*step]==FLUID){
                    for(i=0;i<Q_LBM;i++){
                        nx=x-LATTICE_VELOCITIES[i][0];
                        ny=y-LATTICE_VELOCITIES[i][1];
                        nz=z-LATTICE_VELOCITIES[i][2];
                        
                        stream_field[Q_LBM*(x+y*step+z*step*step)+i]=
                                collide_field[Q_LBM*(nx+ny*step+nz*step*step)+i];
                    }
                }
            }
        }
    }
}

collision.h

#ifndef _COLLISION_H_
#define _COLLISION_H_

/**
 * Carries out the whole local collision process. Computes density and velocity and equilibrium
 * distributions. Carries out BGK update.
 */
void DoCollision(float *collide_field, int *flag_field, float tau, int xlength);

#endif

collision.c

#include "collision.h"
#include "cell_computation.h"
#include "lbm_model.h"
#include "utils.h"


/**
 * Computes the post-collision distribution functions according to the BGK update rule and
 * stores the results again at the same position.
 */
void ComputePostCollisionDistributions(float *current_cell, float tau, const float *const feq){
    int i;
    for(i=0;i<Q_LBM;i++){
        current_cell[i]=current_cell[i]-(current_cell[i]-feq[i])/tau;
        
        /* Probability distribution function can not be less than 0 */
        if (current_cell[i] < 0)
            ERROR("Probability distribution function can not be negative.");
    }
}


void DoCollision(float *collide_field, int *flag_field, float tau, int xlength){
    float density, velocity[3], feq[Q_LBM], *currentCell;
    int x,y,z,step=xlength+2;
    
    for(x=1;x<step-1;x++){
        for(y=1;y<step-1;y++){
            for(z=1;z<step-1;z++){
                currentCell=&collide_field[Q_LBM*(x+y*step+z*step*step)];
                
                ComputeDensity(currentCell,&density);
                ComputeVelocity(currentCell,&density,velocity);
                ComputeFeq(&density,velocity,feq);
                ComputePostCollisionDistributions(currentCell,tau,feq);
            }
        }
    }
}

 

boundary.h

#ifndef _BOUNDARY_H_
#define _BOUNDARY_H_


/**
 * Handles the boundaries in our simulation setup
 */
void TreatBoundary(float *collide_field, int* flag_field, float *wall_velocity, int xlength);

#endif

boundary.c

#include "boundary.h"
#include "lbm_model.h"
#include "cell_computation.h"

/**
 * Finds an inverse probability distribution for the provided lattice index.
 */
int inv(int i){
    return (Q_LBM-1)-i;
}


void TreatBoundary(float *collide_field, int* flag_field, float *wall_velocity, int xlength){
    int x,nx,y,ny,z,nz,i,step=xlength+2;
    float density,dot_prod;
    
    for(x=0;x<step;x++){
        for(y=0;y<step;y++){
            for(z=0;z<step;z++){
                if(flag_field[x+y*step+z*step*step]!=FLUID){
                    for(i=0;i<Q_LBM;i++){
                        nx=x+LATTICE_VELOCITIES[i][0];
                        ny=y+LATTICE_VELOCITIES[i][1];
                        nz=z+LATTICE_VELOCITIES[i][2];

                        /* We don't need the values outside of our extended domain */
                        if(0<nx && nx<step-1 && 0<ny && ny<step-1 && 0<nz && nz<step-1){
                            if (flag_field[x+y*step+z*step*step]==MOVING_WALL){
                                /* Compute density in the neighbour cell */
                                ComputeDensity(&collide_field[Q_LBM*(nx+ny*step+nz*step*step)],&density);
                                /* Compute dot product */
                                dot_prod=LATTICE_VELOCITIES[i][0]*wall_velocity[0]+
                                        LATTICE_VELOCITIES[i][1]*wall_velocity[1]+
                                        LATTICE_VELOCITIES[i][2]*wall_velocity[2];
                                /* Assign the boudary cell value */
                                collide_field[Q_LBM*(x+y*step+z*step*step)+i]=
                                        collide_field[Q_LBM*(nx+ny*step+nz*step*step)+inv(i)]+
                                        2*LATTICE_WEIGHTS[i]*density*C_S_POW2_INV*dot_prod;
                            }else if(flag_field[x+y*step+z*step*step]==NO_SLIP){
                                collide_field[Q_LBM*(x+y*step+z*step*step)+i]=
                                        collide_field[Q_LBM*(nx+ny*step+nz*step*step)+inv(i)];
                            }
                        }
                    }
                }
            }
        }
    }
}

 

cell_computation.h

#ifndef _CELL_COMPUTATIONS_H_
#define _CELL_COMPUTATIONS_H_

/**
 * Computes the density from the particle distribution functions stored at currentCell.
 * currentCell thus denotes the address of the first particle distribution function of the
 * respective cell. The result is stored in density.
 */
void ComputeDensity(const float *const current_cell, float *density);

/**
 * Computes the velocity within currentCell and stores the result in velocity
 */
void ComputeVelocity(const float * const current_cell, const float * const density, float *velocity);

/**
 * Computes the equilibrium distributions for all particle distribution functions of one
 * cell from density and velocity and stores the results in feq.
 */
void ComputeFeq(const float * const density, const float * const velocity, float *feq);

#endif

cell_computation.c

#include "cell_computation.h"
#include "lbm_model.h"
#include "utils.h"

void ComputeDensity(const float *const current_cell, float *density){
    int i; *density=0;
    for(i=0;i<Q_LBM;i++)
        *density+=current_cell[i];

    /* Density should be close to a unit (ρ~1) */
    if((*density-1.0)>EPS)
        ERROR("Density dropped below error tolerance.");
}


void ComputeVelocity(const float * const current_cell, const float * const density, float *velocity){
    int i;
    /* NOTE:Indexes are hard coded to improve program performance */
    velocity[0]=0;
    velocity[1]=0;
    velocity[2]=0;
    
    for(i=0;i<Q_LBM;i++){
        velocity[0]+=current_cell[i]*LATTICE_VELOCITIES[i][0];
        velocity[1]+=current_cell[i]*LATTICE_VELOCITIES[i][1];
        velocity[2]+=current_cell[i]*LATTICE_VELOCITIES[i][2];
    }
    
    velocity[0]/=*density;
    velocity[1]/=*density;
    velocity[2]/=*density;
}


void ComputeFeq(const float * const density, const float * const velocity, float *feq){
    int i;
    float s1, s2, s3; /* summands */
    /* NOTE:Indexes are hard coded to improve program performance */
    for(i=0;i<Q_LBM;i++){
        s1 = LATTICE_VELOCITIES[i][0]*velocity[0]+LATTICE_VELOCITIES[i][1]*velocity[1]+
                LATTICE_VELOCITIES[i][2]*velocity[2];
        s2 = s1*s1;
        s3 = velocity[0]*velocity[0]+velocity[1]*velocity[1]+velocity[2]*velocity[2];
        
        feq[i]=LATTICE_WEIGHTS[i]*(*density)*(1+s1*C_S_POW2_INV+s2*C_S_POW4_INV/2.0-s3*C_S_POW2_INV/2.0);
        
        /* Probability distribution function can not be less than 0 */
        if (feq[i] < 0)
            ERROR("Probability distribution function can not be negative.");
    }
}

visualization.h

#ifndef _VISUALIZATION_H_
#define _VISUALIZATION_H_


/**
 * Writes the density and velocity field (derived from the distributions in collideField)
 * to a file determined by 'filename' and timestep 't'. You can re-use parts of the code
 * from visual.c (VTK output for Navier-Stokes solver) and modify it for 3D datasets.
 */
void WriteVtkOutput(const float * const collide_field, const int * const flag_field,
        const char * filename, unsigned int t, int xlength);


/**
 * Writes the provided stream or collision field to a file with specified identifiers.
 */
void WriteField(const float * const field, const char * filename, unsigned int t, const int xlength,
		const int rank);


/**
 * Writes the provided flag field to a file with specified identifiers.
 */
void writeFlagField(const int * const flagField, const char * filename, const int xlength, const int rank);


/**
 * Prints out the point by point values of the provided field (4D) to the stdout.
 * @param field
 *          linerized 4D array, with (x,y,z,i)=Q*(x+y*(ncell+2)+z*(ncell+2)*(ncell+2))+i
 * @param ncell
 *          number of inner cells, the ones which are there before adding a boundary layer
 */
void PrintField(float *field, int ncell);

#endif

visualization.c

#include <cstdio>

#include "lbm_model.h"
#include "visualization.h"
#include "cell_computation.h"
#include "utils.h"


void write_vtkHeader( FILE *fp, int xlength) {
    if( fp == NULL ){
        char szBuff[80];
        sprintf( szBuff, "Null pointer in write_vtkHeader" );
        ERROR( szBuff );
        return;
    }

    fprintf(fp,"# vtk DataFile Version 2.0\n");
    fprintf(fp,"generated by CFD-lab course output (written by Denys Korzh, Denys Sobchyshak, Ozel Christo Annamanthadoo  \n");
    fprintf(fp,"ASCII\n");
    fprintf(fp,"\n");
    fprintf(fp,"DATASET STRUCTURED_GRID\n");
    fprintf(fp,"DIMENSIONS  %i %i %i \n", xlength, xlength, xlength);
    fprintf(fp,"POINTS %i float\n", xlength*xlength*xlength );
    fprintf(fp,"\n");
}


void write_vtkPointCoordinates( FILE *fp, int xlength) {
    float originX = 0.0, originY = 0.0;
    int i = 0, j = 0, k = 0;

    for(k = 0; k < xlength; k++)
        for(j = 0; j < xlength; j++)
            for(i = 0; i < xlength; i++)
                fprintf(fp, "%f %f %f\n", originX+(i*1.0/xlength), originY+(j*1.0/xlength), 
                        originY+(k*1.0/xlength) );
}


void WriteVtkOutput(const float * const collideField, const int * const flagField, const char * filename, unsigned int t, int xlength) {
    int i, j, k, len = xlength+2; /* lexicographic order "[ Q * ( z*len*len + y*len + x) + i ]" */
    float velocity[3], density;

    char szFileName[80];
    FILE *fp=NULL;
    sprintf( szFileName, "%s.%i.vtk", filename, t );
    fp = fopen( szFileName, "w");
    if( fp == NULL ){
        char szBuff[80];
        sprintf( szBuff, "Failed to open %s", szFileName );
        ERROR( szBuff );
        return;
    }

    write_vtkHeader( fp, xlength);
    write_vtkPointCoordinates( fp, xlength);

    fprintf(fp,"POINT_DATA %i \n", xlength*xlength*xlength );
    fprintf(fp,"\n");
    fprintf(fp, "VECTORS velocity float\n");
    for(k = 1; k < xlength+1; k++) {
        for(j = 1; j < xlength+1; j++) {
            for(i = 1; i < xlength+1; i++) {
                ComputeDensity (&collideField[ 19 * ( k*len*len + j*len + i) ], &density);
                ComputeVelocity(&collideField[ 19 * ( k*len*len + j*len + i) ], &density, velocity);
                fprintf(fp, "%f %f %f\n",  velocity[0], velocity[1], velocity[2]);
            }
        }
    }

    fprintf(fp,"\n");
    fprintf(fp, "SCALARS density float 1 \n");
    fprintf(fp, "LOOKUP_TABLE default \n");
    for(k = 1; k < xlength+1; k++) {
        for(j = 1; j < xlength+1; j++) {
            for(i = 1; i < xlength+1; i++) {
                ComputeDensity (&collideField[ 19 * ( k*len*len + j*len + i) ], &density);
                fprintf(fp, "%f\n",  density);
            }
        }
    }

    if(fclose(fp)){
        char szBuff[80];
        sprintf( szBuff, "Failed to close %s", szFileName );
        ERROR( szBuff );
    }
}


void PrintField(float *field, int ncell){
    int x,y,z,i,step=ncell+2;

    for(x=0;x<step;x++){
        for(y=0;y<step;y++){
            for(z=0;z<step;z++){
                printf("(%d,%d,%d): ",x,y,z);
                for(i=0;i<Q_LBM;i++)
                    printf("%f ",field[Q_LBM*(x+y*step+z*step*step)+i]);
                printf("\n");
            }
        }
    }
}


void WriteField(const float * const field, const char * filename, unsigned int t, const int xlength, const int rank) {
	int x,y,z,i, stepX=xlength+2,stepY=xlength+2,stepZ=xlength+2;

	char szFileName[80];
	FILE *fp=NULL;
	sprintf( szFileName, "%s-rank%i.%i.out", filename, rank, t );
	fp = fopen( szFileName, "w");
	if( fp == NULL ){
		char szBuff[80];
		sprintf( szBuff, "Failed to open %s", szFileName );
		ERROR( szBuff );
		return;
	}

	for(z=0;z<stepZ;z++){
		for(x=0;x<stepX;x++){
			for(y=0;y<stepY;y++){
				fprintf(fp, "(%d,%d,%d): ",x,y,z);
					for(i=0;i<Q_LBM;i++)
					  fprintf(fp, "%f ",field[Q_LBM*(x+y*stepX+z*stepX*stepY)+i]);
				fprintf(fp, "\n");
			}
		}
	}

	if(fclose(fp)){
		char szBuff[80];
		sprintf( szBuff, "Failed to close %s", szFileName );
		ERROR( szBuff );
	}
}


void writeFlagField(const int * const flagField, const char * filename, const int xlength, const int rank) {
  int x,y,z, stepX=xlength+2,stepY=xlength+2,stepZ=xlength+2;

  char szFileName[80];
  FILE *fp=NULL;
  sprintf( szFileName, "%s-rank%i.out", filename, rank );
  fp = fopen( szFileName, "w");
  if( fp == NULL ){
    char szBuff[80];
    sprintf( szBuff, "Failed to open %s", szFileName );
    ERROR( szBuff );
    return;
  }

  fprintf(fp, "(y,z)\\x  ");
  for(x=0; x<stepX; x++) {
      fprintf(fp, "%2i ",x);
  }
  fprintf(fp, "\n-------");
  for(x=0; x<stepX; x++) {
      fprintf(fp, "---");
  }
  fprintf(fp, "\n");

  for(z=0;z<stepZ;z++){
      for(y=stepY-1;y>=0;y--){
          fprintf(fp, "(%2d,%2d): ",y,z);
          for(x=0;x<stepX;x++){
              fprintf(fp, "%2i ",flagField[x+y*stepX+z*stepX*stepY]);
          }
          fprintf(fp, "\n");
      }
      fprintf(fp, "\n");
  }

  if(fclose(fp)){
    char szBuff[80];
    sprintf( szBuff, "Failed to close %s", szFileName );
    ERROR( szBuff );
  }
}

 

Posted on

FOCUS ANDROID DEBUG BRIDGE QUICK TIPS

 To load apps onto your Focus device, we suggest using Android Debug Bridge (adb) or your favorite command line tool. For more information refer to the Android Debug Bridge User Guide.

 

  You can use Terminal from your desktop or from within Android Studio to navigate to and utilize the debug bridge. (NOTE: On a Mac in order to use adb, you must type “./adb” and then the command each time).

For example to show connected devices:

  “adb devices”

  To install an app on a connected device:

“adb install [directory location of apk]”

 If an app is already installed, and you wish to replace it:

“adb install -r [directory location of apk]”

To uninstall an app from a connected device:

“adb uninstall com.CompanyName.AppName”

Even more Tips:

# set and start launcher:
adb shell cmd package set-home-activity com.htc.mobilevr.launcher/com.htc.vr.unity.WVRUnityVRActivity
adb shell am start -n com.htc.mobilevr.launcher/com.htc.vr.unity.WVRUnityVRActivity

 

# reset launcher back to the Android phone launcher:
adb shell cmd package set-home-activity com.android.launcher/com.android.launcher2.Launcher
adb shell am start -n com.android.launcher/com.android.launcher2.Launcher

 

# run an app:
adb shell am start -n <package-name>/<activity-name>
# e.g. adb shell am start -n com.htc.bowshot_3DOF/com.htc.vr.unity.WVRUnityVRActivity

 

# list 3rd party installed packages:
adb shell cmd package list packages -3

 

# find activity of an installed package: (from an apk: aapt dump badging yourapp.apk)
adb shell cmd package resolve-activity –brief <package-name>

 

# uninstall a package e.g. Vysor:
adb uninstall com.koushikdutta.vysor

 

# turn bluetooth on:
adb root
adb shell service call bluetooth_manager 6

# use 8 to turn off

# for wifi use: adb shell svc wifi enable # or disable
adb unroot

 

# power down:
adb shell reboot -p

 

# check battery level:
adb shell “dumpsys battery | grep level”

 

# adb over wifi (setup wifi connection first):
# find ip address:
adb shell “ifconfig wlan0 | grep ‘inet addr:'”

# restart adb over wifi:
adb tcpip 5555

# connect over wifi: (you can now disconnect usb cable)
adb connect <ip address>:5555

# to revert to usb:
adb usb # or just reboot

Posted on

CODE REFERENCE

Unity Native Spline

 

 

 

SPHERE EXPANDER – FROM Fundamental Flower

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using ProceduralToolkit;
using UnityEngine.UI;
using System;
using ArgosTweenSpacePuppy;
using VRTK;
using System.IO;
using UnityEngine.EventSystems;
using System.Linq;

public class Sphere_Expander : MonoBehaviour
{
    public bool bSeparate_Faces = false;
    public bool bDrawTetrahedron = false;

    public int totalVerts = 0;
    public int filteredVerts = 0;

    public int vSort_Count_0;
    public int vSort_Count_1;
    public int vSort_Count_2;

    private ArgosMeshDraft   aMD_TetraTerrean;
    private ArgosMeshDraft   aMD_TT_Sleeve;
    private ArgosMeshDraft[] aMD_tt_Faces = new ArgosMeshDraft[4];
    private GameObject[]     tt_Faces = new GameObject[4];
    private ArgosMeshDraft   aMD_Spheres;
    private ArgosMeshDraft   cylinder_MD;
    private ArgosMeshDraft   scythe_MD;
    private GameObject       scythe_GO;
    public  Color            col_edge_Cylinder;
    public  Color            col_Sphere;
    private GameObject[]     spheres = new GameObject[4];
    private Vector3[]        sphere_pos = new Vector3[4];
    public  GameObject       meshPF;
    public  GameObject       tetra_Vert_GOPF;
    private GameObject[]     tetra_Verts_GO = new GameObject[4];
    private GameObject[]     tetra_Edges    = new GameObject[6];
    private GameObject       TetraTerrien_GO;
    private GameObject       TT_Sleeve_GO;

    public float[] lengths = new float[256];
    public float[] hyperbola_of_R = new float[256];

    public int all_Shape_Depth = 4;

    Vector3[] vC = new Vector3[4];

    public class Vector_Sort
    {
        public Vector3 vert;
        public float dist = 0;
        public int vIDX = 0;
    }

    internal class PlaneHelper
    {
        /// <summary>
        /// Returns a value indicating what side (positive/negative) of a plane a point is
        /// </summary>
        /// <param name="point">The point to check with</param>
        /// <param name="plane">The plane to check against</param>
        /// <returns>Greater than zero if on the positive side, less than zero if on the negative size, 0 otherwise</returns>
        public static float ClassifyPoint(ref Vector3 point, ref Plane plane)
        {
            return point.x * plane.Normal.x + point.y * plane.Normal.y + point.z * plane.Normal.z + plane.D;
        }

        /// <summary>
        /// Returns the perpendicular distance from a point to a plane
        /// </summary>
        /// <param name="point">The point to check</param>
        /// <param name="plane">The place to check</param>
        /// <returns>The perpendicular distance from the point to the plane</returns>
        public static float PerpendicularDistance(ref Vector3 point, ref Plane plane)
        {
            // dist = (ax + by + cz + d) / sqrt(a*a + b*b + c*c)
            return Mathf.Abs((plane.Normal.x * point.x + plane.Normal.y * point.y + plane.Normal.z * point.z)
                                    / Mathf.Sqrt(plane.Normal.x * plane.Normal.x + plane.Normal.y * plane.Normal.y + plane.Normal.z * plane.Normal.z));
        }
    }

    public struct Plane 
    {
        #region Public Fields

        public float D;
        public Vector3 Normal;

        #endregion Public Fields

        #region Constructors

        public Plane(Vector4 value)
            : this(new Vector3(value.x, value.y, value.z), value.w)
        {

        }

        public Plane(Vector3 normal, float d)
        {
            Normal = normal;
            D = d;
        }

        public Plane(Vector3 a, Vector3 b, Vector3 c)
        {
            Vector3 ab = b - a;
            Vector3 ac = c - a;

            Vector3 cross = Vector3.Cross(ab, ac);
            Normal = Vector3.Normalize(cross);
            D = -(Vector3.Dot(Normal, a));
        }

        public Plane(float a, float b, float c, float d)
            : this(new Vector3(a, b, c), d)
        {

        }

        #endregion Constructors

        public override string ToString()
        {
            return string.Format("{{Normal:{0} D:{1}}}", Normal, D);
        }
    }

    public enum TT_TYPE
    {
        SPHERICAL,
        HYPERBOLIC,
        OTHER,
        NONE,
    }
    public TT_TYPE ttType = TT_TYPE.SPHERICAL;

    public bool bSpheresON = true;
    public bool bTetraTerrean_On = false;

    [Range(0, 50f)]
    public float tetra_Radius;
    private float tetra_Radius_Last = 0;

    [Range(0, 100f)]
    public float scale_TT= 100f;
 
    [Range(0, 500f)]
    public float sphere_Radius;
    private float sphere_Radius_Last = 0;

    [Range(0, 0.5f)]
    public float cylinder_Radius;

    private Vector3[] tetra_Verts  = new Vector3[4];
    private Vector3[] face_Centers = new Vector3[4];

    [Space(12)]
    [Header("Hyperbolic Settings")]
    [Space(12)]

    [Range(0, 1f)]
    public float base_Distance;

    [Range(0, 1f)]
    public float tangent_Base;

    [Range(0, 1f)]
    public float tangent_H;

    [Range(0, 1f)]
    public float smidge;

    public GameObject H_Editor_Label;
    public GameObject C_Editor_Label;
    public GameObject B_Editor_Label;
    public GameObject I_Editor_Label;
    public GameObject T1_Editor_Label;
    public GameObject T2_Editor_Label;
    private int frameCount = 1;
    private Vector3 vCHn, vIHn, vT1n, vT2n;
    private Vector3 vIH;

    StreamWriter sWrite;
    public bool bPrint_Trace = false;

    private Vector3[] beeS = new Vector3[4];

    StreamWriter sLineVectors;

    private ArgosMeshDraft aMD_Nurbs_Mesh;
    GameObject pQuad_CARBON_LINE_go;

    private ArgosMeshDraft aMD_PQUAD_TMP;

    public enum RESOLUTION
    {
        RES_LOW,
        RES_MED,
        RES_HIGH,
    }

    public RESOLUTION rESOLUTION = RESOLUTION.RES_LOW;
    private RESOLUTION resolution_Last = RESOLUTION.RES_LOW;

    private int nRESOLUTION = 7;

    private List<Vector3> lst_mesh_v = new List<Vector3>();
    public float[] spline_lengths = new float[256];

    [Range(0, 10f)]
    public float short_Tan_mag = 1;

    [Range(0, 10f)]
    public float long_Tan_mag = 1;

    void Start()
    {
        sLineVectors = new StreamWriter("sLineVectors.txt");

        aMD_TetraTerrean = new ArgosMeshDraft();
        aMD_Spheres = new ArgosMeshDraft();
        aMD_Spheres.Add(MeshDraft.Sphere(1,16,16));
        cylinder_MD = new ArgosMeshDraft();
        aMD_TT_Sleeve = new ArgosMeshDraft();

        scythe_MD = new ArgosMeshDraft();
        scythe_GO = Instantiate(meshPF, transform);
        scythe_GO.name = "SCYTHE_OUTLINE";

        aMD_Nurbs_Mesh = new ArgosMeshDraft();
        pQuad_CARBON_LINE_go = Instantiate(meshPF, transform);
        pQuad_CARBON_LINE_go.name = "pQuad_CARBON_LINE_go";

        pQuad_CARBON_LINE_go.GetComponent<MeshRenderer>().enabled = false;//so we can see the TT

        sWrite = new StreamWriter("Hyperbolic_Params.txt");

        for (int i = 0; i<4; i++)
        {
            spheres[i] = Instantiate(meshPF, transform);
            spheres[i].name = "sphere_" + i.ToString();
            SetColor_Element(col_Sphere, spheres[i]);
            spheres[i].GetComponent<MeshFilter>().mesh = aMD_Spheres.ToMeshInternal();

            tetra_Verts_GO[i] = Instantiate(tetra_Vert_GOPF, transform);
            tetra_Verts_GO[i].name = "TetVERT_" + i.ToString();

            tt_Faces[i] = Instantiate(meshPF, transform);
            tt_Faces[i].name = "TT_Face_" + i.ToString();
            aMD_tt_Faces[i] = new ArgosMeshDraft();
        }
        for(int i = 0; i<6; i++)
        {
            tetra_Edges[i] = Instantiate(meshPF, transform);
            tetra_Edges[i].name = "edge_" + i.ToString();
            SetColor_Element(col_edge_Cylinder, tetra_Edges[i]);
            tetra_Edges[i].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();
        }
        Set_Tetra_Verts(tetra_Radius);
        TetraTerrien_GO = Instantiate(meshPF, transform);
        TetraTerrien_GO.name = "TetraTerrien_GO";
        aMD_TetraTerrean.Clear();
        //Create_TETRA_TERREAN(all_Shape_Depth);
        //TetraTerrien_GO.GetComponent<MeshFilter>().mesh = aMD_TetraTerrean.ToMeshInternal();

        TT_Sleeve_GO = Instantiate(meshPF, transform);
        TT_Sleeve_GO.name = "TT_Sleeve_GO";
        aMD_TT_Sleeve.Clear();

        aMD_PQUAD_TMP = new ArgosMeshDraft();
    }

    /// <image url="$(SolutionDir)\EMB\Quad_Nurb.png" scale="0.15" /> 


    public void PQuad_SetUp_VertsAndTans()
    {
        Vector3 v0;
        Vector3 v1;
        Vector3 v2;
        Vector3 v3;
        Vector3 t01;
        Vector3 t10;
        Vector3 t13;
        Vector3 t31;
        Vector3 t23;
        Vector3 t32;
        Vector3 t20;
        Vector3 t02;

        v0 = tetra_Verts[2];
        v1 = v0;
        v1.y = -v1.y;

        v2 = tetra_Verts[3];
        v3 = v2;
        v3.y = -v3.y;

        //short_Tan_mag;
        //long_Tan_mag;

        Vector3[] vTansNorm = new Vector3[4];

        vTansNorm[0] = -v0.normalized;
        vTansNorm[1] = -v1.normalized;
        vTansNorm[2] = -v2.normalized;
        vTansNorm[3] = -v3.normalized;

        t01 = v0 + vTansNorm[0] * short_Tan_mag;
        t02 = v0 + vTansNorm[0] * long_Tan_mag;

        t10 = v1 + vTansNorm[1] * short_Tan_mag;
        t13 = v1 + vTansNorm[1] * long_Tan_mag;

        t23 = v2 + vTansNorm[2] * short_Tan_mag;
        t20 = v2 + vTansNorm[2] * long_Tan_mag;

        t32 = v3 + vTansNorm[3] * short_Tan_mag;
        t31 = v3 + vTansNorm[3] * long_Tan_mag;

        P_Quad_Generate(v0, v1, v2, v3,
                       t01, t10, t13, t31,
                       t23, t32, t20, t02);
    }


    public void P_Quad_Generate(Vector3 v0, Vector3 v1, Vector3 v2, Vector3 v3,
                       Vector3 t01, Vector3 t10, Vector3 t13, Vector3 t31,
                       Vector3 t23, Vector3 t32, Vector3 t20, Vector3 t02)
    {
        Set_Resolution();

        lst_mesh_v.Clear();
        aMD_Nurbs_Mesh.Clear();

        float u = 0;
        float v = 0;
        float du = 1 / (float)nRESOLUTION;
        float dv = 1 / (float)nRESOLUTION;

        Vector3 t0, t1;
        Vector3 p0, p1;
        Vector3 p;

        for (int i = 0; i < nRESOLUTION + 1; i++)
        {
            t0 = GetPointOnBezierCurve(t02, t01, t10, t13, get_Adjusted(u, ref spline_lengths));
            p0 = GetPointOnBezierCurve(v0, t01, t10, v1, get_Adjusted(u, ref spline_lengths));
            t1 = GetPointOnBezierCurve(t20, t23, t32, t31, get_Adjusted(u, ref spline_lengths));
            p1 = GetPointOnBezierCurve(v2, t23, t32, v3, get_Adjusted(u, ref spline_lengths));

            Compute_Dist(ref spline_lengths, v0, t02, t20, v2);

            v = 0f;
            for (int j = 0; j < nRESOLUTION + 1; j++)
            {
                p = GetPointOnBezierCurve(p0, t0, t1, p1, get_Adjusted(v, ref spline_lengths));
                lst_mesh_v.Add(p);
                v += dv;
            }
            u += du;
        }

        int stride = nRESOLUTION + 1;
        int k = 0;
        for (int i = 0; i < nRESOLUTION; i++)
        {
            for (int j = 0; j < nRESOLUTION; j++)
            {
                aMD_Nurbs_Mesh.Add(MeshDraft.Quad(lst_mesh_v[k], lst_mesh_v[k + 1], lst_mesh_v[k + stride + 1], lst_mesh_v[k + stride]));
                k++;
            }
            k++;
        }

        Quaternion q;
        float angle = 120;

        ArgosMeshDraft amdTmp = new ArgosMeshDraft();

        amdTmp.Copy_MeshDraft(aMD_Nurbs_Mesh);

        amdTmp.Rotate(Quaternion.AngleAxis(angle, Vector3.up));

        aMD_Nurbs_Mesh.Add(amdTmp);

        amdTmp.Rotate(Quaternion.AngleAxis(angle, Vector3.up));

        aMD_Nurbs_Mesh.Add(amdTmp);

        aMD_PQUAD_TMP.Clear();

        subdivide_PQUAD_TRIANGLE(tetra_Verts[2], tetra_Verts[3], tetra_Verts[1], all_Shape_Depth, spheres[2].transform.localPosition, 3, true);

        aMD_Nurbs_Mesh.Add(aMD_PQUAD_TMP);

        aMD_PQUAD_TMP.Rotate(Quaternion.AngleAxis(180, Vector3.forward));

        aMD_Nurbs_Mesh.Add(aMD_PQUAD_TMP);

        pQuad_CARBON_LINE_go.GetComponent<MeshFilter>().mesh = aMD_Nurbs_Mesh.ToMeshInternal();
        pQuad_CARBON_LINE_go.GetComponent<MeshFilter>().mesh.RecalculateNormals(60);

    }

    void subdivide_PQUAD_TRIANGLE(Vector3 v1, Vector3 v2, Vector3 v3, int depth, Vector3 sector_Foc, int nID, bool bFlipNorm)
    {
        Vector3 v12, v23, v31;
        Vector3 v12_n, v23_n, v31_n;

        if (depth == 0)
        {
            if (bFlipNorm)
            {
                addTriangle_UV_Tag_PQUAD(v1, v3, v2, (float)nID);
            }
            else
            {
                //if (nID == 1)
                //{
                //    vTest[i] = v1;
                //    vTest2[i] = v2;
                //    i++;
                //}
                addTriangle_UV_Tag_PQUAD(v3, v2, v1, (float)nID);
            }
            return;
        }

        v12 = (v1 + v2) / 2.0f;
        v23 = (v2 + v3) / 2.0f;
        v31 = (v3 + v1) / 2.0f;

        v12_n = Sphere_Surf(v12, sector_Foc); //sector_Foc
        v23_n = Sphere_Surf(v23, sector_Foc);
        v31_n = Sphere_Surf(v31, sector_Foc);

        //v12_n = v12; //sector_Foc
        //v23_n = v23;
        //v31_n = v31;

        /* recursively subdivide new triangles */
        subdivide_PQUAD_TRIANGLE(v1, v12_n, v31_n, depth - 1, sector_Foc, 1, bFlipNorm);
        subdivide_PQUAD_TRIANGLE(v12_n, v2, v23_n, depth - 1, sector_Foc, 2, bFlipNorm);
        subdivide_PQUAD_TRIANGLE(v31_n, v23_n, v3, depth - 1, sector_Foc, 3, bFlipNorm);
        subdivide_PQUAD_TRIANGLE(v23_n, v31_n, v12_n, depth - 1, sector_Foc, 4, bFlipNorm);
    }

    private void Create_Top_Bottom(int depth)
    {
        //sleeve(tetra_Verts[1], tetra_Verts[2], tetra_Verts[3], spheres[2].transform.localPosition);
        subdivide_CarbonLine(tetra_Verts[1], tetra_Verts[3], tetra_Verts[2], depth, spheres[2].transform.localPosition, 3, false);
        List<Vector_Sort>[] vsortList = new List<Vector_Sort>[3];

        vsortList[0] = new List<Vector_Sort>();
        vsortList[1] = new List<Vector_Sort>();
        vsortList[2] = new List<Vector_Sort>();

        int[] id = new int[] { 0, 1, 2, 0, 2, 3, 0, 3, 1 };
        Vector3 vert;

        for (int j = 0; j < 3; j++)
        {
            Plane p = new Plane(tetra_Verts[id[3 * j]], tetra_Verts[id[3 * j + 1]], tetra_Verts[id[3 * j + 2]]);
            p.D = p.D * 0.999f;

            Vector_Sort vs;

            filteredVerts = 0;

            for (int i = 0; i < aMD_TT_Sleeve.vertices.Count; i++)
            {
                vert = aMD_TT_Sleeve.vertices[i];
                float res = PlaneHelper.ClassifyPoint(ref vert, ref p);

                if (res > 0)
                {
                    filteredVerts++;
                    vs = new Vector_Sort();
                    vs.vert = vert;
                    vs.vIDX = i;
                    vsortList[j].Add(vs);
                }
            }
            RemoveDuplicates(vsortList[j]);
            for (int i = 0; i < vsortList[j].Count; i++)
            {
                vsortList[j][i].dist = (tetra_Verts[j + 1] - vsortList[j][i].vert).magnitude;
            }

            var ordered = from element in vsortList[j]
                          orderby element.dist
                          select element;

            vsortList[j] = ordered.ToList<Vector_Sort>();
        }

        for (int i = 0; i < aMD_TT_Sleeve.vertices.Count; i++)
        {
            aMD_TT_Sleeve.vertices[i] = Sphere_Surf(aMD_TT_Sleeve.vertices[i], spheres[2].transform.localPosition);
        }

        Plane p2 = new Plane(Vector3.up, 0);
        int vertCount = aMD_TT_Sleeve.vertices.Count;

        ArgosMeshDraft amdTemp = new ArgosMeshDraft();

        amdTemp.Add(aMD_TT_Sleeve);
        aMD_TT_Sleeve.FlipTriangles();

        aMD_TT_Sleeve.Add(amdTemp);
        aMD_TT_Sleeve.FlipNormals();
        float d;

        for (int i = 0; i < vertCount; i++)//FLIP
        {
            vert = aMD_TT_Sleeve.vertices[i];
            d = PlaneHelper.ClassifyPoint(ref vert, ref p2);
            aMD_TT_Sleeve.vertices[i] = vert - 2f * d * Vector3.up;
        }

        List<Vector3> vInnerLower = new List<Vector3>();
        List<Vector3> vInnerUpper = new List<Vector3>();

        List<Vector3> vPrint = new List<Vector3>();

        for (int j = 0; j < 3; j++)
        {
            for (int i = 0; i < vsortList[j].Count; i++)
            {
                vert = aMD_TT_Sleeve.vertices[vsortList[j][i].vIDX];
                //vCurr.y = 0;
                d = PlaneHelper.ClassifyPoint(ref vert, ref p2);

                if (j == 0)
                {
                    vPrint.Add(vert - 2f * d * Vector3.up);
                }

                vInnerLower.Add(vert);
                vInnerUpper.Add(vert - 2f * d * Vector3.up);
            }
        }
        aMD_TT_Sleeve.Add(MeshDraft.Band(vInnerUpper, vInnerLower));
        aMD_TT_Sleeve.Add(MeshDraft.Band(vInnerLower, vInnerUpper));

        Generate_Line_List(vPrint);

        vSort_Count_0 = vsortList[0].Count;
        vSort_Count_1 = vsortList[1].Count;
        vSort_Count_2 = vsortList[2].Count;
    }

    private float get_Adjusted(float t, ref float[] lengths)
    {
        int i = 0;
        while (i < 256 && lengths[i] < t)
        {
            i++;
        }
        return (float)i / 256;
    }

    private void Compute_Dist(ref float[] lengths, Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3)
    {
        float t = 0.0f;
        float dt = 1f / 256;

        Vector3 vB = GetPointOnBezierCurve(p0, p1, p2, p3, t);
        Vector3 vB_Last = vB;
        float running_Len = 0;

        for (int i = 0; i < 256; i++)
        {
            vB = GetPointOnBezierCurve(p0, p1, p2, p3, t);
            running_Len += (vB - vB_Last).magnitude;
            lengths[i] = running_Len;
            vB_Last = vB;
            t += dt;
        }
        for (int i = 0; i < 256; i++)
        {
            lengths[i] /= running_Len;
        }
    }

    private void Set_Resolution()
    {
        if (rESOLUTION != resolution_Last)
        {
            if (rESOLUTION == RESOLUTION.RES_LOW)
            {
                nRESOLUTION = 7;
            }
            else if (rESOLUTION == RESOLUTION.RES_MED)
            {
                nRESOLUTION = 18;
            }
            else if (rESOLUTION == RESOLUTION.RES_HIGH)
            {
                nRESOLUTION = 36;
            }
        }
        resolution_Last = rESOLUTION;
    }


    //void OnDrawGizmos()
    //{
    //    Gizmos.DrawIcon(H_Editor_Label.transform.position, "H.png", true);
    //    Gizmos.DrawIcon(C_Editor_Label.transform.position, "C.png", true);
    //    Gizmos.DrawIcon(I_Editor_Label.transform.position, "I.png", true);
    //    Gizmos.DrawIcon(B_Editor_Label.transform.position, "B.png", true);
    //    Gizmos.DrawIcon(T1_Editor_Label.transform.position, "T1.png", true);
    //    Gizmos.DrawIcon(T2_Editor_Label.transform.position, "T2.png", true);
    //    Gizmos.color = Color.green;
    //    Gizmos.DrawLine(B_Editor_Label.transform.position, T1_Editor_Label.transform.position);
    //    Gizmos.DrawLine(H_Editor_Label.transform.position, T2_Editor_Label.transform.position);
    //}

    private void OnApplicationQuit()
    {
        sWrite.Close();
    }

    /// <image url="$(SolutionDir)\EMB\hyperbole2.png" scale="0.35"></image>

    private void build_Reference_Hyperbolic_Spline(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3, Vector3 vI)
    {
        //p4 = H
        float t = 0.0f;
        float dt = 1f / 256;

        Vector3 vB = GetPointOnBezierCurve(p0, p1, p2, p3, t);
        Vector3 vB_Last = vB;
        float running_Len = 0;

        float   vI_Len = vI.magnitude;
        Vector3 vIn    = vI.normalized;

        for (int i = 0; i < 256; i++)
        {
            vB = GetPointOnBezierCurve(p0, p1, p2, p3, t);
            running_Len += (vB - vB_Last).magnitude;
            lengths[i] = running_Len;
            hyperbola_of_R[i] = Vector3.Dot(vB, vIn);
            vB_Last = vB;
            t += dt;
        }
        for (int i = 0; i < 256; i++)
        {
            lengths[i] /= running_Len;
        }
    }

    public void Print_Trace(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3, Vector3 vI)
    {
        for (int i = 0; i < 256; i++)
        {
            sWrite.WriteLine(i.ToString() + "\t " + " -\t " + hyperbola_of_R[i].ToString("F2"));
        }
        sWrite.WriteLine("p0 Base " + "\t " + " -\t " + p0.ToString("F2"));
        sWrite.WriteLine("p1 Tangent_0 " + "\t " + " -\t " + p1.ToString("F2"));
        sWrite.WriteLine("p2 Tangent_1 " + "\t " + " -\t " + p2.ToString("F2"));
        sWrite.WriteLine("p3 H " + "\t " + " -\t " + p3.ToString("F2"));
        sWrite.WriteLine("vI " + "\t " + " -\t " + vI.ToString("F2"));
    }

    private void SetColor_Element(Color col, GameObject go)
    {
        float alpha = col.a;
        float alpha_out = alpha / 3f;

        Color cola = col;
        cola.a = alpha_out;

        go.GetComponent<MeshRenderer>().material.SetColor("_EmissionColor", cola);
        go.GetComponent<MeshRenderer>().material.SetColor("_Color", col);
    }

    public void Set_Tetra_Verts(float radius)
    {
        float tetrahedralAngle = Mathf.PI * -19.471220333f / 180;
        float segmentAngle = Mathf.PI * 2 / 3;
        float currentAngle = 0f;

        Vector3 v = new Vector3(0, radius, 0);
        tetra_Verts[0] = v;
        for (var i = 1; i < 4; i++)
        {
            tetra_Verts[i] = PTUtils.PointOnSphere(radius, currentAngle, tetrahedralAngle);
            currentAngle += segmentAngle;
        }

        face_Centers[0] = (tetra_Verts[0] + tetra_Verts[1] + tetra_Verts[2]) / 3f;
        face_Centers[1] = (tetra_Verts[0] + tetra_Verts[2] + tetra_Verts[3]) / 3f;
        face_Centers[2] = (tetra_Verts[1] + tetra_Verts[2] + tetra_Verts[3]) / 3f;
        face_Centers[3] = (tetra_Verts[0] + tetra_Verts[1] + tetra_Verts[3]) / 3f;

        float edge_length = (tetra_Verts[0] - tetra_Verts[1]).magnitude;
        edge_length = radius * 2f * Mathf.Sqrt(2f) / Mathf.Sqrt(3);

        cylinder_MD.Add(MeshDraft.Cylinder_Z(cylinder_Radius, 5, edge_length));

        //tetra_Edges[0].transform.localPosition = tetra_Verts[0];
        //tetra_Edges[0].transform.localRotation = Quaternion.LookRotation(tetra_Verts[1] - tetra_Verts[0]);
        //tetra_Edges[0].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();

        //tetra_Edges[1].transform.localPosition = tetra_Verts[1];
        //tetra_Edges[1].transform.localRotation = Quaternion.LookRotation(tetra_Verts[2] - tetra_Verts[1]);
        //tetra_Edges[1].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();

        //tetra_Edges[2].transform.localPosition = tetra_Verts[2];
        //tetra_Edges[2].transform.localRotation = Quaternion.LookRotation(tetra_Verts[3] - tetra_Verts[2]);
        //tetra_Edges[2].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();

        //tetra_Edges[3].transform.localPosition = tetra_Verts[3];
        //tetra_Edges[3].transform.localRotation = Quaternion.LookRotation(tetra_Verts[0] - tetra_Verts[3]);
        //tetra_Edges[3].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();

        //tetra_Edges[4].transform.localPosition = tetra_Verts[1];
        //tetra_Edges[4].transform.localRotation = Quaternion.LookRotation(tetra_Verts[3] - tetra_Verts[1]);
        //tetra_Edges[4].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();

        //tetra_Edges[5].transform.localPosition = tetra_Verts[2];
        //tetra_Edges[5].transform.localRotation = Quaternion.LookRotation(tetra_Verts[0] - tetra_Verts[2]);
        //tetra_Edges[5].GetComponent<MeshFilter>().mesh = cylinder_MD.ToMeshInternal();
    }

    private void Create_Scythe(int depth)
    {


        subdivide_CarbonLine(tetra_Verts[1], tetra_Verts[3], tetra_Verts[2], depth, spheres[2].transform.localPosition, 3, false);

    }

    void Update() 
    {
        PQuad_SetUp_VertsAndTans();

        cylinder_MD.Clear();
        scythe_MD.Clear();
        Set_Tetra_Verts(tetra_Radius);

        if (sphere_Radius_Last != sphere_Radius)
        {
            float face_dist;
            float tetra_Height;
            Vector3 face_Norm;
            float a, b;

            aMD_Spheres.Clear();
            aMD_Spheres.Add(MeshDraft.Sphere(sphere_Radius, 36, 32));
            for (int i = 0; i < 4; i++)
            {
                face_dist = face_Centers[i].magnitude;
                tetra_Height = 1.33333f * tetra_Radius;
                a = 2 * Mathf.Sqrt(2) / 3;
                b = Mathf.Sqrt(sphere_Radius * sphere_Radius - tetra_Radius * tetra_Radius * 8f / 9f);
                face_Norm = face_Centers[i].normalized;

                spheres[i].transform.localPosition = face_Norm * (tetra_Radius / 3 + b);
                spheres[i].transform.localRotation = Quaternion.LookRotation(sphere_pos[i]);
                spheres[i].GetComponent<MeshFilter>().mesh = aMD_Spheres.ToMeshInternal();
                tetra_Verts_GO[i].transform.localPosition = tetra_Verts[i];
            }
        }
        sphere_Radius_Last = sphere_Radius;
        tetra_Radius_Last = tetra_Radius;
        if (bSpheresON)
        {
            ShowSpheres(true);
        }
        else
        {
            ShowSpheres(false);
        }

        aMD_TetraTerrean.Clear();

        for (int i = 0; i < 4; i++)
        {
            aMD_tt_Faces[i].Clear();
        }

        Create_TETRA_TERREAN(all_Shape_Depth);
        if (bSeparate_Faces)
        {
            for (int i = 0; i < 4; i++)
            {
                tt_Faces[i].GetComponent<MeshFilter>().mesh = aMD_tt_Faces[i].ToMeshInternal();
            }
        }
        else
        {
            aMD_TetraTerrean.FlipTriangles();
            aMD_TetraTerrean.FlipNormals();

            TetraTerrien_GO.GetComponent<MeshFilter>().mesh = aMD_TetraTerrean.ToMeshInternal();
            TetraTerrien_GO.GetComponent<MeshFilter>().mesh.RecalculateNormals(60);
        }

        //Create_Scythe(all_Shape_Depth);

        //aMD_TT_Sleeve.Clear();
        //Create_TETRA_TERREAN_Sleeve(all_Shape_Depth);
        //TT_Sleeve_GO.GetComponent<MeshFilter>().mesh = aMD_TT_Sleeve.ToMeshInternal();
        //TT_Sleeve_GO.GetComponent<MeshFilter>().mesh.RecalculateNormals(60);

        sphere_Radius_Last = sphere_Radius;
        tetra_Radius_Last  = tetra_Radius;
    }

    private void ShowSpheres(bool bOn)
    {
        for(int i = 0; i<4; i++)
        {
            spheres[i].GetComponent<MeshRenderer>().enabled = bOn;
        }
    }

    public void Calc_Tetra_Verts(float radius)
    {
        float tetrahedralAngle = Mathf.PI * -19.471220333f / 180;
        float segmentAngle = Mathf.PI * 2 / 3;
        float currentAngle = 0f;

        Vector3 v = new Vector3(0, radius, 0);
        sphere_pos[0] = v;
        for (var i = 1; i < 4; i++)
        {
            sphere_pos[i] = PTUtils.PointOnSphere(radius, currentAngle, tetrahedralAngle);
            currentAngle += segmentAngle;
        }
    }

    private Vector3 Sphere_Surf(Vector3 vP0, Vector3 sector_Focus)
    {
        Vector3 l = vP0.normalized;

        Vector3 OminC = -sector_Focus;
        float rsqr = sphere_Radius * sphere_Radius;

        float dotLOminC = Vector3.Dot(l, OminC);
        float Omag = OminC.magnitude;

        float d = -dotLOminC - Mathf.Sqrt(dotLOminC * dotLOminC - Omag * Omag + sphere_Radius * sphere_Radius);

        return l * d;
    }

    private Vector3 Hyperbolic(Vector3 vP0, Vector3 face_Center)
    {
        Vector3 vFCn = face_Center.normalized;
        Vector3 dotV = Vector3.Dot(vP0, vFCn)* vFCn;
        Vector3 radV = vP0 - dotV;
        Vector3 vP0n = vP0.normalized;
        float   radial_mag = radV.magnitude;
        float   perc = (radial_mag / vIH.magnitude)*255;
        float   gVal = 0f;

        //if (perc < 256)
        //{
        gVal = hyperbola_of_R[(int)perc];
        //}

        return  (gVal*vFCn + radV);
    }
     
    void Extrude_1(Vector3 v1, Vector3 v2, Vector3 v3, int nID, Vector3 sector_Foc)
    {
        Vector3 nSect = tetra_Verts[0].normalized;
        float spar = (1f / 9f) * sphere_Radius;

        float d1, d2, d3;

        d1 = Vector3.Dot(v1, nSect);
        d2 = Vector3.Dot(v2, nSect);
        d3 = Vector3.Dot(v3, nSect);

        Vector3 ob1 = (spar - d1) * nSect;
        Vector3 ob2 = (spar - d2) * nSect;
        Vector3 ob3 = (spar - d3) * nSect;

        if (nID == 0)
        {
            addTriangle(v1 + 2 * ob1, v2 + 2 * ob2, v3 + 2 * ob3);
        }
        else
        {
            ob1 = PosCheck(nSect, d1, (d1 / spar), v1 - d1 * nSect);
            ob2 = PosCheck(nSect, d2, (d2 / spar), v2 - d2 * nSect);
            ob3 = PosCheck(nSect, d3, (d3 / spar), v3 - d3 * nSect);
            addTriangle(v1 + ob1, v2 + ob2, v3 + ob3);
        }
    }

    private Vector3 PosCheck(Vector3 v1, float d, float multiplier, Vector3 vHorz)
    {
        return Vector3.zero;
    }

    public static void RemoveDuplicates(List<Vector_Sort> list)
    {
        if (list == null)
        {
            return;
        }
        int i = 1;
        while (i < list.Count)
        {
            int j = 0;
            bool remove = false;
            while (j < i && !remove)
            {
                if (list[i].vert.Equals(list[j].vert))
                {
                    remove = true;
                }
                j++;
            }
            if (remove)
            {
                list.RemoveAt(i);
            }
            else
            {
                i++;
            }
        }
    }

    private void Create_TETRA_TERREAN_Sleeve(int depth)
    {
        //sleeve(tetra_Verts[1], tetra_Verts[2], tetra_Verts[3], spheres[2].transform.localPosition);
        subdivide_CarbonLine(tetra_Verts[1], tetra_Verts[3], tetra_Verts[2], depth, spheres[2].transform.localPosition, 3, false);
        List<Vector_Sort>[] vsortList = new List<Vector_Sort>[3];

        vsortList[0] = new List<Vector_Sort>();
        vsortList[1] = new List<Vector_Sort>();
        vsortList[2] = new List<Vector_Sort>();

        int[] id = new int [] { 0, 1, 2, 0, 2, 3, 0, 3, 1 }; 
        Vector3   vert;

        for (int j = 0; j < 3; j++)
        {
            Plane p = new Plane(tetra_Verts[id[3*j]], tetra_Verts[id[3*j+1]], tetra_Verts[id[3*j+2]]);
            p.D = p.D * 0.999f;
            
            Vector_Sort vs;

            filteredVerts = 0;

            for (int i = 0; i < aMD_TT_Sleeve.vertices.Count; i++)
            {
                vert = aMD_TT_Sleeve.vertices[i];
                float res = PlaneHelper.ClassifyPoint(ref vert, ref p);

                if (res > 0)
                {
                    filteredVerts++;
                    vs = new Vector_Sort();
                    vs.vert = vert;
                    vs.vIDX = i;
                    vsortList[j].Add(vs);
                }
            }
            RemoveDuplicates(vsortList[j]);
            for (int i = 0; i < vsortList[j].Count; i++)
            {
                vsortList[j][i].dist = (tetra_Verts[j+1] - vsortList[j][i].vert).magnitude;
            }

            var ordered = from element in vsortList[j]
                      orderby element.dist
                      select element;

            vsortList[j] = ordered.ToList<Vector_Sort>();
        }

        for(int i = 0; i < aMD_TT_Sleeve.vertices.Count; i++)
        {
            aMD_TT_Sleeve.vertices[i] = Sphere_Surf(aMD_TT_Sleeve.vertices[i], spheres[2].transform.localPosition);
        }

        Plane p2 = new Plane(Vector3.up, 0);
        int vertCount = aMD_TT_Sleeve.vertices.Count;

        ArgosMeshDraft amdTemp = new ArgosMeshDraft();

        amdTemp.Add(aMD_TT_Sleeve);
        aMD_TT_Sleeve.FlipTriangles();
        
        aMD_TT_Sleeve.Add(amdTemp);
        aMD_TT_Sleeve.FlipNormals();
        float d;

        for(int i = 0; i < vertCount; i++)//FLIP
        {
            vert = aMD_TT_Sleeve.vertices[i];
            d = PlaneHelper.ClassifyPoint(ref vert, ref p2);
            aMD_TT_Sleeve.vertices[i] = vert - 2f * d * Vector3.up;
        }

        List<Vector3> vInnerLower = new List<Vector3>();
        List<Vector3> vInnerUpper = new List<Vector3>();

        List<Vector3> vPrint = new List<Vector3>();

        for (int j = 0; j < 3; j++)
        {
            for (int i = 0; i < vsortList[j].Count; i++)
            {
                vert = aMD_TT_Sleeve.vertices[vsortList[j][i].vIDX];
                //vCurr.y = 0;
                d = PlaneHelper.ClassifyPoint(ref vert, ref p2);

                if(j==0)
                {
                    vPrint.Add(vert - 2f * d * Vector3.up);
                }

                vInnerLower.Add(vert);
                vInnerUpper.Add(vert - 2f * d * Vector3.up);
            }
        }
        aMD_TT_Sleeve.Add(MeshDraft.Band(vInnerUpper, vInnerLower));
        aMD_TT_Sleeve.Add(MeshDraft.Band(vInnerLower, vInnerUpper));

        Generate_Line_List(vPrint);

        vSort_Count_0 = vsortList[0].Count;
        vSort_Count_1 = vsortList[1].Count;
        vSort_Count_2 = vsortList[2].Count;
    }


    bool bwritten = false;
    void Generate_Line_List(List<Vector3> vPrint)
    {
        if (!bwritten)
        {
            vPrint.Insert(0, tetra_Verts[1]);
            vPrint.Add(tetra_Verts[2]);
            //vPrint[0] = tetra_Verts[1];
            //vPrint[vPrint.Count - 1] = tetra_Verts[2];

            for (int i = 0; i < vPrint.Count-1; i++)
            {
                sLineVectors.Write("new Vector3(" + vPrint[i].x.ToString("F3") + "f, " + vPrint[i].y.ToString("F3") + "f, " + vPrint[i].z.ToString("F3") + "f), " +
                                   "new Vector3(" + vPrint[i+1].x.ToString("F3") + "f, " + vPrint[i+1].y.ToString("F3") + "f, " + vPrint[i+1].z.ToString("F3") + "f), ");
            }

            Quaternion q;
            float angle = 120f;
            Vector3 axis = Vector3.up;

            q = Quaternion.AngleAxis(angle, axis);
            Vector3 vRota;
            Vector3 vRota2;

            for (int i = 0; i< vPrint.Count-1; i++)
            {
                vRota  = q * vPrint[i];
                vRota2 = q * vPrint[i+1];
                sLineVectors.Write("new Vector3(" + vRota.x.ToString("F3") + "f, " + vRota.y.ToString("F3") + "f, " + vRota.z.ToString("F3") + "f), " +
                    "new Vector3(" + vRota2.x.ToString("F3") + "f, " + vRota2.y.ToString("F3") + "f, " + vRota2.z.ToString("F3") + "f), ");


            }
            angle = 240;
            q = Quaternion.AngleAxis(angle, axis);
            for (int i = 0; i < vPrint.Count - 1; i++)
            {
                vRota = q * vPrint[i];
                vRota2 = q * vPrint[i + 1];
                sLineVectors.Write("new Vector3(" + vRota.x.ToString("F3") + "f, " + vRota.y.ToString("F3") + "f, " + vRota.z.ToString("F3") + "f), " +
                    "new Vector3(" + vRota2.x.ToString("F3") + "f, " + vRota2.y.ToString("F3") + "f, " + vRota2.z.ToString("F3") + "f), ");

            }
            /////////////////////////UPPERS//////////////////////

            List<Vector3> vlUppers = new List<Vector3>();
            axis = tetra_Verts[1].normalized;
            angle = 120;
            q = Quaternion.AngleAxis(angle, axis);
            for (int i = 0; i < vPrint.Count - 1; i++)
            {
                vRota = q * vPrint[i];
                vRota2 = q * vPrint[i + 1];

                vlUppers.Add(vRota);
                vlUppers.Add(vRota2);

                sLineVectors.Write("new Vector3(" + vRota.x.ToString("F3") + "f, " + vRota.y.ToString("F3") + "f, " + vRota.z.ToString("F3") + "f), " +
                    "new Vector3(" + vRota2.x.ToString("F3") + "f, " + vRota2.y.ToString("F3") + "f, " + vRota2.z.ToString("F3") + "f), ");

            }

            axis = Vector3.up;
            angle = 120;
            q = Quaternion.AngleAxis(angle, axis);
            for (int i = 0; i < vlUppers.Count - 1; i++)
            {
                vRota = q * vlUppers[i];
                vRota2 = q * vlUppers[i + 1];
                sLineVectors.Write("new Vector3(" + vRota.x.ToString("F3") + "f, " + vRota.y.ToString("F3") + "f, " + vRota.z.ToString("F3") + "f), " +
                    "new Vector3(" + vRota2.x.ToString("F3") + "f, " + vRota2.y.ToString("F3") + "f, " + vRota2.z.ToString("F3") + "f), ");
            }
            angle = 240;
            q = Quaternion.AngleAxis(angle, axis);
            for (int i = 0; i < vlUppers.Count - 1; i++)
            {
                vRota = q * vlUppers[i];
                vRota2 = q * vlUppers[i + 1];
                sLineVectors.Write("new Vector3(" + vRota.x.ToString("F3") + "f, " + vRota.y.ToString("F3") + "f, " + vRota.z.ToString("F3") + "f), " +
                    "new Vector3(" + vRota2.x.ToString("F3") + "f, " + vRota2.y.ToString("F3") + "f, " + vRota2.z.ToString("F3") + "f), ");
            }
        }
        sLineVectors.Close();
        bwritten = true;
    }

    void sleeve(Vector3 v1, Vector3 v2, Vector3 v3, Vector3 sector_Foc)
    {
        Vector3 vCurr;

        //v1.y = v2.y = v3.y = 0;

        List<Vector3> vInnerLower = new List<Vector3>();
        List<Vector3> vInnerUpper = new List<Vector3>();

        Vector3 vUp = Vector3.up;

        Vector3[] vAlpha = new[] { v1, v2, v3 };
        Vector3[] vOmega = new[] { v2, v3, v1 };

        float t = 0;
        float dt = 1 / 60f;
        float d;
        float el;
        for (int j = 0; j < 3; j++)
        {
            t = 0;
            for (int i = 0; i < 60; i++)
            {
                vCurr = Sphere_Surf(Vector3.Lerp(vAlpha[j], vOmega[j], t), sector_Foc);
                //vCurr.y = 0;
                d  =  vCurr.y;
                el = Mathf.Sqrt(tetra_Radius * tetra_Radius - d * d); 

                //vInnerLower.Add(vCurr + el*vUp);
                //vInnerUpper.Add(vCurr - el*vUp);

                vInnerLower.Add(vCurr);
                vInnerUpper.Add(vCurr - 2f * d * vUp);

                t += dt;
            }
        }
        aMD_TT_Sleeve.Add(MeshDraft.Band(vInnerUpper, vInnerLower));
        aMD_TT_Sleeve.Add(MeshDraft.Band(vInnerLower, vInnerUpper));
    }


    Vector3[] vTest = new Vector3[1024];
    Vector3[] vTest2 = new Vector3[1024];
    int i = 0;
    void subdivide_CarbonLine(Vector3 v1, Vector3 v2, Vector3 v3, int depth, Vector3 sector_Foc, int nID, bool bFlipNorm)
    {
        Vector3 v12, v23, v31;
        Vector3 v12_n, v23_n, v31_n;

        if (depth == 0)
        {
            if (bFlipNorm)
            {
                addTriangle_UV_Tag_CarbonLine(v1, v3, v2, (float)nID);
            }
            else
            {
                //if (nID == 1)
                //{
                //    vTest[i] = v1;
                //    vTest2[i] = v2;
                //    i++;
                //}
                addTriangle_UV_Tag_CarbonLine(v1, v2, v3, (float)nID);
            }
            return;
        }

        v12 = (v1 + v2) / 2.0f;
        v23 = (v2 + v3) / 2.0f;
        v31 = (v3 + v1) / 2.0f;

        //v12_n = Sphere_Surf(v12, sector_Foc); //sector_Foc
        //v23_n = Sphere_Surf(v23, sector_Foc);
        //v31_n = Sphere_Surf(v31, sector_Foc);

        v12_n = v12; //sector_Foc
        v23_n = v23;
        v31_n = v31;

        /* recursively subdivide new triangles */
        subdivide_CarbonLine(v1, v12_n, v31_n, depth - 1, sector_Foc, 1, bFlipNorm);
        subdivide_CarbonLine(v12_n, v2, v23_n, depth - 1, sector_Foc, 2, bFlipNorm);
        subdivide_CarbonLine(v31_n, v23_n, v3, depth - 1, sector_Foc, 3, bFlipNorm);
        subdivide_CarbonLine(v23_n, v31_n, v12_n, depth - 1, sector_Foc, 4, bFlipNorm);
    }

    private void Create_TETRA_TERREAN(int depth)
    {
        float scl = scale_TT / 100f;

        if (ttType == TT_TYPE.SPHERICAL)
        {
            subdivide(tetra_Verts[0], tetra_Verts[1], tetra_Verts[2], depth, spheres[0].transform.localPosition, 0, true);
            subdivide(tetra_Verts[0], tetra_Verts[2], tetra_Verts[3], depth, spheres[1].transform.localPosition, 1, true);
            subdivide(tetra_Verts[0], tetra_Verts[3], tetra_Verts[1], depth, spheres[3].transform.localPosition, 2, true);
            subdivide(tetra_Verts[1], tetra_Verts[3], tetra_Verts[2], depth, spheres[2].transform.localPosition, 3, true);
            //subdivide(tetra_Verts[0], tetra_Verts[1], tetra_Verts[2], depth, spheres[0].transform.localPosition, 0, false);//Inside
            //subdivide(tetra_Verts[0], tetra_Verts[2], tetra_Verts[3], depth, spheres[1].transform.localPosition, 1, false);
            //subdivide(tetra_Verts[0], tetra_Verts[3], tetra_Verts[1], depth, spheres[3].transform.localPosition, 2, false);
        }
        else if(ttType == TT_TYPE.HYPERBOLIC)
        {
            Vector3[] fc = new Vector3[4];

            fc[0] = (tetra_Verts[0] + tetra_Verts[1] + tetra_Verts[2]) / 3f;
            fc[1] = (tetra_Verts[0] + tetra_Verts[2] + tetra_Verts[3]) / 3f;
            fc[2] = (tetra_Verts[1] + tetra_Verts[3] + tetra_Verts[2]) / 3f;
            fc[3] = (tetra_Verts[0] + tetra_Verts[3] + tetra_Verts[1]) / 3f;

            beeS[0] = base_Distance * tetra_Radius * fc[0].normalized;
            beeS[1] = base_Distance * tetra_Radius * fc[1].normalized;
            beeS[2] = base_Distance * tetra_Radius * fc[2].normalized;
            beeS[3] = base_Distance * tetra_Radius * fc[3].normalized;

            subdivide(tetra_Verts[0], tetra_Verts[1], tetra_Verts[2], depth, fc[0], 0, false);
            subdivide(tetra_Verts[0], tetra_Verts[2], tetra_Verts[3], depth, fc[1], 1, false);
            subdivide(tetra_Verts[1], tetra_Verts[3], tetra_Verts[2], depth, fc[2], 2, false);
            subdivide(tetra_Verts[0], tetra_Verts[3], tetra_Verts[1], depth, fc[3], 3, false);

            Scale_Radially();
        }
    }

    private void Scale_Radially()
    {
        float ratio = 0;
        Vector3 v;
        float fidx;
        float len;

        for (int i = 0; i<aMD_TetraTerrean.vertices.Count; i++)
        {  
            fidx  = aMD_TetraTerrean.uv[i].x;

            ratio = aMD_TetraTerrean.vertices[i].magnitude / tetra_Radius;

            aMD_TetraTerrean.vertices[i] *= ((smidge * ratio) + (1 - smidge));
        }
    }

    void subdivide(Vector3 v1, Vector3 v2, Vector3 v3, int depth, Vector3 sector_Foc, int nID, bool bFlipNorm)
    {
        Vector3 v12, v23, v31;
        Vector3 v12_n, v23_n, v31_n;

        if (depth == 0)
        {
            if (ttType == TT_TYPE.SPHERICAL || ttType == TT_TYPE.HYPERBOLIC)
            {
                if (bFlipNorm)
                {
                    addTriangle_UV_Tag(v1, v3, v2, (float)nID);
                }
                else
                {
                    addTriangle_UV_Tag(v1, v2, v3, (float)nID);
                }
            }
            else if (ttType == TT_TYPE.OTHER)
            {
                Extrude_1(v1, v2, v3, nID, sector_Foc);
            }
            return;
        }

        v12 = (v1 + v2) / 2.0f;
        v23 = (v2 + v3) / 2.0f;
        v31 = (v3 + v1) / 2.0f;

        //intrude midpoints
        if (ttType == TT_TYPE.SPHERICAL)
        {
            v12_n = Sphere_Surf(v12, sector_Foc); //sector_Foc
            v23_n = Sphere_Surf(v23, sector_Foc);
            v31_n = Sphere_Surf(v31, sector_Foc);
        }
        else if(ttType == TT_TYPE.HYPERBOLIC)
        {
            v12_n = Hyperbolic(v12, sector_Foc); //sector_Foc
            v23_n = Hyperbolic(v23, sector_Foc);
            v31_n = Hyperbolic(v31, sector_Foc);
        }
        else if (ttType == TT_TYPE.OTHER)
        {
            v12_n = Extrude_1(v12, sector_Foc); //sector_Foc
            v23_n = Extrude_1(v23, sector_Foc);
            v31_n = Extrude_1(v31, sector_Foc);
        }
        else
        {
            v12_n = Vector3.zero;
            v23_n = Vector3.zero;
            v31_n = Vector3.zero;
        }
        /* recursively subdivide new triangles */
        subdivide(v1,    v12_n, v31_n, depth - 1, sector_Foc, nID, bFlipNorm);
        subdivide(v12_n, v2,    v23_n, depth - 1, sector_Foc, nID, bFlipNorm);
        subdivide(v31_n, v23_n, v3,    depth - 1, sector_Foc, nID, bFlipNorm);
        subdivide(v23_n, v31_n, v12_n, depth - 1, sector_Foc, nID, bFlipNorm);
    }

    private Vector3 Extrude_1(Vector3 vP0, Vector3 sector_Focus)
    {
        Vector3 l = vP0.normalized;

        Vector3 OminC = -sector_Focus;
        float rsqr = sphere_Radius * sphere_Radius;

        float dotLOminC = Vector3.Dot(l, OminC);
        float Omag = OminC.magnitude;

        float d = -dotLOminC - Mathf.Sqrt(dotLOminC * dotLOminC - Omag * Omag + sphere_Radius * sphere_Radius);

        return (l * d);
    }

    void addTriangle(Vector3 v0, Vector3 v1, Vector3 v2)
    {
        aMD_TetraTerrean.Add(MeshDraft.Triangle(v0, v1, v2));
    }

    void addTriangle_UV_Tag(Vector3 v0, Vector3 v1, Vector3 v2, float uvTag)
    {
        if (bSeparate_Faces)
        {
            aMD_tt_Faces[(int)uvTag].Add(MeshDraft.Triangl_UV_Tag(v0, v1, v2, uvTag));
        }
        else
        {
            aMD_TetraTerrean.Add(MeshDraft.Triangl_UV_Tag(v0, v1, v2, uvTag));
        }
    }

    void addTriangle_UV_Tag_CarbonLine(Vector3 v0, Vector3 v1, Vector3 v2, float uvTag)
    {
        aMD_TT_Sleeve.Add(MeshDraft.Triangl_UV_Tag(v0, v1, v2, uvTag));
    }

    void addTriangle_Skythe(Vector3 v0, Vector3 v1, Vector3 v2, float uvTag)
    {
        scythe_MD.Add(MeshDraft.Triangl_UV_Tag(v0, v1, v2, uvTag));
    }

    void addTriangle_UV_Tag_PQUAD(Vector3 v0, Vector3 v1, Vector3 v2, float uvTag)
    {
        aMD_PQUAD_TMP.Add(MeshDraft.Triangl_UV_Tag(v0, v1, v2, uvTag));
    }

    public Vector3 GetPointOnBezierCurve(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3, float t)
    {
        float u = 1f - t;
        float t2 = t * t;
        float u2 = u * u;
        float u3 = u2 * u;
        float t3 = t2 * t;

        Vector3 result =
            (u3) * p0 +
            (3f * u2 * t) * p1 +
            (3f * u * t2) * p2 +
            (t3) * p3;

        return result;
    }

    /// <image url="$(SolutionDir)\EMB\hyperbole.png" scale="0.45"></image>
    /// 



}







/// <image url="$(SolutionDir)\EMB\CD.png" scale="0.16803" /> 


//int countDown = 10;
//void Update() // Original
//{
//    H_Editor_Label.transform.position  = tetra_Verts_GO[0].transform.position;
//    C_Editor_Label.transform.position  = transform.position;

//    Vector3 v = (tetra_Verts_GO[0].transform.localPosition + tetra_Verts_GO[1].transform.localPosition + tetra_Verts_GO[2].transform.localPosition) /3;

//    I_Editor_Label.transform.position = v + transform.position;
//    B_Editor_Label.transform.position = transform.position + base_Distance*tetra_Radius*v.normalized;

//    vCHn = (H_Editor_Label.transform.localPosition - C_Editor_Label.transform.localPosition).normalized;
//    vIH  = (H_Editor_Label.transform.localPosition - I_Editor_Label.transform.localPosition);
//    vIHn = vIH.normalized;

//    vT1n = Vector3.Dot(vCHn, vIHn) * vIHn;
//    vT2n = vCHn;

//    T1_Editor_Label.transform.localPosition = (vT1n * tangent_Base * tetra_Radius + B_Editor_Label.transform.localPosition);
//    T2_Editor_Label.transform.localPosition = (H_Editor_Label.transform.localPosition - vT2n * tangent_H * tetra_Radius);

//    if(frameCount-- == 0 || bPrint_Trace)
//    {
//        frameCount = 10;
//        build_Reference_Hyperbolic_Spline(B_Editor_Label.transform.localPosition,  T1_Editor_Label.transform.localPosition, 
//                                          T2_Editor_Label.transform.localPosition, H_Editor_Label.transform.localPosition,
//                                          I_Editor_Label.transform.localPosition);

//        if(--countDown == 0 || bPrint_Trace)
//        {
//            Print_Trace(B_Editor_Label.transform.localPosition, T1_Editor_Label.transform.localPosition,
//                                          T2_Editor_Label.transform.localPosition, H_Editor_Label.transform.localPosition,
//                                          I_Editor_Label.transform.localPosition);             
//        }
//    }
//    bPrint_Trace = false;
//    aMD_Spheres.Clear();
//    aMD_Spheres.Add(MeshDraft.Sphere(sphere_Radius, 36, 32));
//    cylinder_MD.Clear();

//    Set_Tetra_Verts(tetra_Radius);

//    float face_dist;
//    float tetra_Height;
//    Vector3 face_Norm;
//    float a, b;

//    if (sphere_Radius_Last != sphere_Radius)
//    {
//        for (int i = 0; i < 4; i++)
//        {
//            face_dist = face_Centers[i].magnitude;
//            tetra_Height = 1.33333f * tetra_Radius;
//            a = 2 * Mathf.Sqrt(2) / 3;
//            b = Mathf.Sqrt(sphere_Radius * sphere_Radius - tetra_Radius * tetra_Radius * 8f / 9f);
//            face_Norm = face_Centers[i].normalized;

//            spheres[i].transform.localPosition = face_Norm * (tetra_Radius / 3 + b);
//            spheres[i].transform.localRotation = Quaternion.LookRotation(sphere_pos[i]);
//            spheres[i].GetComponent<MeshFilter>().mesh = aMD_Spheres.ToMeshInternal();
//            tetra_Verts_GO[i].transform.localPosition = tetra_Verts[i];
//        }
//    }
//    sphere_Radius_Last = sphere_Radius;
//    tetra_Radius_Last = tetra_Radius;
//    //if ((sphere_Radius_Last != sphere_Radius || tetra_Radius_Last != tetra_Radius) && bTetraTerrean_On)
//    //{
//    aMD_TetraTerrean.Clear();

//    for (int i = 0; i < 4; i++)
//    {
//        aMD_tt_Faces[i].Clear();
//    }

//    Create_TETRA_TERREAN(all_Shape_Depth);
//    if (bSeparate_Faces)
//    {
//        for (int i = 0; i < 4; i++)
//        {
//            tt_Faces[i].GetComponent<MeshFilter>().mesh = aMD_tt_Faces[i].ToMeshInternal();
//        }
//    }
//    else
//    {
//        TetraTerrien_GO.GetComponent<MeshFilter>().mesh = aMD_TetraTerrean.ToMeshInternal();
//        TetraTerrien_GO.GetComponent<MeshFilter>().mesh.RecalculateNormals(60);
//    }

//    aMD_TT_Sleeve.Clear();
//    Create_TETRA_TERREAN_Sleeve();
//    TT_Sleeve_GO.GetComponent<MeshFilter>().mesh = aMD_TT_Sleeve.ToMeshInternal();
//    TT_Sleeve_GO.GetComponent<MeshFilter>().mesh.RecalculateNormals(60);

//    sphere_Radius_Last = sphere_Radius;
//    tetra_Radius_Last  = tetra_Radius;

//    if (bSpheresON)
//    {
//        ShowSpheres(true);
//    }
//    else
//    {
//        ShowSpheres(false);
//    }
//}

 

Posted on

Unexpected Result

    public void OnShow_Differential_Rings()
    {
        if (amd == null)
        {
            amd = new ArgosMeshDraft();
        }
        else
        {
            amd.Clear();
        }

        GetComponent<MeshRenderer>().enabled = true;
        setColor(orbitColor);

        List<Vector3> vInnerLower = new List<Vector3>();
        List<Vector3> vOuterLower = new List<Vector3>();
        List<Vector3> vInnerUpper = new List<Vector3>();
        List<Vector3> vOuterUpper = new List<Vector3>();

        List<Vector3> vDisk = new List<Vector3>();

        float op = Mathf.PI * 2f;
        float delta_theta = op / 180;
        float theta = 0;

        Vector3 vPos_Lower;
        Vector3 vPos_Upper;
        Vector3 vNorm;
        float w_by_2 = cylinder_Width * 0.66666f;//to differentiate from Helix
        float rad = cylinder_Radius;

        Vector3 vZero = (cylinder_Offset + cylinder_Height / 2) * Vector3.up;
        vDisk.Add(vZero);

        for (int i = 0; i < 180; i++)
        {
            vPos_Lower = rad * Mathf.Cos(theta) * Vector3.right + rad * Mathf.Sin(theta) * Vector3.forward;

            float theta_Mod = theta % (Mathf.PI*2f/3f);
            float lout = rad;

            if (theta_Mod < 60)
            {
                lout = (1f - 0.5f * theta_Mod / (Mathf.PI/ 3f));
            }
            else
            {
                lout = (0.5f + 0.5f * (theta_Mod - (Mathf.PI / 3f)) / (Mathf.PI / 3f));
            }
            vPos_Lower *= lout;
            vPos_Upper = vPos_Lower;
            vNorm = vPos_Lower.normalized;
            vPos_Lower += cylinder_Offset * Vector3.up;

            vInnerLower.Add(vPos_Lower - w_by_2 * vNorm);
            vOuterLower.Add(vPos_Lower + w_by_2 * vNorm);

            vPos_Lower += (cylinder_Height / 2) * Vector3.up;
            vDisk.Add(vPos_Lower - w_by_2 * vNorm);

            vPos_Upper += cylinder_Height * Vector3.up + cylinder_Offset * Vector3.up;

            vInnerUpper.Add(vPos_Upper - w_by_2 * vNorm);
            vOuterUpper.Add(vPos_Upper + w_by_2 * vNorm);

            theta += delta_theta;
        }
        amd.Add(MeshDraft.Band(vOuterLower, vInnerLower));
        amd.Add(MeshDraft.Band(vInnerUpper, vOuterUpper));
        amd.Add(MeshDraft.Band(vInnerLower, vInnerUpper));
        amd.Add(MeshDraft.Band(vOuterUpper, vOuterLower));

        if (bAdd_Disk)
        {
            amd.Add(MeshDraft.TriangleFan(vDisk));
        }
        GetComponent<MeshFilter>().mesh = amd.ToMeshInternal();
    }

 

 

Posted on

Code Reference – AllShape

 

Thor Brigsted

using UnityEngine;
using System.Collections;
using System.Collections.Generic;

/// <summary>
/// Helper class that represents a parameterized vertex
/// </summary>
public class Vector3Param {
	///bernstein polynomial packing 
	public List<List<float>> bernPolyPack;

	///Point after applying s,t,u to p0, should result in original point
	public Vector3 p = Vector3.zero;

	///Origin
	public Vector3 p0 = Vector3.zero;

	///Distances along S/T/U axes
	public float s,t,u;

	public Vector3Param()
	{
		s = 0.0f;
		t = 0.0f;
		u = 0.0f;
	}

	public Vector3Param(Vector3Param v)
	{
		s = v.s;
		t = v.t;
		u = v.u;
		p = v.p;
		p0 = v.p0;
	}

};

/// <summary>
/// Free form deformation class
/// 
/// Based off of the paper 'Free-Form Deformation of Solid Geometric Models'[1] this class
/// creates a system of control points that can deform a mesh as if that mesh was embedded
/// in a flexible parallelpiped.
/// 
/// Confused?  Yeah, who uses the term parallelpiped.  The idea is to create a uniformly spaced
/// grid of control points around some mesh.  Each control point has some effect on the mesh, like
/// bone weighting in animation.  The effect each control point has is directly proportional to the
/// total number of control points in the entire grid, each exerts some control.
/// 
/// [1] - http://pages.cpsc.ucalgary.ca/~blob/papers/others/ffd.pdf
/// </summary>
public class FreeFormDeformer : MonoBehaviour {

	/// <summary>
	/// Allow FixedUpdate to modify the mesh.
	/// </summary>
	public bool AllowMeshUpdate = false;

	/// <summary>
	/// Target to be morphed
	/// </summary>
	Mesh MorphTarget = null;

	/// <summary>
	/// Target to be filtered (assumed to contain a meshfilter and valid mesh)
	/// </summary>
	public MeshFilter MorphTargetFilter = null;

	/// <summary>
	/// Update frequency in seconds
	/// </summary>
	public float UpdateFrequency = 1.0f;

	/// <summary>
	/// Game object to represent a control point.  Can be anything really, I suggest spheres.
	/// </summary>
	public GameObject ControlPointPrefab;

	/// <summary>
	/// Local coordinate system
	/// </summary>
	Vector3 S, T, U;

	/// <summary>
	/// Number of controls for S, T, & U respectively.  (L,M, and N MUST be >= 1)
	/// </summary>
	public int L=1, M=1, N=1;

	/// <summary>
	/// Time elapsed since last update
	/// </summary>
	float elapsedTime = 0.0f;

	/// <summary>
	/// Grid of controls points. Stored as 3D grid for easier because width,height, and depth can randomly vary. 
	/// </summary>
	GameObject[, ,] controlPoints;

	/// <summary>
	/// Original vertices from MorphTarget
	/// </summary>
	Vector3[] originalVertices;

	/// <summary>
	/// Current updated vertices for MorphTarget
	/// </summary>
	Vector3[] transformedVertices;


	/// <summary>
	/// Vertex parameters
	/// 
	/// Each vertex is given a set of parameters that will define
	/// its final position based on a local coordinate system.
	/// </summary>
	List<Vector3Param> vertexParams = new List<Vector3Param>();

	void Start () {
		MorphTarget = MorphTargetFilter.mesh ;
		originalVertices = MorphTarget.vertices;
		transformedVertices = new Vector3[originalVertices.Length];

		Parameterize();
	}

	/// <summary>
	/// Calculate a binomial coefficient using the multiplicative formula
	/// </summary>
	float binomialCoeff(int n, int k){
		float total = 1.0f;
		for(int i = 1; i <= k; i++){
			total *= (n - (k - i)) / (float)i;
		}
		return total;
	}

	/// <summary>
	/// Calculate a bernstein polynomial
	/// </summary>
	float bernsteinPoly(int n, int v, float x)
	{
		return binomialCoeff(n,v) * Mathf.Pow(x, (float)v) * Mathf.Pow((float)(1.0f - x), (float)(n - v));
	}

	/// <summary>
	/// Calculate local coordinates
	/// </summary>
	void calculateSTU(Vector3 max, Vector3 min){
		S = new Vector3(max.x - min.x, 0.0f, 0.0f);
		T = new Vector3(0.0f, max.y - min.y, 0.0f);
		U = new Vector3(0.0f, 0.0f, max.z - min.z);
	}

	/// <summary>
	/// Calculate the trivariate bernstein polynomial as described by [1]
	/// 
	/// My method adapts [1] slightly by precalculating the BP coefficients and storing
	/// them in Vector3Param.  When it comes time to extract a world coordinate, 
	/// it's just a matter of summing up multiplications through each polynomial from eq (2).
	/// </summary>
	/// <links>
	///  [1] - Method based on: http://pages.cpsc.ucalgary.ca/~blob/papers/others/ffd.pdf
	/// </links>
	/// <param name="p0">Origin of our coordinate system (where STU meet)</param>
	void calculateTrivariateBernsteinPolynomial(Vector3 p0){
		Vector3 TcU = Vector3.Cross(T, U);
		Vector3 ScU = Vector3.Cross(S, U);
		Vector3 ScT = Vector3.Cross(S, T);

		float TcUdS = Vector3.Dot(TcU, S);
		float ScUdT = Vector3.Dot(ScU, T);
		float ScTdU = Vector3.Dot(ScT, U);

		for (int v = 0; v < originalVertices.Length; v++)
		{
			Vector3 diff = originalVertices[v] - p0;

			Vector3Param tmp = new Vector3Param();
			tmp.s = Vector3.Dot(TcU, diff / TcUdS);
			tmp.t = Vector3.Dot(ScU, diff / ScUdT);
			tmp.u = Vector3.Dot(ScT, diff / ScTdU);
			tmp.p = p0 + (tmp.s * S) + (tmp.t * T) + (tmp.u * U);
			tmp.p0 = p0;
			tmp.bernPolyPack = new List<List<float>>();

			{	// Reserve room for each bernstein polynomial pack.
				tmp.bernPolyPack.Add(new List<float>(L));	//outer bernstein poly
				tmp.bernPolyPack.Add(new List<float>(M));	//middle bernstein poly
				tmp.bernPolyPack.Add(new List<float>(N));	//inner bernstein poly
			}

			{	// Pre-calculate bernstein polynomial expansion.  It only needs to be done once per parameterization
				for (int i = 0; i <= L; i++)
				{
					for (int j = 0; j <= M; j++)
					{
						for (int k = 0; k <= N; k++)
						{
							tmp.bernPolyPack[2].Add(bernsteinPoly(N, k, tmp.u));
						}
						tmp.bernPolyPack[1].Add(bernsteinPoly(M, j, tmp.t));
					}
					tmp.bernPolyPack[0].Add(bernsteinPoly(L, i, tmp.s));
				}
			}
			vertexParams.Add(tmp);
			if (Vector3.Distance(tmp.p, originalVertices[v]) > 0.001f)
			{
				//Debug.Log("Warning, mismatched parameterization");
			}
		}
	}

	/// <summary>
	/// Parameterize MorphTarget's vertices
	/// </summary>
	void Parameterize(){
		Vector3 min = new Vector3(Mathf.Infinity,Mathf.Infinity,Mathf.Infinity);
		Vector3 max = new Vector3(-Mathf.Infinity,-Mathf.Infinity,-Mathf.Infinity);
		foreach(Vector3 v in originalVertices){
			max = Vector3.Max(v,max);
			min = Vector3.Min(v,min);
		}
		calculateSTU(max, min);
		calculateTrivariateBernsteinPolynomial(min);
		createControlPoints(min);
	}


	/// <summary>
	/// Create grid of control points.
	/// </summary>
	void createControlPoints(Vector3 origin){
		controlPoints = new GameObject[L + 1, M + 1, N + 1];
		for(int i = 0; i <= L; i++){
			for(int j = 0; j <= M; j++){
				for(int k = 0; k <= N; k++){
					controlPoints[i, j, k] = createControlPoint(origin, i, j, k);
				}
			}
		}
	}

	/// <summary>
	/// Create a single control point.
	/// </summary>
	GameObject createControlPoint(Vector3 p0, int i, int j, int k)
	{
		Vector3 position = p0 + (i / (float)L * S) + (j / (float)M * T) + (k / (float)N * U);
		return (GameObject)Instantiate(ControlPointPrefab, position, Quaternion.identity);
	}

	/// <summary>
	/// Convert parameterized vertex in to a world coordinate
	/// </summary>
	Vector3 getWorldVector3(Vector3Param r){
		int l = L;
		int m = M;
		int n = N;

		Vector3 tS = Vector3.zero;
		for(int i = 0; i <= l; i++){
			
			Vector3 tM = Vector3.zero;
			for(int j = 0; j <= m; j++){

				Vector3 tK = Vector3.zero;
				for(int k = 0; k <= n; k++){
					tK += r.bernPolyPack[2][k] * controlPoints[i,j,k].transform.position;
				}
				tM += r.bernPolyPack[1][j] * tK;
			}
			tS += r.bernPolyPack[0][i] * tM;
		}
		return tS;
	}


	void UpdateMesh(){
		elapsedTime = 0.0f;

		int idx = 0;
		foreach(Vector3Param vp in vertexParams){
			Vector3 p = getWorldVector3(vp);
			transformedVertices[idx++] = p;
		}

		MorphTarget.vertices = transformedVertices;
		MorphTarget.RecalculateBounds();
		MorphTarget.RecalculateNormals();
		MorphTarget.Optimize();
	}

	void FixedUpdate()
	{
		elapsedTime += Time.fixedDeltaTime;
		if (AllowMeshUpdate)
		{
			if (elapsedTime >= UpdateFrequency) UpdateMesh();
		}
	}
	// Update is called once per frame
	void Update () {
	
	}
}

 

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using VRTK;
using VRTK.GrabAttachMechanics;
using ArgosTweenSpacePuppy;

public class Tetra_MAJOR : MonoBehaviour
{
    public enum State
    {
        FREE,
        ATTACHED,
    }

    public bool bTurn_Off_Collider = false;

    public State state;

    public enum SpokeState
    {
        OPEN,
        OCCUPIED,
    }

    public class Spoke_Tracker
    {
        public SpokeState spokeState = SpokeState.OPEN;
        public GameObject go_occupying = null;
        public float touchTimer = 0;
        public GameObject genSpoke;
    }

    private int id = 0;

    public GameObject[] vertsHT = new GameObject[4]; 

    public Spoke_Tracker[] spoke_Tracker = new Spoke_Tracker[4];
    [Range(0, 10f)]
    public float Radius_of_Action = 7f;

    [Range(0, 1.61803f)]
    public float tDurr = 0.7f;

    private bool bForce_update = false;

    private float tAccum = 0f;

    private float tween_Val = 0f;

    private float touchTimer = 0;

    private Spoke_Tracker attachedSpoke;
    private GameObject attachedTetra_go;

    private Vector3 start_Pos;
    private Vector3 end_Pos;
    private Quaternion qStart;
    private Quaternion qEnd;

    private HMD_Ctrl_Tracking hmd_Ctrl_Tracking;
    private VU_UI_MANAGER VU_UI;
    private ALL_Shape_Instancer all_shape_instancer;
    private bool bTouchGrabEnabled = true;
    private bool bGrabbed = false;

    private VRTK_InteractableObject vrtk_interact;

    public bool bInteractive = true;

    void Start()
    {
        hmd_Ctrl_Tracking = HMD_Ctrl_Tracking.Instance;
        VU_UI = VU_UI_MANAGER.Instance;

        all_shape_instancer = VU_UI.GetComponent<ALL_Shape_Instancer>();

        vrtk_interact = GetComponent<VRTK_InteractableObject>();

        vrtk_interact.InteractableObjectGrabbed   += new InteractableObjectEventHandler(DoObjectGrabbed);
        vrtk_interact.InteractableObjectUngrabbed += new InteractableObjectEventHandler(DoObjectUnGrabbed);

        Vector3[] face_Center = new Vector3[4];

        face_Center[0] = (vertsHT[0].transform.localPosition + vertsHT[1].transform.localPosition + vertsHT[2].transform.localPosition) / 3f;
        face_Center[1] = (vertsHT[1].transform.localPosition + vertsHT[2].transform.localPosition + vertsHT[3].transform.localPosition) / 3f;
        face_Center[2] = (vertsHT[0].transform.localPosition + vertsHT[2].transform.localPosition + vertsHT[3].transform.localPosition) / 3f;
        face_Center[3] = (vertsHT[0].transform.localPosition + vertsHT[1].transform.localPosition + vertsHT[3].transform.localPosition) / 3f;

        float avg_FC = (face_Center[0].magnitude + face_Center[1].magnitude + face_Center[2].magnitude + face_Center[3].magnitude) / 4f;

        Vector3 axis = new Vector3();
        Quaternion q;

        for (int i = 0; i < 4; i++)
        {
            spoke_Tracker[i] = new Spoke_Tracker();
            spoke_Tracker[i].spokeState = SpokeState.OPEN;
            spoke_Tracker[i].genSpoke = new GameObject();

            spoke_Tracker[i].genSpoke.name = "genSpoke_" + i.ToString();

            spoke_Tracker[i].genSpoke.transform.parent = transform;

            spoke_Tracker[i].genSpoke.transform.localPosition = face_Center[i].normalized * 2f*avg_FC;

            axis = vertsHT[i].transform.localPosition - face_Center[i];

            q = Quaternion.AngleAxis(180f, axis);

            spoke_Tracker[i].genSpoke.transform.localRotation = q;
        }


        //if (bInteractive)
        //{
        //    for (int i = 0; i < tet.Length; i++)
        //    {
        //        tet[i].id = i;
        //        tetras.Add(tet[i]);
        //    }
        //}


    }

    private void DoObjectGrabbed(object sender, InteractableObjectEventArgs e)
    {
        bGrabbed = true;

        touchTimer = 1f;

        //if (state == State.ATTACHED)
        //{
        //    //octa_Major.Detach(id);
        //    state = State.FREE;
        //}
    }

    private void DoObjectUnGrabbed(object sender, InteractableObjectEventArgs e)
    {
        bGrabbed = false;
    }

    public void Detach(int id)
    {
        //for (int i = 0; i < spoke_Tracker.Length; i++)
        //{
        //    if (spoke_Tracker[i].go_occupying != null)
        //    {
        //        if (spoke_Tracker[i].go_occupying.GetComponent<Tetra_MAJOR>().id == id)
        //        {
        //            spoke_Tracker[i].touchTimer = 2f;
        //            spoke_Tracker[i].go_occupying = null;
        //            print("OPEN HAPPENED");
        //        }
        //    }
        //}
    }

    void Update()
    {
        touchTimer -= Time.deltaTime;
        if (bGrabbed)
        {
            float dist;
            foreach (GameObject g in all_shape_instancer.getList())
            {
                Spoke_Tracker[] sp = g.GetComponent<Tetra_MAJOR>().spoke_Tracker;

                if (g != gameObject)
                {                
                    for (int i = 0; i < 4; i++)
                    {
                        dist = (sp[i].genSpoke.transform.position - transform.position).magnitude;
   
                        if (dist < Radius_of_Action && sp[i].spokeState == SpokeState.OPEN && touchTimer < 0)
                        {
                            print("Within Range of idx: " + i.ToString());

                            //disconnect
                            GetComponent<VRTK_FixedJointGrabAttach>().StopGrab(false);

                            //sp[i].spokeState = SpokeState.OCCUPIED;
                            sp[i].go_occupying = this.gameObject;
                            //t.state = Tetra_MAJOR.State.ATTACHED;

                            attachedSpoke = sp[i];
                            attachedTetra_go = gameObject;
                            start_Pos = transform.position;
                            end_Pos = sp[i].genSpoke.transform.position;
                            qStart = transform.rotation;
                            qEnd = sp[i].genSpoke.transform.rotation;
                            tAccum = 0;

                            StartCoroutine(Tween_Move_to_Spoke());
                        }
                    }
                }
            }
            //Switch off Touch Grab if in the VU
            //if (!bTouchGrabEnabled && hmd_Ctrl_Tracking.Get_Current_USER_Location() != HMD_Ctrl_Tracking.USER_LOCATION.IN_THE_VU)
            //{
            //    Enable_Touch_Grab();
            //    bTouchGrabEnabled = true;
            //}
            //else if (bTouchGrabEnabled && hmd_Ctrl_Tracking.Get_Current_USER_Location() == HMD_Ctrl_Tracking.USER_LOCATION.IN_THE_VU)
            //{
            //    Disable_Touch_Grab();
            //    bTouchGrabEnabled = false;
            //}
        }
    }

    private void Enable_Touch_Grab()
    {
        //GetComponent<SphereCollider>().enabled = true;
        //foreach (Tetra_MAJOR t in tetras)
        //{
        //    t.GetComponent<SphereCollider>().enabled = true;
        //}
    }

    private void Disable_Touch_Grab()
    {
        //GetComponent<SphereCollider>().enabled = false;
        //foreach (GameObject g in all_shape_instancer.getList())
        //{
        //    g.GetComponent<SphereCollider>().enabled = false;
        //}
    }

    private IEnumerator Tween_Move_to_Spoke()
    {
        int i = 0; //Sanity check for infinite loops
        float linProg = tAccum / tDurr;

        while (i < 180 && (linProg < 1))
        {
            bForce_update = true;
            tAccum += Time.deltaTime;
            linProg = tAccum / tDurr;
            Mathf.Clamp(tAccum, 0, tDurr);

            tween_Val = EaseMethods.QuintEaseOut(tAccum, 0, 1, tDurr);
            Move_To_Spoke(tween_Val);

            yield return true;
            i++;
        }

        tAccum = tDurr;
        tween_Val = 1;

        Move_To_Spoke(tween_Val);
        bForce_update = false;
        StopCoroutine("Tween_State");
    }

    private void Move_To_Spoke(float val)
    {
        transform.position = Vector3.Lerp(start_Pos, end_Pos, val);
        transform.rotation = Quaternion.Slerp(qStart, qEnd, val);
    }
}

 

Posted on

Voronoi Joint

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using VRTK;

//Set to Template option

public class Voronoi_Joint : MonoBehaviour
{
    public enum GEOMETRY
    {
        TETRA,
        OCTA,
    };

    public GEOMETRY geo = GEOMETRY.TETRA;

    public GameObject[] centerRbs = new GameObject[6];//4
    public GameObject[] connected_Joints = new GameObject[6];//4
    private Vector3[] cj_12_EndPoints = new Vector3[12];
    private Vector3[] cj_24_EndPoints = new Vector3[24];

    private Vector3 startPosition;

    //private int[,] vN_Points_to;

    private Six_Splines    six_splines;
    private Twelve_Splines twelve_Splines;

    public bool bPrint_Debug = false;

    public Button_VU bvToggle_Hi_Res;


    /// <image url="$(SolutionDir)\EMB\EMB_Edge_Vert.png" scale="0.36" /> 

    /// <image url="$(SolutionDir)\EMB\EMB_Edge_Vert_2.png" scale="0.1" /> 

    public float radius = 1;
    public float outerforceMag = 10;
    public float centering_Mult = 1;
    public float torque_Coefficient = 1f;
    private VRTK_InteractableObject vrtk_Interact_1;

    public enum Handle_States
    {
        NORMAL,
        TOUCHED,
        GRABBED,
    }

    private Handle_States handle_State = Handle_States.NORMAL;

    public Color handle_Base_Col;
    public Color handle_Touched_Col;
    public Color handle_Grabbed_Col;
    private Color handle_Base_Emissive_Col;

    private VU_UI_MANAGER VU_UI;

    public bool bUsePhysics = true;
    public bool bSetTetraPositions = false;
    public GameObject TetraVoronoi;

    void Start()
    {
        startPosition = transform.position;

        bvToggle_Hi_Res.isON = false;

        vrtk_Interact_1 = GetComponent<VRTK_InteractableObject>();
        VU_UI = VU_UI_MANAGER.Instance;

        if (geo == GEOMETRY.TETRA)
        {
            six_splines = GetComponentInChildren<Six_Splines>();
            six_splines.Set_CJ_0_GO(connected_Joints[0]);//For OutII
        }
        else if(geo == GEOMETRY.OCTA)
        {
            twelve_Splines = GetComponentInChildren<Twelve_Splines>();
        }

        vrtk_Interact_1.InteractableObjectTouched += new InteractableObjectEventHandler(DoObject_1_Touched);
        vrtk_Interact_1.InteractableObjectUntouched += new InteractableObjectEventHandler(DoObject_1_Untouched);

        vrtk_Interact_1.InteractableObjectGrabbed += new InteractableObjectEventHandler(DoObject_1_Grabbed);
        vrtk_Interact_1.InteractableObjectUngrabbed += new InteractableObjectEventHandler(DoObject_1_Ungrabbed);

        if (!bUsePhysics)
        {
            Rigidbody[] rbs = GetComponentsInChildren<Rigidbody>();
            foreach(Rigidbody rb in rbs)
            {
                rb.isKinematic = true;
            }
             
            for (int i = 0; i < centerRbs.Length; i++)
            {
                centerRbs[i].GetComponent<Rigidbody>().isKinematic = true;
                centerRbs[i].GetComponent<ConfigurableJoint>().enablePreprocessing = false;
                connected_Joints[i].GetComponent<Rigidbody>().isKinematic = true;
            }
        }
        if(bSetTetraPositions)
        {
            Voronoi_Joint vj = TetraVoronoi.GetComponent<Voronoi_Joint>();

            Vector3 thisBasePosition = transform.position;

            float edgeLength = (connected_Joints[0].transform.position - connected_Joints[1].transform.position).magnitude;

            Vector3 vfaceCenter = ( connected_Joints[0].transform.localPosition +
                                    connected_Joints[1].transform.localPosition +
                                    connected_Joints[2].transform.localPosition) / 3f;

            Vector3 vDir = vfaceCenter.normalized;

            float tetraHeight = edgeLength * Mathf.Sqrt(0.666666667f);

            TetraVoronoi.transform.position = transform.position + vfaceCenter + vDir * edgeLength * Mathf.Sqrt(6) / 12f;

            Vector3[] vTetraPos = new Vector3[4];

            vTetraPos[0] = connected_Joints[0].transform.position;
            vTetraPos[1] = connected_Joints[1].transform.position;
            vTetraPos[2] = connected_Joints[2].transform.position;
            vTetraPos[3] = thisBasePosition + vfaceCenter + vDir*edgeLength*Mathf.Sqrt(6)/3;

            vj.SetPositions(ref vTetraPos);
        }
    }

    public void SetPositions(ref Vector3[] positions)
    {
        for(int i = 0; i < connected_Joints.Length; i++)
        {
            connected_Joints[i].transform.position = positions[i];
        }
    }

    private void DoObject_1_Touched(object sender, InteractableObjectEventArgs e)
    {
        handle_State = Handle_States.TOUCHED;
        VU_UI.Report_Touched(vrtk_Interact_1.gameObject, e.interactingObject);
    }

    private void DoObject_1_Untouched(object sender, InteractableObjectEventArgs e)
    {
        handle_State = Handle_States.NORMAL;
        VU_UI.Report_UnTouched(vrtk_Interact_1.gameObject, e.interactingObject);
    }

    private void DoObject_1_Grabbed(object sender, InteractableObjectEventArgs e)
    {
        handle_State = Handle_States.GRABBED;
    }

    private void DoObject_1_Ungrabbed(object sender, InteractableObjectEventArgs e)
    {
        handle_State = Handle_States.NORMAL;
    }

    private bool isValid_Vector(Vector3 v)
    {
        if (!float.IsNaN(v.x) && !float.IsNaN(v.y) && !float.IsNaN(v.z))
        {
            return true;
        }
        else
        {
            return false;
        }
    }

    public void Return_To_Start_Position()
    {
        transform.position = startPosition;
    }

    void FixedUpdate()
    {
        if (bUsePhysics)
        {
            //Angular Adjustment
            for (int j = 0; j < connected_Joints.Length; j++)
            {
                Vector3 vA = connected_Joints[j].transform.right;
                GameObject gAlign = connected_Joints[0];

                //subtract y component to get projection on current XZ plane
                float yDot = Vector3.Dot(connected_Joints[j].transform.up, gAlign.transform.up);
                Vector3 yProjected = gAlign.transform.up + yDot * gAlign.transform.up;

                float torque = torque_Coefficient * Vector3.Dot(yProjected, connected_Joints[j].transform.forward);
                Rigidbody rbc = connected_Joints[j].GetComponent<Rigidbody>();
                rbc.AddTorque(rbc.transform.up * torque);
            }

            //Distribute around sphere - Voronoi Action
            foreach (GameObject go in connected_Joints)
            {
                Rigidbody rb = go.GetComponent<Rigidbody>();
                Rigidbody rbInner;
                Vector3 vFrom;
                Vector3 vSum;
                Vector3 vRecip;

                foreach (GameObject gInner in connected_Joints)
                {
                    vSum = Vector3.zero;
                    float r;//Radius from Neighbor

                    if (go != gInner)
                    {
                        rbInner = gInner.GetComponent<Rigidbody>();
                        vFrom = rb.position - rbInner.position;
                        r = vFrom.magnitude;
                        vRecip = (1f / r) * vFrom.normalized;
                        vSum += vRecip;
                    }
                    if (isValid_Vector(vSum))
                    {
                        rb.AddForce(outerforceMag * vSum);
                    }
                }
            }
            //Pull centerRbs into Center of Sphere
            Vector3 vDiff;
            foreach (GameObject gCent in centerRbs)
            {
                Rigidbody rb = gCent.GetComponent<Rigidbody>();
                vDiff = transform.position - rb.position;

                rb.velocity = vDiff / Time.fixedDeltaTime;
            }

            if (bPrint_Debug)
            {
                //for (int i = 0; i < 4; i++)
                //{
                //    for (int j = 0; j < 3; j++)
                //    {
                //        print("VN_[" + i.ToString() + "," + j.ToString() + "] points to " + vN_Points_to[i, j].ToString());
                //    }
                //}
                bPrint_Debug = false;
            }
        }
    }

    Color curr_Color;
    Color curr_dot_Color;
    Color targ_Color;
    Color targ_dot_Color;

    float theta = 390;
    float target_theta = 0;

    public void Animate_Handle_Material()//UI HANDLE
    {
        curr_Color = Color.Lerp(curr_Color, targ_Color, Time.deltaTime * 1.61803f);
        curr_dot_Color = Color.Lerp(curr_dot_Color, targ_dot_Color, Time.deltaTime * 1.61803f);

        if (handle_State == Handle_States.NORMAL)
        {
            //target_theta = 300;
            targ_Color = handle_Base_Col;

        }
        else if (handle_State == Handle_States.TOUCHED)
        {
            //target_theta = 390;
            targ_Color = handle_Touched_Col;
        }
        else if (handle_State == Handle_States.GRABBED)
        {
            //target_theta = 570;
            targ_Color = handle_Grabbed_Col;

        }
        Renderer rend = GetComponent<Renderer>();
        rend.material.SetColor("_Color", curr_Color);
        rend.material.SetColor("_EmissionColor", curr_Color);
    }

    public void Toggle_Hi_Res_From_Voronoi()
    {
        bvToggle_Hi_Res.isON = !bvToggle_Hi_Res.isON;
        six_splines = TetraVoronoi.GetComponentInChildren<Six_Splines>();
        six_splines.Set_Hi_Res(bvToggle_Hi_Res.isON);
        twelve_Splines.Set_Hi_Res(bvToggle_Hi_Res.isON);
    }

    private void Update()
    {
        //Animate_Handle_Material();
        if (geo == GEOMETRY.TETRA)
        {
            cj_12_EndPoints[0] = connected_Joints[0].transform.localPosition;
            cj_12_EndPoints[1] = connected_Joints[1].transform.localPosition;
            cj_12_EndPoints[2] = connected_Joints[0].transform.localPosition;
            cj_12_EndPoints[3] = connected_Joints[2].transform.localPosition;
            cj_12_EndPoints[4] = connected_Joints[0].transform.localPosition;
            cj_12_EndPoints[5] = connected_Joints[3].transform.localPosition;

            cj_12_EndPoints[6] = connected_Joints[1].transform.localPosition;
            cj_12_EndPoints[7] = connected_Joints[2].transform.localPosition;
            cj_12_EndPoints[8] = connected_Joints[1].transform.localPosition;
            cj_12_EndPoints[9] = connected_Joints[3].transform.localPosition;
            cj_12_EndPoints[10] = connected_Joints[2].transform.localPosition;
            cj_12_EndPoints[11] = connected_Joints[3].transform.localPosition;

            six_splines.Generate_Splines(cj_12_EndPoints, transform.position);
        }
        else if(geo == GEOMETRY.OCTA)
        {
            cj_24_EndPoints[0] = connected_Joints[0].transform.localPosition;
            cj_24_EndPoints[1] = connected_Joints[1].transform.localPosition;
            cj_24_EndPoints[2] = connected_Joints[0].transform.localPosition;
            cj_24_EndPoints[3] = connected_Joints[2].transform.localPosition;

            cj_24_EndPoints[4] = connected_Joints[0].transform.localPosition;
            cj_24_EndPoints[5] = connected_Joints[3].transform.localPosition;
            cj_24_EndPoints[6] = connected_Joints[0].transform.localPosition;
            cj_24_EndPoints[7] = connected_Joints[4].transform.localPosition;

            cj_24_EndPoints[8] = connected_Joints[1].transform.localPosition;
            cj_24_EndPoints[9] = connected_Joints[2].transform.localPosition;
            cj_24_EndPoints[10] = connected_Joints[1].transform.localPosition;
            cj_24_EndPoints[11] = connected_Joints[5].transform.localPosition;

            cj_24_EndPoints[12] = connected_Joints[1].transform.localPosition;
            cj_24_EndPoints[13] = connected_Joints[4].transform.localPosition;
            cj_24_EndPoints[14] = connected_Joints[2].transform.localPosition;
            cj_24_EndPoints[15] = connected_Joints[5].transform.localPosition;

            cj_24_EndPoints[16] = connected_Joints[2].transform.localPosition;
            cj_24_EndPoints[17] = connected_Joints[3].transform.localPosition;
            cj_24_EndPoints[18] = connected_Joints[3].transform.localPosition;
            cj_24_EndPoints[19] = connected_Joints[5].transform.localPosition;

            cj_24_EndPoints[20] = connected_Joints[3].transform.localPosition;
            cj_24_EndPoints[21] = connected_Joints[4].transform.localPosition;
            cj_24_EndPoints[22] = connected_Joints[4].transform.localPosition;
            cj_24_EndPoints[23] = connected_Joints[5].transform.localPosition;

            twelve_Splines.Generate_Splines(cj_24_EndPoints, transform.position);
        }
    }
}

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using ProceduralToolkit;
using VRTK;
using UnityEngine.UI;
using VRTK.GrabAttachMechanics;

public class Six_Splines : MonoBehaviour
{
    public Face_Tetra faces;
    public GameObject test_Octahedron;

    public GameObject Horn_Handle_Prefab;
    private GameObject[] tangent_Handles = new GameObject[14];
    [Range(0, 100f)]
    public  float scale_tangent_Handles = 27f;
    private bool  bEdge_Rotation_On = false;
    private float edge_Angular_Rate = 0f;
    public Slider sld_Edge_Angular_Rate;
    private bool bHandles_On = true;
    public Button_VU bvAngular_On;

    private GameObject go_cj_0;

    public Button_VU[] bvEdge_Switches;

    private bool[] bEdge_Switch = new bool[6];

    public Button_VU bvHorns_Handles_On;

    public Slider sld_Tangent;
    public Slider sld_Extend;

    public float tan_Slide_Val = 0;
    public float extend_Slide_Val = 0;

    public Slider_VU sv_Tangent_Adjust;
    public Slider_VU sv_Extend_Adjust;

    public bool bHi_Res_On = false;
    public bool bRun_Once = false;

    public Button_VU bv_Outii; //One Face Opposite Tangents

    public Vector3[] tetra_Edges = new Vector3[12]; //edge 0 -> 0,1 | edge 1 -> 2-3 etc.

    private VRTK_InteractableObject interact;

    private ArgosMeshDraft aMD_Spline_Green;
    private ArgosMeshDraft aMD_Spline_Red;
    private ArgosMeshDraft aMD_Spline_Blue;
    private ArgosMeshDraft aMD_Spline_Yellow;
    private ArgosMeshDraft aMD_Spline_Orange;
    private ArgosMeshDraft aMD_Spline_Violet;

    private ArgosMeshDraft aMD_Spline_Cross;

    public GameObject spline_Green;
    public GameObject spline_Red;
    public GameObject spline_Blue;
    public GameObject spline_Yellow;
    public GameObject spline_Orange;
    public GameObject spline_Violet;

    public GameObject spline_Cross;

    public float[] spline_lens_Green  = new float[256];
    public float[] spline_lens_Red    = new float[256];
    public float[] spline_lens_Blue   = new float[256];
    public float[] spline_lens_Yellow = new float[256];
    public float[] spline_lens_Orange = new float[256];
    public float[] spline_lens_Violet = new float[256];

    public float[] spline_lens_Cross = new float[256];

    private List<Vector3> lst_SplineGreen  = new List<Vector3>();
    private List<Vector3> lst_SplineRed    = new List<Vector3>();
    private List<Vector3> lst_SplineBlue   = new List<Vector3>();
    private List<Vector3> lst_SplineYellow = new List<Vector3>();
    private List<Vector3> lst_SplineOrange = new List<Vector3>();
    private List<Vector3> lst_SplineViolet = new List<Vector3>();

    //private List<Vector3> lst_SplineCross = new List<Vector3>();

    private Vector3 p0_1, p0_2;
    private Vector3 p1_1, p1_2;
    private Vector3 p2_1, p2_2;
    private Vector3 p3_1, p3_2;
    private Vector3 p4_1, p4_2;
    private Vector3 p5_1, p5_2;

    [Range(0, 0.5f)]
    public float radial_Mag = 0.5f;

    [Range(-20, 20f)]
    public float tangent_Mag = 1f;

    [Range(-20, 20f)]
    public float tRad = 0.0f;

    public bool bUpdate = false;
    public bool bTangents_From_Octahedron = false;
    public Twelve_Splines twelve_Splines;

    private VU_UI_MANAGER VU_UI;

    void Start ()
    {
        bv_Outii.isON = false;

        //tan_Slide_Val = sld_Tangent.value;
        //extend_Slide_Val = sld_Extend.value;
        //edge_Angular_Rate = sld_Edge_Angular_Rate.value;

        for (int i = 0; i<6; i++)
        {
            bEdge_Switch[i] = false;
        }

        VU_UI = VU_UI_MANAGER.Instance;

        aMD_Spline_Green  = new ArgosMeshDraft();
        aMD_Spline_Red    = new ArgosMeshDraft();
        aMD_Spline_Blue   = new ArgosMeshDraft();
        aMD_Spline_Yellow = new ArgosMeshDraft();
        aMD_Spline_Orange = new ArgosMeshDraft();
        aMD_Spline_Violet = new ArgosMeshDraft();

        aMD_Spline_Cross = new ArgosMeshDraft();

        //Generate_Splines();

        spline_Green.GetComponent<MeshFilter>().mesh.MarkDynamic(); 
        spline_Red.GetComponent<MeshFilter>().mesh.MarkDynamic();
        spline_Blue.GetComponent<MeshFilter>().mesh.MarkDynamic();
        spline_Yellow.GetComponent<MeshFilter>().mesh.MarkDynamic();
        spline_Orange.GetComponent<MeshFilter>().mesh.MarkDynamic();
        spline_Violet.GetComponent<MeshFilter>().mesh.MarkDynamic();

        spline_Cross.GetComponent<MeshFilter>().mesh.MarkDynamic();

        for (int i = 0; i < 12; i++)
        {
            tangent_Handles[i] = Instantiate(Horn_Handle_Prefab);
            tangent_Handles[i].transform.localPosition = Vector3.zero;
            tangent_Handles[i].transform.localRotation = Quaternion.identity;

            if (!GetComponentInParent<Voronoi_Joint>().bUsePhysics)
            {
                tangent_Handles[i].GetComponent<Rigidbody>().isKinematic = true;
            }
        }
        for (int i = 0; i < 12; i++)
        {
            if (i % 2 == 0)
            {
                tangent_Handles[i].GetComponent<Handle_Tangent_Interact>().Set_Dual(tangent_Handles[i + 1].GetComponent<Handle_Tangent_Interact>());
            }
            else
            {
                tangent_Handles[i].GetComponent<Handle_Tangent_Interact>().Set_Dual(tangent_Handles[i - 1].GetComponent<Handle_Tangent_Interact>());
            }
        }
    }

    public void Set_CJ_0_GO(GameObject cj_0_go)
    {
        go_cj_0 = cj_0_go;
    }

    public void Show_Handles(bool bON)
    {
        for (int i = 0; i < 14; i++)
        {
            tangent_Handles[i].SetActive(bON);
        }
    }

    //public Slider sld_Tangent;
    //public Slider sld_Extend;

    public void onSld_Tangent()
    {
        tan_Slide_Val = sld_Tangent.value;
        VU_UI.Set_Active_Slider(sv_Tangent_Adjust.gameObject);
        UI_Ladder.Instance.fine_Tuner.Active_FTS(sld_Tangent);
    }

    public void onSld_Extend()
    {
        extend_Slide_Val = sld_Extend.value;
        VU_UI.Set_Active_Slider(sv_Extend_Adjust.gameObject);
        UI_Ladder.Instance.fine_Tuner.Active_FTS(sld_Extend);
    }

    public void Set_Hi_Res(bool bON)
    {
        bHi_Res_On = bON;
        faces.Set_Hi_Res(bON);
        bRun_Once = bON;
    }

    public void onBV_Outii()
    {
        bv_Outii.isON = !bv_Outii.isON;
    }

    public void On_Edge_Switch_BV(int i)
    {
        bEdge_Switch[i] = !bEdge_Switch[i];
        if(bvEdge_Switches[0] != null)
        {
            bvEdge_Switches[i].isON = bEdge_Switch[i];
        }      
    }

    public void onSld_Cross_Tangent()
    {
        tan_Slide_Val = sld_Tangent.value;
        VU_UI.Set_Active_Slider(sv_Tangent_Adjust.gameObject);
        UI_Ladder.Instance.fine_Tuner.Active_FTS(sld_Tangent);
    }

    public void onSld_Cross_Extend()
    {
        extend_Slide_Val = sld_Extend.value;
        VU_UI.Set_Active_Slider(sv_Extend_Adjust.gameObject);
        UI_Ladder.Instance.fine_Tuner.Active_FTS(sld_Extend);
    }

    public void On_Edge_Angular_Rotation_sld()
    {
        edge_Angular_Rate = sld_Edge_Angular_Rate.value;
    }

    public void On_Edge_Rotation_Button()
    {
        bEdge_Rotation_On = !bEdge_Rotation_On;
        bvAngular_On.isON = bEdge_Rotation_On;
    }

    public void OnHandles_OnOff()
    {
        bHandles_On = !bHandles_On;
        bvHorns_Handles_On.isON = bHandles_On;

        for(int i = 0; i<12; i++)
        {
            tangent_Handles[i].GetComponent<Handle_Tangent_Interact>().Show(bHandles_On);
        }
    }

    float ang = 0;
    void Update()
    {
        if (bUpdate)
        {
            if (bEdge_Rotation_On)
            {
                ang += edge_Angular_Rate * Time.deltaTime;
            }
            else
            {
                ang = 0;
            }
            for (int i = 0; i < 12; i++)
            {
                tangent_Handles[i].transform.localScale = scale_tangent_Handles * Vector3.one;
                tangent_Handles[i].transform.position   = transform.position + transform.parent.TransformVector(tetra_Edges[i]);
            }
            Vector3 v_i1_i0;
            Quaternion q0, q1;
            Quaternion qZ_0, qZ_1; 
            for (int i = 0; i < 6; i++)
            {
                v_i1_i0 = tangent_Handles[i*2].transform.position - tangent_Handles[i*2+1].transform.position;

                q0 = Quaternion.LookRotation(-v_i1_i0, tangent_Handles[i * 2].transform.position - transform.position);
                q1 = Quaternion.LookRotation(v_i1_i0, tangent_Handles[i * 2 + 1].transform.position - transform.position);

                float direction = 1f;
                if(bEdge_Switch[i])
                {
                    direction = -1f;
                }

                float tanDir = 1f;

                Quaternion q0_bub = Quaternion.identity;
                Quaternion q1_bub = Quaternion.identity;

                qZ_0 = Quaternion.Euler(0, 0, direction*ang);
                qZ_1 = Quaternion.Euler(0, 0, -direction*ang);

                tangent_Handles[i * 2].transform.rotation = q0 * qZ_0 * Quaternion.Euler(new Vector3(tanDir*tan_Slide_Val, 0, 0))*q0_bub;
                tangent_Handles[i * 2 + 1].transform.rotation = q1 * qZ_1 * Quaternion.Euler(new Vector3(tanDir*tan_Slide_Val, 0, 0)) * q1_bub;

                //if (bv_Outii.isON && i > 2)
                //{
                //    if (i == 4)
                //    {

                //        //tanDir = -1f;
                //        Vector3 nYellow0 = Vector3.Cross(go_cj_0.transform.up, tangent_Handles[i * 2].transform.forward).normalized;
                //        Vector3 nYellow1 = Vector3.Cross(tangent_Handles[i * 2 + 1].transform.forward, go_cj_0.transform.up).normalized;

                //        Vector3 nPlanar0 = Vector3.Cross(go_cj_0.transform.up, nYellow0).normalized;
                //        Vector3 nPlanar1 = Vector3.Cross(go_cj_0.transform.up, nYellow1).normalized;

                //        float cosTheta0 = Vector3.Dot(nPlanar0, nYellow0);
                //        float cosTheta1 = Vector3.Dot(nPlanar1, nYellow1);

                //        float ang0_R = Mathf.Acos(cosTheta0);
                //        float ang1_R = Mathf.Acos(cosTheta1);

                //        float convR_to_D = 180f / Mathf.PI;
                //        float ang0_D = ang0_R * convR_to_D;
                //        float ang1_D = ang1_R * convR_to_D;

                //        q0_bub = Quaternion.AngleAxis(ang0_D, nYellow0);
                //        q1_bub = Quaternion.AngleAxis(2 * ang1_D, nYellow1);

                //        tangent_Handles[i * 2].transform.localRotation = q0_bub * tangent_Handles[i * 2].transform.localRotation;
                //        //tangent_Handles[i * 2 + 1].transform.rotation = q1_bub;
                //    }
                //}

                tangent_Handles[i * 2].GetComponent<Handle_Tangent_Interact>().Set_Tangent_Position(extend_Slide_Val);
                tangent_Handles[i * 2 + 1].GetComponent<Handle_Tangent_Interact>().Set_Tangent_Position(extend_Slide_Val);
            }       
        }
    }

    public void Generate_Splines(Vector3[] points, Vector3 vCenter)
    {
        if (!bHi_Res_On || bRun_Once)
        {
            bRun_Once = false;

            aMD_Spline_Green.Clear();
            aMD_Spline_Red.Clear();
            aMD_Spline_Blue.Clear();
            aMD_Spline_Yellow.Clear();
            aMD_Spline_Orange.Clear();
            aMD_Spline_Violet.Clear();

            //aMD_Spline_Cross.Clear();

            lst_SplineGreen.Clear();
            lst_SplineRed.Clear();
            lst_SplineBlue.Clear();
            lst_SplineYellow.Clear();
            lst_SplineOrange.Clear();
            lst_SplineViolet.Clear();

            //lst_SplineCross.Clear();

            for (int i = 0; i < 12; i++)
            {
                tetra_Edges[i] = points[i];
            }

            p0_1 = points[0];//0
            p0_2 = points[1];//1
            p1_1 = points[2];//0
            p1_2 = points[3];//2
            p2_1 = points[4];//0
            p2_2 = points[5];//3

            p3_1 = points[6];//1
            p3_2 = points[7];//2
            p4_1 = points[8];//1
            p4_2 = points[9];//3
            p5_1 = points[10];//2
            p5_2 = points[11];//3

            Vector3 vTn_01, vTn_02, vTn_11, vTn_12, vTn_31, vTn_32;

            if (bTangents_From_Octahedron)
            {
                //0-1
                vTn_01 = twelve_Splines.getTan_Han(0).GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_01 = transform.parent.InverseTransformVector(vTn_01);

                //1-0
                vTn_02 = twelve_Splines.getTan_Han(1).GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_02 = transform.parent.InverseTransformVector(vTn_02);

                //--------------------------
                //0-2
                vTn_11 = twelve_Splines.getTan_Han(2).GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_11 = transform.parent.InverseTransformVector(vTn_11);
                
                //2-0
                vTn_12 = twelve_Splines.getTan_Han(3).GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_12 = transform.parent.InverseTransformVector(vTn_12);

                //---------------------------
                //1-2
                vTn_31 = twelve_Splines.getTan_Han(8).GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_31 = transform.parent.InverseTransformVector(vTn_31);
                
                //2-1
                vTn_32 = twelve_Splines.getTan_Han(9).GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_32 = transform.parent.InverseTransformVector(vTn_32);
            }
            else
            {
                //0-1
                vTn_01 = tangent_Handles[0].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_01 = transform.parent.InverseTransformVector(vTn_01);

                //1-0
                vTn_02 = tangent_Handles[1].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_02 = transform.parent.InverseTransformVector(vTn_02);

                //--------------------------
                //0-2
                vTn_11 = tangent_Handles[2].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_11 = transform.parent.InverseTransformVector(vTn_11);
                //2-0
                vTn_12 = tangent_Handles[3].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_12 = transform.parent.InverseTransformVector(vTn_12);

                //---------------------------
                //1-2
                vTn_31 = tangent_Handles[6].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_31 = transform.parent.InverseTransformVector(vTn_31);
                //2-1
                vTn_32 = tangent_Handles[7].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
                vTn_32 = transform.parent.InverseTransformVector(vTn_32);
            }
            
            //---------------------------
            //0-3
            Vector3 vTn_21 = tangent_Handles[4].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
            vTn_21 = transform.parent.InverseTransformVector(vTn_21);
            //3-0
            Vector3 vTn_22 = tangent_Handles[5].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
            vTn_22 = transform.parent.InverseTransformVector(vTn_22);

            //---------------------------
            //1-3
            Vector3 vTn_41 = tangent_Handles[8].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
            vTn_41 = transform.parent.InverseTransformVector(vTn_41);
            //3-1
            Vector3 vTn_42 = tangent_Handles[9].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
            vTn_42 = transform.parent.InverseTransformVector(vTn_42);

            //---------------------------
            //2-3
            Vector3 vTn_51 = tangent_Handles[10].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
            vTn_51 = transform.parent.InverseTransformVector(vTn_51);
            //3-2
            Vector3 vTn_52 = tangent_Handles[11].GetComponent<Handle_Tangent_Interact>().Get_Tangent_Position() - vCenter;
            vTn_52 = transform.parent.InverseTransformVector(vTn_52);

            //Local Coordinates
            if (bv_Outii.isON)
            {
                //Mirror
                Vector3 vLocal_Up_0 = transform.parent.InverseTransformVector(go_cj_0.transform.up).normalized;

                Vector3 vHandle_Pos = tangent_Handles[6].gameObject.transform.position - vCenter;
                vHandle_Pos = transform.parent.InverseTransformVector(vHandle_Pos);

                Vector3 vTo_Tangent = vTn_31 - vHandle_Pos;

                float dot;
                dot = Vector3.Dot(vTo_Tangent, vLocal_Up_0);
                float mir = -2f * dot;
                vTn_31 += mir * vLocal_Up_0;
                vTn_32 += mir * vLocal_Up_0;
                vTn_41 += mir * vLocal_Up_0;
                vTn_42 += mir * vLocal_Up_0;
                vTn_51 += mir * vLocal_Up_0;
                vTn_52 += mir * vLocal_Up_0;
            }

            //---------------------------
            Compute_Dist_Spline(ref lst_SplineGreen, ref spline_lens_Green, p0_1, vTn_01, vTn_02, p0_2);
            Compute_Dist_Spline(ref lst_SplineRed, ref spline_lens_Red, p1_1, vTn_11, vTn_12, p1_2);
            Compute_Dist_Spline(ref lst_SplineBlue, ref spline_lens_Blue, p2_1, vTn_21, vTn_22, p2_2);

            //---------------------------
            Compute_Dist_Spline(ref lst_SplineYellow, ref spline_lens_Yellow, p3_1, vTn_31, vTn_32, p3_2);
            Compute_Dist_Spline(ref lst_SplineOrange, ref spline_lens_Orange, p4_1, vTn_41, vTn_42, p4_2);
            Compute_Dist_Spline(ref lst_SplineViolet, ref spline_lens_Violet, p5_1, vTn_51, vTn_52, p5_2);

            aMD_Spline_Green.Add(MeshDraft.Along_Spline(radial_Mag, 6, lst_SplineGreen));
            aMD_Spline_Red.Add(MeshDraft.Along_Spline(radial_Mag, 6, lst_SplineRed));
            aMD_Spline_Blue.Add(MeshDraft.Along_Spline(radial_Mag, 6, lst_SplineBlue));
            aMD_Spline_Yellow.Add(MeshDraft.Along_Spline(radial_Mag, 6, lst_SplineYellow));
            aMD_Spline_Orange.Add(MeshDraft.Along_Spline(radial_Mag, 6, lst_SplineOrange));
            aMD_Spline_Violet.Add(MeshDraft.Along_Spline(radial_Mag, 6, lst_SplineViolet));

            spline_Green.GetComponent<MeshFilter>().mesh = aMD_Spline_Green.ToMeshInternal();
            spline_Red.GetComponent<MeshFilter>().mesh = aMD_Spline_Red.ToMeshInternal();
            spline_Blue.GetComponent<MeshFilter>().mesh = aMD_Spline_Blue.ToMeshInternal();
            spline_Yellow.GetComponent<MeshFilter>().mesh = aMD_Spline_Yellow.ToMeshInternal();
            spline_Orange.GetComponent<MeshFilter>().mesh = aMD_Spline_Orange.ToMeshInternal();
            spline_Violet.GetComponent<MeshFilter>().mesh = aMD_Spline_Violet.ToMeshInternal();

            //p0 = p0_1
            //p1 = p1_1
            //p2 = p2_1

            //Tangents
            //  6 & 7 -> tan 1,2 tu1 & tu2
            //
            //  tv1 = 

            //cj_12_EndPoints[0] = connected_Joints[0].transform.localPosition;
            //cj_12_EndPoints[1] = connected_Joints[1].transform.localPosition;
            //cj_12_EndPoints[2] = connected_Joints[0].transform.localPosition;
            //cj_12_EndPoints[3] = connected_Joints[2].transform.localPosition;
            //cj_12_EndPoints[4] = connected_Joints[0].transform.localPosition;
            //cj_12_EndPoints[5] = connected_Joints[3].transform.localPosition;

            //cj_12_EndPoints[6] = connected_Joints[1].transform.localPosition;
            //cj_12_EndPoints[7] = connected_Joints[2].transform.localPosition;
            //cj_12_EndPoints[8] = connected_Joints[1].transform.localPosition;
            //cj_12_EndPoints[9] = connected_Joints[3].transform.localPosition;
            //cj_12_EndPoints[10] = connected_Joints[2].transform.localPosition;
            //cj_12_EndPoints[11] = connected_Joints[3].transform.localPosition;

            /*test_Octahedron.transform.localPosition = */

            faces.Clear_Faces();
            faces.Generate_Face(tetra_Edges[0],  tetra_Edges[3], tetra_Edges[1],  vTn_32, vTn_31,  vTn_12, vTn_02, vTn_11, vTn_01 );
            faces.Generate_Face(tetra_Edges[0],  tetra_Edges[5], tetra_Edges[3],  vTn_52, vTn_51,  vTn_22, vTn_12, vTn_21, vTn_11 );
            faces.Generate_Face(tetra_Edges[0],  tetra_Edges[1], tetra_Edges[5],  vTn_41, vTn_42,  vTn_02, vTn_22, vTn_01, vTn_21 );
            faces.Generate_Face(tetra_Edges[6],  tetra_Edges[7], tetra_Edges[9],  vTn_51, vTn_52,  vTn_32, vTn_42, vTn_31, vTn_41 );

            faces.toMesh_Internal();

            //Generate_Face(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 tu1, Vector3 tu2, Vector3 tv1, Vector3 tv2, Vector3 tv01, Vector3 tv02)
        }
        /// <image url="$(SolutionDir)\EMB\EMB_tetra_Faces.png" scale="0.15" /> 

        //aMD_Spline_Red.Add(MeshDraft.Along_Spline(0.007f, 6, lst_SplineRed));
        //aMD_Spline_Blue.Add(MeshDraft.Along_Spline(0.007f, 6, lst_SplineBlue));
    }

    private void Compute_Spline(ref List<Vector3> spline_lst, ref float[] lengths,  Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3)
    {
        float t = 0.0f;
        float dt = 1f / 72f;

        Vector3 vB = GetPointOnBezierCurve(p0, p1, p2, p3, get_Adjusted(t, ref lengths));

        for (int i = 0; i < 72; i++)
        {
            vB = GetPointOnBezierCurve(p0, p1, p2, p3, get_Adjusted(t, ref lengths));

            spline_lst.Add(vB);
            t += dt;
        }
    }

    private float get_Adjusted(float t, ref float[] lengths)
    {
        int i = 0;
        while (i < 256 && lengths[i] < t)
        {
            i++;
        }
        return (float)i / 256;
    }

    private void Compute_Dist_Spline(ref List<Vector3> spline_lst, ref float[] lengths, Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3)
    {
        float t = 0.0f;
        float dt = 1f / 256;

        Vector3 vB = GetPointOnBezierCurve(p0, p1, p2, p3, t);
        Vector3 vB_Last = vB;
        float running_Len = 0;

        for (int i = 0; i < 256; i++)
        {
            vB = GetPointOnBezierCurve(p0, p1, p2, p3, t);
            running_Len += (vB - vB_Last).magnitude;
            lengths[i] = running_Len;
            vB_Last = vB;
            t += dt;
        }
        for (int i = 0; i < 256; i++)
        {
            lengths[i] /= running_Len;
        }
        Compute_Spline(ref spline_lst, ref lengths, p0, p1, p2, p3);
    }

    public Vector3 GetPointOnBezierCurve(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 p3, float t)
    {
        float u = 1f - t;
        float t2 = t * t;
        float u2 = u * u;
        float u3 = u2 * u;
        float t3 = t2 * t;

        Vector3 result =
            (u3) * p0 +
            (3f * u2 * t) * p1 +
            (3f * u * t2) * p2 +
            (t3) * p3;

        return result;
    }
}