(*有限要素法による解析*)

(*汎関数の計算*)

(*前準備*)

(*各要素内部の電位分布&PHgr;を一次式(a+b*x+c*y)で近似表現する.*)

z={1,x,y};

(*重みベクトルの設定*)

k=Table[{a[i],b[i],c[i]},{i,nElements}];

(*各接点の座標を与える。面積計算を行列式を用いて行うために、先頭に１を付加する。*)

vertices=(Append[{1},#] // Flatten)& /@ Pts;

(*各要素の接点座標のリストを作成する*)

p=(vertices[[#]]&) /@ ele;

(*各要素の面積をリストsとして求める。*)

s=1/2*Det /@ p;

(*各要素の接点電位を用いて、重みベクトルkを計算する。*)

sol=Table[Solve[v/@ele[[i]]==p[[i]].k[[i]],
						k[[i]]],{i,nElements}];

k=k /. Flatten[sol];

(*各要素内の電位phi[i]を計算する。*)

phi=k.z;

(*重心座標の修正*)

Do[center[i*(2*n+1)+1]=center[i*(2*n+1)+3],{i,0,TotalY-1}]

(*全エネルギーを計算する。汎関数を計算する。2次元Laplaceの汎関数*)

(*交互に電荷をおく場合*)

(*2次元場の場合*)

(*Laplace方程式*)

u0=Sum[ipsilon0/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]],
		{i,1,n0Elements}] //Simplify;

(*この表面層のみPosson方程式にしておく。*)

u1=Sum[ ipsilon1/2 * (D[phi[[i]],x]^2 + D[phi[[i]],y]^2 )*s[[i]]
- EQ[[i-n0Elements]] ( phi[[i]] /. {x->center[1][[1]], y->center[1][[2]]}),
		{i,n0Elements+1,n0Elements+n1Elements}];

u2=Sum[ipsilon2/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]],
		{i,n0Elements+n1Elements+1,n0Elements+n1Elements+n2Elements}];

u3=Sum[ipsilon3/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]],
		{i,n0Elements+n1Elements+n2Elements+1,
			n0Elements+n1Elements+n2Elements+n3Elements}];

(*2次元回転場の場合*)

u0=Sum[ipsilon0/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]]*Abs[center[i][[1]]],
{i,1,n0Elements}];

u1=Sum[ipsilon1/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]]*Abs[center[i][[1]]],
{i,n0Elements+1,n0Elements+n1Elements}];

u2=Sum[ipsilon2/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]]*Abs[center[i][[1]]],
{i,n0Elements+n1Elements+1,n0Elements+n1Elements+n2Elements}];

u3=Sum[ipsilon3/2*(D[phi[[i]],x]^2+D[phi[[i]],y]^2)*s[[i]]*Abs[center[i][[1]]],
{i,n0Elements+n1Elements+n2Elements+1,
n0Elements+n1Elements+n2Elements+n3Elements}];

(*叉して全エネルギーにする*)

u=u0+u1+u2+u3 ;

(*境界条件の設定*)

(*条件は「導電層の電位、両端の電位が同じ」*)

(*電極の電位を与える*)

BC1=Thread[v/@Table[i,{i,2*n+1}]->Table[CarrierVoltage,{i,2*n+1}]];

(*表面電位の設定*)

BC5=Thread[v/@Table[(nCto1)(2 n +1)+1+2* i,{i,0,n}]->Table[Layer1Voltage+(CarrierVoltage-Layer1Voltage)i/n,{i,0,n}]];

(*4Layerの電位を与える*)

BC2=Thread[v/@Table[i,
	{i,(TotalY-1)*(2*n+1)+1,(TotalY)*(2*n+1)}]->Table[Layer4Voltage,{i,2*n+1}]];

(*左右の境界が同じ電位であるという条件*)

BC3=Thread[v/@Table[(i+1)*(2*n+1),{i,1,TotalY-2}]->v/@Table[(i+1)*(2*n+1)-2,
			{i,1,TotalY-2}]];

BC4=Thread[v/@Table[i*(2*n+1)+1,{i,1,TotalY-2}]->v/@Table[i*(2*n+1)+3,
			{i,1,TotalY-2}]];

(*全ての境界条件を総合すると*)

BC=Join[BC1,BC2,BC3,BC4];(*, BC5*)

(*汎関数に境界条件を代入する*)

u=u/. BC (*//Simplify*);

(*方程式を解く*)

(*未知電位で変分し、汎関数が極値になる条件を求める。*)

Solve[Flatten[Table[D[u,#]&[v[j*(2*n+1)+i+1]]==0,{i,2*n-1},{j,TotalY-2}]],
		Flatten[Table[v[j*(2*n+1)+i+1],{i,2*n-1},{j,TotalY-2}]] ];

(*計算結果をvに代入する*)

Result0=Flatten[ v/@Flatten[Table[j*(2*n+1)+i+1,{i,2*n-1},{j,TotalY-2}]] /.%];

Aug0=Flatten[Table[j*(2*n+1)+i+1,{i,2*n-1},{j,TotalY-2}]];

Do[v[Aug0[[i]]]=Part[Result0,i],{i,Length[Aug0]}]

(*電極側の電位を設定*)

Do[v[i]=CarrierVoltage,{i,2*n+1}]

(*4Layer側の導電部の電位を代入*)

Do[v[i]=Layer4Voltage,{i,(TotalY-1)*(2*n+1)+1,(TotalY)*(2*n+1)}]

(*両端の電位を等しくなるように代入する*)

Do[v[(i+1)*(2*n+1)]=v[(i+1)*(2*n+1)-2],{i,1,TotalY-2}]

Do[v[i*(2*n+1)+1]=v[i*(2*n+1)+3],{i,1,TotalY-2}]

(*Do[v[(nCto1)(2 n +1)+1+2* i]=Layer1Voltage+(CarrierVoltage-Layer1Voltage)i/n,{i,0,n}]*)