#include "Database.h"
#include <string>

namespace FracTHM{
	using namespace std;

	static bool FileEnd(std::ifstream &ifp) { return ifp.peek() == EOF; };
	static bool UGStrCmp(const char* a, const char* b) { return (_stricmp(a, b) == 0); };
	static bool UGIsDigit(const char a) { return (isdigit(a) || a == '-' || a == '.'); };

	static coord8P getcpgcoord8p(Vertex *XYZden, const int& xp, const int& yp, const int& zp, const int& Nx, const int& Ny, const int& Nz){

		int temp;
		coord8P base8p;
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * xp;
		base8p.p1 = XYZden[temp];
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * xp + 1;
		base8p.p2 = XYZden[temp];			
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * Nx + 2 * xp + 1;
		base8p.p3 = XYZden[temp];
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * Nx + 2 * xp;
		base8p.p4 = XYZden[temp];

		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * xp;
		base8p.p5 = XYZden[temp];
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * xp + 1;
		base8p.p6 = XYZden[temp];				
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * Nx + 2 * xp + 1;
		base8p.p7 = XYZden[temp];
		temp = 2 * zp * 2 * Nx * 2 * Ny + 2 * Nx * 2 * Ny + 2 * yp * 2 * Nx + 2 * Nx + 2 * xp;
		base8p.p8 = XYZden[temp];

		return base8p;
	}
	static double TetraVol(const Vertex& P, const Vertex& a, const Vertex& b, const Vertex& c) {
		double ax = a.Xcood - P.Xcood, ay = a.Ycood - P.Ycood, az = a.Zcood - P.Zcood;
		double bx = b.Xcood - P.Xcood, by = b.Ycood - P.Ycood, bz = b.Zcood - P.Zcood;
		double cx = c.Xcood - P.Xcood, cy = c.Ycood - P.Ycood, cz = c.Zcood - P.Zcood;
		double Volume = fabs((ay*bz - az*by)*cx + (az*bx - ax*bz)*cy + (ax*by - ay*bx)*cz) / (double)6.0;
		return Volume;
	};
	static double TetraVol(const Vertex &p, const Vertex &a, const Vertex &b, const Vertex &c, Vertex &e) {
		e.Xcood = 0.25*(a.Xcood + b.Xcood + c.Xcood + p.Xcood);
		e.Ycood = 0.25*(a.Ycood + b.Ycood + c.Ycood + p.Ycood);
		e.Zcood = 0.25*(a.Zcood + b.Zcood + c.Zcood + p.Zcood);

		double ax = a.Xcood - p.Xcood, ay = a.Ycood - p.Ycood, az = a.Zcood - p.Zcood;
		double bx = b.Xcood - p.Xcood, by = b.Ycood - p.Ycood, bz = b.Zcood - p.Zcood;
		double cx = c.Xcood - p.Xcood, cy = c.Ycood - p.Ycood, cz = c.Zcood - p.Zcood;
		double Volume = fabs((ay*bz - az*by)*cx + (az*bx - ax*bz)*cy + (ax*by - ay*bx)*cz) / (double)6.0;
		return Volume;
	};
	static Vertex getcellcenter(const Vertex grid[8]){
		Vertex c1, c2, c3, c4, c5;  //4572,8457,6257,3247,1245
		double v1 = TetraVol(grid[3], grid[4], grid[6], grid[1], c1);
		double v2 = TetraVol(grid[7], grid[3], grid[4], grid[6], c2);
		double v3 = TetraVol(grid[5], grid[1], grid[4], grid[6], c3);
		double v4 = TetraVol(grid[2], grid[1], grid[3], grid[6], c4);
		double v5 = TetraVol(grid[0], grid[1], grid[3], grid[4], c5);

		Vertex e;
		double v = v1 + v2 + v3 + v4 + v5;

		if (v > InfSmall) {
			e.Xcood = (c1.Xcood* v1 + c2.Xcood* v2 + c3.Xcood* v3 + c4.Xcood* v4 + c5.Xcood* v5) / v;
			e.Ycood = (c1.Ycood* v1 + c2.Ycood* v2 + c3.Ycood* v3 + c4.Ycood* v4 + c5.Ycood* v5) / v;
			e.Zcood = (c1.Zcood* v1 + c2.Zcood* v2 + c3.Zcood* v3 + c4.Zcood* v4 + c5.Zcood* v5) / v;
		}
		else {
			e.Xcood = 0.25*(grid[3].Xcood + grid[4].Xcood + grid[6].Xcood + grid[1].Xcood);
			e.Ycood = 0.25*(grid[3].Ycood + grid[4].Ycood + grid[6].Ycood + grid[1].Ycood);
			e.Zcood = 0.25*(grid[3].Zcood + grid[4].Zcood + grid[6].Zcood + grid[1].Zcood);
		}
		return e;
	}
	static double HexVol(const Vertex cell[8]){
		double v1 = TetraVol(cell[0], cell[1], cell[3], cell[4]);
		double v2 = TetraVol(cell[3], cell[4], cell[5], cell[1]);
		double v3 = TetraVol(cell[3], cell[4], cell[5], cell[1]);
		double v4 = TetraVol(cell[1], cell[2], cell[3], cell[6]);
		double v5 = TetraVol(cell[3], cell[5], cell[6], cell[7]);
		double v6 = TetraVol(cell[3], cell[5], cell[6], cell[1]);
		double volume = v1 + v2 + v3 + v4 + v5 + v6;
		return volume;
	}

	static void skip_sharp(std::ifstream &ifp){

		while (!FileEnd(ifp)){
			char ch;
			ifp >> ch;

			if (ch == '#'){
				string s;
				getline(ifp, s, '\n');
			}
			else if (ch == '-' && ifp.peek() == '-'){
				string s;
				getline(ifp, s, '\n');
			}
			else {
				ifp.seekg(-1, ios_base::cur);
				return;
			}
		}
	};
	static void skip_slash(std::ifstream &ifp){
		while (!FileEnd(ifp)) {
			char ch;
			ifp >> ch;

			if (ch == '#') {
				string s;
				getline(ifp, s, '\n');
			}
			else if (ch == '-' && ifp.peek() == '-') {
				string s;
				getline(ifp, s, '\n');
			}
			else if (ch == '/') {
				string s;
				getline(ifp, s, '\n');
			}
			else {
				ifp.seekg(-1, ios_base::cur);
				return;
			}
		}
	};
	static bool meet_NaN(std::ifstream &ifp) {
		skip_sharp(ifp);
		char ch = ifp.peek();
		if (UGIsDigit(ch)) {
			return false;
		}
		else {
			if (ch == '/') {
				string s;
				getline(ifp, s, '\n');
			}
			return true;
		}
	};

	int CpgtoFrac::readCode(std::ifstream &ifp) {
		ifp >> cbuffer;

		int clen = 0;
		while (cbuffer[clen] != '\0' && cbuffer[clen] != '#') {
			++clen;
		}

		if (cbuffer[clen] != '\0') {
			string comment;
			getline(ifp, comment);
			cbuffer[clen] = '\0';
		}

		if (clen > 1 && cbuffer[clen - 1] == '/') {
			cbuffer[--clen] = '\0';
			ifp.putback('/');
		}

		return clen;
	};
	int CpgtoFrac::readInt(std::ifstream &ifp) {
		skip_slash(ifp);
		readCode(ifp);
		return atoi(cbuffer);
	};
	double CpgtoFrac::readDouble(std::ifstream &ifp){
		skip_slash(ifp);
		readCode(ifp);
		return atof(cbuffer);
	}
	int CpgtoFrac::readValue(double* db, std::ifstream &ifp){
		if (meet_NaN(ifp)) {
			return 0;
		}

		int clen = readCode(ifp);

		int starLoc = -1;
		for (int i = 0; i < clen; i++) {    // Check whether there is a "*"
			if (cbuffer[i] == '*') {
				starLoc = i;
				break;
			}
		}
		if (starLoc != -1) {                         // there is a "*"
			cbuffer[starLoc] = '\0';
			*db = atof(cbuffer + starLoc + 1);    // convert string to double, *db = b
			return atoi(cbuffer);                 // convert string to int, return a
		}
		else {                                   // single number
			*db = atof(cbuffer);
			return 1;
		}
	};
	int CpgtoFrac::readArrayD(double* a, const int length, std::ifstream &ifp){
		int count = 0, pos = 0;

		// Read the following values.
		while (pos < length) {
			count = readValue(a + pos, ifp);
			if (count == 0) break;		                 // reading finished
			for (int i = pos + 1; i < pos + count && i < length; ++i) a[i] = a[pos];
			pos += count;
			if (pos >= length){
				break;
			}
		}

		return pos;
	};

	bool CpgtoFrac::readProps(const char* filename, DataBase& CVFEM_DataBase){	
		std::ifstream ifp;
		ifp.open(filename, ios_base::in | ios_base::binary);
		int Ncell = CVFEM_DataBase.GridNum;
		if (!ifp.fail()){
			while (!FileEnd(ifp)){
				skip_slash(ifp);
				readCode(ifp);

				if (UGStrCmp(cbuffer, "PORO")){
					int nlen = Ncell;
					int pos = readArrayD(CVFEM_DataBase.iPORO_, nlen, ifp);
					for (int i = 0; i < nlen; ++i)
						CVFEM_DataBase.iPORO_[i] = CVFEM_DataBase.iPORO0_[i] = 0.01 * CVFEM_DataBase.iPORO_[i];
				}
				else if (UGStrCmp(cbuffer, "PERMEABILITY")){
					int nlen = Ncell;
					int pos = readArrayD(CVFEM_DataBase.permXX, nlen, ifp);

					for (int i = 0; i < nlen; ++i) {
						CVFEM_DataBase.PermTensor[i].VarXX = CVFEM_DataBase.PermTensor[i].VarYY = CVFEM_DataBase.PermTensor[i].VarZZ= CVFEM_DataBase.permXX[i];
						CVFEM_DataBase.PermTensor[i].VarXY = CVFEM_DataBase.PermTensor[i].VarXZ = CVFEM_DataBase.PermTensor[i].VarYZ = 0.0;

						CVFEM_DataBase.permYY[i] = CVFEM_DataBase.permZZ[i] = CVFEM_DataBase.permXX[i];
						CVFEM_DataBase.permXY[i] = CVFEM_DataBase.permXZ[i] = CVFEM_DataBase.permYZ[i] = 0.0;
					}
				}
				else if (UGStrCmp(cbuffer, "POISSONSRATIO")){
					int nlen = Ncell;
					int pos = readArrayD(CVFEM_DataBase.PoissonRatio, nlen, ifp);
				}
				else if (UGStrCmp(cbuffer, "YOUNG'SMODULUS")){
					int nlen = Ncell;
					int pos = readArrayD(CVFEM_DataBase.YoungsModulus, nlen, ifp);
					for (int i = 0; i < nlen; ++i){
						CVFEM_DataBase.YoungsModulus[i] = 1.0e9 * CVFEM_DataBase.YoungsModulus[i];

						//heat conductivity
						CVFEM_DataBase.heatcdXX[i] = CVFEM_DataBase.heatcdYY[i] = CVFEM_DataBase.heatcdZZ[i] = 3.3;
						CVFEM_DataBase.heatcdXY[i] = CVFEM_DataBase.heatcdXZ[i] = CVFEM_DataBase.heatcdYZ[i] = 0.0;

						//biot coef.
						CVFEM_DataBase.CoefBiot_[i] = 1.0;

						CVFEM_DataBase.RockDensity_[i] = 2000;
						CVFEM_DataBase.alpha_R_[i] = 1.2e-6;
						CVFEM_DataBase.C_R_[i] = 690;
						CVFEM_DataBase.GrainM_[i] = 1.0e12;
						
					}
				}
				else if (UGStrCmp(cbuffer, "SWAT") ||
					UGStrCmp(cbuffer, "MAXSTRESS") ||
					UGStrCmp(cbuffer, "MINSTRESS") ||
					UGStrCmp(cbuffer, "V_STRESS") ||
					UGStrCmp(cbuffer, "KANGLA") ||
					UGStrCmp(cbuffer, "POLIE")){
					double* dTemp = new double[Ncell];
					int nlen = Ncell;
					int pos = readArrayD(dTemp, nlen, ifp);
					delete[] dTemp;
				}
			}
		}
		ifp.close();
		return true;
	}



	bool CpgtoFrac::readCPGrid(const char* filename, DataBase& CVFEM_DataBase){

		int Nx, Ny, Nz, Ncell;
		Nx = Ny = Nz = Ncell = 0;
		bool coordValid = false;
		bool zcornValid = false;
		double* coord = nullptr;
		double* zcorn = nullptr;

		std::ifstream ifp;
		ifp.open(filename, ios_base::in | ios_base::binary);

		if (!ifp.fail()){
			while (!FileEnd(ifp)){
				skip_slash(ifp);
				readCode(ifp);

				if (UGStrCmp(cbuffer, "DIMENS")){
					Nx = readInt(ifp);
					Ny = readInt(ifp);
					Nz = readInt(ifp);
					Ncell = Nx*Ny*Nz;
				}
				else if (UGStrCmp(cbuffer, "COORD")){
					if (Ncell > 0){
						int nlen = 6 * (Nx + 1)*(Ny + 1);
						delete[] coord; coord = new double[nlen];
						int pos = readArrayD(coord, nlen, ifp);
						coordValid = (pos == nlen);
					}
					else {
						break;
					}
				}
				else if (UGStrCmp(cbuffer, "ZCORN")){
					if (Ncell > 0){
						int nlen = 8 * Ncell;
						delete[] zcorn;  zcorn = new double[8 * Ncell];
						int pos = readArrayD(zcorn, nlen, ifp);
						zcornValid = (pos == nlen);
					}
					else {
						break;
					}
				}
			}
		}
		ifp.close();

		if (Ncell > 0 && coordValid && zcornValid){

			Vertex* cellPts = new Vertex[8 * Ncell];
			double t;
			int temp = 0, ipl;
			for (int k = 0; k < 2 * Nz; k++){
				for (int j = 0; j < 2 * Ny; j++){
					for (int i = 0; i < 2 * Nx; i++){
						ipl = ((i + 1) / 2 + (j + 1) / 2 * (Nx + 1)) * 6;

						double z = zcorn[temp];

						t = (z - coord[ipl + 2]) / __max(coord[ipl + 5] - coord[ipl + 2], 1e-10);

						cellPts[temp].Xcood = coord[ipl] + (coord[ipl + 3] - coord[ipl])*t - coord[0];
						cellPts[temp].Ycood = coord[ipl + 1] + (coord[ipl + 4] - coord[ipl + 1])*t - coord[1];
						cellPts[temp].Zcood = -z;
						temp += 1;
					}
				}
			}

			CVFEM_DataBase.Nx = Nx;
			CVFEM_DataBase.Ny = Ny;
			CVFEM_DataBase.Nz = Nz;
			int segLen = (Nx + 1)*(Ny + 1)*(Nz + 1);

			//
//			CVFEM_DataBase.NodeAmt = CVFEM_DataBase.NodeNum = segLen;
//			CVFEM_DataBase.Dx = new double[CVFEM_DataBase.Nx];
//			CVFEM_DataBase.Dy = new double[CVFEM_DataBase.Ny];
//			CVFEM_DataBase.Dz = new double[CVFEM_DataBase.Nz];
//
//			CVFEM_DataBase.ElementAmt = CVFEM_DataBase.GridNum = Ncell;
//			CVFEM_DataBase.T0_I = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.S0_I = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.miu_I = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.fricAng_I = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.UCS = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.a0 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.vm = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.Kni = new double[CVFEM_DataBase.GridNum];
//
//			CVFEM_DataBase.iPORO_ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.iPORO0_ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.iVOLUME_ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.iDEPTH_ = new double[CVFEM_DataBase.GridNum];
//
//			CVFEM_DataBase.YoungsModulus = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.aYM1 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.aYM2 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.PoissonRatio = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.aPR1 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.aPR2 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.ShearM = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.aSM1 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.aSM2 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.CoefBiot_ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.GrainM_ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.C_R_ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.Kdr_ = new Vertex[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.alpha_R_ = new double[CVFEM_DataBase.GridNum];
//			//pAllData->alpha_f = new double[pAllData->GridNum];
//			CVFEM_DataBase.RockDensity_ = new double[CVFEM_DataBase.GridNum];
//
//			CVFEM_DataBase.isFrctureTip = new bool[CVFEM_DataBase.GridNum];
//			memset(CVFEM_DataBase.isFrctureTip, 0x0, CVFEM_DataBase.GridNum*sizeof(bool));
//			CVFEM_DataBase.isiELEFailure = new bool[CVFEM_DataBase.GridNum];
//			memset(CVFEM_DataBase.isiELEFailure, 0x0, CVFEM_DataBase.GridNum*sizeof(bool));
//			CVFEM_DataBase.isNaturalFractured = new bool[CVFEM_DataBase.GridNum];
//			memset(CVFEM_DataBase.isNaturalFractured, 0x0, CVFEM_DataBase.GridNum*sizeof(bool));
//			CVFEM_DataBase.isHydraulicFractured = new bool[CVFEM_DataBase.GridNum];
//			memset(CVFEM_DataBase.isHydraulicFractured, 0x0, CVFEM_DataBase.GridNum*sizeof(bool));
//			CVFEM_DataBase.isPerforationLoc = new int[CVFEM_DataBase.GridNum];
//			memset(CVFEM_DataBase.isPerforationLoc, 0x0, CVFEM_DataBase.GridNum*sizeof(int));
//
//			CVFEM_DataBase.BI_TSUCS1 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.BI_TSUCS2 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.BI_TSUCS3 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.BI_TSUCS4 = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.BI_YMPR = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.BI_YMM = new double[CVFEM_DataBase.GridNum];
//
//			CVFEM_DataBase.permXX = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.permYY = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.permZZ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.permXY = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.permYZ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.permXZ = new double[CVFEM_DataBase.GridNum];
//
//			CVFEM_DataBase.PermTensor = new tensorVar[CVFEM_DataBase.GridNum];
//
//			CVFEM_DataBase.heatcdXX = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.heatcdYY = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.heatcdZZ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.heatcdXY = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.heatcdYZ = new double[CVFEM_DataBase.GridNum];
//			CVFEM_DataBase.heatcdXZ = new double[CVFEM_DataBase.GridNum];
			//

			CVFEM_DataBase.Xcoord = new double[segLen];
			CVFEM_DataBase.Ycoord = new double[segLen];
			CVFEM_DataBase.Zcoord = new double[segLen];
			CVFEM_DataBase.NodeNumSequence = new int[CVFEM_DataBase.GridNum * 8];

			for (int k = 0; k < Nz + 1; k++){
				for (int j = 0; j < Ny + 1; j++){
					for (int i = 0; i < Nx + 1; i++){
						Vertex p;
						int np = 0;
						if (k >= 1){
							if (i >= 1 && j >= 1){
								coord8P pts = getcpgcoord8p(cellPts, i - 1, j - 1, k - 1, Nx, Ny, Nz);
								p = p + pts.p7; np += 1;
							}
							if (i >= 1 && j < Ny){
								coord8P pts = getcpgcoord8p(cellPts, i - 1, j, k - 1, Nx, Ny, Nz);
								p = p + pts.p6; np += 1;
							}
							if (i < Nx && j >= 1){
								coord8P pts = getcpgcoord8p(cellPts, i, j - 1, k - 1, Nx, Ny, Nz);
								p = p + pts.p8; np += 1;
							}
							if (i < Nx && j < Ny){
								coord8P pts = getcpgcoord8p(cellPts, i, j, k - 1, Nx, Ny, Nz);
								p = p + pts.p5; np += 1;
							}
						}
						if (k < Nz){
							if (i >= 1 && j >= 1){
								coord8P pts = getcpgcoord8p(cellPts, i - 1, j - 1, k, Nx, Ny, Nz);
								p = p + pts.p3; np += 1;
							}
							if (i >= 1 && j < Ny){
								coord8P pts = getcpgcoord8p(cellPts, i - 1, j, k, Nx, Ny, Nz);
								p = p + pts.p2; np += 1;
							}
							if (i < Nx && j >= 1){
								coord8P pts = getcpgcoord8p(cellPts, i, j - 1, k, Nx, Ny, Nz);
								p = p + pts.p4; np += 1;
							}
							if (i < Nx && j < Ny){
								coord8P pts = getcpgcoord8p(cellPts, i, j, k, Nx, Ny, Nz);
								p = p + pts.p1; np += 1;
							}
						}

						if (np > 0) p = p / np;

						int segId = i + j*(Nx + 1) + (Nz - k)*(Nx + 1)*(Ny + 1);
						CVFEM_DataBase.Xcoord[segId] = p.Xcood;
						CVFEM_DataBase.Ycoord[segId] = p.Ycood;
						CVFEM_DataBase.Zcoord[segId] = p.Zcood;
					}
				}
			}
			//for (int i = 0; i < segLen; i++) CVFEM_DataBase.NodeNumSequence[i] = i + 1;

			CVFEM_DataBase.cells = new CVFEMcell[Ncell];
			CVFEM_DataBase.iVOLUME_ = new double[Ncell];
			CVFEM_DataBase.iDEPTH_ = new double[Ncell];

			int cellId = 0;
			for (int k = 0; k < Nz; k++){
				for (int j = 0; j < Ny; j++){
					for (int i = 0; i < Nx; i++){
						int _EleNum = EleNum(i, j, k, Nx, Ny);
						int* Pt2NdNum = &CVFEM_DataBase.NodeNumSequence[_EleNum * 8];
						Pt2NdNum[0] = CalNodeNum(i, j, k, Nx, Ny) + 1;
						Pt2NdNum[1] = CalNodeNum(i + 1, j, k, Nx, Ny) + 1;
						Pt2NdNum[2] = CalNodeNum(i + 1, j + 1, k, Nx, Ny) + 1;
						Pt2NdNum[3] = CalNodeNum(i, j + 1, k, Nx, Ny) + 1;
						Pt2NdNum[4] = CalNodeNum(i, j, k + 1, Nx, Ny) + 1;
						Pt2NdNum[5] = CalNodeNum(i + 1, j, k + 1, Nx, Ny) + 1;
						Pt2NdNum[6] = CalNodeNum(i + 1, j + 1, k + 1, Nx, Ny) + 1;
						Pt2NdNum[7] = CalNodeNum(i, j + 1, k + 1, Nx, Ny) + 1;
						//
						
						CVFEMcell& cellPtr = CVFEM_DataBase.cells[cellId];
						//顶点坐标
						int xid[8];
						int yid[8];
						int zid[8];
						xid[0] = i;		yid[0] = j;		zid[0] = k;
						xid[1] = i + 1;	yid[1] = j;		zid[1] = k;
						xid[2] = i + 1;	yid[2] = j + 1;	zid[2] = k;
						xid[3] = i;		yid[3] = j + 1;	zid[3] = k;
						xid[4] = i;		yid[4] = j;		zid[4] = k + 1;
						xid[5] = i + 1; yid[5] = j;		zid[5] = k + 1;
						xid[6] = i + 1; yid[6] = j + 1; zid[6] = k + 1;
						xid[7] = i;		yid[7] = j + 1; zid[7] = k + 1;

						for (int ll = 0; ll < 8; ll++){
							int segId = xid[ll] + yid[ll] * (Nx + 1) + zid[ll] * (Nx + 1)*(Ny + 1);
							cellPtr.corner[ll].Xcood = CVFEM_DataBase.Xcoord[segId];
							cellPtr.corner[ll].Ycood = CVFEM_DataBase.Ycoord[segId];
							cellPtr.corner[ll].Zcood = CVFEM_DataBase.Zcoord[segId];
						}
						//面中心坐标
						cellPtr.facec[0] = (cellPtr.corner[4] + cellPtr.corner[5] + cellPtr.corner[6] + cellPtr.corner[7])*0.25;
						cellPtr.facec[1] = (cellPtr.corner[0] + cellPtr.corner[1] + cellPtr.corner[2] + cellPtr.corner[3])*0.25;
						cellPtr.facec[2] = (cellPtr.corner[0] + cellPtr.corner[1] + cellPtr.corner[4] + cellPtr.corner[5])*0.25;
						cellPtr.facec[3] = (cellPtr.corner[2] + cellPtr.corner[3] + cellPtr.corner[7] + cellPtr.corner[6])*0.25;
						cellPtr.facec[4] = (cellPtr.corner[0] + cellPtr.corner[3] + cellPtr.corner[7] + cellPtr.corner[4])*0.25;
						cellPtr.facec[5] = (cellPtr.corner[1] + cellPtr.corner[2] + cellPtr.corner[6] + cellPtr.corner[5])*0.25;
						//线中心坐标
						cellPtr.linec[0] = (cellPtr.corner[4] + cellPtr.corner[5])*0.5;
						cellPtr.linec[1] = (cellPtr.corner[6] + cellPtr.corner[7])*0.5;
						cellPtr.linec[2] = (cellPtr.corner[4] + cellPtr.corner[7])*0.5;
						cellPtr.linec[3] = (cellPtr.corner[5] + cellPtr.corner[6])*0.5;
						cellPtr.linec[4] = (cellPtr.corner[0] + cellPtr.corner[1])*0.5;
						cellPtr.linec[5] = (cellPtr.corner[2] + cellPtr.corner[3])*0.5;
						cellPtr.linec[6] = (cellPtr.corner[0] + cellPtr.corner[3])*0.5;
						cellPtr.linec[7] = (cellPtr.corner[1] + cellPtr.corner[2])*0.5;
						cellPtr.linec[8] = (cellPtr.corner[0] + cellPtr.corner[4])*0.5;
						cellPtr.linec[9] = (cellPtr.corner[1] + cellPtr.corner[5])*0.5;
						cellPtr.linec[10] = (cellPtr.corner[3] + cellPtr.corner[7])*0.5;
						cellPtr.linec[11] = (cellPtr.corner[2] + cellPtr.corner[6])*0.5;
						//块中心
						cellPtr.blockc = getcellcenter(cellPtr.corner);
						//体积
						CVFEM_DataBase.iVOLUME_[cellId] = HexVol(cellPtr.corner);
						CVFEM_DataBase.iDEPTH_[cellId] = cellPtr.blockc.Zcood;
						cellId++;
					}
				}
			}

			delete[] cellPts;
			delete[] coord;
			delete[] zcorn;
			return true;
		}
		else {
			return false;
		}
	}
}