/*
setupMaterial.c

setup material index
*/

#include "oth.h"
#include "oth_prototype.h"

static void boundingbox(int, const double [], double *, double *, double *, double *, double *, double *);
static void copypbcHx(int, int, int, int);
static void copypbcHy(int, int, int, int);
static void copypbcHz(int, int, int, int);

void setupMaterial(int dofree)
{
	// free (MPI)
	if (dofree) {
		free(mE);
		free(mH);
	}

	// alloc and clear
	mE = (int *)malloc(NE * sizeof(int));
	mH = (int *)malloc(NH * sizeof(int));
	memset(mE, 0, NE * sizeof(int));
	memset(mH, 0, NH * sizeof(int));

	const double eps = EPS * sqrt(
		(Xn[Nx] - Xn[0]) * (Xn[Nx] - Xn[0]) +
		(Yn[Ny] - Yn[0]) * (Yn[Ny] - Yn[0]) +
		(Zn[Nz] - Zn[0]) * (Zn[Nz] - Zn[0]));

	for (int n = 0; n < NGeometry; n++) {
		int shape = Geometry[n].shape;
		int m     = Geometry[n].m;
		double *g = Geometry[n].g;
		int i;

		// bounding box
		int    i1, i2, j1, j2, k1, k2;
		double x1, x2, y1, y2, z1, z2;
		boundingbox(shape, g, &x1, &x2, &y1, &y2, &z1, &z2);

		// Ex
		getspan(Xc, Nx,     x1, x2, &i1, &i2, eps);
		getspan(Yn, Ny + 1, y1, y2, &j1, &j2, eps);
		getspan(Zn, Nz + 1, z1, z2, &k1, &k2, eps);
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (    i = i1; i <= i2; i++) {
		for (int j = j1; j <= j2; j++) {
		for (int k = k1; k <= k2; k++) {
			if (ingeometry(Xc[i], Yn[j], Zn[k], shape, g, eps)) {
				if ((i0Ex - 1 <= i) && (i <= i1Ex + 1) &&
				    (j0Ex - 1 <= j) && (j <= j1Ex + 1) &&
				    (k0Ex - 1 <= k) && (k <= k1Ex + 1)) {  // MPI
					MEX(i, j, k) = m;
				}
			}
		}
		}
		}

		// Ey
		getspan(Xn, Nx + 1, x1, x2, &i1, &i2, eps);
		getspan(Yc, Ny,     y1, y2, &j1, &j2, eps);
		getspan(Zn, Nz + 1, z1, z2, &k1, &k2, eps);
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (    i = i1; i <= i2; i++) {
		for (int j = j1; j <= j2; j++) {
		for (int k = k1; k <= k2; k++) {
			if (ingeometry(Xn[i], Yc[j], Zn[k], shape, g, eps)) {
				if ((i0Ey - 1 <= i) && (i <= i1Ey + 1) &&
				    (j0Ey - 1 <= j) && (j <= j1Ey + 1) &&
				    (k0Ey - 1 <= k) && (k <= k1Ey + 1)) {  // MPI
					MEY(i, j, k) = m;
				}
			}
		}
		}
		}

		// Ez
		getspan(Xn, Nx + 1, x1, x2, &i1, &i2, eps);
		getspan(Yn, Ny + 1, y1, y2, &j1, &j2, eps);
		getspan(Zc, Nz,     z1, z2, &k1, &k2, eps);
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (    i = i1; i <= i2; i++) {
		for (int j = j1; j <= j2; j++) {
		for (int k = k1; k <= k2; k++) {
			if (ingeometry(Xn[i], Yn[j], Zc[k], shape, g, eps)) {
				if ((i0Ez - 1 <= i) && (i <= i1Ez + 1) &&
				    (j0Ez - 1 <= j) && (j <= j1Ez + 1) &&
				    (k0Ez - 1 <= k) && (k <= k1Ez + 1)) {  // MPI
					MEZ(i, j, k) = m;
				}
			}
		}
		}
		}

		// Hx
		getspan(Xn, Nx + 1, x1, x2, &i1, &i2, eps);
		getspan(Yc, Ny,     y1, y2, &j1, &j2, eps);
		getspan(Zc, Nz,     z1, z2, &k1, &k2, eps);
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (    i = i1; i <= i2; i++) {
		for (int j = j1; j <= j2; j++) {
		for (int k = k1; k <= k2; k++) {
			if (ingeometry(Xn[i], Yc[j], Zc[k], shape, g, eps)) {
				if ((i0Hx <= i) && (i <= i1Hx) &&
				    (j0Hx <= j) && (j <= j1Hx) &&
				    (k0Hx <= k) && (k <= k1Hx)) {  // MPI
					MHX(i, j, k) = m;
					copypbcHx(i, j, k, m);
				}
			}
		}
		}
		}

		// Hy
		getspan(Xc, Nx,     x1, x2, &i1, &i2, eps);
		getspan(Yn, Ny + 1, y1, y2, &j1, &j2, eps);
		getspan(Zc, Nz,     z1, z2, &k1, &k2, eps);
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (    i = i1; i <= i2; i++) {
		for (int j = j1; j <= j2; j++) {
		for (int k = k1; k <= k2; k++) {
			if (ingeometry(Xc[i], Yn[j], Zc[k], shape, g, eps)) {
				if ((i0Hy <= i) && (i <= i1Hy) &&
				    (j0Hy <= j) && (j <= j1Hy) &&
				    (k0Hy <= k) && (k <= k1Hy)) {  // MPI
					MHY(i, j, k) = m;
					copypbcHy(i, j, k, m);
				}
			}
		}
		}
		}

		// Hz
		getspan(Xc, Nx,     x1, x2, &i1, &i2, eps);
		getspan(Yc, Ny,     y1, y2, &j1, &j2, eps);
		getspan(Zn, Nz + 1, z1, z2, &k1, &k2, eps);
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (    i = i1; i <= i2; i++) {
		for (int j = j1; j <= j2; j++) {
		for (int k = k1; k <= k2; k++) {
			if (ingeometry(Xc[i], Yc[j], Zn[k], shape, g, eps)) {
				if ((i0Hz <= i) && (i <= i1Hz) &&
				    (j0Hz <= j) && (j <= j1Hz) &&
				    (k0Hz <= k) && (k <= k1Hz)) {  // MPI
					MHZ(i, j, k) = m;
					copypbcHz(i, j, k, m);
				}
			}
		}
		}
		}
	}

	// feed -> PEC
	if (NFeed) {
		for (int n = 0; n < NFeed; n++) {
			int i = Feed[n].i;
			int j = Feed[n].j;
			int k = Feed[n].k;
			if      ((Feed[n].dir == 'X') &&
				(i >= i0Ex) && (i <= i1Ex) &&
				(j >= j0Ex) && (j <= j1Ex) &&
				(k >= k0Ex) && (k <= k1Ex)) {
				MEX(i, j, k) = PEC;
			}
			else if ((Feed[n].dir == 'Y') &&
				(i >= i0Ey) && (i <= i1Ey) &&
				(j >= j0Ey) && (j <= j1Ey) &&
				(k >= k0Ey) && (k <= k1Ey)) {
				MEY(i, j, k) = PEC;
			}
			else if ((Feed[n].dir == 'Z') &&
				(i >= i0Ez) && (i <= i1Ez) &&
				(j >= j0Ez) && (j <= j1Ez) &&
				(k >= k0Ez) && (k <= k1Ez)) {
				MEZ(i, j, k) = PEC;
			}
			//printf("%d %c %d %d %d\n", n, Feed[n].dir, i, j, k);
		}
	}
}


// get bounding box
static void boundingbox(int shape, const double g[],
	double *x1, double *x2, double *y1, double *y2, double *z1, double *z2)
{
	*x1 = *x2 = *y1 = *y2 = *z1 = *z2 = 0;

	// cube / sphere / cylinder
	if ((shape == 1) || (shape == 2) ||
	    (shape == 11) || (shape == 12) || (shape == 13)) {
		*x1 = MIN(g[0], g[1]);
		*x2 = MAX(g[0], g[1]);
		*y1 = MIN(g[2], g[3]);
		*y2 = MAX(g[2], g[3]);
		*z1 = MIN(g[4], g[5]);
		*z2 = MAX(g[4], g[5]);
	}
	// triangle pillar
	else if (shape == 31) {
		*x1 = MIN(g[0], g[1]);
		*x2 = MAX(g[0], g[1]);
		*y1 = MIN(g[2], MIN(g[3], g[4]));
		*y2 = MAX(g[2], MAX(g[3], g[4]));
		*z1 = MIN(g[5], MIN(g[6], g[7]));
		*z2 = MAX(g[5], MAX(g[6], g[7]));
	}
	else if (shape == 32) {
		*y1 = MIN(g[0], g[1]);
		*y2 = MAX(g[0], g[1]);
		*z1 = MIN(g[2], MIN(g[3], g[4]));
		*z2 = MAX(g[2], MAX(g[3], g[4]));
		*x1 = MIN(g[5], MIN(g[6], g[7]));
		*x2 = MAX(g[5], MAX(g[6], g[7]));
	}
	else if (shape == 33) {
		*z1 = MIN(g[0], g[1]);
		*z2 = MAX(g[0], g[1]);
		*x1 = MIN(g[2], MIN(g[3], g[4]));
		*x2 = MAX(g[2], MAX(g[3], g[4]));
		*y1 = MIN(g[5], MIN(g[6], g[7]));
		*y2 = MAX(g[5], MAX(g[6], g[7]));
	}
	// pyramid / cone
	else if ((shape == 41) || (shape == 51)) {
		const double hy = MAX(g[4], g[6]) / 2;
		const double hz = MAX(g[5], g[7]) / 2;
		*x1 = MIN(g[0], g[1]);
		*x2 = MAX(g[0], g[1]);
		*y1 = g[2] - hy;
		*y2 = g[2] + hy;
		*z1 = g[3] - hz;
		*z2 = g[3] + hz;
	}
	else if ((shape == 42) || (shape == 52)) {
		const double hz = MAX(g[4], g[6]) / 2;
		const double hx = MAX(g[5], g[7]) / 2;
		*y1 = MIN(g[0], g[1]);
		*y2 = MAX(g[0], g[1]);
		*z1 = g[2] - hz;
		*z2 = g[2] + hz;
		*x1 = g[3] - hx;
		*x2 = g[3] + hx;
	}
	else if ((shape == 43) || (shape == 53)) {
		const double hx = MAX(g[4], g[6]) / 2;
		const double hy = MAX(g[5], g[7]) / 2;
		*z1 = MIN(g[0], g[1]);
		*z2 = MAX(g[0], g[1]);
		*x1 = g[2] - hx;
		*x2 = g[2] + hx;
		*y1 = g[3] - hy;
		*y2 = g[3] + hy;
	}
}


// copy PBC
static void copypbcHx(int i, int j, int k, int m)
{
	if (Pbc.ypbc) {
		if      ((j == 0     ) && (Ny <= j1Hx)) {
			MHX(i, Ny, k) = m;
		}
		else if ((j == Ny - 1) && (-1 >= j0Hx)) {
			MHX(i, -1, k) = m;
		}
	}
	if (Pbc.zpbc) {
		if      ((k == 0     ) && (Nz <= k1Hx)) {
			MHX(i, i, Nz) = m;
		}
		else if ((k == Nz - 1) && (-1 >= k0Hx)) {
			MHX(i, j, -1) = m;
		}
	}
}
static void copypbcHy(int i, int j, int k, int m)
{
	if (Pbc.zpbc) {
		if      ((k == 0     ) && (Nz <= k1Hy)) {
			MHY(i, j, Nz) = m;
		}
		else if ((k == Nz - 1) && (-1 >= k0Hy)) {
			MHY(i, j, -1) = m;
		}
	}
	if (Pbc.xpbc) {
		if      ((i == 0     ) && (Nx <= i1Hy)) {
			MHY(Nx, j, k) = m;
		}
		else if ((i == Nx - 1) && (-1 >= i0Hy)) {
			MHY(-1, j, k) = m;
		}
	}
}
static void copypbcHz(int i, int j, int k, int m)
{
	if (Pbc.xpbc) {
		if      ((i == 0     ) && (Nx <= i1Hz)) {
			MHZ(Nx, j, k) = m;
		}
		else if ((i == Nx - 1) && (-1 >= i0Hz)) {
			MHZ(-1, j, k) = m;
		}
	}
	if (Pbc.ypbc) {
		if      ((j == 0     ) && (Ny <= j1Hz)) {
			MHZ(i, Ny, k) = m;
		}
		else if ((j == Ny - 1) && (-1 >= j0Hz)) {
			MHZ(i, -1, k) = m;
		}
	}
}

/*
// debug
void debugMaterial(void)
{

	for (int k = k0Ex; k <= k1Ex; k++) {
		printf("Ex k=%d\n", k);
		for (int j = j1Ex; j >= j0Ex; j--) {
			for (int i = i0Ex; i <= i1Ex; i++) {
				printf("%d", MEX(i, j, k));
			}
			printf("\n");
		}
	}
	fflush(stdout);

	for (int k = k0Ey; k <= k1Ey; k++) {
		printf("Ey k=%d\n", k);
		for (int j = j1Ey; j >= j0Ey; j--) {
			for (int i = i0Ey; i <= i1Ey; i++) {
				printf("%d", MEY(i, j, k));
			}
			printf("\n");
		}
	}
	fflush(stdout);

	for (int k = k0Ez; k <= k1Ez; k++) {
		printf("Ez k=%d\n", k);
		for (int j = j1Ez; j >= j0Ez; j--) {
			for (int i = i0Ez; i <= i1Ez; i++) {
				printf("%d", MEZ(i, j, k));
			}
			printf("\n");
		}
	}
	fflush(stdout);

	for (int k = k0Hx; k <= k1Hx; k++) {
		printf("Hx k=%d\n", k);
		for (int j = j1Hx; j >= j0Hx; j--) {
			for (int i = i0Hx; i <= i1Hx; i++) {
				printf("%d", MHX(i, j, k));
			}
			printf("\n");
		}
	}
	fflush(stdout);

	for (int k = k0Hy; k <= k1Hy; k++) {
		printf("Hy k=%d\n", k);
		for (int j = j1Hy; j >= j0Hy; j--) {
			for (int i = i0Hy; i <= i1Hy; i++) {
				printf("%d", MHY(i, j, k));
			}
			printf("\n");
		}
	}
	fflush(stdout);

	for (int k = k0Hz; k <= k1Hz; k++) {
		printf("Hz k=%d\n", k);
		for (int j = j1Hz; j >= j0Hz; j--) {
			for (int i = i0Hz; i <= i1Hz; i++) {
				printf("%d", MHZ(i, j, k));
			}
			printf("\n");
		}
	}
	fflush(stdout);

}
*/
