/**
* @file    simHM.cpp
* @author  Sanbai Li
* @date    2014/9~
* @brief   Deal with different kinds of Boundary Conditions

*
* @version
* Copyright (c) 2019, NEU & SINOPEC.
*/

#include "Constants.h"
#include "SimHM.h"
#include "FEMfunction.h"
#include "MatrixManuplt.h"

namespace KarstVolSys{

	/****************Deal with different boundary conditions***********************/

	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//     1. Displacement B.C.: Use penalty function method to deal with the displacement boundary condition      //
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	void SimHM::HandleDispBC(){
		std::list<Constrain>::iterator DispVecItr;
		for (DispVecItr = AllData->Cnstrlist.begin();
			DispVecItr != AllData->Cnstrlist.end(); ++DispVecItr){
			Constrain Cnstr = *DispVecItr;
			int nodeI = Cnstr.cNd_num;
			//disp[3 * nodeI-3] = 0.0;
			double* stLoc = JacMtr.find_GGD(nodeI, 0);
			stLoc[4 * Cnstr.cDirc - 4] *= PenaltyMulti;  //Krr=a*Krr
			Residual[3 * nodeI + Cnstr.cDirc - 1] = stLoc[4 * Cnstr.cDirc - 4] * Cnstr.Dfm;//RHS=a*Krr*u
		}
	}

	///////////////////////////////////////////////////////////////////////////////
	//     2. Node load B.C.: Add to according location of redidual directly     //
	///////////////////////////////////////////////////////////////////////////////
	void SimHM::HandleLoadBC(){
		std::list<Load>::iterator LoadVecItr;
		for (LoadVecItr = AllData->Ldlist.begin();
			LoadVecItr != AllData->Ldlist.end(); ++LoadVecItr){
			Load Cnstr = *LoadVecItr;
			int nodeI = Cnstr.lNd_num;
			Residual[3 * nodeI + Cnstr.lDirc - 1] += Cnstr.Force;
		}
	}

	///////////////////////////////////////////////////////////////////
	//   3. Face load B.C., be written in subroutine ActiveSolid()   //
	///////////////////////////////////////////////////////////////////

	///////////////////////////////////////////////////////////////////
	//               4. Fix concentration & pressure B.C.            //
	///////////////////////////////////////////////////////////////////
	void SimHM::FixPrsBC(const double dt){
		std::list<fixPre>::iterator itrPreslist;
		for (itrPreslist = AllData->Preslist.begin();
			itrPreslist != AllData->Preslist.end(); ++itrPreslist){
			double perm, area, dist;
			int cellID = (*itrPreslist).pEle_num;
			double pressure = (*itrPreslist).pres;

			(*itrPreslist).perm = PermM2F[cellID];

			perm = (*itrPreslist).perm;
			area = (*itrPreslist).height*filmthick0[cellID];
			dist = (*itrPreslist).dist;

			double* ptRsv = &Residual[geom_bh*geom_size];
			double* stLoc = JacMtr.find_RRD(cellID, 0);
			int RowCtr = JacMtr.RRP.bh;
			double tempD1 = perm * area / dist * dt; 

			//Chemical: C_Ca+2
			if (prs[cellID] > pressure){//Production, outflow
				//fluid flow

				stLoc[0] += tempD1 * uBw(cellID) + tempD1 * duBwdp(cellID)* (prs[cellID] - pressure);// AllData->Viscosity;//JAC
				//stLoc[1] += tempD1 * duBwdT(cellID)* (prs[cellID] - pressure);// AllData->Viscosity;//JAC
				double TempQf = tempD1 * uBw(cellID) * (prs[cellID] - pressure);
				ptRsv[RowCtr * cellID] -= TempQf;// AllData->Viscosity;//Residual
				//cout << "Production:" << TempQf<<endl;
				
				AllData->tempQ_out = 1.0E-12*TempQf / dt; //m3/s
				AllData->C_Ca_out = Temp[cellID];    //mol/m3

				stLoc[2] += 1.0E-12 * tempD1 * uBw(cellID)* Temp[cellID];
				stLoc[3] += 1.0E-12*TempQf;
				ptRsv[RowCtr * cellID + 1] -= 1.0E-12 * TempQf*Temp[cellID];
			}
			else if (prs[cellID] <= pressure && !(*itrPreslist).isPro){ //Injection, inflow
				//fluid flow
				stLoc[0] += tempD1 * uBw(cellID);// AllData->Viscosity;//JAC
				//stLoc[1] += tempD1 * duBwdT(cellID)* (prs[cellID] - pressure);// AllData->Viscosity;//JAC
				double TempQf = tempD1 * uBw(cellID) * (pressure - prs[cellID]);
				ptRsv[RowCtr * cellID] += TempQf;// AllData->Viscosity;//Residual
				//cout << "Injection:" << TempQf << endl;
				AllData->tempQ_in = 1.0E-12*TempQf / dt; //m3/s
				AllData->C_Ca_in = Temp[cellID];    //mol/m3

				stLoc[2] += AllData->Cin * 1.0E-12* tempD1 * uBw(cellID);
				ptRsv[RowCtr * cellID + 1] += AllData->Cin * 1.0E-12*TempQf;

				//stLoc[2] += 1.0E-12* tempD1 * uBw(cellID)* (Temp[cellID] - AllData->Cin);
				//stLoc[3] += 1.0E-12*TempQf;
				//ptRsv[RowCtr * cellID + 1] -= 1.0E-12*TempQf * (Temp[cellID] - AllData->Cin);
			}
			else;
		}
	}

	void SimHM::FixConcBC(const double dt){
		std::list<fixConc>::iterator itrConclist;
		for (itrConclist = AllData->Conclist.begin();
			itrConclist != AllData->Conclist.end(); ++itrConclist){
			double diffu, area, dist;
			int dirc = (*itrConclist).pDirc;
			int cellID = (*itrConclist).pEle_num;
			int I = cellID %AllData->Nx,
				J = (cellID / AllData->Nx) % AllData->Ny,
				K = (cellID / AllData->Nx) / AllData->Ny;
			double Concentration = (*itrConclist).Conc;
			if (dirc == 0){
				diffu = InfLarge; // AllData->heatcdZZ[cellID];
				area = AllData->Dx[I] * AllData->Dy[J];
				dist = 0.5*AllData->Dz[K];
			}
			if (dirc == 1){
				diffu = InfLarge; // AllData->heatcdXX[cellID];
				area = AllData->Dy[J] * AllData->Dz[K];
				dist = 0.5*AllData->Dx[I];
			}
			else if (dirc == 2){
				diffu = InfLarge; // AllData->heatcdYY[cellID];
				area = AllData->Dx[I] * AllData->Dz[K];
				dist = 0.5*AllData->Dy[J];
			}
			else{
				diffu = InfLarge;
				area = 1.0;
				dist = 1.0;
			}
			double* ptRsv = &Residual[geom_bh*geom_size];
			double* stLoc = JacMtr.find_RRD(cellID, 0);
			int RowCtr = JacMtr.RRP.bh;
			stLoc[3] += dt * diffu * area / dist;//JAC
			ptRsv[RowCtr * cellID + 1] -= dt * diffu * area / dist * (Temp[cellID] - Concentration);//Residual

		}
	}

	/////////////////////////////////////////////////////////////
	//                 5. Fix temperature B.C.                 //
	/////////////////////////////////////////////////////////////
	void SimHM::FixTempBC(const double dt)
	{
		std::list<fixTemp>::iterator itrTemplist;
		for (itrTemplist = AllData->Templist.begin();
			itrTemplist != AllData->Templist.end(); ++itrTemplist){
			double condc, area, dist;
			int dirc = (*itrTemplist).pDirc;
			int cellID = (*itrTemplist).pEle_num;
			int I = cellID %AllData->Nx,
				J = (cellID / AllData->Nx) % AllData->Ny,
				K = (cellID / AllData->Nx) / AllData->Ny;
			double temperature = (*itrTemplist).pTemp;
			if (dirc == 1){
				condc = AllData->heatcdXX[cellID];
				area = AllData->Dy[J] * AllData->Dz[K];
				dist = 0.5*AllData->Dx[I];
			}
			else if (dirc == 2){
				condc = AllData->heatcdYY[cellID];
				area = AllData->Dx[I] * AllData->Dz[K];
				dist = 0.5*AllData->Dy[J];
			}
			else {
				condc = AllData->heatcdZZ[cellID];
				area = AllData->Dx[I] * AllData->Dy[J];
				dist = 0.5*AllData->Dz[K];
			}
			double* ptRsv = &Residual[geom_bh*geom_size];
			double* stLoc = JacMtr.find_RRD(cellID, 0);
			int RowCtr = JacMtr.RRP.bh;
			stLoc[3] += dt * condc * area / dist / visc[cellID]*1000;//JAC
			ptRsv[RowCtr * cellID + 1] -= dt * condc * area / dist / visc[cellID] * 1000 * \
				                          (Temp[cellID] - temperature);// AllData->Viscosity;//Residual
		}
	}

	/////////////////////////////////////////////////////////
	//            6. Deal with Well B.C.                   //
	/////////////////////////////////////////////////////////
	void SimHM::HandleWellBC(const int NtCount, const double _simTime, const double dt){
		if (NtCount == 1) IniCurrentwelPara(_simTime);
		std::list<welbc>::iterator itrwelbclist;
		double temp1 = 0.0;
		for (itrwelbclist = AllData->welllist.begin();
			itrwelbclist != AllData->welllist.end(); ++itrwelbclist){
			int iEle = (*itrwelbclist).welloc;
			double wi = -1;
			int Ix, Iy, Iz;
			if (iEle >= AllData->ElementAmt) wi = 10e3;
			else{
				Ix = iEle % AllData->Nx;
				Iy = (iEle / AllData->Nx) % AllData->Ny;
				Iz = iEle / AllData->Nx / AllData->Ny;
				wi = cal_WI(AllData->permXX[iEle], AllData->permYY[iEle], AllData->Dx[Ix], AllData->Dy[Iy], \
					        AllData->Dz[Iz], (*itrwelbclist).CurrentwelPara.Rw, (*itrwelbclist).CurrentwelPara.SF);
			}		
			
			double* ptRsv = &Residual[geom_bh*geom_size];
			double* stLoc = JacMtr.find_RRD(iEle, 0);
			int blkH = JacMtr.RRP.bh;
			if ((*itrwelbclist).CurrentwelPara.weltype == 1){/*1>>BHP production well*/
				if (prs[iEle] >= (*itrwelbclist).CurrentwelPara.bhp){
					double Enthalpy_, rho_, Bw_, WellHeadTemperature = Temp[iEle];
					(*itrwelbclist).CurrentwelPara.Tf = WellHeadTemperature;
					if (AllData->FracingFluid == 0){//Slickwater Fracturing
						rho_ = AllData->H2O_thermProps.GetDensity(prs[iEle], Temp4Chem[iEle]);
						Enthalpy_ = AllData->H2O_thermProps.GetEnthalpy(Temp4Chem[iEle]);
						Bw_ = rhoStd / rho_;
					}
					else if (AllData->FracingFluid == 1){//CO2 Fracturing
						double Temperature_ = WellHeadTemperature + Tt;
						rho_ = AllData->CO2_thermProps.getDensity(Temperature_, (*itrwelbclist).CurrentwelPara.CtrParameter);
						Enthalpy_ = AllData->CO2_thermProps.getEnthalpy(Temperature_, rho_);
						Bw_ = rhoStd_CO2 / rho_;
					}
					else if (AllData->FracingFluid == 2){//N2 Fracturing
						double Temperature_ = WellHeadTemperature + Tt;
						AllData->N2_thermProps.getPropsByP2TIni(Temperature_, (*itrwelbclist).CurrentwelPara.CtrParameter);
						rho_ = AllData->N2_thermProps.density_;
						Enthalpy_ = AllData->N2_thermProps.Enthalpy;
						Bw_ = rhoStd_N2 / rho_;
					}
					
					//fluid flow
					double tempdF = wi * dt;
					double deltaP = prs[iEle] - (*itrwelbclist).CurrentwelPara.bhp;
					stLoc[0] += tempdF * uBw(iEle) + tempdF * duBwdp(iEle) * deltaP;
					double TempQf = tempdF * uBw(iEle) * deltaP;
					(*itrwelbclist).CurrentwelPara.cnsRT_f = 1.0E-12*TempQf*Bw_ / dt;//rsv_state flow rate:m3/s
					ptRsv[blkH*iEle] -= TempQf;
					//cout << "Production:" << TempQf << endl;
					
					//heat flow-advection
					stLoc[2] += WellHeadTemperature * 1.0E-12*(tempdF * uBw(iEle) + tempdF * duBwdp(iEle) * deltaP)*Bw_;
					stLoc[3] += 1.0E-12*TempQf*Bw_;
					double TempQh = WellHeadTemperature * (*itrwelbclist).CurrentwelPara.cnsRT_f*dt;
					(*itrwelbclist).CurrentwelPara.cnsRT_H = WellHeadTemperature;//Ca2+ C:mol/m3
					ptRsv[blkH*iEle + 1] -= TempQh;
				}
			}
			else if ((*itrwelbclist).CurrentwelPara.weltype == 4){/*4>>BHP injection well*/
				if (prs[iEle] <= (*itrwelbclist).CurrentwelPara.bhp){
					double Enthalpy_, rho_, Bw_;
					double WellHeadTemperature, tempVisco, BHP = (*itrwelbclist).CurrentwelPara.bhp;
					WellHeadTemperature = (*itrwelbclist).CurrentwelPara.Tf;
					double deltaP = (*itrwelbclist).CurrentwelPara.bhp - prs[iEle];
					if (AllData->FracingFluid == 0){//Slickwater Fracturing
						rho_ = AllData->H2O_thermProps.GetDensity(BHP, Temp4Chem[iEle]);
						Enthalpy_ = AllData->H2O_thermProps.GetEnthalpy(Temp4Chem[iEle]);
						Bw_ = rhoStd / rho_;
						tempVisco = AllData->H2O_thermProps.claVisco(BHP, Temp4Chem[iEle]);
					}
					else if (AllData->FracingFluid == 1){//CO2 Fracturing
						double Temperature_ = WellHeadTemperature + Tt;
						rho_ = AllData->CO2_thermProps.getDensity(Temperature_, BHP);
						Enthalpy_ = AllData->CO2_thermProps.getEnthalpy(Temperature_, rho_);
						Bw_ = rhoStd_CO2 / rho_;
						tempVisco = AllData->CO2_thermProps.getvisco(Temperature_, rho_);
					}
					else if (AllData->FracingFluid == 2){//N2 Fracturing
						double Temperature_ = WellHeadTemperature + Tt;
						AllData->N2_thermProps.getPropsByP2TIni(Temperature_, BHP);
						rho_ = AllData->N2_thermProps.density_;
						Enthalpy_ = AllData->N2_thermProps.Enthalpy;
						Bw_ = rhoStd_N2 / rho_;
						tempVisco = AllData->N2_thermProps.visco;
					}
					
					//fluid flow
					double tempdF = wi * dt;
					double uBw_ = 1.0 / Bw_ / tempVisco, rhohu_ = rho_*Enthalpy_ / tempVisco;
					stLoc[0] += tempdF * uBw_;
					double TempQf = tempdF * uBw_ * deltaP;
					(*itrwelbclist).CurrentwelPara.cnsRT_f = -1.0E-12*TempQf*Bw_ / dt;//rsv_state flow rate:m3/s
					ptRsv[blkH*iEle] += TempQf;
					//cout << "Injection:" << TempQf << endl;

					//heat flow-advection
					stLoc[2] += WellHeadTemperature * 1.0E-12*tempdF * uBw_*Bw_;
					double TempQh = WellHeadTemperature * (*itrwelbclist).CurrentwelPara.cnsRT_f*dt;
					(*itrwelbclist).CurrentwelPara.cnsRT_H = WellHeadTemperature;//Ca2+ C:mol/m3
					ptRsv[blkH*iEle + 1] += TempQh;
				}
			}
			else if ((*itrwelbclist).CurrentwelPara.weltype == 2 || (*itrwelbclist).CurrentwelPara.weltype == 3){/*2>>Constant injection/production well*/
				double tempD = 1.0e12 * (*itrwelbclist).CurrentwelPara.cnsRT_f * dt; //mD = 1.0e-15 m^2     <-- fluid
				double Enthalpy_, rho_, Bw_, BHP;
				bool isInjWell = (*itrwelbclist).CurrentwelPara.cnsRT_f < 0.0;
				double WellHeadTemperature;
				if (isInjWell) WellHeadTemperature = (*itrwelbclist).CurrentwelPara.Tf;
				else WellHeadTemperature = Temp0[iEle];
								
				//***************************Iterate BHP accorging to given rate***************************//
				if (AllData->FracingFluid == 0){//Slickwater Fracturing
					rho_ = rho_f[iEle];
					Bw_ = rhoStd / rho_;
					BHP = prs[iEle] - tempD * visc[iEle] * Bw_ / wi / dt;   //mD/mP.S = 1.0e-12 N.s <-- heat
					if (isInjWell){
						double BHP_ref, rho_ref, Bw_ref, visc_ref;
						for (int i = 0;; i++){
							visc_ref = AllData->H2O_thermProps.claVisco(BHP/10e6, Temp4Chem[iEle]+Tt);
							rho_ref = AllData->H2O_thermProps.GetDensity(BHP, Temp4Chem[iEle]);
							Bw_ref = rhoStd / rho_ref;
							BHP_ref = prs[iEle] - tempD * visc_ref * Bw_ref / wi / dt;
							if (abs(BHP_ref - BHP) / BHP_ref < 1.0e-4 || i>10){
								BHP = BHP_ref;
								rho_ = rho_ref;
								Bw_ = Bw_ref;								
								break;
							}
							else{
								BHP = BHP_ref;
							}
						}
					} else;
					//Enthalpy_ = (*itrwelbclist).CurrentwelPara.Tf;
				}
				else if (AllData->FracingFluid == 1){//CO2 Fracturing
					rho_ = rho_f[iEle];
					Bw_ = rhoStd_CO2 / rho_;
					BHP = prs[iEle] - tempD * visc[iEle] * Bw_ / wi / dt;   //mD/mP.S = 1.0e-12 N.s <-- heat
					double Temperature_ = WellHeadTemperature + Tt;
					if (isInjWell){
						double BHP_ref, rho_ref, Bw_ref, visc_ref;
						for (int i = 0;; i++){
							visc_ref = AllData->CO2_thermProps.getvisco(WellHeadTemperature, BHP);
							rho_ref = AllData->CO2_thermProps.getDensity(WellHeadTemperature, BHP);
							Bw_ref = rhoStd / rho_ref;
							BHP_ref = prs[iEle] - tempD * visc_ref * Bw_ref / wi / dt;
							if (abs(BHP_ref - BHP) / BHP_ref < 1.0e-4 || i>10){
								BHP = BHP_ref;
								rho_ = rho_ref;
								Bw_ = Bw_ref;
								break;
							}
							else{
								BHP = BHP_ref;
							}
						}
					} else;
					Enthalpy_ = AllData->CO2_thermProps.getEnthalpy(Temperature_, rho_);
				}
				else if (AllData->FracingFluid == 2){//N2 Fracturing
					rho_ = rho_f[iEle];
					Bw_ = rhoStd_N2 / rho_;
					BHP = prs[iEle] - tempD * visc[iEle] * Bw_ / wi / dt;   //mD/mP.S = 1.0e-12 N.s <-- heat
					double Temperature_ = WellHeadTemperature + Tt;
					if (isInjWell){
						double BHP_ref, rho_ref, Bw_ref, visc_ref;
						for (int i = 0;; i++){
							AllData->N2_thermProps.getPropsByP2TIni(WellHeadTemperature, BHP);
							visc_ref = AllData->N2_thermProps.visco;
							rho_ref = AllData->N2_thermProps.density_;
							Bw_ref = rhoStd_N2 / rho_ref;
							BHP_ref = prs[iEle] - tempD * visc_ref * Bw_ref / wi / dt;
							if (abs(BHP_ref - BHP) / BHP_ref < 1.0e-4 || i>10){
								BHP = BHP_ref;
								rho_ = rho_ref;
								Bw_ = Bw_ref;
								break;
							}
							else{
								BHP = BHP_ref;
							}
						}
					} else;
					AllData->N2_thermProps.getPropsByP2TIni(Temperature_, rho_);
					Enthalpy_ = AllData->N2_thermProps.Enthalpy;
				}

				(*itrwelbclist).CurrentwelPara.bhp = BHP;
				//**********************************************************************************************//
				double deltaP = prs[iEle] - (*itrwelbclist).CurrentwelPara.CtrParameter;
				if (isInjWell){ //Negtive value, "-", denotes injection well
					if (true/*BHP >= (*itrwelbclist).CurrentwelPara.CtrParameter*/){
						//Fluid flow
						double TempQf = tempD / Bw_;
						(*itrwelbclist).CurrentwelPara.cnsRT_f = 1.0E-12*TempQf*Bw_ / dt;//rsv_state flow rate:m3/s
						ptRsv[blkH*iEle] -= TempQf;  //"+" or "-" standard volume

						//Heat flow
						double TempQh = WellHeadTemperature * (*itrwelbclist).CurrentwelPara.cnsRT_f*dt;
						//(*itrwelbclist).CurrentwelPara.cnsRT_H = 1.0E-3*TempQh / dt;//STD Ca2+ flux rate:mol/s
						ptRsv[blkH*iEle + 1] -= TempQh; //rho_f*hf*Vf;  "+" or "-" energy
					}
					else{
						//fluid flow
						double tempdF = wi / visc[iEle] * dt;
						stLoc[0] += tempdF / Bw_;
						double TempQf = tempdF / Bw[iEle] * deltaP;
						(*itrwelbclist).CurrentwelPara.cnsRT_f = 1.0E-12*TempQf*Bw_ / dt;//rsv_state flow rate:m3/s
						ptRsv[blkH*iEle] -= TempQf;

						//heat flow-advection
						double tempdH = 1.0e-9 * enthalpy[iEle] * rho_f[iEle] * tempdF;
						stLoc[2] += tempdH;
						double TempQh = tempdH * deltaP;
						(*itrwelbclist).CurrentwelPara.cnsRT_H = 1.0E-3*TempQh / dt;//STD heat rate:KJ/s
						ptRsv[blkH*iEle + 1] -= TempQh;
					}
				}
				else {// Positive value, "+", denotes production well
					if (true/*BHP <= (*itrwelbclist).CurrentwelPara.CtrParameter*/){
						//Fluid flow
						double TempQf = tempD / Bw_;
						(*itrwelbclist).CurrentwelPara.cnsRT_f = 1.0E-12*TempQf*Bw_ / dt;//rsv_state flow rate:m3/s
						ptRsv[blkH*iEle] -= TempQf;

						//Heat flow
						double TempQh = WellHeadTemperature * (*itrwelbclist).CurrentwelPara.cnsRT_f * dt;
						//(*itrwelbclist).CurrentwelPara.cnsRT_H = 1.0E-3*TempQh / dt;//STD Ca2+ flux rate:mol/s
						ptRsv[blkH*iEle + 1] -= TempQh; //rho_f*hf*Vf;  "+" or "-" energy
					}
					else{
						//fluid flow
						double tempdF = wi / visc[iEle] * dt;
						stLoc[0] += tempdF / Bw_;
						double TempQf = tempdF / Bw[iEle] * deltaP;
						(*itrwelbclist).CurrentwelPara.cnsRT_f = 1.0E-12*TempQf*Bw_ / dt;//rsv_state flow rate:m3/s
						ptRsv[blkH*iEle] -= TempQf;

						//heat flow-advection
						double tempdH = 1.0e-9 * rho_ * enthalpy[iEle] * tempdF;
						stLoc[2] += tempdH;
						double TempQh = tempdH * deltaP;
						(*itrwelbclist).CurrentwelPara.cnsRT_H = 1.0E-3*TempQh / dt;//STD heat rate:KJ/s
						ptRsv[blkH*iEle + 1] -= TempQh;
					}
				}
			}

			//else if ((*itrwelbclist).weltype == 3){
			//	double EnthalpyStd = AllData->thermProps.getEnthalpyStd((*itrwelbclist).Tf);
			//	//fluid flow
			//	ptRsv[blkH*iEle] += 1.0e15*(*itrwelbclist).cnsRT_f*dt;   //Residual
			//	//heat flow
			//	ptRsv[blkH*iEle + 1] += 1.0e-12 * rhoStd * EnthalpyStd *(*itrwelbclist).cnsRT_f * dt;  //Residual
			//}

			else cout << "Do not have this kind of well control!";
		}
	}

	void SimHM::IniCurrentwelPara(const double _simTime){
		std::list<welbc>::iterator itrwelbclist;
		double temp1 = 0.0;
		for (itrwelbclist = AllData->welllist.begin();
			itrwelbclist != AllData->welllist.end(); ++itrwelbclist){
			std::list<welPara>::iterator itrwelParaList;
			for (itrwelParaList = (*itrwelbclist).welParaList.begin();
				itrwelParaList != (*itrwelbclist).welParaList.end(); ++itrwelParaList){
				std::list<welPara>::iterator itrwelParaList_ = itrwelParaList;
				++itrwelParaList_;
				bool either = itrwelParaList_ == (*itrwelbclist).welParaList.end(),
					     or = itrwelParaList_ != (*itrwelbclist).welParaList.end() &&
					(*itrwelParaList).workTime <= _simTime && (*itrwelParaList_).workTime > _simTime;
				if (either || or){
					(*itrwelbclist).CurrentwelPara.bhp = (*itrwelParaList).bhp;
					(*itrwelbclist).CurrentwelPara.cnsRT_f = (*itrwelParaList).cnsRT_f;
					(*itrwelbclist).CurrentwelPara.CtrParameter = (*itrwelParaList).CtrParameter;
					(*itrwelbclist).CurrentwelPara.Rw = (*itrwelParaList).Rw;
					(*itrwelbclist).CurrentwelPara.SF = (*itrwelParaList).SF;
					(*itrwelbclist).CurrentwelPara.Tf = (*itrwelParaList).Tf;
					(*itrwelbclist).CurrentwelPara.cnsRT_H = (*itrwelParaList).cnsRT_H;
					(*itrwelbclist).CurrentwelPara.weltype = (*itrwelParaList).weltype;
					(*itrwelbclist).CurrentwelPara.workTime = (*itrwelParaList).workTime;
					break;
				}
			}
		}
	}

	//*****************************Calculate well index***************************************//
	double cal_WI(double kx, double ky, double dx, double dy, double dz, double rw, double sf){
		double tmp = 0.28*pow(pow(ky / kx, 0.5)*dx*dx + pow(kx / ky, 0.5)*dy*dy, 0.5)  \
			         / (pow(ky / kx, 0.25) + pow(kx / ky, 0.25));
		tmp = 2 * PI*pow(kx*ky, 0.5)*dz / (log(tmp / rw) + sf);
		return tmp;
	}

	void SimHM::source_(const double dt, double fsorc){
		std::list<Element>::iterator ElmItr;
		for (ElmItr = AllData->Elmlist.begin(); ElmItr != AllData->Elmlist.end(); ++ElmItr){
			if (ElmItr == AllData->Elmlist.begin()){
				Residual[geom_bh*geom_size] += 1.0e15*fsorc*dt;
			};
		}
	}

	void GenEleProp(const int I, const int J, const int K){
		ofstream OUTPUT("YoungsModulus.dat", ios::out);
		if (!OUTPUT){
			cerr << "open error!" << endl;
			exit(1);
		}
		OUTPUT << "Young's Modulus (Gpa) \n";
		int counter = 0;
		for (int k = 0; k < K; ++k){
			for (int j = 0; j < J; ++j){
				for (int i = 0; i < I; ++i){
					//int nodeNum = i + j*I + k*I*J;		
					if (i>9 && i<15){
						OUTPUT << setiosflags(ios::left) << setw(6) << 10.0;
					}
					else if (i<10) {
						if (k < 3){
							OUTPUT << setiosflags(ios::left) << setw(6) << 30.0;
						}
						else if (k == 27){
							OUTPUT << setiosflags(ios::left) << setw(6) << 20.0;
						}
						else{
							OUTPUT << setiosflags(ios::left) << setw(6) << 15.0;
						}
					}
					else{
						if (k == 0){
							OUTPUT << setiosflags(ios::left) << setw(6) << 30.0;
						}
						else if (k >24){
							OUTPUT << setiosflags(ios::left) << setw(6) << 20.0;
						}
						else{
							OUTPUT << setiosflags(ios::left) << setw(6) << 15.0;
						}
					}
					//////////////////////////////////////////////
					if (counter == 0);
					else if (counter % 20 == 19) OUTPUT << endl;
					++counter;
				}
			}
		}
		OUTPUT.close();
	}

	void DataBase::GenConcBCs(){
		fixConc fixConcList;
		bool isConcUpper = ConcUpper < 0.0;
		bool isConcLower = ConcLower < 0.0;
		bool isConcFront = ConcFront < 0.0;
		bool isConcBack = ConcBack < 0.0;
		bool isConcLeft = ConcLeft < 0.0;
		bool isConcRight = ConcRight < 0.0;

		//***************************** Concentration B.C.************************************//
		for (int k = 0; k < Nz; ++k){
			for (int j = 0; j < Ny; ++j){
				for (int i = 0; i < Nx; ++i){
					int EleNum = i + j*Nx + k*Nx*Ny;
					if (!isConcUpper && k == Nz - 1){
						fixConcList.pEle_num = EleNum;
						fixConcList.pDirc = 3;
						fixConcList.Conc = ConcUpper;
						Conclist.push_back(fixConcList);
					};
					if (!isConcLower && k == 0){
						fixConcList.pEle_num = EleNum;
						fixConcList.pDirc = 3;
						fixConcList.Conc = ConcLower;
						Conclist.push_back(fixConcList);
					};
					if (!isConcFront && j == 0){
						fixConcList.pEle_num = EleNum;
						fixConcList.pDirc = 2;
						fixConcList.Conc = ConcFront;
						Conclist.push_back(fixConcList);
					};
					if (!isConcBack && j == Ny - 1){
						fixConcList.pEle_num = EleNum;
						fixConcList.pDirc = 2;
						fixConcList.Conc = ConcBack;
						Conclist.push_back(fixConcList);
					};
					if (!isConcLeft && i == 0){
						fixConcList.pEle_num = EleNum;
						fixConcList.pDirc = 1;
						fixConcList.Conc = ConcLeft;
						Conclist.push_back(fixConcList);
					};
					if (!isConcRight && i == Nx - 1){
						fixConcList.pEle_num = EleNum;
						fixConcList.pDirc = 1;
						fixConcList.Conc = ConcRight;
						Conclist.push_back(fixConcList);
					};
				}
			}
		}
	}

	void DataBase::GenPresBCs(){
		fixPre fixPreList;
		bool isPresUpper = PresUpper < 0.0;
		bool isPresLower = PresLower < 0.0;
		bool isPresFront = PresFront < 0.0;
		bool isPresBack = PresBack < 0.0;
		bool isPresLeft = PresLeft < 0.0;
		bool isPresRight = PresRight < 0.0;

		//***************************** Pressure B.C.************************************//
		for (int k = 0; k < Nz; ++k){
			for (int j = 0; j < Ny; ++j){
				for (int i = 0; i < Nx; ++i){
					int EleNum = i + j*Nx + k*Nx*Ny;
					
					H2O_thermProps.getPropsByP2TIni(refPres, T_ref);

					if (!isPresUpper && k == Nz - 1){
						fixPreList.pEle_num = EleNum;
						fixPreList.pDirc = 3;
						fixPreList.pres = PresUpper + (iDEPTH_[EleNum] + 0.5*Dz[0]) / H2O_thermProps.Volume *cntg;

						fixPreList.perm = permZZ[EleNum];
						fixPreList.area=Dx[i]*Dy[j];
						fixPreList.dist = 1.0e-15;// 0.5*Dz[k];
						fixPreList.height = (Dx[i] + Dy[j] + Dz[k]) / 3;

						Preslist.push_back(fixPreList);
					};
					if (!isPresLower && k == 0){
						fixPreList.pEle_num = EleNum;
						fixPreList.pDirc = 3;
						fixPreList.pres = PresLower + (iDEPTH_[EleNum] + 0.5*Dz[0]) / H2O_thermProps.Volume *cntg;

						fixPreList.perm = permZZ[EleNum];
						fixPreList.area = Dx[i] * Dy[j];
						fixPreList.dist = 1.0e-15;// 0.5*Dz[k];
						fixPreList.height = (Dx[i] + Dy[j] + Dz[k]) / 3;

						Preslist.push_back(fixPreList);
					};
					if (!isPresFront && j == 0){
						fixPreList.pEle_num = EleNum;
						fixPreList.pDirc = 2;
						fixPreList.pres = PresFront + (iDEPTH_[EleNum] + 0.5*Dz[0]) / H2O_thermProps.Volume *cntg;

						fixPreList.perm = permYY[EleNum];
						fixPreList.area = Dx[i] * Dz[k];
						fixPreList.dist = 1.0e-15;// 0.5*Dy[j];
						fixPreList.height = (Dx[i] + Dy[j] + Dz[k]) / 3;

						Preslist.push_back(fixPreList);
					};
					if (!isPresBack && j == Ny - 1){
						fixPreList.pEle_num = EleNum;
						fixPreList.pDirc = 2;
						fixPreList.pres = PresBack + (iDEPTH_[EleNum] + 0.5*Dz[0]) / H2O_thermProps.Volume *cntg;

						fixPreList.perm = permYY[EleNum];
						fixPreList.area = Dx[i] * Dz[k];
						fixPreList.dist = 1.0e-15;// 0.5*Dy[j];
						fixPreList.height = (Dx[i] + Dy[j] + Dz[k]) / 3;

						Preslist.push_back(fixPreList);
					};
					if (!isPresLeft && i == 0){
						fixPreList.pEle_num = EleNum;
						fixPreList.pDirc = 1;
						fixPreList.pres = PresLeft + (iDEPTH_[EleNum] + 0.5*Dz[0]) / H2O_thermProps.Volume *cntg;

						fixPreList.perm = permXX[EleNum];
						fixPreList.area = Dy[j] * Dz[k];
						fixPreList.dist = 1.0e-15;// 0.5*Dx[i];
						fixPreList.height = (Dx[i] + Dy[j] + Dz[k]) / 3;

						Preslist.push_back(fixPreList);
					};
					if (!isPresRight && i == Nx - 1){
						fixPreList.pEle_num = EleNum;
						fixPreList.pDirc = 1;
						fixPreList.pres = PresRight + (iDEPTH_[EleNum] + 0.5*Dz[0]) / H2O_thermProps.Volume *cntg;

						fixPreList.perm = permXX[EleNum];
						fixPreList.area = Dy[j] * Dz[k];
						fixPreList.dist = 1.0e-15;// 0.5*Dx[i];
						fixPreList.height = (Dx[i] + Dy[j] + Dz[k]) / 3;

						Preslist.push_back(fixPreList);
					};
				}
			}
		}
	}

	void DataBase::GenDisp2LoadBCs(){
		Trac _Trac;
		Constrain _Constrainlist;
		bool isDispUpper = DispUpper < 1.0;
		bool isDispLower = DispLower < 1.0;
		bool isDispFront = DispFront < 1.0;
		bool isDispBack = DispBack < 1.0;
		bool isDispLeft = DispLeft < 1.0;
		bool isDispRight = DispRight < 1.0;

		//*****************************1. disp B.C.************************************//
		for (int k = 0; k < Nz + 1; ++k){
			for (int j = 0; j < Ny + 1; ++j){
				for (int i = 0; i < Nx + 1; ++i){
					int nodeNum = i + j*(Nx + 1) + k*(Nx + 1)*(Ny + 1);					
					if (isDispUpper && k == Nz){
						_Constrainlist.cNd_num = nodeNum;
						_Constrainlist.cDirc = 3;
						_Constrainlist.Dfm = DispUpper;
						Cnstrlist.push_back(_Constrainlist);
					};
					if (isDispLower && k == 0){
						_Constrainlist.cNd_num = nodeNum;
						_Constrainlist.cDirc = 3;
						_Constrainlist.Dfm = DispLower;
						Cnstrlist.push_back(_Constrainlist);
					};
					if (isDispFront && j == 0){
						_Constrainlist.cNd_num = nodeNum;
						_Constrainlist.cDirc = 2;
						_Constrainlist.Dfm = DispFront;
						Cnstrlist.push_back(_Constrainlist);
					};
					if (isDispBack && j == Ny){
						_Constrainlist.cNd_num = nodeNum;
						_Constrainlist.cDirc = 2;
						_Constrainlist.Dfm = DispBack;
						Cnstrlist.push_back(_Constrainlist);
					};
					if (isDispLeft && i == 0){
						_Constrainlist.cNd_num = nodeNum;
						_Constrainlist.cDirc = 1;
						_Constrainlist.Dfm = DispLeft;
						Cnstrlist.push_back(_Constrainlist);
					};
					if (isDispRight && i == Nx){
						_Constrainlist.cNd_num = nodeNum;
						_Constrainlist.cDirc = 1;
						_Constrainlist.Dfm = DispRight;
						Cnstrlist.push_back(_Constrainlist);
					};
				}
			}
		}
		//*****************************2. Traction B.C.************************************//
		double refDepth = iDEPTH_[ElementAmt - 1] - 0.5*Dz[Nz - 1];
		for (int k = 0; k < Nz; ++k){
			for (int j = 0; j < Ny; ++j){
				for (int i = 0; i < Nx; ++i){
					int EleNum = i + j*Nx + k*Nx*Ny;
					if (!isDispUpper && k == Nz-1){
						_Trac.tEle_num = EleNum;
						_Trac.tDirc = 3;
						_Trac.tFaceNum = 0;
						_Trac.traction = -SvTop;
						Traclist.push_back(_Trac);
					};
					if (!isDispLower && k == 0){
						_Trac.tEle_num = EleNum;
						_Trac.tDirc = 3;
						_Trac.tFaceNum = 1;
						_Trac.traction = SvTop;
						Traclist.push_back(_Trac);
					};
					if (!isDispFront && j == 0){
						double py = SyRef - (refDepth - iDEPTH_[EleNum]) * dpydz;
						_Trac.tEle_num = EleNum;
						_Trac.tDirc = 2;
						_Trac.tFaceNum = 2;
						_Trac.traction = py;
						Traclist.push_back(_Trac);
					};
					if (!isDispBack && j == Ny-1){
						double py = SyRef - (refDepth - iDEPTH_[EleNum]) * dpydz;
						_Trac.tEle_num = EleNum;
						_Trac.tDirc = 2;
						_Trac.tFaceNum = 3;
						_Trac.traction = -py;
						Traclist.push_back(_Trac);
					};
					if (!isDispLeft && i == 0){
						double px = SxRef - (refDepth - iDEPTH_[EleNum]) * dpxdz;
						_Trac.tEle_num = EleNum;
						_Trac.tDirc = 1;
						_Trac.tFaceNum = 4;
						_Trac.traction = px;
						Traclist.push_back(_Trac);
					};
					if (!isDispRight && i == Nx-1){
						double px = SxRef - (refDepth - iDEPTH_[EleNum]) * dpxdz;
						_Trac.tEle_num = EleNum;
						_Trac.tDirc = 1;
						_Trac.tFaceNum = 5;
						_Trac.traction = -px;
						Traclist.push_back(_Trac);
					};
				}
			}
		}
    }

	//***************Displacement boundary condition disposal******************// 
	void DataBase::GenDisNode(const int I, const int J, const int K, const double disp, 
			                  const double dpxdz, const double dpydz, 
						      const double pxRef, const double pyRef, const double SvRef)
	{
		//for THM fault problem
		ofstream OUTPUT("DispSeriM.dat", ios::out);
		if (!OUTPUT){
			cerr << "open error!" << endl;
			exit(1);
		}
		OUTPUT << "*****************************disp B.C.************************************\n";
		for (int k = 0; k < K; ++k){
			for (int j = 0; j < J; ++j){
				for (int i = 0; i < I; ++i){
					int nodeNum = i + j*I + k*I*J;
					//OUTPUT << nodeNum << "   " << "   2    " << disp << endl;
					if (i == 0){
						OUTPUT << nodeNum << "   " << "   1    " << disp << endl;
					};
					//if (j == J-1){
					//	OUTPUT << nodeNum << "   " << "   2    " << disp << endl;
					//};
					if (k == 0){
						OUTPUT << nodeNum << "   " << "   3    " << disp << endl;
					};
				}
			}
		}

		//For 5-cluster stage case
		OUTPUT << "*****************************Traction B.C.************************************\n";
		double refDepth = -iDEPTH_[ElementAmt-1] + 0.5*Dz[Nz-1];
		for (int k = 0; k < K - 1; ++k){
			for (int j = 0; j < J - 1; ++j){
				for (int i = 0; i < I - 1; ++i){
					int EleNum = i + j*(I - 1) + k*(I - 1)*(J - 1);	
					//if (i == 0) {
					//	double px = -pxRef - (refDepth + iDEPTH[EleNum]) * dpxdz;
					//	OUTPUT << EleNum << "            1          4           " << -px << endl;
					//};
					if (i == I - 2) {
						double px = -pxRef - (refDepth + iDEPTH_[EleNum]) * dpxdz;
						OUTPUT << EleNum << "            1          5          " << px << endl;
					};
					if (j == 0) {
						double py = pyRef + (refDepth + iDEPTH_[EleNum]) * dpydz;
						OUTPUT << EleNum << "            2          2           "<< py << endl;
					};
					if (j == J - 2) {
						double py = pyRef + (refDepth + iDEPTH_[EleNum]) * dpydz;
						OUTPUT << EleNum << "            2          3          " << -py << endl;
					};
					if (k == K - 2) {
						OUTPUT << EleNum << "            3          0          "<<-SvRef << endl;
					};
				}
			}
		}

		OUTPUT.close();
	}
	//****************************************************************************************//

	void SimHM::Simulate_testR(){
		//*******************************************************************************************
		//rsv test
		strainV = new double[rsv_size]; memset(strainV, 0x0, rsv_size*sizeof(double));
		double ddt = 10.0, sumT = 0.0, real_dt = 0.0;
		bool dx_acpt = true;
		int NtCount = 0;
		Init(AllData->isRestart);
		while (true){
			real_dt += ddt;
			sumT += real_dt;
			for (NtCount = 1;; ++NtCount){

				if (dx_acpt && NtCount > 1){
					--NtCount;
					break;
				}

				if (Jval == NULL){
					int CSR_nnz = JacMtr.RRP.NNZ*JacMtr.RRP.bh*JacMtr.RRP.bw;
					int RowNum = rsv_size*JacMtr.RRP.bh;
					double* resi = &Residual[geom_bh*geom_size];
					Jval = new double[CSR_nnz];
					colind = new int[CSR_nnz];
					rowptr = new int[RowNum + 1];
					setzeor();
					//source_(5.0e-8, 0.01);
					ActiveFlow(real_dt);
					AssembleGRG(real_dt);
					memset(disp, 0x0, geom_bh*geom_size*sizeof(double));
					memset(Residual, 0x0, (geom_size*geom_bh + rsv_size*rsv_bh)*sizeof(double));
					FixPrsBC(real_dt);
					for (int i = geom_bh*geom_size; i < geom_bh*geom_size + rsv_size*rsv_bh; ++i){
						cout << Residual[i + geom_bh*geom_size] << endl;
					}
					JacMtr.RRP.comb_to_CSR(rowptr, colind, JacMtr.RRP.bh, JacMtr.RRP.bw);
					JacMtr.RRP.copy_to_CSR(Jval, rowptr, JacMtr.RRP.bh, JacMtr.RRP.bw);
					directsoler.Init_Pardiso(RowNum, Jval, rowptr, colind, rowptr[0],
						KarstVolSysSOL::Pardiso_USYM, 1, KarstVolSysSOL::Pardiso_NOCGS, 1);
					directsoler.Pardiso_Solve(RowNum, Jval, rowptr, colind, resi);
					for (int i = 0; i < rsv_size; ++i){
						cout << resi[i] << endl;
					}
					dx_acpt = CheckDx(real_dt);
					UpdateXp();
					for (int i = 0; i < rsv_size; ++i){
						cout << prs[i] << endl;
					}
					const double* ptP = &Residual[geom_size*geom_bh];
					//UpdateRsvProp(ptP);
					ofstream OUTPUT("JACRSV.dat", ios::out);
					if (!OUTPUT){ cerr << "open error!" << endl; exit(1); }
					OUTPUT << " >> Chech R-R block Jacbion Matxir!" << endl;
					for (int i = 0, j = 0, k = 0; i < 270 * 270, j < CSR_nnz; ++k, i++)
					{
						if (k == colind[j]){ OUTPUT << setw(10) << Jval[j]; j++; }
						else OUTPUT << setw(10) << " ";
						if (i % 270 == 269){ k = -1; OUTPUT << endl; }
						else;
					}
					OUTPUT.close();
					resi = NULL;
					Jval = NULL;
					colind = NULL;
					rowptr = NULL;
				}
			}
			UpdateX0(1.0);
		}
		//********************************************************************************************
	}

	void SimHM::Simulate_testG(){
		//*******************************************************************************************
		//rsv test
		double ddt = 5.0, sumT = 0.0, real_dt = 0.0;
		bool dx_acpt = true;
		int NtCount = 0;
		Init(AllData->isRestart);
		while (true){
			real_dt += ddt;
			sumT += real_dt;
			for (NtCount = 1;; ++NtCount){

				if (dx_acpt && NtCount > 1){
					--NtCount;
					break;
				}

				if (Jval == NULL){
					int CSR_nnz = JacMtr.GGP.NNZ*JacMtr.GGP.bh*JacMtr.GGP.bw;
					int RowNum = geom_size*JacMtr.GGP.bh;
					double* resi = Residual;
					Jval = new double[CSR_nnz];
					colind = new int[CSR_nnz];
					rowptr = new int[RowNum + 1];
					setzeor();
					//source_(5.0e-8, 0.01);
					ActiveSolid();
					AssembleGRG(real_dt);
					HandleLoadBC(); //Load boundary-->Geom
					HandleDispBC(); //Displacement boundary-->Geom
					for (int i = 0; i < geom_bh*geom_size; ++i){
						cout << Residual[i] << endl;
					}
					JacMtr.GGP.comb_to_CSR(rowptr, colind, JacMtr.GGP.bh, JacMtr.GGP.bw);
					JacMtr.GGP.copy_to_CSR(Jval, rowptr, JacMtr.GGP.bh, JacMtr.GGP.bw);
					directsoler.Init_Pardiso(RowNum, Jval, rowptr, colind, rowptr[0],
						KarstVolSysSOL::Pardiso_SSYM, 1, KarstVolSysSOL::Pardiso_NOCGS, 1);
					directsoler.Pardiso_Solve(RowNum, Jval, rowptr, colind, resi);
					for (int i = 0; i < geom_bh*geom_size; ++i){
						cout << resi[i] << endl;
					}

					calStrainStress();

					dx_acpt = CheckDx(real_dt);
					UpdateXp();

					ofstream OUTPUT("JACGeom.dat", ios::out);
					if (!OUTPUT){ cerr << "open error!" << endl; exit(1); }
					OUTPUT << " >> Chech G-G block Jacbion Matxir!" << endl;
					for (int i = 0, j = 0, k = 0; i < geom_bh*geom_size*geom_bh*geom_size, j < CSR_nnz; ++k, i++)
					{
						if (k == colind[j]){ OUTPUT << setw(10) << Jval[j]; j++; }
						else OUTPUT << setw(10) << " ";
						if (i % geom_bh*geom_size == geom_bh*geom_size - 1){ k = -1; OUTPUT << endl; }
						else;
					}
					OUTPUT.close();
					resi = NULL;
					Jval = NULL;
					colind = NULL;
					rowptr = NULL;
				}
			}
			UpdateX0(1.0);
		}
	}

	void SimHM::Simulate_testGR(){
		//*******************************************************************************************
		//rsv test
		double ddt = 5.0, sumT = 0.0, real_dt = 0.0;
		bool dx_acpt = true;
		int NtCount = 0;
		Init(AllData->isRestart);
		while (true){
			real_dt += ddt;
			sumT += real_dt;
			for (NtCount = 1;; ++NtCount){

				if (dx_acpt && NtCount > 1){
					--NtCount;
					break;
				}

				if (Jval == NULL){
					int CSR_nnz = JacMtr.Get_ACTnnz(),
						countGeo = geom_size * geom_bh,
						countRsv = rsv_size  * rsv_bh,
						RowNum = countGeo + countRsv;
					double* resi = Residual;
					Jval = new double[CSR_nnz];
					colind = new int[CSR_nnz];
					rowptr = new int[RowNum + 1];
					setzeor();
					ActiveFlow(real_dt); //R-R
					ActiveSolid();  //G-G
					AssembleGRG(real_dt);  //G-R & R-G
					FixPrsBC(real_dt);   //Constant pressure boundary-->Rsv
					HandleLoadBC(); //Load boundary-->Geom
					HandleDispBC(); //Displacement boundary-->Geom

					for (int i = 0; i < geom_bh*geom_size + rsv_size*rsv_bh; ++i){
						cout << Residual[i] << endl;
					}

					JacMtr.comb_to_CSR(rowptr, colind);
					JacMtr.copy_to_CSR(Jval, rowptr);
					directsoler.Init_Pardiso(RowNum, Jval, rowptr, colind, rowptr[0],
						KarstVolSysSOL::Pardiso_USYM, 1, KarstVolSysSOL::Pardiso_NOCGS, 1);
					directsoler.Pardiso_Solve(RowNum, Jval, rowptr, colind, resi);
					dx_acpt = CheckDx(real_dt);
					calStrainStress();
					dx_acpt = CheckDx(real_dt);
					UpdateXp();
					const double* ptPrs = &Residual[geom_size*geom_bh];
					//UpdateRsvProp(ptPrs);
					resi = NULL;
					Jval = NULL;
					colind = NULL;
					rowptr = NULL;
				}
			}
			UpdateX0(1.0);
		}
	}
}