钻孔孔弹性应力分析benchmark

钻孔孔弹性应力分析benchmark

基于 Cheng.H.D(1988)Poroelastic Responseof a Borehole in a Non-hydrostatic StressField
1.为方便理论（解析）解求解，Cheng.H.D&E. DETOURNAY*将钻孔泄压问题分成3个子问题，分别推导解析解，再利用叠加原理，将三个解相加得到完整解。其中，子问题1和3不考虑渗流作用，任何一本弹性力学书上都有Lame解的推导过程，在此不再赘述。
2.针对子问题2，即渗流作用对孔口应力场的影响，是孔隙线弹性区别传统线弹性独有的现象，我们称之为“负压传递（Negative pressure transfer）”。简单说，就是有一个压应力区像波一样，短时间内由孔口开始径向荡开。其压缩区域越来越大， 且区域越来越大，一定时间后求解域全部为压应力。
3.该问题的有限元解国内外已经做出，所得结论与现象一致，在此不再赘述。本博客提供体积元（FVM）求解方案，包括完整的代码以及数据结果。有限体积元基函数具有能够完全反映单元内参数变化的优势，因此计算流体力学问题首选该方法，但在固体力学上应用不多见。
4.本构与控制方程
该渗流应力耦合问题本质上是去求解一个椭圆型方程和一个抛物线型方程。通过两个耦合项体现流体和孔隙骨架相互影响。加上适合的定解条件即可确定结果。
5.求解结果展示
图1为孔隙压力随时间的分布
图2为Cheng.H.D的数值结果（仿真目标）
图3为切向应力场随时间的分布
图4为Cheng.H.D的数值结果

可以看到体积元（FVM）成功地对该问题进行了计算。值得一提的是相对于传统有限元方案，我们提出的FVM 方法对划分网格质量要求没有那么高，且求解收敛性好很多（更容易收敛），同时计算速度更快。根据实际测试，传统有限元计算时间为3分21秒，而本文提出的FVM计算时间仅为47秒。可以想见，针对未来更加复杂的模型，其优势会更加明显。
6.以下为计算源代码

#include <iostream>
#include <math.h>
#include <string>
#include <vector>

using namespace std;

class coefficientsAssembly
{
   
public:
	// Class variables
	vector<vector<double>> coefficientsMatrix;
	vector<vector<int>> boundaryConditionType;
	int Nu, Nv, NP, NPM;
	vector<vector<int>> uDisplacementFVIndex;
	vector<vector<int>> vDisplacementFVIndex;
	vector<vector<int>> pressureFVIndex;
	vector<vector<int>> uDisplacementFVCoordinates;
	vector<vector<int>> vDisplacementFVCoordinates;
	vector<vector<int>> pressureFVCoordinates;
	vector<vector<int>> horizontalFacesStatus;
	vector<vector<int>> verticalFacesStatus;
	string gridType;
	string interpScheme;
	vector<double> sparseCoefficientsRow;
	vector<double> sparseCoefficientsColumn;
	vector<double> sparseCoefficientsValue;

	// Class functions
	void resizeLinearProblem();
	int getUDisplacementFVPosition(int,int);
	int getVDisplacementFVPosition(int,int);
	int getPressureFVPosition(int,int);
	int getMacroPressureFVPosition(int,int);
	void assemblyCoefficientsMatrix(double,double,double,double,double,double,double,double,double,
		double,double);
	void assemblyXMomentum(double,double,double,double,double);
	void assemblyYMomentum(double,double,double,double,double);
	void assemblyContinuity(double,double,double,double,double,double,double,double,double);
	void addUDisplacementToXMomentum(double,double,double,double);
	void addVDisplacementToXMomentum(double,double,double,double);
	void addPressureToXMomentum(double,double);
	void addBCToXMomentum(double,double,double);
	void addUDisplacementToYMomentum(double,double,double,double);
	void addVDisplacementToYMomentum(double,double,double,double);
	void addPressureToYMomentum(double,double);
	void addBCToYMomentum(double,double,double);
	void addTransientToContinuity(double,double,double,double);
	void addFluidFlowToContinuity(double,double,double,double,double,double,double,double);
	void addDisplacementToContinuity(double,double,double,double,double,double);
	void addBCToContinuity();
	void addStaggeredVDisplacementToXMomentum(double,double,double,double);
	void addStaggeredPressureToXMomentum(double,double);
	void addStaggeredUDisplacementToYMomentum(double,double,double,double);
	void addStaggeredPressureToYMomentum(double,double);
	void addStaggeredFluidFlowToContinuity(double,double,double,double);
	void addStaggeredDisplacementToContinuity(double,double,double,double);
	void addCollocatedFluidFlowToContinuity(double,double,double,double);
	void addCDSVDisplacementToXMomentum(double,double,double,double);
	void addCDSPressureToXMomentum(double,double);
	void addCDSUDisplacementToYMomentum(double,double,double,double);
	void addCDSPressureToYMomentum(double,double);
	void addCDSDisplacementToContinuity(double,double,double,double);
	void addDirichletBCToXMomentum(double,double,double,int);
	void addDirichletBCToYMomentum(double,double,double,int);
	void addDirichletBCToContinuity(int);
	void add1DPISFluidFlowToContinuity(double,double,double,double,double,double);
	void addI2DPISFluidFlowToContinuity(double,double,double,double,double);
	void addI2DPISDisplacementToContinuity(double,double,double,double);
	void addC2DPISFluidFlowToContinuity(double,double,double,double,double,double);
	void addC2DPISDisplacementToContinuity(double,double,double,double,double,double);
	void assemblySparseMatrix(vector<vector<double>>);
	void assemblyMandelCoefficientsMatrix(double,double,double,double,double);
	void addMandelRigidMotion();
	void increaseMandelCoefficientsMatrixSize();
	void addMandelStaggeredStressToVDisplacement(double);
	void addMandelStaggeredStress(double,double,double,double,double);
	void addMandelCollocatedStressToVDisplacement(double);
	void addMandelCollocatedStress(double,double,double,double,double);
	void addStripfootBC(int,double,double);
	void assemblyDoublePorosityMatrix(double,double,double,double,double,double,double,double,double,double,double,double,double,double,double);
	void increaseMacroPorosityCoefficientsMatrixSize();
	void addMacroPressureToXMomentum(double,double);
	void addStaggeredMacroPressureToXMomentum(double,double);
	void addCDSMacroPressureToXMomentum(double,double);
	void addMacroPressureToYMomentum(double,double);
	void addStaggeredMacroPressureToYMomentum(double,double);
	void addCDSMacroPressureToYMomentum(double,double);
	void addMacroTransientToContinuity(double,double,double,double,double);
	void addMacroDisplacementToContinuity(double,double,double,double,double,double);
	void addStaggeredMacroDisplacementToContinuity(double,double,double,double);
	void addCDSMacroDisplacementToContinuity(double,double,double,double);
	void addI2DPISMacroDisplacementToContinuity(double,double,double,double);
	void addMacroFluidFlowToContinuity(double,double,double,double,double,double,double,double,
		double);
	void addStaggeredMacroFluidFlowToContinuity(double,double,double,double);
	void addCollocatedMacroFluidFlowToContinuity(double,double,double,double);
	void addI2DPISFluidFlowToMicroContinuity(double,double,double,double,double,double);
	void addLeaktoContinuity(double);
	void addMacroBCToContinuity();
	void addMacroDirichletBCToContinuity(int);
	void assignFakePressure(double,double);
	void addMacroStripfootBC(int,double,double);

	// Constructor
	coefficientsAssembly(vector<vector<int>>,int,int,int,vector<vector<int>>,vector<vector<int>>,
		vector<vector<int>>,vector<vector<int>>,vector<vector<int>>,vector<vector<int>>,
		vector<vector<int>>,vector<vector<int>>,string,string);

	// Destructor
	~coefficientsAssembly();
};

coefficientsAssembly::coefficientsAssembly(vector<vector<int>> bcType, 
	int numberOfActiveUDisplacementFV, int numberOfActiveVDisplacementFV,
	int numberOfActivePressureFV, vector<vector<int>> idU, vector<vector<int>> idV,
	vector<vector<int>> idP, vector<vector<int>> cooU, vector<vector<int>> cooV,
	vector<vector<int>> cooP, vector<vector<int>> horFStatus, vector<vector<int>> verFStatus,
	string myGridType, string myInterpScheme)
{
   
	boundaryConditionType=bcType;
	Nu=numberOfActiveUDisplacementFV;
	Nv=numberOfActiveVDisplacementFV;
	NP=numberOfActivePressureFV;
	NPM=NP;
	uDisplacementFVIndex=idU;
	vDisplacementFVIndex=idV;
	pressureFVIndex=idP;
	uDisplacementFVCoordinates=cooU;
	vDisplacementFVCoordinates=cooV;
	pressureFVCoordinates=cooP;
	horizontalFacesStatus=horFStatus;
	verticalFacesStatus=verFStatus;
	gridType=myGridType;
	interpScheme=myInterpScheme;

	resizeLinearProblem();
}

coefficientsAssembly::~coefficientsAssembly(){
   }

void coefficientsAssembly::resizeLinearProblem()
{
   
	int rowNo, colNo;

	// Resize coefficientsMatrix
	rowNo=NP+Nv+Nu;
	colNo=rowNo;
	coefficientsMatrix.resize(rowNo);
	for(int i=0; i<rowNo; i++)
	{
   
		coefficientsMatrix[i].resize(colNo);
		for(int j=0; j<colNo; j++)
		{
   
			coefficientsMatrix[i][j]=0;
		}
	}

	return;
}

int coefficientsAssembly::getUDisplacementFVPosition(int x, int y)
{
   
	int uDisplacementFVPosition;

	uDisplacementFVPosition=uDisplacementFVIndex[x][y]-1;

	return uDisplacementFVPosition;
}

int coefficientsAssembly::getVDisplacementFVPosition(int x, int y)
{
   
	int vDisplacementFVPosition;

	vDisplacementFVPosition=vDisplacementFVIndex[x][y]+Nu-1;

	return vDisplacementFVPosition;
}

int coefficientsAssembly::getPressureFVPosition(int x, int y)
{
   
	int pressureFVPosition;

	pressureFVPosition=pressureFVIndex[x][y]+Nu+Nv-1;

	return pressureFVPosition;
}

int coefficientsAssembly::getMacroPressureFVPosition(int x, int y)
{
   
	int pressureFVPosition;

	pressureFVPosition=pressureFVIndex[x][y]+Nu+Nv+NP-1;

	return pressureFVPosition;
}

void coefficientsAssembly::assemblyCoefficientsMatrix(double dx, double dy, double dt, double G,
	double lambda, double alpha, double K, double mu_f, double Q, double rho, double g)
{
   
	assemblyXMomentum(dx,dy,G,lambda,alpha);
	assemblyYMomentum(dx,dy,G,lambda,alpha);
	assemblyContinuity(dx,dy,dt,alpha,K,mu_f,Q,G,lambda);
	assemblySparseMatrix(coefficientsMatrix);

	return;
}

void coefficientsAssembly::assemblyXMomentum(double dx, double dy, double G, double lambda,
	double alpha)
{
   	
	addUDisplacementToXMomentum(dx,dy,G,lambda);
	addVDisplacementToXMomentum(dx,dy,G,lambda);
	addPressureToXMomentum(dy,alpha);
	addBCToXMomentum(dx,dy,G);

	return;
}

void coefficientsAssembly::assemblyYMomentum(double dx, double dy, double G, double lambda,
	double alpha)
{
   	
	addUDisplacementToYMomentum(dx,dy,G,lambda);
	addVDisplacementToYMomentum(dx,dy,G,lambda);
	addPressureToYMomentum(dx,alpha);
	addBCToYMomentum(dx,dy,G);

	return;
}

void coefficientsAssembly::assemblyContinuity(double dx, double dy, double dt, double alpha,
	double K, double mu_f, double Q, double G, double lambda)
{
   
	addTransientToContinuity(dx,dy,dt,Q);
	addFluidFlowToContinuity(dx,dy,dt,K,mu_f,alpha,G,lambda);
	addDisplacementToContinuity(dx,dy,dt,alpha,G,lambda);
	addBCToContinuity();

	return;
}

void coefficientsAssembly::addUDisplacementToXMomentum(double dx, double dy, double G,
	double lambda)
{
   
	int u_P, u_E, u_W, u_N, u_S;
	int FVCounter;
	int i, j;
	double value=1;

	for(FVCounter=0; FVCounter<Nu; FVCounter++)
	{
   
		i=uDisplacementFVCoordinates[FVCounter][0]-1;
		j=uDisplacementFVCoordinates[FVCounter][1]-1;

		u_P=getUDisplacementFVPosition(i,j);

		if(i==0) // Northern border
		{
   
			u_S=getUDisplacementFVPosition(i+1,j);

			if(j==0 || j==uDisplacementFVIndex[0].size()-1)	value=0.5;

			coefficientsMatrix[u_P][u_S]-=G*(dx/dy)*value;
			coefficientsMatrix[u_P][u_P]+=G*(dx/dy)*value;

			value=1;
		}
		else if(i==uDisplacementFVIndex.size()-1) // Southern border
		{
   
			u_N=getUDisplacementFVPosition(i-1,j);

			if(j==0 || j==uDisplacementFVIndex[0].size()-1) value=0.5;

			coefficientsMatrix[u_P][u_N]-=G*(dx/dy)*value;
			coefficientsMatrix[u_P][u_P]+=G*(dx/dy)*value;

			value=1;
		}
		else
		{
   
			u_N=getUDisplacementFVPosition(i-1,j);
			u_S=getUDisplacementFVPosition(i+1,j);

			if(j==0 || j==uDisplacementFVIndex[0].size()-1) value=0.5;

			coefficientsMatrix[u_P][u_N]-=G*(dx/dy)*value;
			coefficientsMatrix[u_P][u_S]-=G*(dx/dy)*value;
			coefficientsMatrix[u_P][u_P]+=2*G*(dx/dy)*value;

			value=1;
		}

		if(j==0) // Western border
		{
   
			u_E=getUDisplacementFVPosition(i,j+1);

			if(i==0 || i==uDisplacementFVIndex.size()-1) if(gridType!="staggered") value=0.5;

			coefficientsMatrix[u_P][u_E]-=(2*G+lambda)*(dy/dx)*value;
			coefficientsMatrix[u_P][u_P]+=(2*G+lambda)*(dy/dx)*value;

			value=1;
		}
		else if(j==uDisplacementFVIndex[0].size()-1) // Eastern border
		{
   
			u_W=getUDisplacementFVPosition(i,j-1);

			if(i==0 || i==uDisplacementFVIndex.size()-1) if(gridType!="staggered") value=0.5;

			coefficientsMatrix[u_P][u_W]-=(2*G+lambda)*(dy/dx)*value;
			coefficientsMatrix[u_P][u_P]+=(2*G+lambda)*(dy/dx)*value;

			value=1;
		}
		else
		{
   
			u_E=getUDisplacementFVPosition(i,j+1);
			u_W=getUDisplacementFVPosition(i,j-1);

			if(i==0 || i==uDisplacementFVIndex.size()-1) if(gridType!="staggered") value=0.5;

			coefficientsMatrix[u_P][u_E]-=(2*G+lambda)*(dy/dx)*value;
			coefficientsMatrix[u_P][u_W]-=(2*G+lambda)*(dy/dx)*value;
			coefficientsMatrix[u_P][u_P]+=2*(2*G+lambda)*(dy/dx)*value;

			value=1;
		}
	}

	return;
}

void coefficientsAssembly::addVDisplacementToXMomentum(double dx, double dy, double G,
	double lambda)
{
   
	if(gridType=="staggered") addStaggeredVDisplacementToXMomentum(dx,dy,G,lambda);
	else if(gridType=="collocated")
	{
   
		if(interpScheme=="CDS") addCDSVDisplacementToXMomentum(dx,dy,G,lambda);
		if(interpScheme=="1DPIS") addCDSVDisplacementToXMomentum(dx,dy,G,lambda);
		if(interpScheme=="I2DPIS") addCDSVDisplacementToXMomentum(dx,dy,G,lambda);
		if(interpScheme=="C2DPIS") addCDSVDisplacementToXMomentum(dx,dy,G,lambda);
	}

	return;
}

void coefficientsAssembly::addPressureToXMomentum(double dy, double alpha)
{
   
	if(gridType=="staggered") addStaggeredPressureToXMomentum(dy,alpha);
	else if(gridType=="collocated")
	{
   
		if(interpScheme=="CDS") addCDSPressureToXMomentum(dy,alpha);
		if(interpScheme=="1DPIS") addCDSPressureToXMomentum(dy,alpha);
		if(interpScheme=="I2DPIS") addCDSPressureToXMomentum(dy,alpha);
		if(interpScheme=="C2DPIS") addCDSPressureToXMomentum(dy,alpha);
	}

	return;
}

void coefficientsAssembly::addBCToXMomentum(double dx, double dy, double G)
{
   
	addDirichletBCToXMomentum(dx,dy,G,0);
	addDirichletBCToXMomentum(dx,dy,G,2);
	addDirichletBCToXMomentum(dx,dy,G,1);
	addDirichletBCToXMomentum(dx,dy,G,3);

	return;
}

void coefficientsAssembly::addUDisplacementToYMomentum(double dx, double dy, double G,
	double lambda)
{
   
	if(gridType=="staggered") addStaggeredUDisplacementToYMomentum(dx,dy,G,lambda);
	else if(gridType=="collocated")
	{
   
		if(interpScheme=="CDS") addCDSUDisplacementToYMomentum(dx,dy,G,lambda);
		if(interpScheme=="1DPIS") addCDSUDisplacementToYMomentum(dx,dy,G,lambda);
		if(interpScheme=="I2DPIS") addCDSUDisplacementToYMomentum(dx,dy,G,lambda);
		if(interpScheme=="C2DPIS") addCDSUDisplacementToYMomentum(dx,dy,G,lambda);
	}

	return;
}

void coefficientsAssembly::addVDisplacementToYMomentum(double dx, double dy, double G,
	double lambda)
{
   
	int v_P, v_E, v_W, v_N, v_S;
	int FVCounter;
	int i, j;
	double value=1;

	for(FVCounter=0; FVCounter<Nv; FVCounter++)
	{
   
		i=vDisplacementFVCoordinates[FVCounter][0]-1;
		j=vDisplacementFVCoordinates[FVCounter][1]-1;

		v_P=getVDisplacementFVPosition(i,j);

		if(i==0) // Northern border
		{
   
			v_S=getVDisplacementFVPosition(i+1,j);

			if(j==0 || j==vDisplacementFVIndex[0].size()-1) if(gridType!="staggered") value=0.5;

			coefficientsMatrix[v_P][v_S]-=(2*G+lambda)*(dx/dy)*value;
			coefficientsMatrix[v_P][v_P]+=(2*G+lambda)*(dx/dy)*value;

			value=1;
		}
		else if(i==vDisplacementFVIndex.size()-1) // Southern border
		{
   
			v_N=getVDisplacementFVPosition(i-1,j);

			if(j==0 || j==vDisplacementFVIndex[0].size()-1) if(gridType!="staggered") value=0.5;

			coefficientsMatrix[v_P][v_N]-=(2*G+lambda)*(dx/dy)*value;
			coefficientsMatrix[v_P][v_P]+=(2*G+lambda)*(dx/dy)*value;

			value=1;
		}
		else
		{
   
			v_N=getVDisplacementFVPosition(i-1,j);
			v_S=getVDisplacementFVPosition(i+1,j);

			if(j==0 || j==vDisplacementFVIndex[0].size()-1) if(gridType!="staggered") value=0.5;

			coefficientsMatrix[v_P][v_N]-=(2*G+lambda)*(dx/dy)*value;
			coefficientsMatrix[v_P][v_S]-=(2*G+lambda)*(dx/dy)*value;
			coefficientsMatrix[v_P][v_P]+=2*(2*G+lambda)*(dx/dy)*value;

			value=1;
		}

		if(j==0) // Western border
		{
   
			v_E=getVDisplacementFVPosition(i,j+1);

			if(i==0 || i==vDisplacementFVIndex.size()-1) value=0.5;

			coefficientsMatrix[v_P][v_E]-=G*(dy/dx)*value;
			coefficientsMatrix[v_P][v_P]+=G*(dy/dx)*value;

			value=1;
		}
		else if(j==vDisplacementFVIndex[0].size()-1) // Eastern border
		{
   
			v_W=getVDisplacementFVPosition(i,j-1);

			if(i==0 || i==vDisplacementFVIndex.size()-1) value=0.5;

			coefficientsMatrix[v_P][v_W]-=G*(dy/dx)*value;
			coefficientsMatrix[v_P][v_P]+=G*(dy/dx)*value;

			value=1;
		}
		else
		{
   
			v_E=getVDisplacementFVPosition(i,j+1);
			v_W=getVDisplacementFVPosition(i,j-1);

			if(i==0 || i==vDisplacementFVIndex.size()-1) value=0.5;

			coefficientsMatrix[v_P][v_E]-=G*(dy/dx)*value;
			coefficientsMatrix[v_P][v_W]-=G*(dy/dx)*value;
			coefficientsMatrix[v_P][v_P]+=2*G*(dy/dx)*value;

			value=1;
		}
	}

	return;
}

void coefficientsAssembly::addPressureToYMomentum(double dx, double alpha)
{
   
	if(gridType=="staggered") addStaggeredPressureToYMomentum(dx,alpha);
	else if(gridType=="collocated")
	{
   
		if(interpScheme=="CDS") addCDSPressureToYMomentum(dx,alpha);
		if(interpScheme=="1DPIS") addCDSPressureToYMomentum(dx,alpha);
		if(interpScheme=="I2DPIS") addCDSPressureToYMomentum(dx,alpha);
		if(interpScheme=="C2DPIS") addCDSPressureToYMomentum(dx,alpha);
	}

	return;
}

void coefficientsAssembly::addBCToYMomentum(double dx, double dy, double G)
{
   
	addDirichletBCToYMomentum(dx,dy,G,1);
	addDirichletBCToYMomentum(dx,dy,G,3);
	addDirichletBCToYMomentum(dx,dy,G,0);
	addDirichletBCToYMomentum(dx,dy,G,2);

	return;
}

void coefficientsAssembly::addTransientToContinuity(double dx, double dy, double dt, double Q)
{
   	
	int P_P;
	int FVCounter;
	int i, j;
	double Mp=(1/Q)*(dx*dy/dt);
	int borderCounter=0;

	for(FVCounter=0; FVCounter<NP; FVCounter++)
	{
   
		i=pressureFVCoordinates[FVCounter][0]-1;
		j=pressureFVCoordinates[FVCounter][1]-1;

		P_P=getPressureFVPosition(i,j);

		if(gridType=="staggered")
		{
   
			coefficientsMatrix[P_P][P_P]+=Mp;
		}
		else if(gridType=="collocated")
		{
   
			if(i==0 || i==pressureFVIndex.size()-1) borderCounter++;
			if(j==0 || j==pressureFVIndex[0].size()-1) borderCounter++;

			coefficientsMatrix[P_P][P_P]+=Mp/pow(2,borderCounter);

			borderCounter=0;
		}
	}
}

void coefficientsAssembly::addFluidFlowToContinuity(double dx, double dy, double dt, double K,
	double mu_f, double alpha, double G, double lambda)
{
   
	if(gridType=="staggered") addStaggeredFluidFlowToContinuity(dx,dy,K,mu_f);
	else if(gridType=="collocated") 
	{
   
		addCollocatedFluidFlowToContinuity(dx,dy,K,mu_f);

		if(interpScheme=="1DPIS") add1DPISFluidFlowToContinuity(dx,dy,dt,alpha,G,lambda);
		if(interpScheme=="I2DPIS") addI2DPISFluidFlowToContinuity(dx,dy,dt,alpha,G);
		if(interpScheme=="C2DPIS") addC2DPISFluidFlowToContinuity(dx,dy,dt,alpha,G,lambda);
	}

	return;
}

void coefficientsAssembly::addDisplacementToContinuity(double dx, double dy, double dt,
	double alpha, double G, double lambda)
{
   
	if(gridType=="staggered") addStaggeredDisplacementToContinuity(dx,dy,dt,alpha);
	else if(gridType=="collocated")
	{
   
		if(interpScheme=="CDS") addCDSDisplacementToContinuity(dx,dy,dt,alpha);
		if(interpScheme=="1DPIS") addCDSDisplacementToContinuity(dx,dy,dt,alpha);
		if(interpScheme=="I2DPIS") addI2DPISDisplacementToContinuity(dx,dy,dt,alpha);
		if(interpScheme=="C2DPIS") addC2DPISDisplacementToContinuity(dx,dy,dt,alpha,G,lambda);
	}
	
	return;
}

void coefficientsAssembly::addBCToContinuity()
{
   
	if(gridType=="collocated")
		for(int counter=0; counter<4; counter++)
			addDirichletBCToContinuity(counter);

	return;
}

void coefficientsAssembly::addStaggeredVDisplacementToXMomentum(double dx, double dy, double G,
	double lambda)
{
   
	int u_P;
	int v_P, v_W, v_S, v_SW;
	int FVCounter;
	int i, j;

	for(FVCounter=0; FVCounter<Nu; FVCounter++)
	{
   
		i=uDisplacementFVCoordinates[FVCounter][0]-1;
		j=uDisplacementFVCoordinates[FVCounter][1]-1;

		u_P=getUDisplacementFVPosition(i,j);

		if(j==0) // Western border
		{
   
			v_P=getVDisplacementFVPosition(i,j);
			v_S=getVDisplacementFVPosition(i+1,j);

			coefficientsMatrix[u_P][v_P]-=lambda;
			coefficientsMatrix[u_P][v_S]-=-lambda;
		}
		else if(j==uDisplacementFVIndex[0].size()-1) // Eastern border
		{
   
			v_W=getVDisplacementFVPosition(i,j-1);
			v_SW=getVDisplacementFVPosition(i+1,j-1);

			coefficientsMatrix[u_P][v_W]-=-lambda;
			coefficientsMatrix[u_P][v_SW]-=lambda;
		}
		else
		{
   
			if(i==0) // Northern border
			{
   
				v_P=getVDisplacementFVPosition(i,j);
				v_W=getVDisplacementFVPosition(i,j-1);
				v_S=getVDisplacementFVPosition(i+1,j);
				v_SW=getVDisplacementFVPosition(i+1,j-1);

				coefficientsMatrix[u_P][v_P]-=lambda;
				coefficientsMatrix[u_P][v_W]-=-lambda;
				coefficientsMatrix[u_P][v_S]-=-G-lambda;
				coefficientsMatrix[u_P][v_SW]-=G+lambda;
			}
			else if(i==uDisplacementFVIndex.size()-1) // Southern border
			{
   
				v_P=getVDisplacementFVPosition(i,j);
				v_W=getVDisplacementFVPosition(i,j-1);
				v_S=getVDisplacementFVPosition(i+1,j);
				v_SW=getVDisplacementFVPosition(i+1,j-1);

				coefficientsMatrix[u_P][v_P]-=G+lambda;
				coefficientsMatrix[u_P][v_W]-=-G-lambda;
				coefficientsMatrix[u_P][v_S]-=-lambda;
				coefficientsMatrix[u_P][v_SW]-=lambda;
			}	
			else
			{
   
				v_P=getVDisplacementFVPosition(i,j);
				v_W=getVDisplacementFVPosition(i,j-1);
				v_S=getVDisplacementFVPosition(i+1,j);
				v_SW=getVDisplacementFVPosition(i+1,j-1);

				coefficientsMatrix[u_P][v_P]-=G+lambda;
				coefficientsMatrix[u_P][v_W]-=-G-lambda;
				coefficientsMatrix[u_P][v_S]-=-G-lambda;
				coefficientsMatrix[u_P][v_SW]-=G+lambda;
			}
		}
	}

	return;
}

void coefficientsAssembly::addStaggeredPressureToXMomentum(double dy, double alpha)
{
   
	int u_P;
	int P_P, P_W;
	int FVCounter;
	int i, j;

	for(FVCounter=0; FVCounter<Nu; FVCounter++)
	{
   
		i=uDisplacementFVCoordinates[FVCounter][0]-1;
		j=uDisplacementFVCoordinates[FVCounter][1]-1;

		u_P=getUDisplacementFVPosition(i,j);

		if(j==0) // FV on the western border
		{
   
			P_P=getPressureFVPosition(i,j);

			coefficientsMatrix[u_P][P_P]-=-alpha*dy;
		}
		else if(j==uDisplacementFVIndex[0].size()-1) // FV on the eastern border
		{
   
			P_W=getPressureFVPosition(i,j-1);

			coefficientsMatrix[u_P][P_W]-=alpha*dy;
		}
		else // FV not on the western or eastern border
		{
   
			P_P=getPressureFVPosition(i,j);
			P_W=getPressureFVPosition(i,j-1);

			coefficientsMatrix[u_P][P_P]-=-alpha*dy;
			coefficientsMatrix[u_P][P_W]-=alpha*dy;
		}
	}

	return;
}

void coefficientsAssembly::addStaggeredUDisplacementToYMomentum(double dx, double dy, double G,
	double lambda)
{
   
	int u_P, u_E, u_N, u_NE;
	int v_P;
	int FVCounter;
	int i, j;

	for(FVCounter=0; FVCounter<Nv; FVCounter++)
	{
   
		i=vDisplacementFVCoordinates[FVCounter][0]-1;
		j=vDisplacementFVCoordinates[FVCounter][1]-1;

		v_P=getVDisplacementFVPosition(i,j);

		if(i==0) // Northern border
		{
   
			u_P=getUDisplacementFVPosition(i,j);
			u_E=getUDisplacementFVPosition(i,j+1);

			coefficientsMatrix[v_P][u_P]-=lambda;
			coefficientsMatrix[v_P][u_E]-=-lambda;
		}
		else if(i==vDisplacementFVIndex.size()-1) // Southern border
		{
   
			u_N=getUDisplacementFVPosition(i-1,j);
			u_NE=getUDisplacementFVPosition(i-1,j+1);

			coefficientsMatrix[v_P][u_N]-=-lambda;
			coefficientsMatrix[v_P][u_NE]-=lambda;
		}
		else
		{
   
			if(j==0) // Western border
			{
   
				u_P=getUDisplacementFVPosition(i,j);
				u_E=getUDisplacementFVPosition(i,j+1);
				u_N=getUDisplacementFVPosition(i-1,j);
				u_NE=getUDisplacementFVPosition(i-1,j+1);

				coefficientsMatrix[v_P][u_P]-=lambda;
				coefficientsMatrix[v_P][u_E]-=-G-lambda;
				coefficientsMatrix[v_P][u_N]-=-lambda;
				coefficientsMatrix[v_P][u_NE]-=G+lambda;
			}
			else if(j==vDisplacementFVIndex[0].size()-1) // Eastern border
			{
   
				u_P=getUDisplacementFVPosition(i,j);
				u_E=getUDisplacementFVPosition(i,j+1);
				u_N=getUDisplacementFVPosition(i-1,j);
				u_NE=getUDisplacementFVPosition(i-1,j+1);

				coefficientsMatrix[v_P][u_P]-=G+lambda;
				coefficientsMatrix[v_P][u_E]-=-lambda;
				coefficientsMatrix[v_P][u_N]-=-G-lambda;
				coefficientsMatrix[v_P][u_NE]-=lambda;
			}
			else
			{
   
				u_P=getUDisplacementFVPosition(i,j);
				u_E=getUDisplacementFVPosition(i,j+1);
				u_N=getUDisplacementFVPosition(i-1,j);
				u_NE=getUDisplacementFVPosition(i-1,j+1);

				coefficientsMatrix[v_P][u_P]-=G+lambda;
				coefficientsMatrix[v_P][u_E]-=-G-lambda;
				coefficientsMatrix[v_P][u_N]-=-G-lambda;
				coefficientsMatrix[v_P][u_NE]-=G+lambda;
			}
		}
	}

	return;
}

void coefficientsAssembly::addStaggeredPressureToYMomentum(double dx, double alpha)
{
   
	int v_P;
	int P_P, P_N;
	int FVCounter;
	int i, j;
	
	for(FVCounter=0; FVCounter<Nv; FVCounter++)
	{
   
		i=vDisplacementFVCoordinates[FVCounter][0]-1;
		j=vDisplacementFVCoordinates[FVCounter][1]-1;

		v_P=getVDisplacementFVPosition(i,j);

		if(i==0) // FV on the northern border
		{
   
			P_P=getPressureFVPosition(i,j);

			coefficientsMatrix[v_P][P_P]-=alpha*dx;
		}
		else if(i==vDisplacementFVIndex.size()-1) // FV on the southern border
		{
   
			P_N=getPressureFVPosition(i-1,j);

			coefficientsMatrix[v_P][P_N]-=-alpha*dx;
		}
		else // FV not on the northern or southern border
		{
   
			P_P=getPressureFVPosition(i,j);
			P_N=getPressureFVPosition(i-1,j);

			coefficientsMatrix[v_P][P_P]-=alpha*dx;
			coefficientsMatrix[v_P][P_N]-=-alpha*dx;
		}
	}

	return;
}

void coefficientsAssembly::addStaggeredFluidFlowToContinuity(double dx, double dy, double K,
	double mu_f)
{
   
	int P_P, P_E, P_W, P_N, P_S;
	int FVCounter;
	int i, j;
	int bcType;

	for(FVCounter=0; FVCounter<NP; FVCounter++)
	{
   
		i=pressureFVCoordinates[FVCounter][0]-1;
		j=pressureFVCoordinates[FVCounter][1]-1;

		P_P=getPressureFVPosition(i,j);

		if(i==0) // Northern border
		{
   
			P_S=getPressureFVPosition(i+1,j);

			coefficientsMatrix[P_P][P_S]-=(K/mu_f)*(dx/dy);
			coefficientsMatrix[P_P][P_P]+=(K/mu_f)*(dx/dy);

			bcType=boundaryConditionType[0][2];
			if(bcType==1) coefficientsMatrix[P_P][P_P]+=2*(K/mu_f)*(dx/dy);
		}
		else if(i==pressureFVIndex.size()-1) // Southern border
		{
   
			P_N=getPressureFVPosition(i-1,j);

			coefficientsMatrix[P_P][P_N]-=(K/mu_f)*(dx/dy);
			coefficientsMatrix[P_P][P_P]+=(K/mu_f)*(dx/dy);

			bcType=boundaryConditionType[2][2];
			if(bcType==1) coefficientsMatrix[P_P][P_P]+=2*(K/mu_f)*(dx/dy);
		}
		else
		{
   
			P_N=getPressureFVPosition(i-1,j);
			P_S=getPressureFVPosition(i+1,j);

			coefficientsMatrix[P_P][P_N]-=(K/mu_f)*(dx/dy);
			coefficientsMatrix[P_P][P_S]-=(K/mu_f)*(dx/dy);
			coefficientsMatrix[P_P][P_P]+=2*(K/mu_f)*(dx/dy);
		}

		if(j==0) // Western border
		{
   
			P_E=getPressureFVPosition(i,j+1);

			coefficientsMatrix[P_P][P_E]-=(K/mu_f)*(dy/dx);
			coefficientsMatrix[P_P][P_P]+=(K/mu_f)*(dy/dx);

			bcType=boundaryConditionType[1][2];
			if(bcType==1) coefficientsMatrix[P_P][P_P]+=2*(K/mu_f)*(dy/dx);
		}
		else if(j==pressureFVIndex[0].size()-1) // Eastern border
		{
   
			P_W=getPressureFVPosition(i,j-1);

			coefficientsMatrix[P_P][P_W]-=(K/mu_f)*(dy/dx);
			coefficientsMatrix[P_P][P_P]+=(K/mu_f)*(dy/dx);

			bcType=boundaryConditionType[3][2];
			if(bcType==1) coefficientsMatrix[P_P][P_P]+=2*(K/mu_f)*(dy/dx);
		}
		else
		{
   
			P_E=getPressureFVPosition(i,j+1);
			P_W=getPressureFVPosition(i,j-1);

			coefficientsMatrix[P_P][P_E]-=(K/mu_f)*(dy/dx);
			coefficientsMatrix[P_P][P_W]-=(K/mu_f)*(dy/dx);
			coefficientsMatrix[P_P][P_P]+=2*(K/mu_f)*(dy/dx);
		}
	}

	return;
}

void coefficientsAssembly::addStaggeredDisplacementToContinuity(double dx, double dy, double dt,
	double alpha)
{
   
	int u_P, u_E;
	int v_P, v_S;
	int P_P;
	int FVCounter;
	int i, j;

	for(FVCounter=0; FVCounter<NP; FVCounter++)
	{
   
		i=pressureFVCoordinates[FVCounter][0]-1;
		j=pressureFVCoordinates[FVCounter][1]-1;

		u_P=getUDisplacementFVPosition(i,j);
		u_E=getUDisplacementFVPosition(i,j+1);
		v_P=getVDisplacementFVPosition(i,j);
		v_S=getVDisplacementFVPosition(i+1,j);