/**
* @file    GeneHeterogeneity.h
* @author  Haoyu Tang
* @date    2016/08
* @brief   Generate fields of parameters/properties such as mechanical props and flow props for Hydraulic Fracturing
* @version
* FracTHM
*/

#include "GeneHeterogeneity.h"

namespace FracTHM{
	using namespace std;
	Hetero_field::Hetero_field(bool isGeneNewHeterField, bool isRestart, const int nx, const int ny, const int nz,
		                       double* Dx, double* Dy, double* Dz, int*NDx, int* NDy, int* NDz){
		AllocateArray(nx, ny, nz, Dx, Dy, Dz, NDx, NDy, NDz);
		GenerateRandom(isGeneNewHeterField, isRestart);
		PostGenerateRandom();
		InsituStress();
		MANNIE2();
		EmpiricalRelation();
		StanfordKozenyCarmen();
		MohrCoulomb();
		BrittlenessIndex();
	}

	void Hetero_field::AllocateArray(const int _nx, const int _ny, const int _nz,
		                             double* _Dx, double* _Dy, double* _Dz, int* _NDx, int* _NDy, int* _NDz){
		nx = _nx;
		ny = _ny;
		nz = _nz;
		NDx = new int[nx];
		NDy = new int[ny];
		NDz = new int[nz];
		Dx = new double[nx];
		Dy = new double[ny];
		Dz = new double[nz];
		
		int i, j, k;
		Nx = 0; Ny = 0; Nz = 0;

		for (k = 0; k < nz; k++) {
			Dz[k] = _Dz[k];
			NDz[k] = _NDz[k];
			Nz += NDz[k];
		}
		for (j = 0; j < ny; j++) {
			Dy[j] = _Dy[j];
			NDy[j] = _NDy[j];
			Ny += NDy[j];
		}
		for (i = 0; i < nx; i++) {
			Dx[i] = _Dx[i];
			NDx[i] = _NDx[i];
			Nx += NDx[i];
		}
		
		Dx_min = 100.0;
		Dy_min = 100.0;
		Dz_min = 100.0;


		for (k = 0; k<nz; k++)
		if (Dz[k] < Dz_min) Dz_min = Dz[k];
		for (j = 0; j<ny; j++)
		if (Dy[j] < Dy_min) Dy_min = Dy[j];
		for (i = 0; i < nx; i++)
		if (Dx[i] < Dx_min) Dx_min = Dx[i];


		x_max = 0;
		y_max = 0;
		z_max = 0;
		for (k = 0; k < nz; k++)
			z_max += NDz[k] * (int)(Dz[k] / Dz_min);
		for (j = 0; j<ny; j++)
			y_max += NDy[j] * (int)(Dy[j] / Dy_min);
		for (i = 0; i < nx; i++)
			x_max += NDx[i] * (int)(Dx[i] / Dx_min);
        
		//////////////////////////////////////////////
		size = Nx*Ny*Nz;
		Vp = new double[size];
		Vs = new double[size];
		phi = new double[size];
		perm = new double[size];
		T0 = new double[size];
		C = new double[size];
		miu_i = new double[size];
		Ev = new double[size];
		Eh = new double[size];
		niuv = new double[size];
		niuh = new double[size];
		niuhv = new double[size];
		UCS = new double[size];
		Pp = new double[size];
		Sv = new double[size];
		Sh = new double[size];
		SH = new double[size];
		C11 = new double[size];
		C12 = new double[size];
		C13 = new double[size];
		C33 = new double[size];
		C44 = new double[size];
		C66 = new double[size];
		K = new double[size];
		

		BI_TSUCS1 = new double[size];
		BI_TSUCS2 = new double[size];
		BI_TSUCS3 = new double[size];
		BI_TSUCS4 = new double[size];
		BI_YMPR = new double[size];
		BI_YMM = new double[size];

		R = new double[x_max*y_max*z_max];
	}

Hetero_field::~Hetero_field(){
	delete[] Vp; Vp = NULL;
	delete[] Vs; Vs = NULL;
	delete[] phi; phi = NULL;
	delete[] perm; perm = NULL;
	delete[] T0; T0 = NULL;
	delete[] C; C = NULL;
	delete[] miu_i; miu_i = NULL;
	delete[] Ev; Ev = NULL;
	delete[] Eh; Eh = NULL;
	delete[] niuv; niuv = NULL;
	delete[] niuh; niuh = NULL;
	delete[] niuhv; niuhv = NULL;
	delete[] UCS; UCS = NULL;
	delete[] Pp; Pp = NULL;
	delete[] Sv; Sv = NULL;
	delete[] Sh; Sh = NULL;
	delete[] SH; SH = NULL;
	delete[] C11; C11 = NULL;
	delete[] C12; C12 = NULL;
	delete[] C13; C13 = NULL;
	delete[] C33; C33 = NULL;
	delete[] C44; C44 = NULL;
	delete[] C66; C66 = NULL;
	delete[] K; K = NULL;
	delete[] R; R = NULL;

	delete[] BI_TSUCS1; BI_TSUCS1 = NULL;
	delete[] BI_TSUCS2; BI_TSUCS2 = NULL;
	delete[] BI_TSUCS3; BI_TSUCS3 = NULL;
	delete[] BI_TSUCS4; BI_TSUCS4 = NULL;
	delete[] BI_YMPR; BI_YMPR = NULL;
	delete[] BI_YMM; BI_YMM = NULL;

	delete[] NDx; NDx = NULL;
	delete[] NDy; NDy = NULL;
	delete[] NDz; NDz = NULL;
	delete[] Dx; Dx = NULL;
	delete[] Dy; Dy = NULL;
	delete[] Dz; Dz = NULL;
}

void Hetero_field::PostGenerateRandom()
{
	int i, j, k;

	int a, b, c;
	a = 0;
	b = 0;
	c = 0;
	for (k = 0; k<z_max; k++)
	for (j = 0; j<y_max; j++)
	for (i = 0; i<x_max; i++)
	{
		if (R[i + j*x_max + k*x_max*y_max] != 0)
		{
			K[a + b*Nx + c*Nx*Ny] = R[i + j*x_max + k*x_max*y_max];
			Vp[a + b*Nx + c*Nx*Ny] = 3600 + 150 * K[a + b*Nx + c*Nx*Ny];
			Vs[a + b*Nx + c*Nx*Ny] = Vp[a + b*Nx + c*Nx*Ny] / 1.88 + 133;
			a++;
			if (a == Nx){
				b++; a = 0;
			}
			if (b == Ny){
				c++; b = 0;
			}
		}
		//printf("%lf",Vp[i + j*Nx + k*Nx*Ny]);
	}

}

void Hetero_field::GenerateRandom(bool isGeneNewHeterField, bool isRestart)
{
	if (isRestart || !isGeneNewHeterField){
		cout << "\n*************************************************************************************\n";
		cout << " We are using RESTART option OR, Do not need generate NEW heterogeneious fields. SO,\n";
		cout << "        The data file \"case1.dat\" used before should be placed in this folder \n";
		cout << "*************************************************************************************\n";
		Sleep(5 * 1000);
	}
	else{
		//*******Write the parameter file for sgsim.exe
		ofstream OUTPUT("OK", ios::out);
		if (!OUTPUT){
			cerr << "open error!" << endl;
			exit(1);
		}
		srand((int)time(NULL));
		char Char_i[8] = { '0' };
		for (int i = 0; i < 8; i++){
			Char_i[i] = rand() % 10 + 48;
		}
		OUTPUT << "                  Parameters for SGSIM\n"
			<< "                  ********************\n"
			<< "\n"
			<< "START OF PARAMETERS:\n"
			<< "nodata                        \\file with conditioning data\n"
			<< "0  0  0  0  0  0              \\columns for X,Y,Z,vr,wt\n"
			<< "-10.0       10.0              \\trimming limits\n"
			<< "0                             \transform the data (0=no, 1=yes)\n"
			<< "sgsim.trn                     \\file for output trans table\n"
			<< "0                             \\consider ref. dist (0=no, 1=yes)\n"
			<< "histsmth.out                  \\file with ref. dist distribution\n"
			<< "1  3                          \\columns for vr and wt in the file\n"
			<< "0.0001   100000.0             \\zmin,zmax(tail extrapolation)\n"
			<< "1       0.0001                \\lower tail option, parameter \n"
			<< "4       10000.0               \\upper tail option, parameter\n"
			<< "1                             \\debugging level: 0,1,2,3\n"
			<< "case1.dbg                     \\file for debugging output\n"
			<< "case1.dat                     \\file for simulation output\n"
			<< "1                             \\number of realizations to generate\n"
			<< left << setw(4) << x_max << "   0.00 20                \\nx,xmn,xsiz\n"
			<< left << setw(4) << y_max << "   0.00 20                \\ny, ymn, ysiz\n"
			<< left << setw(4) << z_max << "    0.0 10                \\nz,zmn,zsiz\n"
			<< Char_i << "                      \\random number seed\n"  //01812163//<< "01812163                      \\random number seed\n"
			<< "0     16                      \\min and max original data for sim\n"
			<< "20                            \\number of simulated nodes to use\n"
			<< "1                             \\assign data to nodes (0=no, 1=yes) \n"
			<< "0     3                       \\multiple grid search (0=no, 1=yes),num \n"
			<< "16                            \\maximum data per octant (0=not used) \n"
			<< "3000.0   3000.0    300.0              \\maximum search radii (hmax,hmin,vert)\n"
			<< "0.0   0.0   0.0               \\angles for search ellipsoid \n"
			<< "0     0.60   1.0              \\ktype: 0=SK,1=OK,2=LVM,3=EXDR,4=COLC \n"
			<< "nofile                        \\file with LVM, EXDR, or COLC variable\n"
			<< "0                             \\column for secondary variable\n"
			<< "1    0.0                      \\nst, nugget effect\n"
			<< "2    1.0  0.0   0.0   0.0     \\it,cc(variance),ang1,ang2,ang3  \n"
			<< "          600.0   600.0   100.0     \\a_hmax, a_hmin, a_vert \n"
			<< " \n"
			<< " \n"
			<< "Note: it = 2 --> exponential covariance model \n"
			<< "      a_hmax = 6.0: range of correlation scale is equivalent to integral scale of 2.0\n"
			<< "                   scale\n"
			<< " \n"
			<< "value = ln K \n"
			<< " loop\n"
			<< "\n"
			<< "k = 1, Nz\n"
			<< "j = 1, Ny\n"
			<< "i = 1, Nx\n\n";
		OUTPUT.close();
		//////////////////////////////////////////////////////////////////////////////////
		WinExec("sgsim.exe", SW_SHOW);
		Sleep(5 * 100);
		cout << "\nPlease input \"OK\" AND then press \"Enter\" to continue! \n";
		Sleep(8 * 1000);

		//////////////////////////////////////////////////////
		int round = 0;										//
		while (round<60){ //Wait 2 min totally here			//
			ofstream OpenFile("case1.dat", ios::in); 		//
			if (OpenFile) {									//
				OpenFile.close();							//
				break;										//
			}												// Wait the subroutine WinExec() running successfully 
			else {											//
				OpenFile.close();							//
				round++;									//
				Sleep(5 * 1000);  //Wait 5 sec				//
			}												//
		}													//
		//////////////////////////////////////////////////////

	}

	int i, j, k, l, t;
	int sum, sum1;
	ifstream ifp;
	char    cbuffer[100];
	int temp = 0;
	ifp.open("case1.dat", ios::in);
	if (!ifp){
		cerr << "open error!" << endl;
		exit(1);
	}
	ifp.getline(cbuffer, 1000, '\n');
	ifp.getline(cbuffer, 1000, '\n');
	ifp.getline(cbuffer, 1000, '\n');
	for (k = 0; k<z_max; k++)
	for (j = 0; j<y_max; j++)
	for (i = 0; i < x_max; i++){
		int count = i + j*x_max + k*x_max*y_max;
		double tempD;
		ifp >> cbuffer;
		tempD = atof(cbuffer);
		if (abs(tempD) < 1.0e-3)
			R[count] = 1.0e-3; 
		else
			R[count] = tempD;
	}
	ifp.close();

	sum = 0;
	for (l = 0; l < nx; l++)
	{
		sum1 = 0;
		for (k = 0; k < z_max; k++)
		for (j = 0; j < y_max; j++)
		{
			for (i = 0; i < NDx[l]; i++)
			{
				sum1 = (int)(Dx[l] / Dx_min)*i;
				for (t = 1; t < (int)(Dx[l] / Dx_min); t++)
				{
					R[(sum + sum1) + j * x_max + k*x_max*y_max] += R[(sum + sum1 + t) + j*x_max + k*x_max*y_max];
					R[(sum + sum1 + t) + j*x_max + k*x_max*y_max] = 0;
				}
				R[(sum + sum1) + j * x_max + k*x_max*y_max] = R[(sum + sum1) + j * x_max + k*x_max*y_max] / ((int)(Dx[l] / Dx_min));
			}
		}
		sum += NDx[l] * ((int)(Dx[l] / Dx_min));
	}
	sum = 0;
	for (l = 0; l < ny; l++)
	{
		sum1 = 0;
		for (k = 0; k < z_max; k++)
		for (j = 0; j < x_max; j++)
		{
			for (i = 0; i < NDy[l]; i++)
			{
				sum1 = (int)(Dy[l] / Dy_min)*i;
				for (t = 1; t < (int)(Dy[l] / Dy_min); t++)
				{
					R[(sum + sum1)*x_max + j + k*x_max*y_max] += R[(sum + sum1 + t)*x_max + j + k*x_max*y_max];
					R[(sum + sum1 + t)*x_max + j + k*x_max*y_max] = 0;
				}
				R[(sum + sum1) * x_max + j + k*x_max*y_max] = R[(sum + sum1)* x_max + j + k*x_max*y_max] / ((int)(Dy[l] / Dy_min));
			}
		}
		sum += NDy[l] * ((int)(Dy[l] / Dy_min));
	}

	sum = 0;
	for (l = 0; l < nz; l++)
	{
		sum1 = 0;
		for (k = 0; k < x_max; k++)
		for (j = 0; j < y_max; j++)
		{
			for (i = 0; i < NDz[l]; i++)
			{
				sum1 = (int)(Dz[l] / Dz_min)*i;
				for (t = 1; t < (int)(Dz[l] / Dz_min); t++)
				{
					R[(sum + sum1)*x_max*y_max + j * x_max + k] += R[(sum + sum1 + t)*x_max*y_max + j*x_max + k];
					R[(sum + sum1 + t)*x_max*y_max + j*x_max + k] = 0;
				}
				R[(sum + sum1)*x_max*y_max + j * x_max + k] = R[(sum + sum1)*x_max*y_max + j * x_max + k] / ((int)(Dz[l] / Dz_min));
			}
		}
		sum += NDz[l] * ((int)(Dz[l] / Dz_min));
	}
}

void Hetero_field::InsituStress()
{
	double d1 = 6000,  d2 = 7000,  Pp_d = 0.48,  SH_d = 1.1,  Sh_d = 0.73,  Sv_d = 0.65;
	int i, j, k;
	for (k = 0; k < Nz; k++)
	{
		for (j = 0; j < Ny; j++)
		{
			for (i = 0; i < Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				Pp[iEle] = Pp_d*((d2 - d1) / Nz*(k + 0.5) + d1);
				SH[iEle] = SH_d*((d2 - d1) / Nz*(k + 0.5) + d1);
				Sh[iEle] = Sh_d*((d2 - d1) / Nz*(k + 0.5) + d1);
				Sv[iEle] = Sv_d*((d2 - d1) / Nz*(k + 0.5) + d1);
			}

		}

	}
}

void Hetero_field::MANNIE2()
{
	int i, j, k;
	double k1, k2, k3;
	double rou = 2.6;

	k1 = 1.12;
	k2 = 0.76;
	k3 = 0.9283;
	for (k = 0; k<Nz; k++)
	{
		for (j = 0; j<Ny; j++)
		{
			for (i = 0; i<Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
		
				C33[iEle] = rou*Vp[iEle] * Vp[iEle];
				C44[iEle] = rou*Vs[iEle] * Vs[iEle];
				C66[iEle] = (2 * k1*C44[iEle] / C33[iEle] + 1 - k1 - k3) / (2 * k1 / C33[iEle] - k3 / C44[iEle]);
				C11[iEle] = k1*(2 * C66[iEle] - 2 * C44[iEle] + C33[iEle]);
				C12[iEle] = C11[iEle] - 2 * C66[iEle];
				C13[iEle] = k2*C12[iEle];
				
				Ev[iEle] = C33[iEle] - 2 * C13[iEle] * C13[iEle] / (C11[iEle] + C12[iEle]);
				Eh[iEle] = C11[iEle] + (C13[iEle] * C13[iEle] * (C12[iEle] - C11[iEle]) + C12[iEle] * (C13[iEle] * C13[iEle] - C12[iEle] * C33[iEle])) / (C33[iEle] * C11[iEle] - C13[iEle] * C13[iEle]);
				niuv[iEle] = C13[iEle] / (C11[iEle] + C12[iEle]);
				niuh[iEle] = (C33[iEle] * C12[iEle] - C13[iEle] * C13[iEle]) / (C33[iEle] * C11[iEle] - C13[iEle] * C13[iEle]);
				niuhv[iEle] = Eh[iEle] / Ev[iEle] * niuv[iEle];
			}

		}

	}
}

void Hetero_field::EmpiricalRelation()

{
	int i, j, k;
	double k1;
	double multiplier1, multiplier2, multiplier3;
	double K_max, K_min;
	multiplier1 = 0.5*0.08;
	multiplier2 = 0.6;
	K_max = 0;
	K_min = 1;

	for (k = 0; k < Nz; k++)
	{
		for (j = 0; j < Ny; j++)
		{
			for (i = 0; i < Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				phi[iEle] = 2.85 * 1000000000 * pow(Vp[iEle], -2.95);
				UCS[iEle] = multiplier1*0.03*exp(0.0025*Vp[iEle]);

				C[iEle] = multiplier2*50.88*exp(-21.28*phi[iEle]);

				if (K_max < K[iEle]) K_max = K[iEle];
				if (K_min > K[iEle]) K_min = K[iEle];

			}
		}
	}
	for (k = 0; k < Nz; k++)
	for (j = 0; j < Ny; j++)
	for (i = 0; i < Nx; i++)
	{
		int iEle = i + j*Nx + k*Nx*Ny;
		UCS[iEle] = UCS[iEle] / 2 + 20;

		T0[iEle] = 3*(((0.125 - 0.1)*(K[iEle] - K_min) / (K_max - K_min) + 0.1)*UCS[iEle])-6.0;
	}
}

void Hetero_field::KozenyCarmen()

{
	int i, j, k;
	double B;
	B = 1;

	for (k = 0; k < Nz; k++)
	for (j = 0; j < Ny; j++)
	for (i = 0; i < Nx; i++){
		int iEle = i + j*Nx + k*Nx*Ny;
		perm[iEle] = B*phi[iEle] * phi[iEle] * phi[iEle] / ((1 - phi[iEle])*(1 - phi[iEle]));
	}
}

void Hetero_field::StanfordKozenyCarmen()
{
	int i, j, k;
	double k0, phis, phish, rouc, rouq, T, Ts, Tsh, Tpore, S, Ss, Ssh;
	double clay;
	k0 = 0.0000000001;
	phis = 0.32; phish = 0.25;
	rouq = 2.568, rouc = 2.77;
	Ts = 1.5; Tsh = 10;
	Ss = 0.15; Ssh = 50;
	for (k = 0; k<Nz; k++)
	{
		for (j = 0; j<Ny; j++)
		{
			for (i = 0; i<Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				clay = phi[iEle] / phish;
				if (clay<phis)
				{
					clay = (phis - phi[iEle]) / (1 - phish);
					Tpore = clay / phis*(Tsh - 1) + 1;
					T = Ts*Tpore;
					S = Ss + Ssh*clay;
				}
				else
				{
					T = Tsh*(1 + (Ts - 1) / (phis - 1)*(clay - 1));
					S = Ss*(1 - clay) / (1 - phis) + Ssh*clay;
				}
				perm[iEle] = phi[iEle] * phi[iEle] * phi[iEle] / (k0*T*T*S*S);
			}

		}

	}
}

void Hetero_field::MohrCoulomb()
{
	int i, j, k;
	for (k = 0; k < Nz; k++)
	for (j = 0; j < Ny; j++)
	for (i = 0; i < Nx; i++){
		int iEle = i + j*Nx + k*Nx*Ny;
		//miu_i[iEle] = tan(2 * (atan(UCS[iEle] / (2 * C[iEle])) - PI / 4));
		miu_i[iEle] = tan(2 * (atan(UCS[iEle] / (2 * C[iEle])) - PI / 4)) - 1;
	}
	for (k = 0; k < Nz; k++)
	for (j = 0; j < Ny; j++)
	for (i = 0; i < Nx; i++){
		int iEle = i + j*Nx + k*Nx*Ny;
		C[iEle] = C[iEle] / 1.5 + 10; 
	}
}

void Hetero_field::BrittlenessIndex()
{
	int i, j, k;
	double E_max = 0, niu_max = 0, E_min = 200000000, niu_min = 1;
	double BI_TSUCS1_min = 1.0e8, BI_TSUCS2_min = 1.0e8, BI_TSUCS3_min = 1.0e8, BI_TSUCS4_min = 1.0e8, BI_YMPR_min=1.0e8;
	double BI_TSUCS1_max = 0, BI_TSUCS2_max = 0, BI_TSUCS3_max = 0, BI_TSUCS4_max = 0, BI_YMPR_max = 0;
	for (k = 0; k<Nz; k++)
	{
		for (j = 0; j<Ny; j++)
		{
			for (i = 0; i<Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				if (Ev[iEle]>E_max) E_max = Ev[iEle];
				if (Eh[iEle]>E_max) E_max = Eh[iEle];
				if (Ev[iEle]<E_min) E_min = Ev[iEle];
				if (Eh[iEle]<E_min) E_min = Eh[iEle];
				if (niuv[iEle]>niu_max) niu_max = niuv[iEle];
				if (niuh[iEle]> niu_max) niu_max = niuh[iEle];
				if (niuhv[iEle]> niu_max) niu_max = niuhv[iEle];
				if (niuv[iEle]< niu_min) niu_min = niuv[iEle];
				if (niuh[iEle]< niu_min) niu_min = niuh[iEle];
				if (niuhv[iEle]< niu_min) niu_min = niuhv[iEle];
			}
		}
	}
	for (k = 0; k < Nz; k++)
	{
		for (j = 0; j < Ny; j++)
		{
			for (i = 0; i < Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				BI_TSUCS1[iEle] = UCS[iEle] / T0[iEle];
				BI_TSUCS2[iEle] = (UCS[iEle] - T0[iEle]) / (UCS[iEle] + T0[iEle]);
				BI_TSUCS3[iEle] = UCS[iEle] * T0[iEle] / 2;
				BI_TSUCS4[iEle] = pow(UCS[iEle] * T0[iEle], 0.5) / 2;
				BI_YMPR[iEle] = ((Ev[iEle] - E_min) / (E_max - E_min) + (niuv[iEle] - niu_max) / (niu_min - niu_max)) / 2;			
			}
		}
	}
	for (k = 0; k < Nz; k++)
	{
		for (j = 0; j < Ny; j++)
		{
			for (i = 0; i < Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				if (BI_TSUCS1[iEle]>=BI_TSUCS1_max) BI_TSUCS1_max = BI_TSUCS1[iEle];
				if (BI_TSUCS1[iEle]<=BI_TSUCS1_min) BI_TSUCS1_min = BI_TSUCS1[iEle];

				if (BI_TSUCS2[iEle]>=BI_TSUCS2_max) BI_TSUCS2_max = BI_TSUCS2[iEle];
				if (BI_TSUCS2[iEle]<=BI_TSUCS2_min) BI_TSUCS2_min = BI_TSUCS2[iEle];

				if (BI_TSUCS3[iEle]>=BI_TSUCS3_max) BI_TSUCS3_max = BI_TSUCS3[iEle];
				if (BI_TSUCS3[iEle]<=BI_TSUCS3_min) BI_TSUCS3_min = BI_TSUCS3[iEle];

				if (BI_TSUCS4[iEle]>=BI_TSUCS4_max) BI_TSUCS4_max = BI_TSUCS4[iEle];
				if (BI_TSUCS4[iEle]<=BI_TSUCS4_min) BI_TSUCS4_min = BI_TSUCS4[iEle];

				if (BI_YMPR[iEle]>=BI_YMPR_max) BI_YMPR_max = BI_YMPR[iEle];
				if (BI_YMPR[iEle]<=BI_YMPR_min) BI_YMPR_min = BI_YMPR[iEle];
			}
		}
	}
	
	double temD1, temD2, temD3, temD4, temD5;
	temD1 = BI_TSUCS1_max - BI_TSUCS1_min;
	temD2 = BI_TSUCS2_max - BI_TSUCS2_min;
	temD3 = BI_TSUCS3_max - BI_TSUCS3_min;
	temD4 = BI_TSUCS4_max - BI_TSUCS4_min;
	temD5 = BI_YMPR_max - BI_YMPR_min;

	for (k = 0; k < Nz; k++)
	{
		for (j = 0; j < Ny; j++)
		{
			for (i = 0; i < Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				BI_TSUCS1[iEle] = (BI_TSUCS1[iEle] - BI_TSUCS1_min) / temD1;
				BI_TSUCS2[iEle] = (BI_TSUCS2[iEle] - BI_TSUCS2_min) / temD2;
				BI_TSUCS3[iEle] = (BI_TSUCS3[iEle] - BI_TSUCS3_min) / temD3;
				BI_TSUCS4[iEle] = (BI_TSUCS4[iEle] - BI_TSUCS4_min) / temD4;
				BI_YMPR[iEle] = (BI_YMPR[iEle] - BI_YMPR_min) / temD5;
			}
		}
	}
	
}

void Hetero_field::Output()
{
	int i, j, k;
	FILE *fp0, *fp1, *fp2, *fp3, *fp4, *fp5, *fp6, *fp7, *fp8, *fp9, *fp10, *fp11, *fp12;
	fopen_s(&fp0, "Vp.dat", "w");
	fopen_s(&fp1, "Vs.dat", "w");
	fopen_s(&fp2, "porosity.dat", "w");
	fopen_s(&fp3, "permeability.dat", "w");
	fopen_s(&fp4, "tensile strength.dat", "w");
	fopen_s(&fp5, "cohesive strength.dat", "w");
	fopen_s(&fp6, "friction coefficient.dat", "w");
	fopen_s(&fp7, "vertical young's modulus.dat", "w");
	fopen_s(&fp8, "horizontal young's modulus.dat", "w");
	fopen_s(&fp9, "v-v poisson's ratio.dat", "w");
	fopen_s(&fp10, "h-h poisson's ratio.dat", "w");
	fopen_s(&fp11, "h-v poisson's ratio.dat", "w");
	fopen_s(&fp12, "uniaixal compressive strength.dat", "w");
	for (k = 0; k < Nz; k++)
	{
		for (j = 0; j < Ny; j++)
		{
			for (i = 0; i < Nx; i++)
			{
				int iEle = i + j*Nx + k*Nx*Ny;
				fprintf(fp0, "%lf\n", Vp[iEle]);
				fprintf(fp1, "%lf\n", Vs[iEle]);
				fprintf(fp2, "%lf\n", phi[iEle]);
				fprintf(fp3, "%lf\n", perm[iEle]);
				fprintf(fp4, "%lf\n", T0[iEle]);
				fprintf(fp5, "%lf\n", C[iEle]);
				fprintf(fp6, "%lf\n", miu_i[iEle]);
				fprintf(fp7, "%lf\n", Ev[iEle]);
				fprintf(fp8, "%lf\n", Eh[iEle]);
				fprintf(fp9, "%lf\n", niuv[iEle]);
				fprintf(fp10, "%lf\n", niuh[iEle]);
				fprintf(fp11, "%lf\n", niuhv[iEle]);
				fprintf(fp12, "%lf\n", UCS[iEle]);
			}
		}
	}
}
}//namespace FracTHM