/**
* @file    UpdateX_Vars_IntactRev.cpp
* @author  Sanbai Li
* @date    2015/4
* @brief   Update variables during iterative or time step in intact reservoirs.
*
* @version
* Copyright (c) 2018, NEU & SINOPEC.
*/

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

namespace KarstVolSys{
	//Update temperatrue, pore pressure, displacement, porosity, permeability, 
	//       stiffness, strain, stress and thermal properties
	bool SimHM::UpdateX(double dt){
		int t_start = clock(), t_end;
		std::list<Element>::iterator iteELE;
		double* ptPrss = &Residual[geom_size*geom_bh];
		for (iteELE = AllData->Elmlist.begin(); iteELE != AllData->Elmlist.end(); ++iteELE){//1. for grid block
			Element& ele = *iteELE;
			int iEle = ele.Ele_num;
			//-----------------------------------------//
			//1. Update Pore pressure & Temperature
			delta_prs[iEle] = ptPrss[2 * iEle];
			delta_Temp[iEle] = ptPrss[2 * iEle + 1];

			prs[iEle] += delta_prs[iEle];
			Temp[iEle] += delta_Temp[iEle];				

			if (prs[iEle] > 2.00e8){
				cout << "Get trouble in MATRIX block iEle = " << iEle << "regarding Pressure >2.0e8.  Note!!!\n";
				//ptPrss[2 * iEle] = 0.0;
				//prs[iEle] -= delta_prs[iEle];
				return false;
			}
			AllData->Ca_CO2_H2O_therm_Props.getPropsByP2TIni(AllData->isOpenSys, AllData->T_ref, P_CO2[iEle], filmthick[iEle]);
			//AllData->Ca_CO2_H2O_therm_Props.getPropsByP2TIni(AllData->isOpenSys, Temp4Chem[iEle], P_CO2[iEle], filmthick[iEle]);
			Ceq[iEle] = AllData->Ca_CO2_H2O_therm_Props.Ceq;

			if (Temp[iEle] < 0.0){
				//cout << "Get trouble in MATRIX block iEle = " << iEle << " regarding C_Ca2+ < 0.0.  Note!!!\n";
				ptPrss[2 * iEle + 1] = 0.0;
				Temp[iEle] -= delta_Temp[iEle];
			}
			else if (Temp[iEle] > Ceq[iEle]){
				//cout << "Get trouble in MATRIX block iEle = " << iEle << " regarding C_Ca2+ > Ceq  Note!!!\n";
				ptPrss[2 * iEle + 1] = Ceq[iEle];
				Temp[iEle] -= delta_Temp[iEle];
			}
			else;

			//2. Update thermal properties
			UpdateProps(iEle, prs[iEle], Temp4Chem[iEle]);
			//----------------------------------------- //
			//3. Calculate element stress locating at center position
			//4. Calculate the principal stress & strain
			double* D = new double[36];
			double* resi = new double[24];
			if (ele.FrcForm || ele.isNaturalFractured){ //Stiffness matrix for both Matrix and Fractures
				if (AllData->hasNFs)
					UpdatePerm2FormMtrDinNFR(ele, D);
				else
					UpdatePerm2FormMtrD(ele, D);
			}
			else{  //Stiffness matrix for Matrix part
				calMtrD_Stiffness(ele, D);
			}
			FindDispVct(resi, ele); //Find local node numbering from global node numbering
			calEleStrainStress(D, &ele, resi);
			delete[] resi;
			delete[] D;
		}
		for (int iFace = AllData->ElementAmt; iFace<rsv_size; ++iFace){//2. for embedded faces
			//1. Update Pore pressure & Temperature
			delta_prs[iFace] = ptPrss[2 * iFace];
			delta_Temp[iFace] = ptPrss[2 * iFace + 1];

			prs[iFace] += delta_prs[iFace];
			Temp[iFace] += delta_Temp[iFace];

			if (prs[iFace] > 2.00e8){
				cout << "Get trouble in fracture-sgment block iFace = " << iFace << "regarding Pressure >2.0e8.  Note!!!\n";
				ptPrss[2 * iFace] = 0.0;
				prs[iFace] -= delta_prs[iFace];
			}
			AllData->Ca_CO2_H2O_therm_Props.getPropsByP2TIni(AllData->isOpenSys, AllData->T_ref, P_CO2[iFace], filmthick[iFace]);
			//AllData->Ca_CO2_H2O_therm_Props.getPropsByP2TIni(AllData->isOpenSys, Temp4Chem[iFace], P_CO2[iFace], filmthick[iFace]);
			Ceq[iFace] = AllData->Ca_CO2_H2O_therm_Props.Ceq;

			if (Temp[iFace] < 0.0){
				//cout << "Get trouble in fracture-sgment block iFace = " << iFace << " regarding C_Ca2+ < 0.0.  Note!!!\n";
				ptPrss[2 * iFace + 1] = 0.0;
				Temp[iFace] -= delta_Temp[iFace];
			}
			else if (Temp[iFace] > Ceq[iFace]){
				//cout << "Get trouble in fracture-sgment block iFace = " << iFace << " regarding C_Ca2+ > Ceq  Note!!!\n";
				ptPrss[2 * iFace + 1] = Ceq[iFace];
				Temp[iFace] -= delta_Temp[iFace];
			}
			else;

			//2. Update thermal properties
			UpdateProps(iFace, prs[iFace], Temp4Chem[iFace]);

			//cout << "dP[" << iFace << "]=" << delta_prs[iFace] << ";  dC[" << iFace << "]=" << delta_Temp[iFace] << endl;
		}

		//whether update MESH or not
		if (0) AllData->UpdataMESH(Residual);

		//7. Update Node displacement
		for (int i = 0; i < geom_size*geom_bh; ++i){
			delta_disp[i] = Residual[i];
			disp[i] += delta_disp[i];
		}

		//Update porosity, permeability of matrix and fracture, fracture aperture, Ceq
		for (int i = 0; i < rsv_size; ++i){//update porosity, pore size, interficial area, and permeability			
			if (i < AllData->ElementAmt){//update matrix props
				//double F_Ca = AllData->Ca_CO2_H2O_therm_Props.CalFluxRate(filmthick[i], Temp[i]);
				double ceq = AllData->Ca_CO2_H2O_therm_Props.Ceq;
				if (Temp[i]>ceq)
					Temp[i] = ceq;
				else if (Temp[i] < 0)
					Temp[i] = 0.0;
				double rhoRock = AllData->Lightdata->RockDensity_;
				double disPower = AllData->Lightdata->dissPower;
				//filmthick[i] = filmthick0[i] + 0.10009*F_Ca / rhoRock*dt; //M_CaCO3 = 0.1 kg/mol
				
				double Sh = 2.0 + 0.6*3.0 / 7.0*pow(visc[i] / rho_f[i] / AllData->heatcdXX[i], _1to3);
				double kc = 0.5 * Sh * AllData->heatcdXX[i] / rp[i];
				double F_Ca = kc*(ceq - Temp[i]);
				iPORO[i] = iPORO0[i] + F_Ca * av[i] * disPower / rhoRock * dt;

				double a = iPORO[i] / iPORO0[i];
				double b = iPORO[i] * (1.0 - iPORO0[i]) / iPORO0[i] / (1.0 - iPORO[i]);
				av[i] = (a / b) * av0[i];
				rp[i] = b * rp0[i];
				AllData->UpdatePerm(i, a*b*b);

				if (AllData->permXX[i] > 1.0e10) AllData->permXX[i] = 1.0e10;
				if (AllData->permYY[i] > 1.0e10) AllData->permYY[i] = 1.0e10;
				if (AllData->permZZ[i] > 1.0e10) AllData->permZZ[i] = 1.0e10;
				if (AllData->permXY[i] > 1.0e10) AllData->permXY[i] = 1.0e10;
				if (AllData->permYZ[i] > 1.0e10) AllData->permYZ[i] = 1.0e10;
				if (AllData->permXZ[i] > 1.0e10) AllData->permXZ[i] = 1.0e10;
			}
			else{//update frature aperture and permeability
				double F_Ca = AllData->Ca_CO2_H2O_therm_Props.CalFluxRate(filmthick[i], Temp[i]);
				double ceq = AllData->Ca_CO2_H2O_therm_Props.Ceq;
				if (Temp[i]>ceq)
					Temp[i] = ceq;
				else if (Temp[i] < 0)
					Temp[i] = 0.0;
				double rhoRock = AllData->Lightdata->RockDensity_;
				filmthick[i] = filmthick0[i] + 0.10009*F_Ca / ReactionFct / rhoRock*dt*1000;
				//if (!AllData->isGenHeterProps && filmthick0[i]>1.0e-8)
				//	filmthick[i] = filmthick0[i] + 0.10009*F_Ca / rhoRock*dt;
				//else
				//	filmthick[i] = filmthick0[i] + 0.00010009*F_Ca / rhoRock*dt;
			}
		}

		t_end = clock();
		UpdatePerm2ReconJac_time += t_end - t_start;

		return true;
	}

	/* 1. Form continuum/pseudo-continnum elastic/platic matrix; 2. Update permeability tensor */
	void SimHM::UpdatePerm2FormMtrD(Element& ele, double* ElasticMtrD){
		int iELE = ele.Ele_num;
		/*@ Initialize block perm. for intact rock*/
		AllData->IntactRockPerm(iELE);
		/*@ Initialize block D for intact rock*/
		IntactRockCompliance_F(ele, ElasticMtrD);
		//==============================================================================================//
		if (ele.FrcForm == 0); //0. intact rock
		else if (ele.FrcForm == 1){//1. only one tensile fracture block
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			calTran(ele.dipAng[0], ele.azimuthAng[0], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[0] = getSfTensile(iELE, ele.NormVct[0]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[0], ele.sf_T[0], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0; //ele.sf_T;
			//Update fracture width
			ele.wf[0] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[0], Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix	
		}
		else if (ele.FrcForm == 2){//2. secondary tensile fractured block
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			//*******************************************first tensile fracture*******************************//
			calTran(ele.dipAng[0], ele.azimuthAng[0], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[0] = getSfTensile(iELE, ele.NormVct[0]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[0], ele.sf_T[0], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[0] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[0], Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//

			//*******************************************second tensile fracture*******************************//
			calTran(ele.dipAng[1], ele.azimuthAng[1], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability:
			ele.sf_T[1] = getSfTensile(iELE, ele.NormVct[1]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[1], ele.sf_T[1], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[1] = Wf;
			//calculate elastic matrix
			double F_1[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[1], Tran, F_1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//
		}
		else if (ele.FrcForm == 3){//3. only one shear conjugate fracture set
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			calTran(ele.dipAng[2], ele.azimuthAng[2], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[2], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[3], ele.azimuthAng[3], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[3], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
		}
		else if (ele.FrcForm == 4){//4. has two shear conjugate fracture sets
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			//*******************************************first shear fracture*******************************//
			calTran(ele.dipAng[2], ele.azimuthAng[2], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability
			calPermTensor(1.0, ele.NormVct[2], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[3], ele.azimuthAng[3], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability
			calPermTensor(1.0, ele.NormVct[3], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//***********************************************************************************************//

			//*******************************************second shear fracture*******************************//
			calTran(ele.dipAng[4], ele.azimuthAng[4], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability
			calPermTensor(1.0, ele.NormVct[4], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[3] = Wf;
			//calculate elastic matrix
			double F_1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[5], ele.azimuthAng[5], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[5], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;//Sf_;
			//Update fracture width
			ele.wf[3] = Wf;
			//calculate elastic matrix
			double F1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//
		}
		else if (ele.FrcForm == 5){//5. one tensile + one shear fracture sets
			//*******************************************tensile fracture*******************************//
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			calTran(ele.dipAng[0], ele.azimuthAng[0], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[0] = getSfTensile(iELE, ele.NormVct[0]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[0], ele.sf_T[0], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[0] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[0], Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//**********************************************************************************************************//

			//*******************************************shear fracture*******************************//
			calTran(ele.dipAng[2], ele.azimuthAng[2], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[2], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F_1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[3], ele.azimuthAng[3], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[3], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
		}
		else if (ele.FrcForm == 6){//6. two tensile + one shear fracture sets
			//*******************************************1-tensile fracture*******************************//
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			calTran(ele.dipAng[0], ele.azimuthAng[0], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[0] = getSfTensile(iELE, ele.NormVct[0]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[0], ele.sf_T[0], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;//ele.sf_T;
			//Update fracture width
			ele.wf[0] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[0], Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//**********************************************************************************************************//

			//*******************************************2-tensile fracture*******************************//
			calTran(ele.dipAng[1], ele.azimuthAng[1], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[1] = getSfTensile(iELE, ele.NormVct[1]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[1], ele.sf_T[1], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[1] = Wf;
			//calculate elastic matrix
			double F_1[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[1], Tran, F_1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//********************************************************************************************//

			//*******************************************shear fracture*******************************//
			calTran(ele.dipAng[2], ele.azimuthAng[2], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[2], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F_2[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_2, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[3], ele.azimuthAng[3], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[3], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*******************************************tensile fracture*******************************//
		}
		else if (ele.FrcForm == 7){//7. one tensile + two shear fracture sets
			//*******************************************tensile fracture*******************************//
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			calTran(ele.dipAng[0], ele.azimuthAng[0], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[0] = getSfTensile(iELE, ele.NormVct[0]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[0], ele.sf_T[0], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[0] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[0], Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//**********************************************************************************************************//

			//*******************************************1-shear fracture*******************************//
			calTran(ele.dipAng[2], ele.azimuthAng[2], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[2], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F_1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[3], ele.azimuthAng[3], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[3], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//

			//************************************************2-shear fracture*********************************//
			calTran(ele.dipAng[4], ele.azimuthAng[4], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[4], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[3] = Wf;
			//calculate elastic matrix
			double F_2[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_2, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[5], ele.azimuthAng[5], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[5], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[3] = Wf;
			//calculate elastic matrix
			double F2[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F2, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//
		}
		else if (ele.FrcForm == 8){//8. two tensile + two shear fracture sets
			//*******************************************1-tensile fracture*******************************//
			double PlaneStress[3], Tran[18];
			double Wf, Kn;
			tensorVar Kij;
			calTran(ele.dipAng[0], ele.azimuthAng[0], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[0] = getSfTensile(iELE, ele.NormVct[0]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[0], ele.sf_T[0], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[0] = Wf;
			//calculate elastic matrix
			double F_[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[0], Tran, F_, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//**********************************************************************************************************//

			//*******************************************2-tensile fracture*******************************//
			calTran(ele.dipAng[1], ele.azimuthAng[1], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);//get fracture aperture & stiffness
			//Update permeability:
			ele.sf_T[1] = getSfTensile(iELE, ele.NormVct[1]); /* ele.NormVct[0]:the normal vector of least principle stress */
			calPermTensor(1.0, ele.NormVct[1], ele.sf_T[1], Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[0]) / 1.0;// ele.sf_T;
			//Update fracture width
			ele.wf[1] = Wf;
			//calculate elastic matrix
			double F_1[3] = { 1 / Kn, 1 / yita / Kn, 1 / yita / Kn };
			calTFT(ele.sf_T[1], Tran, F_1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//********************************************************************************************//

			//*******************************************1-shear fracture*******************************//
			calTran(ele.dipAng[2], ele.azimuthAng[2], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[2], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F_2[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_2, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[3], ele.azimuthAng[3], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[3], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[2] = Wf;
			//calculate elastic matrix
			double F1[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F1, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//

			//************************************************2-shear fracture*********************************//
			calTran(ele.dipAng[4], ele.azimuthAng[4], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[4], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[1]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[3] = Wf;
			//calculate elastic matrix
			double F_3[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F_3, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix

			//////////////////////////////////////////////////////////////////////////////////////////////////
			calTran(ele.dipAng[5], ele.azimuthAng[5], Tran);//get the transformation matrix
			calPlaneStress(Tran, Tstress[iELE], PlaneStress);
			getKnWf(PlaneStress[0], Kni_, vj_, vjmax_, Wf, Kn);
			//Update permeability
			calPermTensor(1.0, ele.NormVct[5], Sf_, Wf, Kij);
			AllData->UpdatePerm(iELE, Kij);//Update the permeability tensor associated with fractured blcok
			//Update porosity
			//AllData->iPORO[iELE] += (Wf - ele.wf[2]) / 1.0;// Sf_;
			//Update fracture width
			ele.wf[3] = Wf;
			//calculate elastic matrix
			double F2[3] = { 1 / Kn, yita / Kn, yita / Kn };
			calTFT(Sf_, Tran, F2, ElasticMtrD);//add the contributions of single fracture sets on compiance matrix
			//*************************************************************************************************//
		}
		else;
		InverseF6X6(ElasticMtrD);//Martix ivcerse operation: Cijkl==>Dijkl
	}
}