/**
* @file    MPFAL.cpp
* @author  Sanbai Li
* @date    2014/09
* @brief   Calculate transmissibility matrix of MPFA L-method
*
* @version
* FracTHM
*/

#include "TranMPFA.h"
#include "UNCONG_INL.h"

namespace FracTHM{

	void Vtx2Vct(Vertex Vtx, double* Vct){
		Vct[0] = Vtx.Xcood;
		Vct[1] = Vtx.Ycood;
		Vct[2] = Vtx.Zcood;
	}

	void Vct2Vtx(double* Vct, Vertex Vtx){
		Vtx.Xcood = Vct[0];
		Vtx.Ycood = Vct[1];
		Vtx.Zcood = Vct[2];
	}

	double Rijk(Vertex P, Vertex X1, Vertex X2, Vertex X3, double Tj){
		Vertex niu;
		niu = (X1 - P) ^ (X2 - P);
		double r = niu*(X3 - P) / Tj;
		return r;
	}

	Vertex NormalVctRadm4side(Vertex Y1, Vertex Y2, Vertex Y3, Vertex YS){//anticlockwise: Y3->Y1->Y2->YS->Y3
		Vertex tempVct;
		/*Note that when using overloaded operators,
		  the bracket "()" must be used to guarantee the priority */
		tempVct = 0.25*(((Y2 - Y1) ^ (Y3 - Y1)) + ((YS - Y2) ^ (Y1 - Y2)) +
			            ((Y3 - YS) ^ (Y2 - YS)) + ((Y1 - Y3) ^ (YS - Y3)));
		return tempVct;
	}

	double OMEGA_ijk(double* tK, Vertex P, Vertex X1, Vertex X2,
		                         Vertex Y1, Vertex Y2, Vertex Y3, Vertex YS, double Tj){
		double omega, niuD[3], nVctD[3], temp[3];
		Vertex nVct, niu;

		nVct = NormalVctRadm4side(Y1, Y2, Y3, YS);

		niu = (X1 - P) ^ (X2 - P);
		Vtx2Vct(nVct, nVctD);
		Vtx2Vct(niu, niuD);
		MatrA_MatrB<double>(1, 3, 3, nVctD, tK, temp);
		omega = (temp[0] * niuD[0] + temp[1] * niuD[1] + temp[2] * niuD[2]) / Tj;
		return omega;
	}

	double CalArea4(Vertex Y1, Vertex Y2, Vertex Y3, Vertex YS){
		Vertex nVct = NormalVctRadm4side(Y1, Y2, Y3, YS);
		return VctNorm(nVct);
	}

	void inverseA(double *mtxA){
		double tempM[9];
		double DetmtxA = mtxA[0] * (mtxA[4] * mtxA[8] - mtxA[5] * mtxA[7]) -
			             mtxA[1] * (mtxA[3] * mtxA[8] - mtxA[5] * mtxA[6]) +
			             mtxA[2] * (mtxA[3] * mtxA[7] - mtxA[4] * mtxA[6]);
		tempM[0] = mtxA[4] * mtxA[8] - mtxA[5] * mtxA[7];
		tempM[1] = mtxA[2] * mtxA[7] - mtxA[1] * mtxA[8];
		tempM[2] = mtxA[1] * mtxA[5] - mtxA[2] * mtxA[4];
		tempM[3] = mtxA[5] * mtxA[6] - mtxA[3] * mtxA[8];
		tempM[4] = mtxA[0] * mtxA[8] - mtxA[2] * mtxA[6];
		tempM[5] = mtxA[2] * mtxA[3] - mtxA[0] * mtxA[5];
		tempM[6] = mtxA[3] * mtxA[7] - mtxA[4] * mtxA[6];
		tempM[7] = mtxA[1] * mtxA[6] - mtxA[0] * mtxA[7];
		tempM[8] = mtxA[0] * mtxA[4] - mtxA[1] * mtxA[3];

		for (int i = 0; i < 9; ++i) mtxA[i] = tempM[i] / DetmtxA;
	}

	void SequenceMap(int flag, const int _fclnSeq[7], const int dir, const int steps,
		             int fclnSeq[7], int IJKplus[9])
	{
		//IJKplus[9] contain the correct item
		
		/*fclnSeq[7]={facec1, facec2, facec3, linec1, ..., linec3 (, linec4)}
		****************************************************************************************
		**  For type A fclnSeq[7]={ Block1.face1, Block1.face2, Block1.face3,                 **
		**                          Block1.line1, Block1.line2, Block1.line3}                 **
		**  For type B fclnSeq[7]={ Block2.face1, Block2.face2, Block2.face3,                 **
		**                          Block2.line1, Block2.line2, Block2.line3}				  **
		**  For type C fclnSeq[7]={ Block1.face1, Block2.face2, Block1.face3,                 **
		**                          Block1.line1, Block1.line2, Block1.line3, Block2.line4}   **
		**  For type D fclnSeq[7]={ Block1.face1, Block1.face2, Block2.face3,                 **
		**                          Block1.line1, Block1.line2, Block2.line3, Block1.line4}   **
		****************************************************************************************
		*/
		/* See present ELEMENT (usually number 1) as the reference (I,J,K) index
		   then enumerate ELEMENTS, shearing its 6 faces, IJK increments as follwos */
		const int IJKplus_[6][3] =   { { 0, 0, 1 }, { 0, 0, -1 }, { 0, -1, 0},
		                               { 0, 1, 0 }, { -1, 0, 0 }, { 1,  0, 0} };
		
		const int FacecSeq[3][10] = {{ 4, 5,
			                           0, 3, 1, 2, 0, 3, 1, 2 },
		                             { 2, 3,
			                           0, 5, 1, 4, 0, 5, 1, 4 },
                                     { 0, 1,
			                           2, 5, 3, 4, 2, 5, 3, 4 }};	

		const int linecSeq[3][3][8] =
		{
			{ { 0, 1,   5,  4, 0,  1,  5,  4 },
		      { 3, 11,  7,  9, 3, 11,  7,  9 },
		      { 2, 10,  6,  8, 2, 10,  6,  8 } },
            { { 2,  3,  7,  6, 2,  3,  7,  6}, 
		      { 0,  9,  4,  8, 0,  9,  4,  8}, 
		      { 1, 11,  5, 10, 1, 11,  5, 10} },
            { { 9, 11, 10,  8, 9, 11, 10,  8}, 
		      { 4,  7,  5,  6, 4,  7,  5,  6},
		      { 0,  3,  1,  2, 0,  3,  1,  2} }
		};

		int count_ = 0;
		for (int i = 0; i < 6 + int(flag>2); ++i)
		{
			if (i < 3){    //deal with facecs
				for (int j = 0; j < 6; ++j)
				{
					if (_fclnSeq[i] == FacecSeq[dir - 1][j]) {
						if (j < 2)
							fclnSeq[i] = FacecSeq[dir - 1][j];
						else
							fclnSeq[i] = FacecSeq[dir - 1][j + steps];

						//Generate IJK index for relevant 4 cells
						for (int l = 0; l < 3; ++l)
							IJKplus[count_ + l] = IJKplus_[fclnSeq[i]][l];
						count_ += 3;
						break;  //if find then jump out the loop
					}else;
				}
			}
			else{         //deal with linecs
				bool isfind=false;
				for (int j = 0; j < 3; ++j)
				{
					for (int k = 0; k < 4; ++k){
						if (_fclnSeq[i] == linecSeq[dir - 1][j][k]) {
							fclnSeq[i] = linecSeq[dir - 1][j][k + steps];
							isfind = true;
							break;
						}
					}
					if (isfind)break;
					else;
				}
			}
		}

		/*
		  The IJK index increment refer to block(I,J,K), usually cell 1, however some face center, such as of L-stencil B, C,and D,
		  not belong to cell 1, but to cell 2, thus we may add an additional increment of cell 2 referring to cell 1.
		*/
		if (flag == 1);
		else if (flag == 2){
			for (int i = 0; i < 3; ++i)
			{
				IJKplus[3 + i] += IJKplus[i];
				IJKplus[6 + i] += IJKplus[i];
			}
		}
		else if (flag == 3){
			for (int i = 0; i < 3; ++i)
			{
				IJKplus[3 + i] += IJKplus[i];
			}
		}
		else if (flag == 4){
			for (int i = 0; i < 3; ++i)
			{
				IJKplus[6 + i] += IJKplus[i];
			}
		}
	}

	//**************************************  L-stencil Type  A    ***********************************
	double LtypeA(const int flSeq[7], 
		          const CVFEMcell& B1, const CVFEMcell& B2, 
				  const CVFEMcell& B3, const CVFEMcell& B5, const Vertex& SP,
				  double* tK1, double* tK2, double* tK3, double* tK5, double* Ta){

		Vertex facec1 = B1.facec[flSeq[0]];
		Vertex facec2 = B1.facec[flSeq[1]];
		Vertex facec3 = B1.facec[flSeq[2]];
		Vertex linec4 = B1.linec[flSeq[3]];
		Vertex linec5 = B1.linec[flSeq[4]];
		Vertex linec6 = B1.linec[flSeq[5]];
		
		double mtxA[9], mtxB[12], mtxC[9], mtxD[3];
	
		double Tetr1, Tetr2, Tetr3, Tetr5;

		Tetr1 = TetraVol(B1.blockc, facec1, facec2, facec3);
		Tetr2 = TetraVol(B2.blockc, facec1, linec4, linec5);
		Tetr3 = TetraVol(B3.blockc, facec2, linec4, linec6);
		Tetr5 = TetraVol(B5.blockc, facec3, linec5, linec6);

		//Calculate omega_ijk parameters
		double w111, w112, w113, w121, w122, w123, w211, w212, w213, 
			   w231, w232, w233, w311, w312, w313, w351, w352, w353;

		w111 = OMEGA_ijk(tK1, B1.blockc, facec2, facec3, facec1, linec4, linec5, SP, 6 * Tetr1);
		w112 = OMEGA_ijk(tK1, B1.blockc, facec3, facec1, facec1, linec4, linec5, SP, 6 * Tetr1);
		w113 = OMEGA_ijk(tK1, B1.blockc, facec1, facec2, facec1, linec4, linec5, SP, 6 * Tetr1);

		w211 = OMEGA_ijk(tK1, B1.blockc, facec2, facec3, facec2, linec6, linec4, SP, 6 * Tetr1);
		w212 = OMEGA_ijk(tK1, B1.blockc, facec3, facec1, facec2, linec6, linec4, SP, 6 * Tetr1);
		w213 = OMEGA_ijk(tK1, B1.blockc, facec1, facec2, facec2, linec6, linec4, SP, 6 * Tetr1);

		w311 = OMEGA_ijk(tK1, B1.blockc, facec2, facec3, facec3, linec5, linec6, SP, 6 * Tetr1);
		w312 = OMEGA_ijk(tK1, B1.blockc, facec3, facec1, facec3, linec5, linec6, SP, 6 * Tetr1);
		w313 = OMEGA_ijk(tK1, B1.blockc, facec1, facec2, facec3, linec5, linec6, SP, 6 * Tetr1);

		w121 = OMEGA_ijk(tK2, B2.blockc, facec1, linec5, facec1, linec4, linec5, SP, 6 * Tetr2);
		w122 = OMEGA_ijk(tK2, B2.blockc, linec5, linec4, facec1, linec4, linec5, SP, 6 * Tetr2);
		w123 = OMEGA_ijk(tK2, B2.blockc, linec4, facec1, facec1, linec4, linec5, SP, 6 * Tetr2);

		w231 = OMEGA_ijk(tK3, B3.blockc, linec4, linec6, facec2, linec6, linec4, SP, 6 * Tetr3);
		w232 = OMEGA_ijk(tK3, B3.blockc, linec6, facec2, facec2, linec6, linec4, SP, 6 * Tetr3);
		w233 = OMEGA_ijk(tK3, B3.blockc, facec2, linec4, facec2, linec6, linec4, SP, 6 * Tetr3);

		w351 = OMEGA_ijk(tK5, B5.blockc, linec5, facec3, facec3, linec5, linec6, SP, 6 * Tetr5);
		w352 = OMEGA_ijk(tK5, B5.blockc, facec3, linec6, facec3, linec5, linec6, SP, 6 * Tetr5);
		w353 = OMEGA_ijk(tK5, B5.blockc, linec6, linec5, facec3, linec5, linec6, SP, 6 * Tetr5);

		//Calculate r_ijk parameters
		double r111, r121, r131, r112, r122, r132, r113, r123, r133;
		r111 = Rijk(B1.blockc, facec2, facec3, linec4, 6 * Tetr1);
		r112 = Rijk(B1.blockc, facec2, facec3, linec5, 6 * Tetr1);
		r113 = Rijk(B1.blockc, facec2, facec3, linec6, 6 * Tetr1);
		r121 = Rijk(B1.blockc, facec3, facec1, linec4, 6 * Tetr1);
		r122 = Rijk(B1.blockc, facec3, facec1, linec5, 6 * Tetr1);
		r123 = Rijk(B1.blockc, facec3, facec1, linec6, 6 * Tetr1);
		r131 = Rijk(B1.blockc, facec1, facec2, linec4, 6 * Tetr1);
		r132 = Rijk(B1.blockc, facec1, facec2, linec5, 6 * Tetr1);
		r133 = Rijk(B1.blockc, facec1, facec2, linec6, 6 * Tetr1);

		//Calculate matrix A, B, C, and D
		mtxA[0] = w111 - w122      - w121*r111  - w123*r112;
		mtxA[1] = w112 - w121*r121 - w123*r122;
		mtxA[2] = w113 - w121*r131 - w123*r132;
		mtxA[3] = w211 - w232*r111 - w233*r113;
		mtxA[4] = w212 - w231      - w232*r121 - w233*r123;
		mtxA[5] = w213 - w232*r131 - w233*r133;
		mtxA[6] = w311 - w351*r113 - w352*r112;
		mtxA[7] = w312 - w351*r123 - w352*r122;
		mtxA[8] = w313 - w353      - w351*r133 - w352*r132;

		mtxB[0 ] = w111 + w112 + w113 + w121*(1 - r111 - r121 - r131) + w123*(1 - r112 - r122 - r132);
		mtxB[1] = -w121 - w122 - w123;
		mtxB[2 ] = 0.0;
		mtxB[3 ] = 0.0;
		mtxB[4 ] = w211 + w212 + w213 + w232*(1 - r111 - r121 - r131) + w233*(1 - r113 - r123 - r133);
		mtxB[5 ] = 0.0;
		mtxB[6] = -w231 - w232 - w233;
		mtxB[7 ] = 0.0;
		mtxB[8 ] = w311 + w312 + w313 + w351*(1 - r113 - r123 - r133) + w352*(1 - r112 - r122 - r132);
		mtxB[9 ] = 0.0;
		mtxB[10] = 0.0;
		mtxB[11] = -w351 - w352 - w353;

		mtxC[0] = -w111;
		mtxC[1] = -w112;
		mtxC[2] = -w113;
		mtxC[3] = -w211;
		mtxC[4] = -w212;
		mtxC[5] = -w213;
		mtxC[6] = -w311;
		mtxC[7] = -w312;
		mtxC[8] = -w313;

		mtxD[0] = w111 + w112 + w113;//mtxD[0]
		mtxD[1] = w211 + w212 + w213;//mtxD[4]
		mtxD[2] = w311 + w312 + w313;//mtxD[8]
		//for (int i = 0; i < 12; ++i){if (i % 4 != 0)mtxD[i] = 0.0;}
		    
		//Calculate the transmissibility matrix Ta for L-stencil(a)
		double CA[9], CA_B[12];
		SVD_Inv3(mtxA);                   //inverseA(mtxA);
		FracTHMSOL::MM3(CA, mtxC, mtxA);   //MatrA_MatrB(3, 3, 3, mtxC, mtxA, CA  );
		FracTHMSOL::MM334(CA_B, CA, mtxB); // MatrA_MatrB(3, 3, 4, CA,   mtxB, CA_B);
		//plus_Mtr(3, 4, CA_B, mtxD);
		CA_B[0] += mtxD[0];
		CA_B[4] += mtxD[1];
		CA_B[8] += mtxD[2];
		copy_Mtr(3, 4, CA_B, Ta);
		return CalArea4(SP, linec4, linec5, facec1);
	}
	//************************************************************************************************

	//**************************************  L-stencil Type  B    ***********************************
	double LtypeB(const int flSeq[7],
		          const CVFEMcell& B1, const CVFEMcell& B2, 
				  const CVFEMcell& B4, const CVFEMcell& B6, const Vertex& SP,
				  double* tK1, double* tK2, double* tK4, double* tK6, double* Tb){

		Vertex facec1 = B1.facec[flSeq[0]];
		Vertex facec2 = B2.facec[flSeq[1]];
		Vertex facec3 = B2.facec[flSeq[2]];
		Vertex linec4 = B1.linec[flSeq[3]];
		Vertex linec5 = B1.linec[flSeq[4]];
		Vertex linec6 = B2.linec[flSeq[5]];
		
		double mtxA[9], mtxB[12], mtxC[9], mtxD[3];

		double Tetr1, Tetr2, Tetr4, Tetr6;

		Tetr1 = TetraVol(B1.blockc, facec1, linec4, linec5);
		Tetr2 = TetraVol(B2.blockc, facec1, facec2, facec3);
		Tetr4 = TetraVol(B4.blockc, linec4, linec6, facec2);
		Tetr6 = TetraVol(B6.blockc, linec6, facec3, linec5);

		//Calculate omega_ijk parameters
		double w111, w112, w113, w121, w122, w123, w221, w222, w223,
			   w241, w242, w243, w321, w322, w323, w361, w362, w363;

		w111 = OMEGA_ijk(tK1, B1.blockc, linec5, facec1, facec1, linec4, linec5, SP, 6 * Tetr1);
		w112 = OMEGA_ijk(tK1, B1.blockc, linec4, linec5, facec1, linec4, linec5, SP, 6 * Tetr1);
		w113 = OMEGA_ijk(tK1, B1.blockc, facec1, linec4, facec1, linec4, linec5, SP, 6 * Tetr1);

		w121 = OMEGA_ijk(tK2, B2.blockc, facec3, facec2, facec1, linec4, linec5, SP, 6 * Tetr2);
		w122 = OMEGA_ijk(tK2, B2.blockc, facec1, facec3, facec1, linec4, linec5, SP, 6 * Tetr2);
		w123 = OMEGA_ijk(tK2, B2.blockc, facec2, facec1, facec1, linec4, linec5, SP, 6 * Tetr2);

		w221 = OMEGA_ijk(tK2, B2.blockc, facec3, facec2, facec2, linec4, linec6, SP, 6 * Tetr2);
		w222 = OMEGA_ijk(tK2, B2.blockc, facec1, facec3, facec2, linec4, linec6, SP, 6 * Tetr2);
		w223 = OMEGA_ijk(tK2, B2.blockc, facec2, facec1, facec2, linec4, linec6, SP, 6 * Tetr2);

		w241 = OMEGA_ijk(tK4, B4.blockc, linec6, linec4, facec2, linec4, linec6, SP, 6 * Tetr4);
		w242 = OMEGA_ijk(tK4, B4.blockc, facec2, linec6, facec2, linec4, linec6, SP, 6 * Tetr4);
		w243 = OMEGA_ijk(tK4, B4.blockc, linec4, facec2, facec2, linec4, linec6, SP, 6 * Tetr4);

		w321 = OMEGA_ijk(tK2, B2.blockc, facec3, facec2, facec3, linec6, linec5, SP, 6 * Tetr2);
		w322 = OMEGA_ijk(tK2, B2.blockc, facec1, facec3, facec3, linec6, linec5, SP, 6 * Tetr2);
		w323 = OMEGA_ijk(tK2, B2.blockc, facec2, facec1, facec3, linec6, linec5, SP, 6 * Tetr2);

		w361 = OMEGA_ijk(tK6, B6.blockc, facec3, linec5, facec3, linec6, linec5, SP, 6 * Tetr6);
		w362 = OMEGA_ijk(tK6, B6.blockc, linec6, facec3, facec3, linec6, linec5, SP, 6 * Tetr6);
		w363 = OMEGA_ijk(tK6, B6.blockc, linec5, linec6, facec3, linec6, linec5, SP, 6 * Tetr6);

		//Calculate r_ijk parameters
		double r211, r221, r231, r212, r222, r232, r213, r223, r233;
		r211 = Rijk(B2.blockc, facec3, facec2, linec4, 6 * Tetr2);
		r212 = Rijk(B2.blockc, facec3, facec2, linec5, 6 * Tetr2);
		r213 = Rijk(B2.blockc, facec3, facec2, linec6, 6 * Tetr2);
		r221 = Rijk(B2.blockc, facec1, facec3, linec4, 6 * Tetr2);
		r222 = Rijk(B2.blockc, facec1, facec3, linec5, 6 * Tetr2);
		r223 = Rijk(B2.blockc, facec1, facec3, linec6, 6 * Tetr2);
		r231 = Rijk(B2.blockc, facec2, facec1, linec4, 6 * Tetr2);
		r232 = Rijk(B2.blockc, facec2, facec1, linec5, 6 * Tetr2);
		r233 = Rijk(B2.blockc, facec2, facec1, linec6, 6 * Tetr2);

		//Calculate matrix A, B, C, and D
		mtxA[0] = -w121 + w112      + w111*r211 + w113*r212;
		mtxA[1] = -w122 + w111*r221 + w113*r222;
		mtxA[2] = -w123 + w111*r231 + w113*r232;
		mtxA[3] =  w221 - w242*r211 - w243*r213;
		mtxA[4] =  w222 - w241      - w242*r221 - w243*r223;
		mtxA[5] =  w223 - w242*r231 - w243*r233;
		mtxA[6] =  w321 - w361*r213 - w362*r212;
		mtxA[7] =  w322 - w361*r223 - w362*r222;
		mtxA[8] =  w323 - w363      - w361*r233 - w362*r232;

		mtxB[0] =  w111 + w112 + w113;
		mtxB[1] = -w121 - w122 - w123 - w111*(1 - r211 - r221 - r231)
			                          - w113*(1 - r212 - r222 - r232);
		mtxB[2] = 0.0;
		mtxB[3] = 0.0;
		mtxB[4] = 0.0;
		mtxB[5] = w221 + w222 + w223 + w242*(1 - r211 - r221 - r231)
			                         + w243*(1 - r213 - r223 - r233);
		mtxB[6] = -w241 - w242 - w243;
		mtxB[7] = 0.0;
		mtxB[8] = 0.0;
		mtxB[9] = w321 + w322 + w323 + w361*(1 - r213 - r223 - r233)
			                         + w362*(1 - r212 - r222 - r232);
		mtxB[10] = 0.0;
		mtxB[11] = -w361 - w362 - w363;

		mtxC[0] = -w121;
		mtxC[1] = -w122;
		mtxC[2] = -w123;
		mtxC[3] = -w221;
		mtxC[4] = -w222;
		mtxC[5] = -w223;
		mtxC[6] = -w321;
		mtxC[7] = -w322;
		mtxC[8] = -w323;

		mtxD[0] = w121 + w122 + w123;//mtxD[1]
		mtxD[1] = w221 + w222 + w223;//mtxD[5]
		mtxD[2] = w321 + w322 + w323;//mtxD[9]
		/*for (int i = 0; i < 12; ++i){
			if (i % 4 != 1)mtxD[i] = 0.0;
		}*/

		//Calculate the transmissibility matrix Tb for L-stencil(b)
		double CA[9], CA_B[12];
		SVD_Inv3(mtxA);                   //inverseA(mtxA);
		FracTHMSOL::MM3(CA, mtxC, mtxA);   //MatrA_MatrB(3, 3, 3, mtxC, mtxA, CA  );
		FracTHMSOL::MM334(CA_B, CA, mtxB); // MatrA_MatrB(3, 3, 4, CA,   mtxB, CA_B);
		//plus_Mtr(3, 4, CA_B, mtxD);
		CA_B[1] += mtxD[0];
		CA_B[5] += mtxD[1];
		CA_B[9] += mtxD[2];
		copy_Mtr(3, 4, CA_B, Tb);
		return CalArea4(SP, linec4, linec5, facec1);
	}
	//************************************************************************************************

	//**************************************  L-stencil Type  C    ***********************************
	double LtypeC(const int flSeq[7],
		          const CVFEMcell& B1, const CVFEMcell& B2, 
				  const CVFEMcell& B4, const CVFEMcell& B5, const Vertex& SP,
				  double* tK1, double* tK2, double* tK4, double* tK5, double* Tc){

		Vertex facec1 = B1.facec[flSeq[0]];
		Vertex facec2 = B2.facec[flSeq[1]];
		Vertex facec3 = B1.facec[flSeq[2]];
		Vertex linec4 = B1.linec[flSeq[3]];
		Vertex linec5 = B1.linec[flSeq[4]];
		Vertex linec6 = B1.linec[flSeq[5]];
		Vertex linec7 = B2.linec[flSeq[6]];
		
		double mtxA[9], mtxB[12], mtxC[9], mtxD[12];

		double Tetr1, Tetr2, Tetr4, Tetr5;

		Tetr1 = TetraVol(B1.blockc, facec1, facec3, linec4);
		Tetr2 = TetraVol(B2.blockc, facec1, facec2, linec5);
		Tetr4 = TetraVol(B4.blockc, facec2, linec7, linec4);
		Tetr5 = TetraVol(B5.blockc, facec3, linec5, linec6);

		//Calculate omega_ijk parameters
		double w111, w112, w113, w121, w122, w123, w221, w222, w223,
			   w241, w242, w243, w311, w312, w313, w351, w352, w353;

		w111 = OMEGA_ijk(tK1, B1.blockc, linec4, facec3, facec1, linec4, linec5, SP, 6 * Tetr1);
		w112 = OMEGA_ijk(tK1, B1.blockc, facec1, linec4, facec1, linec4, linec5, SP, 6 * Tetr1);
		w113 = OMEGA_ijk(tK1, B1.blockc, facec3, facec1, facec1, linec4, linec5, SP, 6 * Tetr1);

		w121 = OMEGA_ijk(tK2, B2.blockc, linec5, facec2, facec1, linec4, linec5, SP, 6 * Tetr2);
		w122 = OMEGA_ijk(tK2, B2.blockc, facec1, linec5, facec1, linec4, linec5, SP, 6 * Tetr2);
		w123 = OMEGA_ijk(tK2, B2.blockc, facec2, facec1, facec1, linec4, linec5, SP, 6 * Tetr2);

		w221 = OMEGA_ijk(tK2, B2.blockc, linec5, facec2, facec2, linec4, linec7, SP, 6 * Tetr2);
		w222 = OMEGA_ijk(tK2, B2.blockc, facec1, linec5, facec2, linec4, linec7, SP, 6 * Tetr2);
		w223 = OMEGA_ijk(tK2, B2.blockc, facec2, facec1, facec2, linec4, linec7, SP, 6 * Tetr2);

		w241 = OMEGA_ijk(tK4, B4.blockc, linec7, linec4, facec2, linec4, linec7, SP, 6 * Tetr4);
		w242 = OMEGA_ijk(tK4, B4.blockc, facec2, linec7, facec2, linec4, linec7, SP, 6 * Tetr4);
		w243 = OMEGA_ijk(tK4, B4.blockc, linec4, facec2, facec2, linec4, linec7, SP, 6 * Tetr4);

		w311 = OMEGA_ijk(tK1, B1.blockc, linec4, facec3, facec3, linec5, linec6, SP, 6 * Tetr1);
		w312 = OMEGA_ijk(tK1, B1.blockc, facec1, linec4, facec3, linec5, linec6, SP, 6 * Tetr1);
		w313 = OMEGA_ijk(tK1, B1.blockc, facec3, facec1, facec3, linec5, linec6, SP, 6 * Tetr1);

		w351 = OMEGA_ijk(tK5, B5.blockc, linec6, linec5, facec3, linec5, linec6, SP, 6 * Tetr5);
		w352 = OMEGA_ijk(tK5, B5.blockc, facec3, linec6, facec3, linec5, linec6, SP, 6 * Tetr5);
		w353 = OMEGA_ijk(tK5, B5.blockc, linec5, facec3, facec3, linec5, linec6, SP, 6 * Tetr5);

		//Calculate r_ijk parameters
		double r112, r122, r132, r113, r123, r133, 
			   r211, r221, r231, r214, r224, r234;
		r112 = Rijk(B1.blockc, linec4, facec3, linec5, 6 * Tetr1);
		r122 = Rijk(B1.blockc, facec1, linec4, linec5, 6 * Tetr1);
		r132 = Rijk(B1.blockc, facec3, facec1, linec5, 6 * Tetr1);
		r113 = Rijk(B1.blockc, linec4, facec3, linec6, 6 * Tetr1);
		r123 = Rijk(B1.blockc, facec1, linec4, linec6, 6 * Tetr1);
		r133 = Rijk(B1.blockc, facec3, facec1, linec6, 6 * Tetr1);
		r211 = Rijk(B2.blockc, linec5, facec2, linec4, 6 * Tetr2);
		r221 = Rijk(B2.blockc, facec1, linec5, linec4, 6 * Tetr2);
		r231 = Rijk(B2.blockc, facec2, facec1, linec4, 6 * Tetr2);
		r214 = Rijk(B2.blockc, linec5, facec2, linec7, 6 * Tetr2);
		r224 = Rijk(B2.blockc, facec1, linec5, linec7, 6 * Tetr2);
		r234 = Rijk(B2.blockc, facec2, facec1, linec7, 6 * Tetr2);

		//Calculate matrix A, B, C, and D
		double cst = 1 / (1 - r231*r132);
		mtxA[0] = w111 - w121 + cst*w113*(r231*r112 + r211) - cst*w123*(r112+r132*r211);
		mtxA[1] = cst*w113*r221 - w122 - cst*w123*r132*r221;
		mtxA[2] = w112 + cst*w113*r231*r122 - cst*w123*r122;
		mtxA[3] = w221 - w243*r214 + cst*(w223 - w243*r234)*(r112 + r132*r211)
			      - cst*w242*(r231*r112 + r211);
		mtxA[4] = w222 - w243*r224 - w241 + cst*(w223 - w243*r234)*r132*r221-cst*w242*r221;
		mtxA[5] = cst*r122*(w223 - w243*r234 - w242*r231);
		mtxA[6] = w311 - w353*r113 + cst*(w313 - w353*r133)*(r231*r112 + r211)
			      - cst*w352*(r112 + r132*r211);
		mtxA[7] = cst*r221*(w313 - w353*r133 - w352*r132);
		mtxA[8] = w312 - w351 - w353*r123 + cst*(w313 - w353*r133)*r231*r122 - cst*w352*r122;

		mtxB[0] = cst*(w113*r231 - w123)*(r112 + r122 + r132 - 1) + w111 + w112 + w113;
		mtxB[1] = cst*(w113 - w123*r132)*(r211 + r221 + r231 - 1) - w121 - w122 - w123;
		mtxB[2] = 0.0;
		mtxB[3] = 0.0;
		mtxB[4] = cst*(r112 + r122 + r132 - 1)*(w223 - w243*r234 - w242*r231);
		mtxB[5] = w221 + w222 + w223 - w243*(r214 + r224 + r234 - 1) 
			      + cst*(r211 + r221 + r231 - 1)*(w223*r132 - w243*r234*r132 - w242);
		mtxB[6] = -w241 - w242 - w243;
		mtxB[7] = 0.0;
		mtxB[8] = cst*(r112 + r122 + r132 - 1)*(w313*r231 - w353*r133*r231 - w352)
			      + w311 + w312 + w313 - w353*(r113 + r123 + r133 - 1);
		mtxB[9] = cst*(r211 + r221 + r231 - 1)*(w313 - w353*r133 - w352*r132);
		mtxB[10] = 0.0;
		mtxB[11] = -w351 - w352 - w353;

		mtxC[0] = -w111 - cst*w113*r231*r112 - cst*w113*r211;
		mtxC[1] = -cst*w113*r221;
		mtxC[2] = -w112 - cst*w113*r231*r122;
		mtxC[3] = -w221 - cst*w223*r112 - cst*w223*r132*r211;
		mtxC[4] = -w222 - cst*w223*r132*r221;
		mtxC[5] = -cst*w223*r122;
		mtxC[6] = -w311 - cst*w313*r231*r112 - cst*w313*r211;
		mtxC[7] = -cst*w313*r221;
		mtxC[8] = -w312 - cst*w313*r231*r122;

		mtxD[0] = cst*w113*r231*(r112 + r122 + r132 - 1) + w111 + w112 + w113;
		mtxD[1] = cst*w113*(r211 + r221 + r231 - 1);
		mtxD[4] = cst*w223*(r112 + r122 + r132 - 1);
		mtxD[5] = w221 + w222 + w223 + cst*w223*r132*(r211 + r221 + r231 - 1);
		mtxD[8] = cst*w313*r231*(r112 + r122 + r132 - 1) + w311 + w312 + w313;
		mtxD[9] = cst*w313*(r211 + r221 + r231 - 1);
		for (int i = 0; i < 12; ++i){if (i % 4 == 2 || i % 4 == 3)mtxD[i] = 0.0;}

		//Calculate the transmissibility matrix Tc for L-stencil(c)
		double CA[9], CA_B[12];
		SVD_Inv3(mtxA);                   //inverseA(mtxA);
		FracTHMSOL::MM3(CA, mtxC, mtxA);   //MatrA_MatrB(3, 3, 3, mtxC, mtxA, CA  );
		FracTHMSOL::MM334(CA_B, CA, mtxB); // MatrA_MatrB(3, 3, 4, CA,   mtxB, CA_B);
		plus_Mtr(3, 4, CA_B, mtxD);
		copy_Mtr(3, 4, mtxD, Tc);
		return CalArea4(SP, linec4, linec5, facec1);
	}
	//************************************************************************************************

	//**************************************  L-stencil Type  D    ***********************************
	double LtypeD(const int flSeq[7],
		          const CVFEMcell& B1, const CVFEMcell& B2, 
				  const CVFEMcell& B3, const CVFEMcell& B6, const Vertex& SP,
				  double* tK1, double* tK2, double* tK3, double* tK6, double* Td){

		Vertex facec1 = B1.facec[flSeq[0]];
		Vertex facec2 = B1.facec[flSeq[1]];
		Vertex facec3 = B2.facec[flSeq[2]];
		Vertex linec4 = B1.linec[flSeq[3]];
		Vertex linec5 = B1.linec[flSeq[4]];
		Vertex linec6 = B2.linec[flSeq[5]];
		Vertex linec7 = B1.linec[flSeq[6]];
		
		double mtxA[9], mtxB[12], mtxC[9], mtxD[12];

		double Tetr1, Tetr2, Tetr3, Tetr6;

		Tetr1 = TetraVol(B1.blockc, facec1, facec2, linec5);
		Tetr2 = TetraVol(B2.blockc, facec1, facec3, linec4);
		Tetr3 = TetraVol(B3.blockc, facec2, linec4, linec7);
		Tetr6 = TetraVol(B6.blockc, facec3, linec5, linec6);

		//Calculate omega_ijk parameters
		double w111, w112, w113, w121, w122, w123, w211, w212, w213,
			   w231, w232, w233, w321, w322, w323, w361, w362, w363;

		w111 = OMEGA_ijk(tK1, B1.blockc, facec2, linec5, facec1, linec4, linec5, SP, 6 * Tetr1);
		w112 = OMEGA_ijk(tK1, B1.blockc, linec5, facec1, facec1, linec4, linec5, SP, 6 * Tetr1);
		w113 = OMEGA_ijk(tK1, B1.blockc, facec1, facec2, facec1, linec4, linec5, SP, 6 * Tetr1);

		w121 = OMEGA_ijk(tK2, B2.blockc, facec3, linec4, facec1, linec4, linec5, SP, 6 * Tetr2);
		w122 = OMEGA_ijk(tK2, B2.blockc, linec4, facec1, facec1, linec4, linec5, SP, 6 * Tetr2);
		w123 = OMEGA_ijk(tK2, B2.blockc, facec1, facec3, facec1, linec4, linec5, SP, 6 * Tetr2);

		w211 = OMEGA_ijk(tK1, B1.blockc, facec2, linec5, facec2, linec7, linec4, SP, 6 * Tetr1);
		w212 = OMEGA_ijk(tK1, B1.blockc, linec5, facec1, facec2, linec7, linec4, SP, 6 * Tetr1);
		w213 = OMEGA_ijk(tK1, B1.blockc, facec1, facec2, facec2, linec7, linec4, SP, 6 * Tetr1);
															
		w231 = OMEGA_ijk(tK3, B3.blockc, linec4, linec7, facec2, linec7, linec4, SP, 6 * Tetr3);
		w232 = OMEGA_ijk(tK3, B3.blockc, linec7, facec2, facec2, linec7, linec4, SP, 6 * Tetr3);
		w233 = OMEGA_ijk(tK3, B3.blockc, facec2, linec4, facec2, linec7, linec4, SP, 6 * Tetr3);

		w361 = OMEGA_ijk(tK6, B6.blockc, linec5, linec6, facec3, linec6, linec5, SP, 6 * Tetr6);
		w362 = OMEGA_ijk(tK6, B6.blockc, linec6, facec3, facec3, linec6, linec5, SP, 6 * Tetr6);
		w363 = OMEGA_ijk(tK6, B6.blockc, facec3, linec5, facec3, linec6, linec5, SP, 6 * Tetr6);

		w321 = OMEGA_ijk(tK2, B2.blockc, facec3, linec4, facec3, linec6, linec5, SP, 6 * Tetr2);
		w322 = OMEGA_ijk(tK2, B2.blockc, linec4, facec1, facec3, linec6, linec5, SP, 6 * Tetr2);
		w323 = OMEGA_ijk(tK2, B2.blockc, facec1, facec3, facec3, linec6, linec5, SP, 6 * Tetr2);

		//Calculate r_ijk parameters
		double r111, r121, r131, r114, r124, r134,
			   r212, r222, r232, r213, r223, r233;
		r111 = Rijk(B1.blockc, facec2, linec5, linec4, 6 * Tetr1);
		r121 = Rijk(B1.blockc, linec5, facec1, linec4, 6 * Tetr1);
		r131 = Rijk(B1.blockc, facec1, facec2, linec4, 6 * Tetr1);
		r114 = Rijk(B1.blockc, facec2, linec5, linec7, 6 * Tetr1);
		r124 = Rijk(B1.blockc, linec5, facec1, linec7, 6 * Tetr1);
		r134 = Rijk(B1.blockc, facec1, facec2, linec7, 6 * Tetr1);
		r212 = Rijk(B2.blockc, facec3, linec4, linec5, 6 * Tetr2);
		r222 = Rijk(B2.blockc, linec4, facec1, linec5, 6 * Tetr2);
		r232 = Rijk(B2.blockc, facec1, facec3, linec5, 6 * Tetr2);
		r213 = Rijk(B2.blockc, facec3, linec4, linec6, 6 * Tetr2);
		r223 = Rijk(B2.blockc, linec4, facec1, linec6, 6 * Tetr2);
		r233 = Rijk(B2.blockc, facec1, facec3, linec6, 6 * Tetr2);

		//Calculate matrix A, B, C, and D
		double cst = 1 / (1 - r232*r131);
		mtxA[0] = w111 - w121 + cst*w113*(r232*r111 + r212) - cst*w123*(r111 + r131*r212);
		mtxA[1] = cst*w113*r232*r121 + w112 - cst*w123*r121;
		mtxA[2] = -w122 - cst*w123*r131*r222 + cst*w113*r222;
		mtxA[3] = w211 - w233*r114 + cst*(w213 - w233*r134)*(r212 + r232*r111)
			       - cst*w232*(r111+r131*r212);
		mtxA[4] = w212 - w233*r124 - w231 + cst*(w213 - w233*r134)*r232*r121 - cst*w232*r121;
		mtxA[5] = cst*r222*(w213 - w233*r134 - w232*r131);
		mtxA[6] = w321 - w363*r213 + cst*(w323 - w363*r233)*(r111 + r212*r131)
			      - cst*w362*(r212 + r232*r111);
		mtxA[7] = cst*r121*(w323 - w363*r233 - w362*r232);
		mtxA[8] = w322 - w361 - w363*r223 + cst*(w323 - w363*r233)*r131*r222 - cst*w362*r222;

		mtxB[0] = cst*(w113*r232 - w123)*(r111 + r121 + r131 - 1) + w111 + w112 + w113;
		mtxB[1] = cst*(w113 - w123*r131)*(r212 + r222 + r232 - 1) - w121 - w122 - w123;
		mtxB[2] = 0.0;
		mtxB[3] = 0.0;
		mtxB[4] = w211 + w212 + w213 - w233*(r114 + r124 + r134 - 1) 
			      + cst*(r111 + r121 + r131 - 1)*(w213*r232 - w233*r134*r232 - w232);
		mtxB[5] = cst*(r212 + r222 + r232 - 1)*(w213 - w233*r134 - w232*r131);
		mtxB[6] = -w231 - w232 - w233;
		mtxB[7] = 0.0;
		mtxB[8] = cst*(r111 + r121 + r131 - 1)*(w323 - w363*r233 - w362*r232);
		mtxB[9] = cst*(r212 + r222 + r232 - 1)*(w323*r131 - w363*r233*r131 - w362)
			      + w321 + w322 + w323 - w363*(r213 + r223 + r233 - 1);
		mtxB[10] = 0.0;
		mtxB[11] = -w361 - w362 - w363;

		mtxC[0] = -w111 - cst*w113*r232*r111 - cst*w113*r212;
		mtxC[1] = -w112 - cst*w113*r232*r121;
		mtxC[2] = -cst*w113*r222;
		mtxC[3] = -w211 - cst*w213*r212 - cst*w213*r232*r111;
		mtxC[4] = -w212 - cst*w213*r232*r121;
		mtxC[5] = -cst*w213*r222;
		mtxC[6] = -w321 - cst*w323*r131*r212 - cst*w323*r111;
		mtxC[7] = -cst*w323*r121;
		mtxC[8] = -w322 - cst*w323*r131*r222;

		mtxD[0] = cst*w113*r232*(r111 + r121 + r131 - 1) + w111 + w112 + w113;
		mtxD[1] = cst*w113*(r212 + r222 + r232 - 1);
		mtxD[4] = w211 + w212 + w213 + cst*w213*r232*(r111 + r121 + r131 - 1);
		mtxD[5] = cst*w213*(r212 + r222 + r232 - 1);
		mtxD[8] = cst*w323*(r111 + r121 + r131 - 1);
		mtxD[9] = cst*w323*r131*(r212 + r222 + r232 - 1) + w321 + w322 + w323;
		for (int i = 0; i < 12; ++i){if (i % 4 == 2 || i % 4 == 3)mtxD[i] = 0.0;}

		//Calculate the transmissibility matrix Td for L-stencil(d)
		double CA[9], CA_B[12];
		SVD_Inv3(mtxA);                   //inverseA(mtxA);
		FracTHMSOL::MM3(CA, mtxC, mtxA);   //MatrA_MatrB(3, 3, 3, mtxC, mtxA, CA  );
		FracTHMSOL::MM334(CA_B, CA, mtxB); // MatrA_MatrB(3, 3, 4, CA,   mtxB, CA_B);
		plus_Mtr(3, 4, CA_B, mtxD);
		copy_Mtr(3, 4, mtxD, Td);
		return CalArea4(SP, linec4, linec5, facec1);
	}
	//************************************************************************************************
}