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

/*
bool CEqnSystem_Solid2D::SetDomain_Field(unsigned int id_field_base, Fem::Field::CFieldWorld& world)
{
	{	// ì¸óÕÉtÉBÅ[ÉãÉhÇÃç¿ïWêflì_ÉZÉOÉÅÉìÉgÇÃdofÇ™2ǩǫǧǩÉ`ÉFÉbÉNÇ∑ÇÈ
//		unsigned int id_field_base = world.GetFieldBaseID();
		assert( world.IsIdField(id_field_base) );
		const CField& field_base = world.GetField(id_field_base);
		assert( field_base.GetNDimCoord() == 2 );
	}
	m_IdFieldDisp = world.MakeField_FieldElemDim(id_field_base, 2,
        VECTOR2,VALUE|VELOCITY|ACCELERATION,CORNER);
	const CIDConvEAMshCad conv = world.GetIDConverter(id_field_base);
	{
		m_aEqn.clear();
		const CField& field = world.GetField(m_IdFieldDisp);
		const std::vector<unsigned int>& aIdEA = field.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			assert( world.IsIdEA(id_ea) );
			const CElemAry& ea = world.GetEA(id_ea);
			if( ea.ElemType() != TRI && ea.ElemType() != QUAD ) continue;
			m_aEqn.push_back( CEqn_Solid2D(id_ea,m_IdFieldDisp) );
		}
	}
	this->ClearLinearSystemPreconditioner();
	this->m_aIdFixField.clear();
	return true;
}
*/
bool CEqnSystem_Solid2D::SetDomain_FieldEA(unsigned int id_field_base, unsigned int id_ea, 
                                           Fem::Field::CFieldWorld& world)
{
	{	// ì¸óÕÉtÉBÅ[ÉãÉhÇÃç¿ïWêflì_ÉZÉOÉÅÉìÉgÇÃdofÇ™2ǩǫǧǩÉ`ÉFÉbÉNÇ∑ÇÈ
//		unsigned int id_field_base = world.GetFieldBaseID();
		assert( world.IsIdField(id_field_base) );
		const CField& field_base = world.GetField(id_field_base);
		assert( field_base.GetNDimCoord() == 2 );
	}
	m_IdFieldDisp = world.MakeField_FieldElemAry(id_field_base, id_ea,
        VECTOR2,VALUE|VELOCITY|ACCELERATION,CORNER);
	{
		m_aEqn.clear();
		const CField& field = world.GetField(m_IdFieldDisp);
		const std::vector<unsigned int>& aIdEA = field.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			assert( world.IsIdEA(id_ea) );
			const CElemAry& ea = world.GetEA(id_ea);
			if( ea.ElemType() != TRI && ea.ElemType() != QUAD ) continue;
			m_aEqn.push_back( CEqn_Solid2D(id_ea,m_IdFieldDisp) );
		}
	}
	this->ClearLinearSystemPreconditioner();
	this->m_aIdFixField.clear();
	return true;
}

double CEqnSystem_Fluid2D::MakeLinearSystem(const Fem::Field::CFieldWorld& world)
{
    if( pLS==0 || pPrec==0 ) this->InitializeLinearSystem(world);
    // 連立一次方程式を作る
    pLS->InitializeMarge();	// 連立一次方程式を初期化する(0クリア)

    if( m_IsStationary ) {
        for(unsigned int ieqn=0; ieqn<m_aEqn.size(); ieqn++) {
            m_aEqn[ieqn].AddLinSys(*pLS,world);
        }
    }
    else {
        for(unsigned int ieqn=0; ieqn<m_aEqn.size(); ieqn++) {
            //            std::cout << "CEqnSystem_Fluid2D::MakeLinearSystem " << ieqn << " " << m_aEqn[ieqn].GetIdEA() << std::endl;
            m_aEqn[ieqn].AddLinSys_NewmarkBetaAPrime(m_dt,m_gamma_newmark,m_beta_newmark,true,
                    *pLS, world);
        }
    }
    if( world.IsIdField( this->m_id_force ) )
    {
        assert( world.IsIdField(this->m_id_force) );
        const Fem::Field::CField& ff = world.GetField(this->m_id_force);
        const Fem::Field::CNodeAry::CNodeSeg& nsf_v = ff.GetNodeSeg(Fem::Field::CORNER,true, world,Fem::Field::VELOCITY);
        const unsigned int nno = nsf_v.Size();
        assert( nsf_v.Length() == 2 );
        ////////////////
        assert( world.IsIdField(this->m_id_velo) );
        const Fem::Field::CField& fv = world.GetField(this->m_id_velo);
        const Fem::Field::CNodeAry::CNodeSeg& nsv_v = fv.GetNodeSeg(Fem::Field::CORNER,true, world,Fem::Field::VELOCITY);
        assert( nsv_v.Size() == nno );
        assert( nsv_v.Length() == 2 );
        ////////////////
        MatVec::CVector_Blk& vec_res = pLS->GetResidual(this->m_id_velo,Fem::Field::CORNER,world);
        assert( vec_res.NBlk() == nno );
        assert( vec_res.Len() >= 2 );    // Combineのときは2以上
        double nres_add = 0;
        for(unsigned int ino=0; ino<nno; ino++)
        {
            double force[2];
            nsf_v.GetValue(ino,force);
            vec_res.AddValue(ino,0,force[0]);
            vec_res.AddValue(ino,1,force[1]);
            nres_add += force[0]*force[0]+force[1]*force[1];
        }
        //        std::cout << res << " " << nres_add << std::endl;
    }
    const double norm_res = pLS->FinalizeMarge();
    pPrec->SetValue((*pLS).m_ls);
    return norm_res;
}

bool Fem::Eqn::AddLinSys_AdvectionDiffusion_NonStatic_Newmark(
	Fem::Ls::CLinearSystem_SaveDiaM_Newmark& ls,
	double rho, double myu, double source,
	const Fem::Field::CFieldWorld& world,
	unsigned int id_field_val, unsigned int id_field_velo, 
	unsigned int id_ea )
{
	if( !world.IsIdField(id_field_val) ) return false;
	const CField& field_val = world.GetField(id_field_val);

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

	if( field_val.GetFieldType() != SCALAR ) return false;
	if( field_velo.GetFieldType() != VECTOR2 ) return false;

	if( id_ea != 0 ){
		Fem::Field::INTERPOLATION_TYPE intp_type_val  = field_val.GetInterpolationType(id_ea,world);
		Fem::Field::INTERPOLATION_TYPE intp_type_velo = field_velo.GetInterpolationType(id_ea,world);
		if( intp_type_val == TRI11 && intp_type_velo == TRI11 ){
			AddLinSys_AdvectionDiffusion_NonStatic_Newmark_P1P1(
				rho,myu,source,
				ls,
				id_field_val,id_field_velo,world,
				id_ea);
		}
		else{
			assert(0);
		}
	}
	else{
		const std::vector<unsigned int> aIdEA = field_val.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			bool res = Fem::Eqn::AddLinSys_AdvectionDiffusion_NonStatic_Newmark(
				ls,
				rho, myu, source,
				world,
				id_field_val, id_field_velo, 
				id_ea );
			if( !res ) return false;
		}
		return true;
	}



	return true;
}

void CDrawerFace::EnableUVMap(bool is_uv_map, const Fem::Field::CFieldWorld& world)
{
  if( (pUVArray != 0 ) == is_uv_map ){ return; }
  if( is_uv_map ){
    const unsigned int nnode = this->m_vertex_ary.NPoin();
//    const unsigned int ndim  = this->m_vertex_ary.NDim();
    delete[] pUVArray;
    pUVArray = new double [nnode*2];
    if( !world.IsIdField(m_id_field) ) return;
    const Fem::Field::CField& field = world.GetField(m_id_field);
//    unsigned int id_na_c_co = field.GetNodeSegInNodeAry(CORNER).id_na_co;
//    assert( world.IsIdNA(id_na_c_co) );
//    const Fem::Field::CNodeAry& na_c_co = world.GetNA(id_na_c_co);
//    assert( field.IsNodeSeg(CORNER,false,world) );
    const Fem::Field::CNodeAry::CNodeSeg& ns_c_co = field.GetNodeSeg(CORNER,false,world);
//    const unsigned int ndim = ns_c_co.Length();      
//    assert( ndim >= 2 );
    for(unsigned int ino=0;ino<ns_c_co.Size();ino++){
      double c[3]; ns_c_co.GetValue(ino,c);
      pUVArray[ino*2+0] = c[0]*tex_scale;
      pUVArray[ino*2+1] = c[1]*tex_scale;
    }
  }
  else{
    delete[] pUVArray;
    pUVArray = 0;
  }  
}

// add pattern to the diagonal
bool Ls::CLinearSystem_RigidField2::AddPattern_Field(const unsigned int id_field, const Fem::Field::CFieldWorld& world)
{
	if( !world.IsIdField(id_field) ) return false;
  const Fem::Field::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(); }
	}

	const unsigned int nlen_value = field.GetNLenValue();

	int ilss_add;
	{	// Cornerブロックを作る
    unsigned int id_na_val = field.GetNodeSegInNodeAry(Fem::Field::CORNER).id_na_va;
		if( id_na_val != 0 ){
			assert( world.IsIdNA(id_na_val) );
      const Fem::Field::CNodeAry& na = world.GetNA(id_na_val);
      assert( m_ls.GetNLinSysSeg() == this->m_aSegRF.size() );
      ilss_add = m_ls.GetNLinSysSeg();
      this->m_aSegRF.push_back( CLinSysSegRF(id_field_parent,Fem::Field::CORNER) );
      int ires = m_ls.AddLinSysSeg( na.Size(), field.GetNLenValue() );
      assert( ires == ilss_add );
		}
		else{ ilss_add = -1; }
	}
	////////////////////////////////
	const std::vector<unsigned int> aIdEA = field.GetAryIdEA();
  assert( aIdEA.size() > 0 );
	for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
		const unsigned int id_ea = aIdEA[iiea];
    const Fem::Field::CElemAry& ea = world.GetEA(id_ea);
		// CORNER節点について
    if( field.GetIdElemSeg(id_ea,Fem::Field::CORNER,true,world) != 0 ){
			assert( world.IsIdEA(id_ea) );
      const unsigned int id_es_c = field.GetIdElemSeg(id_ea,Fem::Field::CORNER,true,world);
			assert( ea.IsSegID(id_es_c) );
      {
        Com::CIndexedArray crs;
        ea.MakePattern_FEM(id_es_c,crs);
        m_ls.AddMat_Dia(ilss_add, crs );			// cc行列を作る
      }
      if( field.GetIdElemSeg(id_ea,Fem::Field::BUBBLE,true,world) != 0 ){	// CORNER-BUBBLE
        assert(0);
			}
      if( field.GetIdElemSeg(id_ea,Fem::Field::EDGE,true,world) != 0 ){	// CONRER-EDGE
        assert(0);
			}
		}
		// EDGE節点について
    if( field.GetIdElemSeg(id_ea,Fem::Field::EDGE,true,world) != 0 ){
      assert(0);
		}
		// BUBBLE節点について
    if( field.GetIdElemSeg(id_ea,Fem::Field::BUBBLE,true,world) != 0 ){
      assert(0);
		}
	}
	return true;
}

bool Fem::Field::CFieldValueSetter::ExecuteValue
(double cur_time, Fem::Field::CFieldWorld& world)
{
  if( id_field_gradient_ != 0 ){
    SetFieldValue_Gradient(id_field_, world, id_field_gradient_);
  }
  if( !world.IsIdField(id_field_) ){ return false; }  
  Fem::Field::CField& field = world.GetField(id_field_);
  const unsigned int nlen = field.GetNLenValue();
  for(unsigned int ilen=0;ilen<nlen;ilen++){
    if(      aValueFieldDof_[ilen+nlen*0].itype == 1 ){
      SetFieldValue_Constant(id_field_,ilen,VALUE,world,aValueFieldDof_[ilen+nlen*0].val);      
    }
    else if( aValueFieldDof_[ilen+nlen*0].itype == 2 ){
      SetFieldValue_MathExp(id_field_,ilen,VALUE,world,aValueFieldDof_[ilen+nlen*0].math_exp,cur_time);
    }
  }
  for(unsigned int ilen=0;ilen<nlen;ilen++){
    if(      aValueFieldDof_[ilen+nlen*1].itype == 1 ){
      SetFieldValue_Constant(id_field_,ilen,VELOCITY,world,aValueFieldDof_[ilen+nlen*1].val);  
    }
    else if( aValueFieldDof_[ilen+nlen*1].itype == 2 ){
      SetFieldValue_MathExp(id_field_,ilen,VELOCITY,world,aValueFieldDof_[ilen+nlen*1].math_exp,cur_time);
    }
  }
  for(unsigned int ilen=0;ilen<nlen;ilen++){
    if(      aValueFieldDof_[ilen+nlen*2].itype == 1 ){
      SetFieldValue_Constant(id_field_,ilen,ACCELERATION,world,aValueFieldDof_[ilen+nlen*2].val);      
    }
    else if( aValueFieldDof_[ilen+nlen*2].itype == 2 ){
      SetFieldValue_MathExp(id_field_,ilen,ACCELERATION,world,aValueFieldDof_[ilen+nlen*2].math_exp,cur_time);
    }
  }    
  return true;
}

void Update_FrictionalContact
(const CContactTarget3D& ct, double stiff_n, double stiff_f, double myu_s, double myu_k,
 double offset,
 unsigned int id_field_disp,
 Fem::Field::CFieldWorld& world,
 std::vector<CFrictionPoint>& aFrictionPoint )
{	
	assert( world.IsIdField(id_field_disp) );
	const Fem::Field::CField& field_disp = world.GetField(id_field_disp);
	assert( field_disp.GetFieldType() == Fem::Field::VECTOR3 );
	////////////////
	
	const unsigned int ndim = 3;
	const CNodeAry::CNodeSeg& ns_co      = field_disp.GetNodeSeg(  CORNER,false,world,VALUE);
	const CNodeAry::CNodeSeg& ns_udisp   = field_disp.GetNodeSeg(  CORNER,true, world,VALUE);
	const CNodeAry::CNodeSeg& ns_vdisp   = field_disp.GetNodeSeg(  CORNER,true, world,VELOCITY);
	
	assert( aFrictionPoint.size() == ns_co.Size() );
	for(unsigned int inode=0;inode<ns_co.Size();inode++)
	{
		double Co[ndim];	ns_co.GetValue(   inode,Co);
		double ud[ndim];	ns_udisp.GetValue(inode,ud);
		double uv[ndim];	ns_vdisp.GetValue(inode,uv);
		
		double co[3] = { Co[0]+ud[0], Co[1]+ud[1], Co[2]+ud[2] };
		
		double n1[3];
		const double pd1 = ct.Projection(co[0],co[1],co[2],n1)+offset;
		
		CFrictionPoint& fp = aFrictionPoint[inode];
		const double pd0 = fp.pd;

    if( fp.is_pin ) continue;
		
		if( pd1 < 0 ){} // not contact
		else if( pd0 < 0 ){	// h1>0 && h0<0	: contact in next step
			//			std::cout << "contact next time step" << std::endl;
			// put anchor at projected point
			for(unsigned int i=0;i<3;i++){ fp.aloc[i] = co[i]+pd1*n1[i]; }
		}
		else{	// h1>0 && h0>0
			// update the anchor that spring give force equal to dynamic friction direction to the velocity
			if( fp.itype_contact == 2 ){	
				// update anchor
				double v_t[3];
				{	// tangent direction
					const double t = Com::Dot3D(n1,uv);
					for(unsigned int i=0;i<3;i++){ v_t[i] = uv[i] - t*n1[i]; }
				}
				const double len_vt = Com::Length3D(v_t);
				const double invlen_vt = 1.0/len_vt;
				for(unsigned int i=0;i<3;i++){ v_t[i] *= invlen_vt; }
				const double dist = pd1*stiff_n*myu_k/stiff_f;
				for(unsigned int i=0;i<3;i++){ fp.aloc[i] = co[i]-v_t[i]*dist; }
			}
		}
	}
}

Fem::Field::CFieldValueSetter::CFieldValueSetter
(unsigned int id_field, Fem::Field::CFieldWorld& world) : id_field_(0), id_field_gradient_(0)
{
  if( !world.IsIdField(id_field) ){ return; }
  Fem::Field::CField& field = world.GetField(id_field);
  const unsigned int nlen = field.GetNLenValue();
  id_field_ = id_field;
  aValueFieldDof_.resize(nlen*3);
}

double CLinearSystem_Save::MakeResidual(const Fem::Field::CFieldWorld& world)
{
	const unsigned int nseg = this->m_aSegField.size();
	if( nseg == 0 )	return 0.0;

	////////////////////////////////
	// updateに値をセットする
	for(unsigned  int ilss=0;ilss<nseg;ilss++){
		const CLinSysSeg_Field& ls = this->m_aSegField[ilss];
		CVector_Blk* update = this->m_ls.m_Update[ilss];
		unsigned int ilen1 = 0;
		{
			const Fem::Field::CField& field = world.GetField(ls.id_field);
			ilen1 = field.GetNLenValue();
			unsigned int id_na_val = field.GetNodeSegInNodeAry(CORNER).id_na_va;
			unsigned int id_ns_val = field.GetNodeSegInNodeAry(CORNER).id_ns_va;
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			na.GetValueFromNodeSegment(id_ns_val,*update);
		}
		if( world.IsIdField(ls.id_field2) ){
			const Fem::Field::CField& field = world.GetField(ls.id_field2);
			unsigned int id_na_val = field.GetNodeSegInNodeAry(CORNER).id_na_va;
			unsigned int id_ns_val = field.GetNodeSegInNodeAry(CORNER).id_ns_va;
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			na.GetValueFromNodeSegment(id_ns_val,*update,ilen1);
		}
	}

	for(unsigned int iseg=0;iseg<nseg;iseg++){
		(*m_ls.m_Residual[iseg]) = (*m_Force[iseg]);
	}

	for(unsigned int iseg=0;iseg<nseg;iseg++){
	for(unsigned int jseg=0;jseg<nseg;jseg++){
		if( m_Matrix_Boundary[iseg][jseg] == 0 ) continue;
		m_Matrix_Boundary[iseg][jseg]->MatVec( -1.0, (*m_ls.m_Update[jseg]), 1.0, (*m_ls.m_Residual[iseg]), true );
	}
	}
	return sqrt(this->DOT(-1,-1));
}

bool CZLinearSystem_GeneralEigen::SetVector_fromField(int iv1,
	unsigned int id_field, const Fem::Field::CFieldWorld& world, Fem::Field::FIELD_DERIVATION_TYPE fdt )
{
	if( !world.IsIdField(id_field) ){ assert(0); return false; }
	const CField& field = world.GetField(id_field);
	if( field.GetFieldType() != ZSCALAR ) return false;
	
	const unsigned int nseg = this->m_aSeg.size();
	if( nseg == 0 )	return true;

	std::vector< CZVector_Blk* >* paVec1 = 0;
    if(      iv1 >= 0 && iv1 < (int)this->GetTmpVectorArySize() ) paVec1 = &m_TmpVectorArray[iv1];
	else if( iv1 == -1 ) paVec1 = &this->m_Residual;
	else if( iv1 == -2 ) paVec1 = &this->m_Update;
	else assert(0);


	// コーナー節点について値を更新
	{
/*		const int ilss0 = this->FindIndexArray_Seg(id_field,CORNER,world);
		assert( ilss0 != -1 );
		std::cout << ilss0 << std::endl;*/
        unsigned int ilss0 = 0;
		assert( ilss0 < (*paVec1).size() );
		CZVector_Blk* pVec = (*paVec1)[ilss0];
		const CField& field = world.GetField(id_field);
		unsigned int id_na = field.GetNodeSegInNodeAry(CORNER).id_na_va;
		const CNodeAry& na = world.GetNA(id_na);
		const CNodeAry::CNodeSeg& ns = field.GetNodeSeg(CORNER,true,world,fdt);
		const unsigned int nblk = na.Size();
		const unsigned int nlen = ns.Length();
		assert( nlen == 2 );
		double* val = new double [nlen];
		for(unsigned int iblk=0;iblk<nblk;iblk++){
			ns.GetValue(iblk,val);
			pVec->SetValue(iblk,0, Com::Complex(val[0],val[1]) );
		}
		delete[] val;
	}
	return true;
}

bool CEqn_Fluid3D::SetDomain(unsigned int id_base, Fem::Field::CFieldWorld& world)
{
    bool is_mixed = false;
    Field::ELEM_TYPE elem_type = (Field::ELEM_TYPE)0;
    {
        assert( world.IsIdField(id_base) );
        const CField& field_base = world.GetField(id_base);
        const std::vector<unsigned int>& aIdEA = field_base.GetAryIdEA();
        for(unsigned int iiea=0; iiea<aIdEA.size(); iiea++) {
            unsigned int id_ea = aIdEA[iiea];
            CElemAry& ea = world.GetEA(id_ea);
            if( ea.ElemType() == TET ) {
                if( elem_type == 0 ) {
                    elem_type = TET;
                }
                else if( elem_type == HEX ) {
                    is_mixed = true;
                    break;
                }
            }
            else if( ea.ElemType() == HEX ) {
                if( elem_type == 0 ) {
                    elem_type = HEX;
                }
                else if( elem_type == TET ) {
                    is_mixed = true;
                    break;
                }
            }
        }
    }
    if( is_mixed ) {
        assert(0);
        std::cout << "Not Implimented" << std::endl;
        m_id_press = 0;
        m_id_velo = 0;
        return false;
    }
    if( elem_type == TET ) {
        m_id_press = world.MakeField_FieldElemDim(id_base,3,
                     Fem::Field::SCALAR, VELOCITY|ACCELERATION,CORNER);
    }
    else if( elem_type == HEX ) {
        m_id_press = world.MakeField_FieldElemDim(id_base,3,
                     Fem::Field::SCALAR, VELOCITY|ACCELERATION,BUBBLE);
    }
    m_id_velo  = world.MakeField_FieldElemDim(id_base,3,
                 Fem::Field::VECTOR3,VELOCITY|ACCELERATION,CORNER);
    this->ClearLinearSystemPreconditioner();
    return true;
}

bool CEdgeTextureColor::Update(const Fem::Field::CFieldWorld& world)
{
	if( !world.IsIdField(this->id_field_velo) ){ return false; }
	if( !world.IsIdField(this->id_part_field_velo) ){ return false; }
	if( !world.IsIdEA(this->id_ea) ){ return false; }
	const Fem::Field::CElemAry& ea = world.GetEA(id_ea);
	assert( ea.ElemType() == Fem::Field::LINE );
	if( nelem != ea.Size() ){
		nelem  = ea.Size();
		if( m_aXYVeloElem != 0 ){ delete[] m_aXYVeloElem; }
		if( m_aXYElem     != 0 ){ delete[] m_aXYElem; }
		m_aXYVeloElem = new double [nelem*4];
		m_aXYElem     = new double [nelem*4];
	}
	const Fem::Field::CField& fv = world.GetField(this->id_part_field_velo);
//	const Fem::Field::CField::CNodeSegInNodeAry& nans = fv.GetNodeSegInNodeAry(Fem::Field::CORNER);
	const Fem::Field::CNodeAry::CNodeSeg& ns_v = fv.GetNodeSeg(CORNER,true, world,VELOCITY);
	const Fem::Field::CNodeAry::CNodeSeg& ns_c = fv.GetNodeSeg(CORNER,false,world,VELOCITY);
	const Fem::Field::CElemAry::CElemSeg& es_v = fv.GetElemSeg(id_ea,CORNER,true, world);
	const Fem::Field::CElemAry::CElemSeg& es_c = fv.GetElemSeg(id_ea,CORNER,false,world);
	assert( es_v.Length() == 2 );
	assert( es_c.Length() == 2 );
	for(unsigned int ielem=0;ielem<nelem;ielem++)
	{
		unsigned int noes[2];
		es_c.GetNodes(ielem,noes);
		double co[2][2];
		ns_c.GetValue(noes[0],co[0]);
		ns_c.GetValue(noes[1],co[1]);
		////////////////
		es_v.GetNodes(ielem,noes);
		double velo[2][2];
		ns_v.GetValue(noes[0],velo[0]);	
		ns_v.GetValue(noes[1],velo[1]);
		////////////////
		m_aXYElem[ielem*4+0] = co[0][0];
		m_aXYElem[ielem*4+1] = co[0][1];
		m_aXYElem[ielem*4+2] = co[1][0];
		m_aXYElem[ielem*4+3] = co[1][1];
		////////////////
		const double r = velo_scale;
		m_aXYVeloElem[ielem*4+0] = co[0][0] + r*velo[0][0];	
		m_aXYVeloElem[ielem*4+1] = co[0][1] + r*velo[0][1];
		m_aXYVeloElem[ielem*4+2] = co[1][0] + r*velo[1][0];
		m_aXYVeloElem[ielem*4+3] = co[1][1] + r*velo[1][1];
	}
	return true;
}

bool CEqnSystem_Scalar2D::SetDomain_FieldElemAry(unsigned int id_base, unsigned int id_ea, Fem::Field::CFieldWorld& world)
{
	m_IdFieldVal = world.MakeField_FieldElemAry(id_base,id_ea,SCALAR,VELOCITY|VALUE,CORNER);
	{
		m_aEqn.clear();
		const CField& field = world.GetField(m_IdFieldVal);
		const std::vector<unsigned int>& aIdEA = field.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			m_aEqn.push_back( CEqn_Scalar2D(id_ea,m_IdFieldVal) );
		}
	}
	this->ClearLinearSystemPreconditioner();
	this->m_aIdFixField.clear();
	return true;
}

bool CEqnSystem_Solid2D::UpdateDomain_Field(unsigned int id_base, Fem::Field::CFieldWorld& world)
{
	m_IdFieldDisp  = world.MakeField_FieldElemDim(id_base,2,
		Fem::Field::VECTOR2,VALUE|VELOCITY|ACCELERATION,CORNER);
    assert( world.IsIdField(m_IdFieldDisp) );

	{	// ìØÇ∂óvëfîzóÒIDÇéùǬï˚íˆéÆÇ™Ç Ç¡ÇΩÇÁÅCǪÇÍÇégǧÅDÇ»ÇØÇÍÇŒêVãKÇ…í«â¡
		std::vector<CEqn_Solid2D> aEqn_old = m_aEqn;
		m_aEqn.clear();
		const CField& field = world.GetField(m_IdFieldDisp);
		const std::vector<unsigned int>& aIdEA = field.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			unsigned int ieqn0=0;
			for(;ieqn0<aEqn_old.size();ieqn0++){
				if( aEqn_old[ieqn0].GetIdEA() == id_ea ){ 
					const unsigned int ieqn1 = m_aEqn.size();
					m_aEqn.push_back( aEqn_old[ieqn0] );
					m_aEqn[ieqn1].SetIdFieldDisp(m_IdFieldDisp);
					break; 
				}
			}
			if( ieqn0 != aEqn_old.size() ){ continue; }
            CEqn_Solid2D eqn1(id_ea,m_IdFieldDisp);
			eqn1.SetYoungPoisson(m_young_back, m_poisson_back, m_is_plane_stress_back);
			eqn1.SetRho(m_rho_back);
			eqn1.SetGeometricalNonlinear(m_is_geom_nonlin_back);
            m_aEqn.push_back( eqn1 );
		}
	}

    std::cout << "Size Eqn : " << m_aEqn.size() << std::endl;
    for(unsigned int ieqn=0;ieqn<m_aEqn.size();ieqn++){
        const CEqn_Solid2D& eqn = m_aEqn[ieqn];
		double young, poisson;
		eqn.GetYoungPoisson(young,poisson);
//        std::cout << ieqn << " " << eqn.GetIdEA() << " " << young << " " << poisson << " " << std::endl;
//		m_aEqn[ieqn].SetGravitation(0,-0.0);
    }
	

	if( !world.IsIdField(m_IdFieldDisp) ) return false;
	this->ClearLinearSystemPreconditioner();
	this->m_aIdFixField.clear();
	return true;
}

// Dynamic
bool Fem::Eqn::AddLinSys_StVenant2D_NonStatic_BackwardEular
(double dt, 
 Fem::Eqn::ILinearSystem_Eqn& ls,
 double lambda, double myu,
 double rho, double g_x, double g_y,
 const Fem::Field::CFieldWorld& world,
 unsigned int id_field_disp, 
 const MatVec::CVector_Blk& velo_pre,
 bool is_initial,
 unsigned int id_ea )
{
	const CField& field_disp = world.GetField(id_field_disp);
	if( field_disp.GetFieldType() != VECTOR2 ) return false;
	
	if( id_ea != 0 ){
		if( field_disp.GetInterpolationType(id_ea,world) == TRI11 ){
			AddLinSys_StVenant2D_NonStatic_BackwardEular_P1
			(dt,
			 ls,
			 lambda, myu,
			 rho, g_x, g_y,
			 id_field_disp,world,
			 velo_pre,
			 is_initial,
			 id_ea);
		}
		else{ assert(0); }
	}
	else{
		const std::vector<unsigned int>& aIdEA = field_disp.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			bool res = Fem::Eqn::AddLinSys_StVenant2D_NonStatic_BackwardEular
			(dt,
			 ls,
			 lambda, myu, rho, g_x, g_y,
			 world, id_field_disp, 
			 velo_pre,
			 is_initial, 
			 id_ea );
			if( !res ) return false;
		}
		return true;
	}
	return true;
}