C++ CFieldWorld::GetEA方法代码示例

static void BoundaryCondition
(unsigned int id_field, const ELSEG_TYPE& elseg_type,  
 MatVec::CBCFlag& bc_flag, const CFieldWorld& world,  
 unsigned int ioffset=0)
{
  if( !world.IsIdField(id_field) ) return;
  const Fem::Field::CField& field = world.GetField(id_field);  
  {	// Assert
    const unsigned int len = bc_flag.LenBlk();
    assert( ioffset < len );
  }
  const unsigned int nlen = field.GetNLenValue();
  
  const std::vector<unsigned int>& aIdEA = field.GetAryIdEA();
  for(unsigned int iea=0;iea<aIdEA.size();iea++){
    unsigned int id_ea = aIdEA[iea];
    const CElemAry& ea = world.GetEA(id_ea);
    unsigned int noes[256];
    if( elseg_type == CORNER && field.GetIdElemSeg(id_ea,CORNER,true,world) != 0 ){
      const Fem::Field::CElemAry::CElemSeg& es = field.GetElemSeg(id_ea,CORNER,true,world);
      unsigned int nnoes = es.Length();
      for(unsigned int ielem=0;ielem<ea.Size();ielem++){
        es.GetNodes(ielem,noes);
        for(unsigned int inoes=0;inoes<nnoes;inoes++){
          for(unsigned int ilen=0;ilen<nlen;ilen++){
            bc_flag.SetBC(noes[inoes],ilen+ioffset);
          }
        }
      }
    }
    if( elseg_type == BUBBLE && field.GetIdElemSeg(id_ea,BUBBLE,true,world) != 0 ){
      const Fem::Field::CElemAry::CElemSeg& es = field.GetElemSeg(id_ea,BUBBLE,true,world);
      unsigned int nnoes = es.Length();
      for(unsigned int ielem=0;ielem<ea.Size();ielem++){
        es.GetNodes(ielem,noes);
        for(unsigned int inoes=0;inoes<nnoes;inoes++){
          for(unsigned int ilen=0;ilen<nlen;ilen++){
            bc_flag.SetBC(noes[inoes],ilen+ioffset);
          }
        }
      }
    }
    if( elseg_type == EDGE   && field.GetIdElemSeg(id_ea,EDGE,true,world) != 0 ){
      const Fem::Field::CElemAry::CElemSeg& es = field.GetElemSeg(id_ea,EDGE,true,world);
      unsigned int nnoes = es.Length();
      for(unsigned int ielem=0;ielem<ea.Size();ielem++){
        es.GetNodes(ielem,noes);
        for(unsigned int inoes=0;inoes<nnoes;inoes++){
          for(unsigned int ilen=0;ilen<nlen;ilen++){
            bc_flag.SetBC(noes[inoes],ilen+ioffset);
          }
        }
      }
    }
  }
}

static bool AddLinSys_AdvectionDiffusion_Static_P1P1(
		double myu, double source, 
		CLinearSystem_Field& ls, 
		const unsigned int id_field_val, const unsigned int id_field_velo, const CFieldWorld& world, 
		unsigned int id_ea )
{
//	std::cout << "Advection Diffusion Static 2D Triangle 3-point 1st order" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TRI );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& val_field = world.GetField(id_field_val);

	if( !world.IsIdField(id_field_velo) ) return false;
	const CField& field_velo = world.GetField(id_field_velo);

	// 角節点の節点配列ID
	unsigned int id_na_c_co = val_field.GetNodeSegInNodeAry(CORNER).id_na_co;
	unsigned int id_ns_c_co = val_field.GetNodeSegInNodeAry(CORNER).id_ns_co;
	unsigned int id_na_c_val = val_field.GetNodeSegInNodeAry(CORNER).id_na_va;
	unsigned int id_ns_c_val = val_field.GetNodeSegInNodeAry(CORNER).id_ns_va;

	unsigned int id_na_c_velo = field_velo.GetNodeSegInNodeAry(CORNER).id_na_va;
	unsigned int id_ns_c_velo = field_velo.GetNodeSegInNodeAry(CORNER).id_ns_ve;
	assert( id_na_c_co != 0 && id_ns_c_co != 0 );
	assert( id_na_c_val != 0 && id_ns_c_val != 0 );
	assert( id_na_c_velo != 0 );
	assert( id_ns_c_velo != 0 );

	const unsigned int nno = 3;
	const unsigned int ndim = 2;

	const CElemAry::CElemSeg& es_c_val = field_velo.GetElemSeg(id_ea,CORNER,true,world);

	unsigned int no_c[nno];	// 要素節点の全体節点番号

	double value_c[nno];		// 要素節点の値
	double coord_c[nno][ndim];	// 要素節点の座標
	double velo_c[nno][ndim];
	
	double emat[nno][nno];	// 要素剛性行列
	double eres_c[nno];	// 要素節点等価内力、外力、残差ベクトル
				
	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(  id_field_val,CORNER,world);
	CVector_Blk&    res_c  = ls.GetResidual(id_field_val,CORNER,world);

	const CNodeAry& na_c_val = world.GetNA(id_na_c_val);
	const CNodeAry::CNodeSeg& ns_c_val = na_c_val.GetSeg(id_ns_c_val);
	const CNodeAry& na_c_velo = world.GetNA(id_na_c_velo);
	const CNodeAry::CNodeSeg& ns_c_velo = na_c_velo.GetSeg(id_ns_c_velo);
	const CNodeAry& na_c_co = world.GetNA(id_na_c_co);
	const CNodeAry::CNodeSeg& ns_c_co = na_c_co.GetSeg(id_ns_c_co);

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		// 要素配列から要素セグメントの節点番号を取り出す
		es_c_val.GetNodes(ielem,no_c);
		// 節点の値を取って来る
		for(unsigned int inoes=0;inoes<nno;inoes++){
			ns_c_co.GetValue(no_c[inoes],coord_c[inoes]);
			ns_c_val.GetValue(no_c[inoes],&value_c[inoes]);
			ns_c_velo.GetValue(no_c[inoes],velo_c[inoes]);
		}

		////////////////////////////////////////////////////////////////

		// 面積を求める
		const double area = TriArea(coord_c[0],coord_c[1],coord_c[2]);
		// 形状関数の微分を求める			
		double dldx[nno][ndim];	// 形状関数のxy微分
		double const_term[nno];	// 形状関数の定数項
		TriDlDx(dldx,const_term,coord_c[0],coord_c[1],coord_c[2]);

		// 要素剛性行列を作る
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			emat[ino][jno] = myu*area*(dldx[ino][0]*dldx[jno][0]+dldx[ino][1]*dldx[jno][1]);
		}
		}
		{
			const double dtmp1 = area*0.08333333333333333333333;
			for(unsigned int ino=0;ino<nno;ino++){
				const double dtmp_0 = dtmp1*(velo_c[0][0]+velo_c[1][0]+velo_c[2][0]+velo_c[ino][0]);
				const double dtmp_1 = dtmp1*(velo_c[0][1]+velo_c[1][1]+velo_c[2][1]+velo_c[ino][1]);
				for(unsigned int jno=0;jno<nno;jno++){
					emat[ino][jno] += dldx[jno][0]*dtmp_0+dldx[jno][1]*dtmp_1;
				}
			}
		}

		// Calc Stabilization Parameter
		double tau;
		{
			const double velo_ave[2] = { 
				(velo_c[0][0]+velo_c[1][0]+velo_c[2][0])/3.0, 
				(velo_c[0][1]+velo_c[1][1]+velo_c[2][1])/3.0 };
			const double norm_v = sqrt(velo_ave[0]*velo_ave[0]+velo_ave[1]*velo_ave[1]);
			const double velo_dir[2] = { velo_ave[0]/norm_v, velo_ave[1]/norm_v };
//.........这里部分代码省略.........

static bool AddLinearSystem_Diffusion_P1b(
		double rho, double alpha, double source,
		double gamma, double dt,
		CLinearSystem_Field& ls, 
		unsigned int id_field_val, const CFieldWorld& world,
		unsigned int id_ea )
{	
//	std::cout << "Diffusion2D Tri P1b" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TRI );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

	const CElemAry::CElemSeg& es_c = field_val.GetElemSeg(id_ea,CORNER,true,world);
	const CElemAry::CElemSeg& es_b = field_val.GetElemSeg(id_ea,BUBBLE,true,world);

	const unsigned int nno_c = 3;
	const unsigned int nno_b = 1;
	const unsigned int ndim = 2;

	unsigned int no_c[nno_c];
	unsigned int no_b;

	double val_c[nno_c], val_b;
	double vval_c[nno_c], vval_b;
	double coord_c[nno_c][ndim];
				
	double dldx[nno_c][ndim];
	double const_term[nno_c];

	double eCmat_cc[nno_c][nno_c], eCmat_cb[nno_c], eCmat_bc[nno_c], eCmat_bb;
	double eMmat_cc[nno_c][nno_c], eMmat_cb[nno_c], eMmat_bc[nno_c], eMmat_bb;
	double eqf_out_c[nno_c], eqf_out_b;
	double eqf_in_c[nno_c], eqf_in_b;

	double emat_cc[nno_c][nno_c], emat_cb[nno_c], emat_bc[nno_c], emat_bb;
	double eres_c[nno_c], eres_b;	// 要素節点等価内力、外力、残差ベクトル

	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(id_field_val,CORNER, world);
	CMatDia_BlkCrs& mat_bb = ls.GetMatrix(id_field_val,BUBBLE, world);
	CMat_BlkCrs&    mat_cb = ls.GetMatrix(id_field_val,CORNER, id_field_val, BUBBLE, world);
	CMat_BlkCrs&    mat_bc = ls.GetMatrix(id_field_val,BUBBLE, id_field_val, CORNER, world);
	////////////////
	CVector_Blk& res_c = ls.GetResidual(id_field_val,CORNER, world);
	CVector_Blk& res_b = ls.GetResidual(id_field_val,BUBBLE, world);

	const CNodeAry::CNodeSeg& ns_c_val = field_val.GetNodeSeg(CORNER,true,world,VALUE);//na_c_val.GetSeg(id_ns_c_val);
	const CNodeAry::CNodeSeg& ns_c_vval = field_val.GetNodeSeg(CORNER,true,world,VELOCITY);//na_c_val.GetSeg(id_ns_c_vval);
	const CNodeAry::CNodeSeg& ns_b_val = field_val.GetNodeSeg(BUBBLE,true,world,VALUE);//na_b_val.GetSeg(id_ns_b_val);
	const CNodeAry::CNodeSeg& ns_b_vval = field_val.GetNodeSeg(BUBBLE,true,world,VELOCITY);//na_b_val.GetSeg(id_ns_b_vval);
	const CNodeAry::CNodeSeg& ns_c_co = field_val.GetNodeSeg(CORNER,false,world,VALUE);//na_c_val.GetSeg(id_ns_c_co);

	for(unsigned int ielem=0;ielem<ea.Size();ielem++){
		// 要素配列から節点セグメントの節点番号を取り出す
		es_c.GetNodes(ielem,no_c);
		es_b.GetNodes(ielem,&no_b);
		// 節点の値を取ってくる
		for(unsigned int inoes=0;inoes<nno_c;inoes++){
			ns_c_co.GetValue(no_c[inoes],coord_c[inoes]);
			ns_c_val.GetValue(no_c[inoes],&val_c[inoes]);
			ns_c_vval.GetValue(no_c[inoes],&vval_c[inoes]);
		}
		ns_b_val.GetValue(no_b,&val_b);
		ns_b_vval.GetValue(no_b,&vval_b);

		// 面積を求める
		const double area = TriArea(coord_c[0],coord_c[1],coord_c[2]);
		// 形状関数の微分を求める
		TriDlDx(dldx,const_term,coord_c[0],coord_c[1],coord_c[2]);
		{	// 要素剛性行列を作る
			double vc_b[4];
			vc_b[0] = 1.0/3.0; vc_b[1] = 1.0/3.0; vc_b[2] = 1.0/3.0; vc_b[3] = 27.0;
			const double tmp_val1 = vc_b[3]*vc_b[3]*alpha*area/180.0*( 
				dldx[0][0]*dldx[0][0]+dldx[0][1]*dldx[0][1]+
				dldx[1][0]*dldx[1][0]+dldx[1][1]*dldx[1][1]+
				dldx[2][0]*dldx[2][0]+dldx[2][1]*dldx[2][1] );
            for(unsigned int ino_c=0;ino_c<nno_c;ino_c++){
            for(unsigned int jno_c=0;jno_c<nno_c;jno_c++){
				eCmat_cc[ino_c][jno_c] 
					= alpha*area*(dldx[ino_c][0]*dldx[jno_c][0]+dldx[ino_c][1]*dldx[jno_c][1])
					+vc_b[ino_c]*vc_b[jno_c]*tmp_val1;
			}
			}
			for(unsigned int ino_c=0;ino_c<nno_c;ino_c++){
				const double tmp1 = -1.0*vc_b[ino_c]*tmp_val1;
				eCmat_cb[ino_c] = tmp1;
				eCmat_bc[ino_c] = tmp1;
			}
			eCmat_bb = tmp_val1;
			Set_RhoTri_CB_Scalar(eMmat_cc,eMmat_cb,eMmat_bc,eMmat_bb,  area, dldx,vc_b, rho);
		}
		// 要素外力ベクトルを求める
        for(unsigned int ino_c=0;ino_c<nno_c;ino_c++){
			eqf_out_c[ino_c] = source*area*11.0/60.0;
		}
		eqf_out_b = source*area*27.0/60.0;

//.........这里部分代码省略.........

static bool AddLinearSystem_Diffusion3D_P1(
		double rho, double alpha, double source,
		CLinearSystem_SaveDiaM_Newmark& ls, 
		unsigned int id_field_val, const CFieldWorld& world,
		unsigned int id_ea )
{
//	std::cout << "Diffusion3D Tet savemat" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TET );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

	const double gamma = ls.GetGamma();
	const double dt = ls.GetDt();

	const CElemAry::CElemSeg& es_c = field_val.GetElemSeg(id_ea,CORNER,true,world);

	const unsigned int nno = 4;
	const unsigned int ndim = 3;

	unsigned int no_c[nno];	// 要素節点の全体節点番号

	double val_c[nno];		// 要素節点の値
	double vval_c[nno];		// 要素節点の値
	double coord_c[nno][ndim];	// 要素節点の座標
				
	double emat[nno][nno];
	double eCmat[nno][nno];	// 要素剛性行列
	double eMmat[nno][nno];	// 要素剛性行列
	double eqf_out_c[nno];	// 要素節点等価内力、外力、残差ベクトル
				
	CMatDia_BlkCrs& mat_cc  = ls.GetMatrix(id_field_val,CORNER,world);
	CVector_Blk&    force_c = ls.GetForce( id_field_val,CORNER,world);
	
	CMat_BlkCrs& mat_cc_bound = ls.GetMatrix_Boundary(id_field_val,CORNER,  id_field_val,CORNER,  world);
	CDiaMat_Blk& mat_mass     = ls.GetDiaMassMatrix(  id_field_val,CORNER,world);

	const CNodeAry::CNodeSeg& ns_c_val = field_val.GetNodeSeg(CORNER,true,world,VALUE);	
	const CNodeAry::CNodeSeg& ns_c_vval = field_val.GetNodeSeg(CORNER,true,world,VELOCITY);
	const CNodeAry::CNodeSeg& ns_c_co = field_val.GetNodeSeg(CORNER,false,world,VALUE);	

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		// 要素配列から要素セグメントの節点番号を取り出す
		es_c.GetNodes(ielem,no_c);
		// 節点の値を取って来る
		for(unsigned int ino=0;ino<nno;ino++){
			ns_c_co.GetValue(no_c[ino],coord_c[ino]);
			ns_c_val.GetValue(no_c[ino],&val_c[ino]);
			ns_c_vval.GetValue(no_c[ino],&vval_c[ino]);
		}

		////////////////////////////////////////////////////////////////

		// 面積を求める
		const double vol = TetVolume(coord_c[0],coord_c[1],coord_c[2],coord_c[3]);
		// 形状関数の微分を求める
		double dldx[nno][ndim];	// 形状関数のxy微分
		double const_term[nno];	// 形状関数の定数項
		TetDlDx(dldx,const_term,coord_c[0],coord_c[1],coord_c[2],coord_c[3]);
		// 要素剛性行列を作る
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			eCmat[ino][jno] = alpha*vol*( dldx[ino][0]*dldx[jno][0]+dldx[ino][1]*dldx[jno][1]+dldx[ino][2]*dldx[jno][2]);
		}
		}
		{
			const double dtmp1 = rho*vol*0.05;
			for(unsigned int ino=0;ino<nno;ino++){
				for(unsigned int jno=0;jno<nno;jno++){
					eMmat[ino][jno] = dtmp1;
				}
				eMmat[ino][ino] += dtmp1;
			}
		}
		// 要素節点等価外力ベクトルを求める
		for(unsigned int ino=0;ino<nno;ino++){
			eqf_out_c[ino] = source*vol*0.25;
		}

		////////////////////////////////////////////////////////////////

		{	// 要素係数行列を求める
			double dtmp1 = gamma*dt;
			for(unsigned int i=0;i<nno*nno;i++){
				(&emat[0][0])[i] = (&eMmat[0][0])[i]+dtmp1*(&eCmat[0][0])[i];
			}
		}

		// 剛性行列にマージする
		mat_cc      .Mearge(nno,no_c,nno,no_c,  1, &emat[0][0]);
		mat_cc_bound.Mearge(nno,no_c,nno,no_c,  1,&eCmat[0][0]);
		for(unsigned int ino=0;ino<nno;ino++){
			mat_mass.Mearge( no_c[ino], 1, &eMmat[0][0] );
        }

		// 残差ベクトルにマージする
//.........这里部分代码省略.........

static bool AddLinearSystem_Poisson3D_Q1(
		double alpha, double source, 
        Fem::Ls::CLinearSystem_Field& ls, 
		const unsigned int id_field_val, const CFieldWorld& world,
		const unsigned int id_ea)
{
//	std::cout << "Poisson3D Hex" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == HEX );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

	const CElemAry::CElemSeg& es_c_co = field_val.GetElemSeg(id_ea,CORNER,false,world);
	const CElemAry::CElemSeg& es_c_va = field_val.GetElemSeg(id_ea,CORNER,true, world);

	unsigned int num_integral = 1;
	const unsigned int nInt = NIntLineGauss[num_integral];
	const double (*Gauss)[2] = LineGauss[num_integral];

	const unsigned int nno = 8;
	const unsigned int ndim = 3;

	unsigned int no_c_co[nno];	// 要素節点の全体節点番号
	unsigned int no_c_va[nno];	// 要素節点の全体節点番号

	double value_c[nno];		// 要素節点の値
	double coord_c[nno][ndim];	// 要素節点の座標
				
	double dndx[nno][ndim];	// 形状関数のxy微分
	double an_c[nno];

	double emat[nno][nno];	// 要素剛性行列
	double eres_c[nno];	// 要素節点等価内力、外力、残差ベクトル
				
	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(  id_field_val,CORNER,world);
	CVector_Blk&    res_c  = ls.GetResidual(id_field_val,CORNER,world);

	const CNodeAry::CNodeSeg& ns_c_val = field_val.GetNodeSeg(CORNER,true, world);
	const CNodeAry::CNodeSeg& ns_c_co  = field_val.GetNodeSeg(CORNER,false,world);

	for(unsigned int ielem=0;ielem<ea.Size();ielem++){
		// 要素配列から要素セグメントの節点番号を取り出す
		es_c_co.GetNodes(ielem,no_c_co);
		for(unsigned int ino=0;ino<nno;ino++){ 
			ns_c_co.GetValue(no_c_co[ino],coord_c[ino]);
		}
		es_c_va.GetNodes(ielem,no_c_va);
		for(unsigned int ino=0;ino<nno;ino++){ 
			ns_c_val.GetValue(no_c_va[ino],&value_c[ino]);
		}
		
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			emat[ino][jno] = 0.0;
		}
		}
		for(unsigned int ino=0;ino<nno;ino++){
			eres_c[ino] = 0.0;
		}
		double vol = 0.0;
		for(unsigned int ir1=0;ir1<nInt;ir1++){
		for(unsigned int ir2=0;ir2<nInt;ir2++){
		for(unsigned int ir3=0;ir3<nInt;ir3++){
			const double r1 = Gauss[ir1][0];
			const double r2 = Gauss[ir2][0];
			const double r3 = Gauss[ir3][0];
			double detjac;
			ShapeFunc_Hex8(r1,r2,r3,coord_c,detjac,dndx,an_c);
			const double detwei = detjac*Gauss[ir1][1]*Gauss[ir2][1]*Gauss[ir3][1];
			vol += detwei;
			for(unsigned int ino=0;ino<nno;ino++){
			for(unsigned int jno=0;jno<nno;jno++){
				emat[ino][jno] += alpha*detwei*(dndx[ino][0]*dndx[jno][0]+dndx[ino][1]*dndx[jno][1]+dndx[ino][2]*dndx[jno][2]);
			}
			}
			// 要素節点等価外力ベクトルを積n分毎に足し合わせる
			for(unsigned int ino=0;ino<nno;ino++){
				eres_c[ino] += detwei*source*an_c[ino];
			}
		}
		}
		}

		// 要素節点等価内力ベクトルを求める
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			eres_c[ino] -= emat[ino][jno]*value_c[jno];
		}	
		}
		// 要素剛性行列にマージする
		mat_cc.Mearge(nno,no_c_va,nno,no_c_va,1,&emat[0][0]);
		// 残差ベクトルにマージする
		for(unsigned int ino=0;ino<nno;ino++){
			res_c.AddValue( no_c_va[ino],0,eres_c[ino]);
		}
	}
	return true;
//.........这里部分代码省略.........

bool AddLinSys_StVenant2D_NonStatic_BackwardEular_P1
(double dt,
 Fem::Eqn::ILinearSystem_Eqn& ls,
 double lambda, double myu,
 double rho, double g_x, double g_y,
 const unsigned int id_field_disp, const CFieldWorld& world, 
 const MatVec::CVector_Blk& velo_pre,
 bool is_initial,
 const unsigned int id_ea )
{
	//	std::cout << "St.Venant2D Triangle 3-point 1st order" << std::endl;
	
	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TRI );
	
	const CField& field_disp = world.GetField(id_field_disp);
	
	const CElemAry::CElemSeg& es_co = field_disp.GetElemSeg(id_ea,CORNER,false,world);
	const CElemAry::CElemSeg& es_va = field_disp.GetElemSeg(id_ea,CORNER,true,world);
	
	const unsigned int nno = 3;
	const unsigned int ndim = 2;
	
	unsigned int no[nno];	// elment node nuber 2 grobal node number
	
	double eKmat[nno][nno][ndim][ndim];	// element stiffness matrix
	double eMmat[nno][nno][ndim][ndim];	// element mass matrix		
	double eforce_in[nno][ndim];		// element residual vector		
	double eforce_ex[nno][ndim];		// element external force vector
	
	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(id_field_disp,CORNER,world);	// matrix
	CVector_Blk&     res_c = ls.GetResidual(id_field_disp,CORNER,world);// residual vector
	
	const CNodeAry::CNodeSeg& ns_c_val  = field_disp.GetNodeSeg(CORNER,true,world,VALUE);
	const CNodeAry::CNodeSeg& ns_c_velo = field_disp.GetNodeSeg(CORNER,true,world,VELOCITY);
	const CNodeAry::CNodeSeg& ns_c_co   = field_disp.GetNodeSeg(CORNER,false,world);
	assert( ns_c_val.Length() == ndim );
	assert( ns_c_co.Length()  == ndim );
	
	double g[2] = { g_x, g_y };
	
	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		es_co.GetNodes(ielem,no);	// get global node number of coordinate
		double ecoords[nno][ndim];	// element node coordinate
		for(unsigned int ino=0;ino<nno;ino++){
			ns_c_co.GetValue(  no[ino], ecoords[ino]);
		}
		
		es_va.GetNodes(ielem,no);	// get global node nubmer of value
		double edisp[nno][ndim];	// displacement 
		double evelo[nno][ndim];	// velocity
		for(unsigned int ino=0;ino<nno;ino++){
			ns_c_val.GetValue( no[ino], edisp[ino]);
			ns_c_velo.GetValue(no[ino], evelo[ino]);
		}
		
		////////////////////////////////
		
		// set 0 to elemnet stiffness matrix, element residual
		for(unsigned int i=0;i<nno*nno*ndim*ndim;i++){ *(&eKmat[0][0][0][0]+i) = 0.0; }
		for(unsigned int i=0;i<nno*nno*ndim*ndim;i++){ *(&eMmat[0][0][0][0]+i) = 0.0; }
		for(unsigned int i=0;i<         nno*ndim;i++){ *(&eforce_ex[0][0]  +i) = 0.0; }
		for(unsigned int i=0;i<         nno*ndim;i++){ *(&eforce_in[0][0]  +i) = 0.0; }
		
		const double area = TriArea(ecoords[0],ecoords[1],ecoords[2]);
		
		double dldx[nno][ndim];		// spacial derivative of shape function
		double zero_order_term[nno];	// constant term of shape function
		TriDlDx(dldx, zero_order_term,   ecoords[0], ecoords[1], ecoords[2]);
		
		{	// calc stiffness matrix
            double dudx[ndim][ndim] = { { 0.0, 0.0}, {0.0, 0.0} };
			for(unsigned int ino=0;ino<nno;ino++){
				for(unsigned int idim=0;idim<ndim;idim++){
					for(unsigned int jdim=0;jdim<ndim;jdim++){
						dudx[idim][jdim] += edisp[ino][idim]*dldx[ino][jdim];
					}
				}
			}
			AddElemMatFin_StVenant2D( area, myu,lambda, nno,dudx,dldx,  eKmat[0],eforce_in );
		}
		{	// calc mass matrix
			const double tmp1 = rho*area/12.0;
			for(unsigned int ino=0;ino<nno;ino++){
				for(unsigned int jno=0;jno<nno;jno++){
					for(unsigned int idim=0;idim<ndim;idim++){
						eMmat[ino][jno][idim][idim] += tmp1;
					}
				}
				for(unsigned int idim=0;idim<ndim;idim++){
					eMmat[ino][ino][idim][idim] += tmp1;
				}
			}
		}
		// call external force
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int idim=0;idim<ndim;idim++){
			eforce_ex[ino][idim] += area*rho*g[idim]*0.33333333333333333333;
//.........这里部分代码省略.........

bool AddLinSys_StVenant2D_Static_P1(
        Fem::Eqn::ILinearSystem_Eqn& ls,
		double lambda, double myu,
		double  rho, double g_x, double g_y,
		const unsigned int id_field_disp, const CFieldWorld& world, 
		const unsigned int id_ea)
{
	
//	std::cout << "St.Venant2D Static Tri3point" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TRI );
	const CField& field_disp = world.GetField(id_field_disp);

	const CElemAry::CElemSeg& es_co = field_disp.GetElemSeg(id_ea,CORNER,false,world);
	const CElemAry::CElemSeg& es_va = field_disp.GetElemSeg(id_ea,CORNER,true, world);

	const unsigned int nnoes = 3;
	const unsigned int ndim = 2;

	unsigned int noes[nnoes];	// 要素内の節点の節点番号

	double emat[nnoes][nnoes][ndim][ndim];	// 要素剛性行列
	double eforce_ex[nnoes][ndim];		// 要素内外力ベクトル
	double eforce_in[nnoes][ndim];		// 要素内内力ベクトル
	double eres[nnoes][ndim];		// 要素内残差ベクトル

	double ecoords[nnoes][ndim];		// 要素節点座標
	double edisp[  nnoes][ndim];		// 要素節点変位

	double dldx[nnoes][ndim];		// 形状関数の空間微分
	double zero_order_term[nnoes];	// 形状関数の定数項
				
	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(  id_field_disp,CORNER,world);	// 要素剛性行列(コーナ-コーナー)
	CVector_Blk&     res_c = ls.GetResidual(id_field_disp,CORNER,world);	// 要素残差ベクトル(コーナー)

	const CNodeAry::CNodeSeg& ns_c_val = field_disp.GetNodeSeg(CORNER,true,world,VALUE);
	const CNodeAry::CNodeSeg& ns_c_co  = field_disp.GetNodeSeg(CORNER,false,world);
/*
	unsigned int num_integral = 0;
	const unsigned int nInt = NIntTriGauss[num_integral];
	const double (*Gauss)[3] = TriGauss[num_integral];
*/
	double g[2] = { g_x, g_y };

	for(unsigned int ielem=0;ielem<ea.Size();ielem++){
		// 要素の節点番号を取ってくる
		es_co.GetNodes(ielem,noes);
		// 節点の座標、値を取ってくる
		for(unsigned int ino=0;ino<nnoes;ino++){
			ns_c_co.GetValue(noes[ino],ecoords[ino]);
		}
		// 要素の節点番号を取ってくる
		es_va.GetNodes(ielem,noes);
		// 節点の座標、値を取ってくる
		for(unsigned int ino=0;ino<nnoes;ino++){
			ns_c_val.GetValue(noes[ino],edisp[ino]);
		}

		////////////////////////////////

		// 要素剛性行列、残差を０で初期化
		for(unsigned int i=0;i<nnoes*nnoes*ndim*ndim;i++){ *(&emat[0][0][0][0]+i) = 0.0; }
		for(unsigned int i=0;i<           nnoes*ndim;i++){ *(&eforce_ex[0][0] +i) = 0.0; }
		for(unsigned int i=0;i<           nnoes*ndim;i++){ *(&eforce_in[0][0] +i) = 0.0; }

		// 面積を求める
		const double area = TriArea(ecoords[0],ecoords[1],ecoords[2]);

		// 形状関数のｘｙ微分を求める
		TriDlDx(dldx, zero_order_term,   ecoords[0], ecoords[1], ecoords[2]);

        double dudx[ndim][ndim] = { { 0.0, 0.0}, {0.0, 0.0} };
		for(unsigned int ino=0;ino<nnoes;ino++){
			for(unsigned int idim=0;idim<ndim;idim++){
			for(unsigned int jdim=0;jdim<ndim;jdim++){
				dudx[idim][jdim] += edisp[ino][idim]*dldx[ino][jdim];
			}
			}
		}

		AddElemMatFin_StVenant2D( area, myu,lambda, nnoes,dudx,dldx,  emat[0],eforce_in );

		// 外力ベクトルを求める
		for(unsigned int ino=0;ino<nnoes;ino++){
		for(unsigned int idim=0;idim<ndim;idim++){
			eforce_ex[ino][idim] += area*rho*g[idim]/nnoes;
		}
		}

		////////////////////////////////

		// 要素内残差ベクトルを求める
		for(unsigned int ino=0;ino<nnoes;ino++){
		for(unsigned int idim=0;idim<ndim;idim++){
			eres[ino][idim] = eforce_ex[ino][idim] - eforce_in[ino][idim];
		}
		}

//.........这里部分代码省略.........

bool AddLinSys_StVenant3D_NonStatic_NewmarkBeta_P1
(double gamma, double beta, double dt,
 Fem::Eqn::ILinearSystem_Eqn& ls,
 double lambda, double myu,
 double rho, double g_x, double g_y, double g_z,
 const unsigned int id_field_disp, const CFieldWorld& world, 
 bool is_initial,
 const unsigned int id_ea)
{
//	std::cout << "St.Venant3D Tet1st" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TET );

	const CField& field_disp = world.GetField(id_field_disp);

	const CElemAry::CElemSeg& es_c = field_disp.GetElemSeg(id_ea,CORNER,true,world);

	const unsigned int nnoes = 4;
	const unsigned int ndim = 3;

	double emat[nnoes][nnoes][ndim][ndim];	// coefficient element matrix
	double eMmat[nnoes][nnoes][ndim][ndim];	// mass element matrix
	double eKmat[nnoes][nnoes][ndim][ndim];	// stiffness element matrix

	double eforce_in[nnoes][ndim];		// 要素内内力ベクトル
	double eres[nnoes][ndim];		// 要素内残差ベクトル

	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(  id_field_disp,CORNER,world);// 要素剛性行列(コーナ-コーナー)
	CVector_Blk&     res_c = ls.GetResidual(id_field_disp,CORNER,world);// 要素残差ベクトル(コーナー)

	const CNodeAry::CNodeSeg& ns_c_val  = field_disp.GetNodeSeg(CORNER,true,world,VALUE);
	const CNodeAry::CNodeSeg& ns_c_velo = field_disp.GetNodeSeg(CORNER,true,world,VELOCITY);
	const CNodeAry::CNodeSeg& ns_c_acc  = field_disp.GetNodeSeg(CORNER,true,world,ACCELERATION);
	const CNodeAry::CNodeSeg& ns_c_co  = field_disp.GetNodeSeg(CORNER,false,world);
	assert( ns_c_val.Length() == ndim );
	assert( ns_c_co.Length()  == ndim );

	const double g[ndim] = { g_x, g_y, g_z };

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
    {
	    unsigned int noes[nnoes];
		es_c.GetNodes(ielem,noes);
	    double ecoords[nnoes][ndim];		// coordinate 
	    double edisp[  nnoes][ndim];		// displacement 
	    double evelo[  nnoes][ndim];		// velocity
	    double eacc[  nnoes][ndim];			// acceleration
		for(unsigned int ino=0;ino<nnoes;ino++){
			ns_c_co.GetValue(  noes[ino], ecoords[ino]);
			ns_c_val.GetValue( noes[ino], edisp[ino]);
			ns_c_velo.GetValue(noes[ino], evelo[ino]);
			ns_c_acc.GetValue( noes[ino], eacc[ ino]);
		}

		////////////////////////////////

		const double vol = TetVolume(ecoords[0],ecoords[1],ecoords[2],ecoords[3]);

	    double dldx[nnoes][ndim];		// derivative of shape function
	    double zero_order_term[nnoes];	// constant term of shape function
		TetDlDx(dldx, zero_order_term,  ecoords[0],ecoords[1],ecoords[2],ecoords[3]);
		
		{	// calc stiffness matrix
			double dudx[ndim][ndim] = { {0.0,0.0,0.0}, {0.0,0.0,0.0}, {0.0,0.0,0.0} };
			for(unsigned int ino=0;ino<nnoes;ino++){
			for(unsigned int idim=0;idim<ndim;idim++){
			for(unsigned int jdim=0;jdim<ndim;jdim++){
				dudx[idim][jdim] += edisp[ino][idim]*dldx[ino][jdim];
			}
			}
			}
			SetElemMatFin_StVenant3D( vol, myu,lambda, nnoes,dudx,dldx,  eKmat[0],eforce_in );
		}

		{	// calc mass matrix
		    for(unsigned int i=0;i<nnoes*nnoes*ndim*ndim;i++){ *(&eMmat[0][0][0][0]+i) = 0.0; }
			const double tmp1 = rho*vol*0.05;
			for(unsigned int ino=0;ino<nnoes;ino++){
			for(unsigned int jno=0;jno<nnoes;jno++){
				eMmat[ino][jno][0][0] += tmp1;
				eMmat[ino][jno][1][1] += tmp1;
				eMmat[ino][jno][2][2] += tmp1;
            }
            }
			for(unsigned int ino=0;ino<nnoes;ino++){
				eMmat[ino][ino][0][0] += tmp1;
				eMmat[ino][ino][1][1] += tmp1;
				eMmat[ino][ino][2][2] += tmp1;
			}
		}

		////////////////////////////////

		for(unsigned int i=0;i<nnoes*nnoes*ndim*ndim;i++){
			(&emat[0][0][0][0])[i] = (&eMmat[0][0][0][0])[i] + beta*dt*dt*(&eKmat[0][0][0][0])[i];
		}
		for(unsigned int ino=0;ino<nnoes;ino++){
		for(unsigned int idim=0;idim<ndim;idim++){
//.........这里部分代码省略.........

static bool AddLinearSystem_Wave2D_P1(
		double rho, double alpha, double source, 
		CLinearSystem_Eigen& ls, 
		const unsigned int id_field_val, const CFieldWorld& world,
		const unsigned int id_ea )
{
//	std::cout << "Wave2D TriP1" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TRI );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

	const CElemAry::CElemSeg& es_c_co = field_val.GetElemSeg(id_ea,CORNER,false,world);

	const unsigned int nno = 3;
	const unsigned int ndim = 2;

	unsigned int no_c[nno];	// 要素節点の全体節点番号
	double coord_c[nno][ndim];	// 要素節点の座標
				
	double eKmat[nno][nno];	// 要素剛性行列
	double eMmat[nno];	// 要素節点等価内力、外力、残差ベクトル
				
	CMatDia_BlkCrs& Kmat_cc = ls.GetMatrix(       id_field_val,CORNER,world);
	CDiaMat_Blk&    Mmat_cc = ls.GetDiaMassMatrix(id_field_val,CORNER,world);

	const CNodeAry::CNodeSeg& ns_c_co  = field_val.GetNodeSeg(CORNER,false,world);

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		// 要素配列から要素セグメントの節点番号を取り出す
		es_c_co.GetNodes(ielem,no_c);
		for(unsigned int inoes=0;inoes<nno;inoes++){
			ns_c_co.GetValue(no_c[inoes],coord_c[inoes]);
		}
		// 面積を求める
		const double area = TriArea(coord_c[0],coord_c[1],coord_c[2]);
		// 形状関数の微分を求める
		double dldx[nno][ndim];	// 形状関数のxy微分
		double const_term[nno];	// 形状関数の定数項
		TriDlDx(dldx,const_term,coord_c[0],coord_c[1],coord_c[2]);
		// 要素剛性行列を作る
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			eKmat[ino][jno] = alpha*area*(dldx[ino][0]*dldx[jno][0]+dldx[ino][1]*dldx[jno][1]);
		}
		}
		// 要素節点等価外力ベクトルを求める
		for(unsigned int ino=0;ino<nno;ino++){
			eMmat[ino] = rho*area*0.33333333333333333;
		}
		// 要素剛性行列にマージする
		Kmat_cc.Mearge(nno,no_c,nno,no_c,1,&eKmat[0][0]);
		// 残差ベクトルにマージする
		for(unsigned int ino=0;ino<nno;ino++){
			Mmat_cc.Mearge(no_c[ino],1,&eMmat[ino]);
		}
	}
	return true;
}

bool AddLinearSystem_Poisson2D_P1(
		double alpha, double source, 
		CLinearSystem_Save& ls, 
		unsigned int id_field_val, const CFieldWorld& world,
		unsigned int id_ea)
{
//	std::cout << "Poisson2D Tri P1" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TRI );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

	const CElemAry::CElemSeg& es_c_va = field_val.GetElemSeg(id_ea,CORNER,true, world);
	const CElemAry::CElemSeg& es_c_co = field_val.GetElemSeg(id_ea,CORNER,false,world);

	const unsigned int nno = 3;
	const unsigned int ndim = 2;

	unsigned int no_c[nno];	// 要素節点の全体節点番号

	double value_c[nno];		// 要素節点の値
	double coord_c[nno][ndim];	// 要素節点の座標
				
	double emat[nno][nno];	// 要素剛性行列
	double eqf_out_c[nno];	// 要素節点等価内力、外力、残差ベクトル
				
	CMatDia_BlkCrs&      mat_cc  = ls.GetMatrix(id_field_val,CORNER,world);
	CVector_Blk&         force_c = ls.GetForce(id_field_val,CORNER,world);
	CMat_BlkCrs& mat_cc_boundary = ls.GetMatrix_Boundary(id_field_val,CORNER,id_field_val,CORNER,world);

	const CNodeAry::CNodeSeg& ns_c_val = field_val.GetNodeSeg(CORNER,true, world);//na_c_val.GetSeg(id_ns_c_val);
	const CNodeAry::CNodeSeg& ns_c_co  = field_val.GetNodeSeg(CORNER,false,world);//na_c_co.GetSeg(id_ns_c_co);

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		// 要素配列から要素セグメントの節点番号を取り出す
		es_c_co.GetNodes(ielem,no_c);
		// 節点の値を取って来る
		for(unsigned int ino=0;ino<nno;ino++){
			ns_c_co.GetValue(no_c[ino],coord_c[ino]);
		}
		es_c_va.GetNodes(ielem,no_c);
		// 節点の値を取って来る
		for(unsigned int ino=0;ino<nno;ino++){
			ns_c_val.GetValue(no_c[ino],&value_c[ino]);
		}

		// 面積を求める
		const double area = TriArea(coord_c[0],coord_c[1],coord_c[2]);
		// 形状関数の微分を求める
		double dldx[nno][ndim];	// 形状関数のxy微分
		double const_term[nno];	// 形状関数の定数項
		TriDlDx(dldx,const_term,coord_c[0],coord_c[1],coord_c[2]);
		// 要素剛性行列を作る
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			emat[ino][jno] = alpha*area*(dldx[ino][0]*dldx[jno][0]+dldx[ino][1]*dldx[jno][1]);
		}
		}
		// 要素節点等価外力ベクトルを求める
		for(unsigned int ino=0;ino<nno;ino++){
			eqf_out_c[ino] = source*area*0.333333333333333333;
		}
		// 要素剛性行列にマージする
		mat_cc.Mearge(         nno,no_c,nno,no_c,1,&emat[0][0]);
		mat_cc_boundary.Mearge(nno,no_c,nno,no_c,1,&emat[0][0]);
		// 残差ベクトルにマージする
		for(unsigned int ino=0;ino<nno;ino++){
			force_c.AddValue( no_c[ino],0,eqf_out_c[ino]);
		}
	}
	return true;
}

bool AddLinSys_StVenant3D_Static_Q1
(Fem::Eqn::ILinearSystem_Eqn& ls,
 double lambda, double myu,
 double  rho, double g_x, double g_y, double g_z,
 const unsigned int id_field_disp, const CFieldWorld& world, 
 const unsigned int	id_ea)
{
//	std::cout << "StVenant3D Hex 8-point 1st order" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == HEX );

	assert( world.IsIdField(id_field_disp) );
	const CField& field_disp = world.GetField(id_field_disp);

	const CElemAry::CElemSeg& es_c = field_disp.GetElemSeg(id_ea,CORNER,true,world);

	unsigned int num_integral = 1;
	const unsigned int nInt = NIntLineGauss[num_integral];
	const double (*Gauss)[2] = LineGauss[num_integral];
	double detjac, detwei;

	const unsigned int nnoes = 8;	assert( nnoes == es_c.Length() );
	const unsigned int ndim = 3;

	unsigned int noes[nnoes];	// 要素内の節点の節点番号

	double emat[nnoes][nnoes][ndim][ndim];	// 要素剛性行列
	double eforce_ex[nnoes][ndim];		// 要素内外力ベクトル
	double eforce_in[nnoes][ndim];		// 要素内内力ベクトル
	double eres[nnoes][ndim];		// 要素内残差ベクトル

	double ecoords[nnoes][ndim];		// 要素節点座標
	double edisp[  nnoes][ndim];		// 要素節点変位

	double dndx[nnoes][ndim];		// 形状関数の空間微分
	double an[nnoes];
				
	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(  id_field_disp,CORNER,world);// 要素剛性行列(コーナ-コーナー)
	CVector_Blk&     res_c = ls.GetResidual(id_field_disp,CORNER,world);// 要素残差ベクトル(コーナー)

	const CNodeAry::CNodeSeg& ns_c_val = field_disp.GetNodeSeg(CORNER,true,world,VALUE);//.GetSeg(id_ns_c_val);
	const CNodeAry::CNodeSeg& ns_c_co  = field_disp.GetNodeSeg(CORNER,false,world,VALUE);//na_c_co.GetSeg(id_ns_c_co);

//	double g[ndim] = { g_x, g_y, g_z };

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		// 要素の節点番号を取ってくる
		es_c.GetNodes(ielem,noes);
		// 節点の座標、値を取ってくる
		for(unsigned int ino=0;ino<nnoes;ino++){
			ns_c_co.GetValue(noes[ino],ecoords[ino]);
			ns_c_val.GetValue(noes[ino],edisp[ino]);
		}

		////////////////////////////////

		// 要素剛性行列、残差を０で初期化
		for(unsigned int i=0;i<nnoes*nnoes*ndim*ndim;i++){ *(&emat[0][0][0][0]+i) = 0.0; }
		for(unsigned int i=0;i<           nnoes*ndim;i++){ *(&eforce_ex[0][0] +i) = 0.0; }
		for(unsigned int i=0;i<           nnoes*ndim;i++){ *(&eforce_in[0][0] +i) = 0.0; }

		double vol = 0.0;
		for(unsigned int ir1=0;ir1<nInt;ir1++){
		for(unsigned int ir2=0;ir2<nInt;ir2++){
		for(unsigned int ir3=0;ir3<nInt;ir3++){
			const double r1 = Gauss[ir1][0];
			const double r2 = Gauss[ir2][0];
			const double r3 = Gauss[ir3][0];
			ShapeFunc_Hex8(r1,r2,r3,ecoords,detjac,dndx,an);
			detwei = detjac*Gauss[ir1][1]*Gauss[ir2][1]*Gauss[ir3][1];
			vol += detwei;
			{	// 要素剛性行列，要素内力ベクトルを作る
				double dudx[ndim][ndim] = { {0.0,0.0,0.0}, {0.0,0.0,0.0}, {0.0,0.0,0.0} };
				for(unsigned int ino=0;ino<nnoes;ino++){
					for(unsigned int idim=0;idim<ndim;idim++){
					for(unsigned int jdim=0;jdim<ndim;jdim++){
						dudx[idim][jdim] += edisp[ino][idim]*dndx[ino][jdim];
					}
					}
				}	
				SetElemMatFin_StVenant3D( vol,  myu,lambda,  nnoes,dudx,dndx,  emat[0],eforce_in );
			}
			// 要素節点等価外力ベクトルを積n分毎に足し合わせる
			for(unsigned int ino=0;ino<nnoes;ino++){
				eforce_ex[ino][0] += detwei*rho*g_x;
				eforce_ex[ino][1] += detwei*rho*g_y;
				eforce_ex[ino][2] += detwei*rho*g_z;
			}
		}
		}
		}

		////////////////////////////////

		// 要素内残差ベクトルを求める
		for(unsigned int ino=0;ino<nnoes;ino++){
		for(unsigned int idim=0;idim<ndim;idim++){
//.........这里部分代码省略.........

// 対角にパターンを追加
bool CZLinearSystem::AddPattern_Field(const unsigned int id_field, const CFieldWorld& world)
{
	if( !world.IsIdField(id_field) ) return false;
	const CField& field = world.GetField(id_field);
	unsigned int id_field_parent;
	{
		if( field.GetIDFieldParent() == 0 ){ id_field_parent = id_field; }
		else{ id_field_parent = field.GetIDFieldParent(); }
	}

	unsigned int nlen_value;
	{	// 複素数の場合は２で割る
		const unsigned int nlen = field.GetNLenValue();
		assert( nlen % 2 == 0 );
		nlen_value = nlen / 2;
	}

	int ESType2iLS[3];
	{	// Bubbleブロックを作る
		unsigned int id_na_val = field.GetNodeSegInNodeAry(BUBBLE).id_na_va;
		if( id_na_val != 0 ){
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			CLinSysSeg seg;
			{
				seg.id_field = id_field_parent; seg.node_config = BUBBLE;
				seg.len=nlen_value; seg.nnode=na.Size();
			}
			ESType2iLS[2] = this->AddLinSysSeg(seg);
		}
		else{ ESType2iLS[2] = -1; }
	}
	{	// Edgeブロックを作る
		unsigned int id_na_val = field.GetNodeSegInNodeAry(EDGE).id_na_va;
		if( id_na_val != 0 ){
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			CLinSysSeg seg;
			{
				seg.id_field = id_field_parent; seg.node_config = EDGE;
				seg.len=nlen_value; seg.nnode=na.Size();
			}
			ESType2iLS[1] = this->AddLinSysSeg(seg);
		}
		else{ ESType2iLS[1] = -1; }
	}
	{	// Cornerブロックを作る
		unsigned int id_na_val = field.GetNodeSegInNodeAry(CORNER).id_na_va;
		if( id_na_val != 0 ){
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			CLinSysSeg seg;
			{
				seg.id_field = id_field_parent; seg.node_config = CORNER;
				seg.len=nlen_value; seg.nnode=na.Size();
			}
			ESType2iLS[0] = this->AddLinSysSeg(seg);
		}
		else{ ESType2iLS[0] = -1; }
	}
	////////////////////////////////
	const std::vector<unsigned int> aIdEA = field.GetAryIdEA();
	if( aIdEA.size() == 0 ){ // 剛体モードのための行列
		unsigned int ils0 = ESType2iLS[0];
		if( m_Matrix_Dia[ils0] == 0 ){
			m_Matrix_Dia[ils0] = new CZMatDia_BlkCrs(1, nlen_value);
		}
		return true;
	}

	for(unsigned int iiea=0;iiea<aIdEA.size();iiea++)
	{
		const unsigned int id_ea = aIdEA[iiea];
		const CElemAry& ea = world.GetEA(id_ea);
		// CORNER節点について
		if( field.GetIdElemSeg(id_ea,CORNER,true,world) != 0 ){
			assert( world.IsIdEA(id_ea) );
			const unsigned int id_es_c = field.GetIdElemSeg(id_ea,CORNER,true,world);
			assert( ea.IsSegID(id_es_c) );
			const unsigned int ils0 = ESType2iLS[0];
			this->AddMat_Dia(ils0, ea, id_es_c );			// cc行列を作る
			if( field.GetIdElemSeg(id_ea,BUBBLE,true,world) != 0 ){	// CORNER-BUBBLE
				const unsigned int id_es_b = field.GetIdElemSeg(id_ea,BUBBLE,true,world);
				assert( ea.IsSegID(id_es_b) );
				const unsigned int ils1 = ESType2iLS[2];
				Com::CIndexedArray crs;
				ea.MakePattern_FEM(id_es_c,id_es_b,crs);
				assert( crs.CheckValid() );
				this->AddMat_NonDia(ils0,ils1, crs);		// cb行列を作る
				const unsigned int nnode1 = m_aSeg[ils1].nnode;
				Com::CIndexedArray crs_inv;
				crs_inv.SetTranspose(nnode1,crs);
				this->AddMat_NonDia(ils1,ils0, crs_inv);	// bc行列を作る
			}
			if( field.GetIdElemSeg(id_ea,EDGE,true,world) != 0 ){	// CONRER-EDGE
				const unsigned int id_es_e = field.GetIdElemSeg(id_ea,EDGE,true,world);
				assert( ea.IsSegID(id_es_e) );
				const unsigned int ils1 = ESType2iLS[1];
				Com::CIndexedArray crs;
//.........这里部分代码省略.........

static bool AddLinearSystem_Poisson3D_P1(
		double alpha, double source, 
        Fem::Ls::CLinearSystem_Field& ls, 
		unsigned int id_field_val, const CFieldWorld& world,
		unsigned int id_ea)
{
//	std::cout << "Poisson3D Tet 4-point 1st order" << std::endl;

	assert( world.IsIdEA(id_ea) );
	const CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == TET );

	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

	const CElemAry::CElemSeg& es_c = field_val.GetElemSeg(id_ea,CORNER,false,world);

	const unsigned int nno = 4;
	const unsigned int ndim = 3;

	unsigned int no_c[nno];	// 要素節点の全体節点番号

	double value_c[nno];		// 要素節点の値
	double coord_c[nno][ndim];	// 要素節点の座標
			
	double emat[nno][nno];	// 要素剛性行列
	double eres_c[nno];	// 要素節点等価内力、外力、残差ベクトル
				
	CMatDia_BlkCrs& mat_cc = ls.GetMatrix(  id_field_val,CORNER,world);
	CVector_Blk&    res_c  = ls.GetResidual(id_field_val,CORNER,world);

	const CNodeAry::CNodeSeg& ns_c_val = field_val.GetNodeSeg(CORNER,true,world);
	const CNodeAry::CNodeSeg& ns_c_co  = field_val.GetNodeSeg(CORNER,false,world);

	for(unsigned int ielem=0;ielem<ea.Size();ielem++)
	{
		// 要素配列から要素セグメントの節点番号を取り出す
		es_c.GetNodes(ielem,no_c);
		// 節点の値を取って来る
		for(unsigned int inoes=0;inoes<nno;inoes++){
			ns_c_co.GetValue(no_c[inoes],coord_c[inoes]);
			ns_c_val.GetValue(no_c[inoes],&value_c[inoes]);
		}

		// 面積を求める
		const double vol = TetVolume(coord_c[0],coord_c[1],coord_c[2], coord_c[3]);
		// 形状関数の微分を求める	
		double dldx[nno][ndim];	// 形状関数のxy微分
		double const_term[nno];	// 形状関数の定数項
		TetDlDx(dldx,const_term,coord_c[0],coord_c[1],coord_c[2], coord_c[3]);
		// 要素剛性行列を作る
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			emat[ino][jno] = alpha*vol*(dldx[ino][0]*dldx[jno][0]+dldx[ino][1]*dldx[jno][1]+dldx[ino][2]*dldx[jno][2]);
		}
		}
		// 要素節点等価外力ベクトルを求める
		for(unsigned int ino=0;ino<nno;ino++){
			eres_c[ino] = source*vol*0.25;
		}
		// 要素節点等価内力ベクトルを求める
		for(unsigned int ino=0;ino<nno;ino++){
		for(unsigned int jno=0;jno<nno;jno++){
			eres_c[ino] -= emat[ino][jno]*value_c[jno];
		}	
		}
		// 要素剛性行列にマージする
		mat_cc.Mearge(nno,no_c,nno,no_c,1,&emat[0][0]);
		// 残差ベクトルにマージする
		for(unsigned int inoes=0;inoes<nno;inoes++){
			res_c.AddValue( no_c[inoes],0,eres_c[inoes]);
		}
	}
	return true;
}