(*プリ処理*)

(*使用変数の計算*)

(*Y方向の分割数を与える*)

TotalY=nCto1+n1to2+n2to3+n3to4+1;

(*各層の要素の総数を与える*)

n0Elements=2*n*2*nCto1;

n1Elements=2*n*2*n1to2;

n2Elements=2*n*2*n2to3;

n3Elements=2*n*2*n3to4;

nElements=n0Elements+n1Elements+n2Elements+n3Elements;

(*接点の座標の計算*)

(*計算の前処理*)

(*X軸上の位置から、電極表面上のYを計算する関数を設定*)

(* XtoY[x_]=0;  (* in Flat Case*) *)

(* XtoY[x_]=N[( Abs[CarrierRadius^2 - x^2] )^(1/2),3]; (*in Circle Case*) *)

(*電極表面のY座標を計算する*)

YCPos=Rationalize[ XtoY /@ XPos ,10^-2];

YCPos[[1]]=YCPos[[3]];

(*計算*)

(*電極表面上の位置を設定*)

OnCarrierPtr=Table[{Part[XPos,i],Part[YCPos,i]}, (* 1/10 is Gimick*)
			{i,0+1,2*n+1}];

(*1Layer表面上の位置を設定*)

On1LayerPtr=Table[{Part[XPos,i],CarrierRadius+1/1000},
			{i,0+1,2*n+1}];

(*2Layer表面上の位置を設定*)

On2LayerPtr=Table[{Part[XPos,i],CarrierRadius+Height1Layer},
			{i,0+1,2*n+1}];

(*3Layer表面上の位置を設定*)

On3LayerPtr=Table[{Part[XPos,i],
					CarrierRadius+Height1Layer+Height2Layer},
			{i,0+1,2*n+1}];

(*4Layer表面上の位置を設定*)

On4LayerPtr=Table[{Part[XPos,i],
		CarrierRadius+Height1Layer+Height2Layer+Height3Layer},
			{i,0+1,2*n+1}];

(*電極と1Layer間の接点の位置を設定*)

Cto1Ptr=Flatten[
Table[((nCto1-i)*OnCarrierPtr+i*On1LayerPtr)/nCto1,
					{i,nCto1-1}]
		,1];

(*1Layerの接点の位置を設定*)

Ptr12=Flatten[
Table[((n1to2-i)*On1LayerPtr+i*On2LayerPtr)/n1to2,
					{i,n1to2-1}]
		,1];

(*2Layerの接点の位置を設定*)

Ptr23=Flatten[
Table[((n2to3-i)*On2LayerPtr+i*On3LayerPtr)/n2to3,
					{i,n2to3-1}]
		,1];

(*3Layerの接点の位置を設定*)

Ptr34=Flatten[
Table[((n3to4-i)*On3LayerPtr+i*On4LayerPtr)/n3to4,
					{i,0,n3to4}]
		,1];

(*接点をリスト形式で整理する*)

Pts=
Join[
OnCarrierPtr,
Cto1Ptr,
On1LayerPtr,
Ptr12,
On2LayerPtr,
Ptr23,
Ptr34
];

(*要素の総数を与える*)

nPts=Length[Pts];

(*接点の表示*)

(*接点を描く*)

g1=PlanarGraphPlot[Pts,
		DefaultFont->{"BoldCourier",22.},PlotRange->All]; 

(*誘電体層の境界を描く*) (*あとで、電極境界も付け叉す*)

g0=Graphics[{Thickness[0.01],
   Line[{ {-EvaluateX,CarrierRadius},{EvaluateX,CarrierRadius} }],
   Line[{ {-EvaluateX,CarrierRadius+Height1Layer},
						{EvaluateX,CarrierRadius+Height1Layer} }],
   Line[{ {-EvaluateX,CarrierRadius+Height1Layer+Height2Layer},
				{EvaluateX,CarrierRadius+Height1Layer+Height2Layer} }],
   Line[{ {-EvaluateX,CarrierRadius+Height1Layer+Height2Layer+Height3Layer},
		{EvaluateX,CarrierRadius+Height1Layer+Height2Layer+Height3Layer} }]
}];

Show[g0,Axes->True];

(*接点と誘電体層の境界を描く*)

Show[g1,g0,Axes->True]

(*要素の計算*)

(*各要素の接点番号をリストで設定*)

(*まずは電極と1Layer間の要素*)

ele0un=Flatten[
Table[{1+i+j*(2*n+1),2+i+j*(2*n+1),2+i+(j+1)*(2*n+1)},
						{j,0,nCto1-1},{i,0,2*n-1}]
		,1];

ele0up=Flatten[
Table[{1+i+j*(2*n+1),2+i+(j+1)*(2*n+1),1+i+(j+1)*(2*n+1)},
						{j,0,nCto1-1},{i,0,2*n-1}]
		,1];

usedElementsN=(2*n+1)*nCto1;

(*1Layerの要素*)

ele1un=Flatten[
Table[{usedElementsN+1+i+j*(2*n+1),usedElementsN+2+i+j*(2*n+1),
			usedElementsN+2+i+(j+1)*(2*n+1)},
						{j,0,n1to2-1},{i,0,2*n-1}]
		,1];

ele1up=Flatten[
Table[{usedElementsN+1+i+j*(2*n+1),usedElementsN+2+i+(j+1)*(2*n+1),
		usedElementsN+1+i+(j+1)*(2*n+1)},
						{j,0,n1to2-1},{i,0,2*n-1}]
		,1];

usedElementsN=(2*n+1)*nCto1+(2*n+1)*n1to2;

(*2Layerの要素*)

ele2un=Flatten[
Table[{usedElementsN+1+i+j*(2*n+1),usedElementsN+2+i+j*(2*n+1),
			usedElementsN+2+i+(j+1)*(2*n+1)},
						{j,0,n2to3-1},{i,0,2*n-1}]
		,1];

ele2up=Flatten[
Table[{usedElementsN+1+i+j*(2*n+1),usedElementsN+2+i+(j+1)*(2*n+1),
				usedElementsN+1+i+(j+1)*(2*n+1)},
						{j,0,n2to3-1},{i,0,2*n-1}]
		,1];

usedElementsN=(2*n+1)*nCto1+(2*n+1)*n1to2+(2*n+1)*n2to3;

(*3Layerの要素*)

ele3un=Flatten[
Table[{usedElementsN+1+i+j*(2*n+1),usedElementsN+2+i+j*(2*n+1),
			usedElementsN+2+i+(j+1)*(2*n+1)},
						{j,0,n3to4-1},{i,0,2*n-1}]
		,1];

ele3up=Flatten[
Table[{usedElementsN+1+i+j*(2*n+1),usedElementsN+2+i+(j+1)*(2*n+1),
					usedElementsN+1+i+(j+1)*(2*n+1)},
						{j,0,n3to4-1},{i,0,2*n-1}]
		,1];

(*すべての要素を一つのリストにする*)

ele=
Join[
ele0un,
ele0up,
ele1un,
ele1up,
ele2un,
ele2up,
ele3un,
ele3up
];

(*各要素の重心を計算する*)

center[i_]:=Apply[Plus,Pts[[ele[[i]]]]]/3

(*設定した要素の表示*)

(*設定した要素の重心を表示する*)

gCenter=Show[
	Graphics[
		Table[Text[i,center[i]],
			{i,nElements}]],
		DefaultFont->{"Times",10},Axes->True
		];

(*総図を描く、メモリエラーが起きやすい？*)

Show[g0,g1,gCenter,PlotRange->All]