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

bool CLinearSystem_SaveDiaM_NewmarkBeta::UpdateValueOfField(
		unsigned int id_field, Fem::Field::CFieldWorld& world, Fem::Field::FIELD_DERIVATION_TYPE fdt )
{
//	std::cout << "CLinearSystem_SaveDiaM_NewmarkBeta::UpdateValueOfField" << std::endl;
	if( !world.IsIdField(id_field) ) return false;
	CField& field = world.GetField(id_field);

	assert( fdt == ACCELERATION );
	{
		unsigned int id_na_val = field.GetNodeSegInNodeAry(CORNER).id_na_va;
		unsigned int id_ns_a = field.GetNodeSegInNodeAry(CORNER).id_ns_ac;
		if( id_na_val != 0 ){
			const unsigned int ilss0 = this->FindIndexArray_Seg(id_field,CORNER,world);
			const CLinSysSeg_Field& lss = this->m_aSegField[ilss0];
			CVector_Blk* pUpdate = this->m_ls.m_Update[ilss0];
			assert( pUpdate != 0 );
			CNodeAry& na = world.GetNA(id_na_val);
			const unsigned int nblk = na.Size();
			CNodeAry::CNodeSeg& ns_u = field.GetNodeSeg(CORNER,true,world,VALUE);
			CNodeAry::CNodeSeg& ns_v = field.GetNodeSeg(CORNER,true,world,VELOCITY);
			CNodeAry::CNodeSeg& ns_a = field.GetNodeSeg(CORNER,true,world,ACCELERATION);
			assert( ns_v.Length() == ns_u.Length() );
			const unsigned int nlen = ns_u.Length();
			const CBCFlag& bc_flag = m_ls.GetBCFlag(ilss0);//this->m_ls.m_BCFlag[ilss0];
			assert( nblk == bc_flag.NBlk() );
			double* velo0 = new double [nlen];
			double* acc0  = new double [nlen];
			if( lss.id_field == id_field ){
				for(unsigned int iblk=0;iblk<nblk;iblk++){
					ns_v.GetValue(iblk,velo0);
					ns_a.GetValue(iblk, acc0);
					for(unsigned int ilen=0;ilen<nlen;ilen++){
						if( bc_flag.GetBCFlag(iblk,ilen) == 0 ){	// 固定境界条件が入ってなかったら
							const double acc1 = pUpdate->GetValue(iblk,ilen);	// 更新後の速度
							// 値を更新
							ns_u.AddValue(iblk,ilen,dt*velo0[ilen]+dt*dt*(0.5-beta)*acc0[ilen]+dt*dt*beta*acc1);
							ns_v.AddValue(iblk,ilen,dt*(1-gamma)*acc0[ilen]+gamma*dt*acc1);
							ns_a.AddValue(iblk,ilen,acc1);
						}
						else{	// 固定境界条件が入ってなかったら
							// Updateに固定境界条件つき速度を入れる
							pUpdate->SetValue(iblk,ilen,acc0[ilen]);
						}
					}
				}
				na.SetValueToNodeSegment(id_ns_a,*pUpdate);
			}
			else{
				assert( lss.id_field2 == id_field );
				unsigned int nlen0;
				{
					const unsigned int id_field0 = lss.id_field;
					const CField& field = world.GetField(id_field0);
					nlen0 = field.GetNLenValue();
				}
				for(unsigned int iblk=0;iblk<nblk;iblk++){
					ns_v.GetValue(iblk,velo0);
					ns_a.GetValue(iblk, acc0);
					for(unsigned int ilen=0;ilen<nlen;ilen++){
						if( bc_flag.GetBCFlag(iblk,ilen+nlen0) == 0 ){	// 固定境界条件が入ってなかったら
							const double acc1 = pUpdate->GetValue(iblk,ilen+nlen0);	// 更新後の速度
							// 値を更新
							ns_u.AddValue(iblk,ilen,dt*velo0[ilen]+dt*dt*(0.5-beta)*acc0[ilen]+dt*dt*beta*acc1);
							ns_v.AddValue(iblk,ilen,dt*(1-gamma)*acc0[ilen]+gamma*dt*acc1);
							ns_a.AddValue(iblk,ilen,acc1);
						}
						else{	// 固定境界条件が入ってなかったら
							// Updateに固定境界条件つき速度を入れる
							pUpdate->SetValue(iblk,ilen+nlen0,acc0[ilen]);
						}
					}
				}
				na.SetValueToNodeSegment(id_ns_a,*pUpdate,nlen0);
			}
			delete[] velo0;
			delete[] acc0;
			// 加速度の更新
		}
	}
	return true;
}

bool CDrawerFace::Update
(const Fem::Field::CFieldWorld& world)
{
	const Fem::Field::CField& field = world.GetField(m_id_field);
	// set the vertex array
	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();
	const unsigned int npoin_co = na_c_co.Size();
	unsigned int npoin;	
	{ // calc the number of point
		if( m_is_draw_nsv ){
			unsigned int id_na_c_val = field.GetNodeSegInNodeAry(CORNER).id_na_va;
			assert( world.IsIdNA(id_na_c_val) );
			const Fem::Field::CNodeAry& na_c_val = world.GetNA(id_na_c_val);
			npoin = na_c_val.Size();
		}
		else{ npoin = npoin_co; }
	}
//	std::cout << m_vertex_ary.NPoin() << " " << npoin << std::endl;
	assert( m_vertex_ary.NPoin() == npoin );

	if( !m_isnt_value_disp ){	// including displacement
		unsigned int id_na_c_val = field.GetNodeSegInNodeAry(CORNER).id_na_va;
		assert( world.IsIdNA(id_na_c_val) );
		const Fem::Field::CNodeAry& na_c_val = world.GetNA(id_na_c_val);
		const Fem::Field::CNodeAry::CNodeSeg& ns_c_val = field.GetNodeSeg(CORNER,true,world,VALUE|VELOCITY|ACCELERATION);
		if( ndim == 2 && (field.GetFieldType()==SCALAR||field.GetFieldType()==ZSCALAR) )	// 垂直方向の変位として捉える
		{
			assert( m_vertex_ary.NDim() == 3 );
			double coord[2], value[2];
			for(unsigned int ipoin=0;ipoin<npoin;ipoin++){
				unsigned int ipoin_co = field.GetMapVal2Co(ipoin);
				assert( ipoin_co < npoin_co );
				ns_c_val.GetValue(ipoin,value);
				ns_c_co.GetValue(ipoin_co,coord);
				this->m_vertex_ary.pVertexArray[ipoin*3+0] = coord[0];
				this->m_vertex_ary.pVertexArray[ipoin*3+1] = coord[1];
				this->m_vertex_ary.pVertexArray[ipoin*3+2] = value[0];
			}
		}
		else{
			assert( m_vertex_ary.NDim() == ndim );
			assert( ndim == ns_c_val.Length() ); // Coordの次元とValueの次元が合ってなければならない
			assert( npoin_co == na_c_co.Size() );
			assert( na_c_val.Size() == npoin );
			assert( ndim <= 3 );
			double coord[3], value[3];
			for(unsigned int ipoin=0;ipoin<npoin;ipoin++){
				unsigned int ipoin_co = field.GetMapVal2Co(ipoin);
				assert( ipoin_co < npoin_co );
				ns_c_val.GetValue(ipoin,value);
				ns_c_co.GetValue(ipoin_co,coord);
				for(unsigned int idim=0;idim<ndim;idim++){
					this->m_vertex_ary.pVertexArray[ipoin*ndim+idim] = coord[idim]+value[idim];
				}
			}
		}
	}
	else{		
		assert( m_vertex_ary.NDim() == ndim );
		for(unsigned int ipoin=0;ipoin<npoin;ipoin++){
			unsigned int ipoin_co = field.GetMapVal2Co(ipoin);
			assert( ipoin_co < npoin_co );
			double* pval = &this->m_vertex_ary.pVertexArray[ipoin*ndim];
			ns_c_co.GetValue(ipoin_co,pval);
		}
	}

	////////////////////////////////////////////////
	// make color
  
	if( world.IsIdField(id_field_val) )
	{
		const Fem::Field::CField& field_val = world.GetField(id_field_val);
		Fem::Field::FIELD_DERIVATION_TYPE fdt;
		{	// 描画する値の微分タイプ(できるだけVALUEを選択)
			unsigned int fdt_all = field_val.GetFieldDerivativeType();
			if(      fdt_all & VALUE        ){ fdt = VALUE; }
			else if( fdt_all & VELOCITY     ){ fdt = VELOCITY; }
			else if( fdt_all & ACCELERATION ){ fdt = ACCELERATION; }
			else{ assert(0); }
		}
		if( !color_map->IsMinMaxFix() ){	// 値の最大値最小値を求める
			double min_val, max_val;
			field_val.GetMinMaxValue(min_val,max_val,world,0,fdt);
			color_map->SetMinMax(min_val,max_val);
		}
		unsigned int id_na_c_val = field_val.GetNodeSegInNodeAry(CORNER).id_na_va;
		unsigned int id_na_b_val = field_val.GetNodeSegInNodeAry(BUBBLE).id_na_va;
		if(      world.IsIdNA(id_na_c_val) ){
			const CNodeAry::CNodeSeg& ns_v = field_val.GetNodeSeg(CORNER,true,world,fdt);
			double val[10];
			if( pColorArray == 0 ){ pColorArray = new float [npoin*4]; }
			for(unsigned int ipoin=0;ipoin<npoin;ipoin++){
				ns_v.GetValue(ipoin,val);
        color_map->GetColor(&pColorArray[ipoin*4],val[0]);
//.........这里部分代码省略.........

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

	// 外力項をセット
	for(unsigned int iseg=0;iseg<nseg;iseg++){
		(*m_ls.m_Residual[iseg]) = (*m_Force[iseg]);
	}
	
	////////////////////////////////
	// updateにKを掛けるための値をセットする
	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;
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			unsigned int id_ns_u = field.GetNodeSegInNodeAry(CORNER).id_ns_va;
			unsigned int id_ns_v = field.GetNodeSegInNodeAry(CORNER).id_ns_ve;
			unsigned int id_ns_a = field.GetNodeSegInNodeAry(CORNER).id_ns_ac;
			na.GetValueFromNodeSegment(                 id_ns_u,*update);
			na.AddValueFromNodeSegment(dt,              id_ns_v,*update);
			na.AddValueFromNodeSegment(dt*dt*(0.5-beta),id_ns_a,*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;
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			unsigned int id_ns_u = field.GetNodeSegInNodeAry(CORNER).id_ns_va;
			unsigned int id_ns_v = field.GetNodeSegInNodeAry(CORNER).id_ns_ve;
			unsigned int id_ns_a = field.GetNodeSegInNodeAry(CORNER).id_ns_ac;
			na.GetValueFromNodeSegment(                 id_ns_u,*update,ilen1);
			na.AddValueFromNodeSegment(dt,              id_ns_v,*update,ilen1);
			na.AddValueFromNodeSegment(dt*dt*(0.5-beta),id_ns_a,*update,ilen1);
		}
	}

	// 行列ベクトル積により残差を計算する
	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 );
	}
	}
/*
	////////////////////////////////
	// updateにMを掛けるため加速度値をセットする
	for(unsigned int ilss=0;ilss<nseg;ilss++){
		const CLinSysSeg_Field& ls = this->m_aSegField[ilss];
		CVector_Blk* update = this->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;
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			unsigned int id_ns_a = field.GetNodeSegInNodeAry(CORNER).id_ns_ac;
			na.GetValueFromNodeSegment(id_ns_a,*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;
			assert( world.IsIdNA(id_na_val) );
			const CNodeAry& na = world.GetNA(id_na_val);
			unsigned int id_ns_a = field.GetNodeSegInNodeAry(CORNER).id_ns_ac;
			na.GetValueFromNodeSegment(id_ns_a,*update,ilen1);
		}
	}

	// 行列ベクトル積により残差を計算する
	for(unsigned int iseg=0;iseg<nseg;iseg++){
		m_DiaMassMatrix[iseg]->MatVec( -1.0, *m_Update[iseg], 1.0, *m_Residual[iseg] );
	}
*/
	// 残差への境界条件のセット
	for(unsigned int iseg=0;iseg<nseg;iseg++){
//		m_ls.m_BCFlag[iseg]->SetZeroToBCDof(*m_ls.m_Residual[iseg]);
		m_ls.GetBCFlag(iseg).SetZeroToBCDof(*m_ls.m_Residual[iseg]);
	}

	return sqrt(this->DOT(-1,-1));
}

bool CEqnSystem_Fluid2D::UpdateDomain_Field(unsigned int id_base, Fem::Field::CFieldWorld& world)
{
    m_id_press = world.MakeField_FieldElemDim(id_base,2,
                 Fem::Field::SCALAR, VELOCITY|ACCELERATION,CORNER);
    assert( world.IsIdField(m_id_press) );
    //    std::cout << "press " << m_id_press << std::endl;

    if( m_IsntInterpolationBubble) {
        std::cout << "corner intp" << std::endl;
        m_id_velo  = world.MakeField_FieldElemDim(id_base,2,
                     Fem::Field::VECTOR2,VELOCITY|ACCELERATION,CORNER);
        m_IsntCombine = false;
    }
    else {
        std::cout << "bubble intp" << std::endl;
        m_id_velo  = world.MakeField_FieldElemDim(id_base,2,
                     Fem::Field::VECTOR2,VELOCITY|ACCELERATION,CORNER|BUBBLE);
        m_IsntCombine = true;
    }
    //    std::cout << "velo : " << m_id_velo << std::endl;
    assert( world.IsIdField(m_id_velo) );

    {   // 同じ要素配列IDを持つ方程式があったら，それを使う．なければ新規に追加
        std::vector<CEqn_Fluid2D> aEqn_old = m_aEqn;
        m_aEqn.clear();
        const CField& field = world.GetField(m_id_velo);
        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].SetIdFieldVelocity(m_id_velo);
                    m_aEqn[ieqn1].SetIdFieldPressure(m_id_press);
                    break;
                }
            }
            if( ieqn0 != aEqn_old.size() ) {
                continue;
            }
            CEqn_Fluid2D eqn1(id_ea,m_id_velo,m_id_press);
            eqn1.SetRho(m_rho_back);
            eqn1.SetMyu(m_myu_back);
            if( this->m_is_stokes_back ) {
                eqn1.SetStokes();
            }
            else {
                eqn1.SetNavierStokes();
            }
            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_Fluid2D& eqn = m_aEqn[ieqn];
     std::cout << ieqn << " " << eqn.GetMyu() << " " << eqn.GetRho() << " " << eqn.IsNavierStokes() << std::endl;
     }
     */
    if( !world.IsIdField(m_id_velo) ) return false;
    this->ClearLinearSystemPreconditioner();
    return true;
}

bool CEqn_Fluid3D::AddFixField(const unsigned int id_field, Fem::Field::CFieldWorld& world, int idof)
{
    if( !world.IsIdField( id_field ) ) return false;
    m_aIdFixField.push_back( std::make_pair(id_field,idof) );
    return true;
}

void Fem::Field::CFieldValueSetter::SetGradient
(unsigned int id_field_from, Fem::Field::CFieldWorld& world)
{
  if( !world.IsIdField(id_field_from) ){ return; }
  this->id_field_gradient_ = id_field_from;
}

bool AddLinSys_FrictionalContact_Penalty_NonStatic_Sensitivity
(Fem::Ls::CLinearSystem_Field& ls, 
 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 )
{
	if( !world.IsIdField(id_field_disp) ) return false;
	const Fem::Field::CField& field_disp = world.GetField(id_field_disp);
	if( field_disp.GetFieldType() != Fem::Field::VECTOR3 ) return false;
	////////////////
	
	MatVec::CMatDia_BlkCrs& pmat_dd = ls.GetMatrix(id_field_disp,  CORNER,world);
	MatVec::CVector_Blk& res_d = ls.GetResidual(id_field_disp,  CORNER,world);   
	
	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 co[3] = { Co[0]+ud[0], Co[1]+ud[1], Co[2]+ud[2] };
    CFrictionPoint& fp = aFrictionPoint[inode];    
    double n0[3];
    const double pd = ct.Projection(co[0],co[1],co[2], n0)+offset;
    fp.pd = pd;
    
    if( pd < 0 ){
      fp.itype_contact = 0;
      continue;
    }
    double eKmat[3][3];
    for(unsigned int i=0;i<3;i++){
      for(unsigned int j=0;j<3;j++){
        eKmat[i][j] = stiff_n*n0[i]*n0[j];
      }
    }
    double eres_d[3];
    eres_d[0] = stiff_n*n0[0]*pd;
    eres_d[1] = stiff_n*n0[1]*pd;
    eres_d[2] = stiff_n*n0[2]*pd;
        
    for(unsigned int i=0;i<3;i++){
      for(unsigned int j=0;j<3;j++){
        eKmat[i][j] += -n0[i]*n0[j]*stiff_f;
      }
      eKmat[i][i] += stiff_f;      
    }
    double ap_t[3] = { co[0]-fp.aloc[0], co[1]-fp.aloc[1], co[2]-fp.aloc[2] };    
    for(unsigned int i=0;i<3;i++){ eres_d[i] += -stiff_f*ap_t[i]; }
    pmat_dd.Mearge(1,&inode,  1,&inode,     9,  &eKmat[0][0]);
    res_d.AddValue(inode,0,eres_d[0]);
    res_d.AddValue(inode,1,eres_d[1]);
    res_d.AddValue(inode,2,eres_d[2]);    
  }  
  return true;
}

//! set constant value to the field
bool Fem::Field::SetFieldValue_Constant
(unsigned int id_field_to, unsigned int idofns, Fem::Field::FIELD_DERIVATION_TYPE fdt, 
 Fem::Field::CFieldWorld& world,
 double val)
{
  if( !world.IsIdField(id_field_to) ){ return false; }
  Fem::Field::CField& field = world.GetField(id_field_to);
  
	if( idofns >= field.GetNLenValue() ) return false;
	if( !(field.GetFieldDerivativeType() & fdt) ) return false;
  
  const std::vector<unsigned int>& aIdEA = field.GetAryIdEA();
  
  if( field.GetNodeSegInNodeAry(CORNER).id_na_va != 0 ){
    const CField::CNodeSegInNodeAry& nsna = field.GetNodeSegInNodeAry(CORNER);      
    assert( world.IsIdNA(nsna.id_na_va) );
    CNodeAry& na = world.GetNA(nsna.id_na_va);			
    unsigned int id_ns = 0;	// target node segment;
    {
      if(      fdt & VALUE        ) id_ns = nsna.id_ns_va;
      else if( fdt & VELOCITY     ) id_ns = nsna.id_ns_ve;
      else if( fdt & ACCELERATION ) id_ns = nsna.id_ns_ac;
      else{ assert(0); }
      assert( na.IsSegID(id_ns) );
      const CNodeAry::CNodeSeg& ns = na.GetSeg(id_ns);
      assert( idofns < ns.Length() );
    }
    /*
     if( m_aElemIntp.size() == 0 || !this->IsPartial() ){	// 剛体の場合
     const unsigned int nnode = na.Size();
     CNodeAry::CNodeSeg& ns = na.GetSeg(id_ns);
     for(unsigned int inode=0;inode<nnode;inode++){
     ns.SetValue(inode,idofns,const_val);
     }
     }*/
    {
      for(unsigned int iei=0;iei<aIdEA.size();iei++){
        const unsigned int id_ea = aIdEA[iei];
        assert( world.IsIdEA(id_ea) );
        const CElemAry& ea = world.GetEA(id_ea);
        unsigned int id_es_c_va = field.GetIdElemSeg(id_ea,CORNER,true,world);
        assert( ea.IsSegID(id_es_c_va) );
        na.SetValueToNodeSegment(ea,id_es_c_va,id_ns,idofns,val);
      }
    }
  }
  /*		if( m_na_e.id_na_va != 0 ){
   assert( world.IsIdNA(m_na_e.id_na_va) );
   CNodeAry& na = world.GetNA(m_na_e.id_na_va);
   const CNodeAry::CNodeSeg& ns_val = na.GetSeg(m_na_e.id_ns_va);
   for(unsigned int iei=0;iei<m_aElemIntp.size();iei++){
   unsigned int id_ea = m_aElemIntp[iei].id_ea;
   assert( world.IsIdEA(id_ea) );
   const CElemAry& ea = world.GetEA(id_ea);
   assert( ea.IsSegID( m_aElemIntp[iei].id_es_e_va ) );
   na.SetValueToNodeSegment(ea,m_aElemIntp[iei].id_es_e_va,
   m_na_e.id_ns_va,idofns,const_val);
   }
   }*/
  if( field.GetNodeSegInNodeAry(BUBBLE).id_na_va != 0 ){
    const CField::CNodeSegInNodeAry& nsna = field.GetNodeSegInNodeAry(BUBBLE);
    assert( world.IsIdNA(nsna.id_na_va) );
    CNodeAry& na = world.GetNA(nsna.id_na_va);
    unsigned int id_ns = 0;	// target node segment;
    {
      if(      fdt & VALUE        ) id_ns = nsna.id_ns_va;
      else if( fdt & VELOCITY     ) id_ns = nsna.id_ns_ve;
      else if( fdt & ACCELERATION ) id_ns = nsna.id_ns_ac;
      else{ assert(0); }
      assert( na.IsSegID(id_ns) );
      const CNodeAry::CNodeSeg& ns = na.GetSeg(id_ns);
      assert( idofns < ns.Length() );
    }
    //			const CNodeAry::CNodeSeg& ns_val = na.GetSeg(id_ns);
    for(unsigned int iea=0;iea<aIdEA.size();iea++){
      const unsigned int id_ea = aIdEA[iea];
      assert( world.IsIdEA(id_ea) );
      const CElemAry& ea = world.GetEA(id_ea);
      unsigned int id_es_b_va = field.GetIdElemSeg(id_ea,BUBBLE,true,world);
      assert( ea.IsSegID(id_es_b_va) );
      na.SetValueToNodeSegment(ea,id_es_b_va,id_ns,idofns,val);
    }
  }  
  return true;
}

// ÉXÉJÉâÅ[ÇÃâûóÕëäìñílÇí«â¡Ç∑ÇÈÅDǪÇÃǧÇøÉÇÅ[ÉhÇǬÇØÇÈ
// mode = 0 : É~Å[É[ÉX
// mode = 1 : ç≈ëÂâûóÕ
// äÙâΩäwìIÇ»îÒê¸å`ê´Ççló∂Ç∑ÇÈǩǫǧǩÇÕÅCEqnObjÇå©ÇƒåàÇflÇÈÇÊǧDžÇ∑ÇÈ
bool CEqnSystem_Solid2D::SetEquivStressValue(unsigned int id_field_str, Fem::Field::CFieldWorld& world)
{	
	if( !world.IsIdField(id_field_str) ) return false;
	Fem::Field::CField& field_str = world.GetField(id_field_str);
	if( field_str.GetFieldType() != SCALAR ) return false;

	if( !world.IsIdField(m_IdFieldDisp) ) return false;
	Fem::Field::CField& field_dis = world.GetField(m_IdFieldDisp);

	const std::vector<unsigned int>& aIdEA_from = field_dis.GetAryIdEA();
	const std::vector<unsigned int>& aIdEA_to   = field_str.GetAryIdEA();
	if( aIdEA_from.size() != aIdEA_to.size() ) return false;

	const unsigned int niea = aIdEA_from.size();
	for(unsigned int iiea=0;iiea<niea;iiea++)
	{
		Fem::Field::INTERPOLATION_TYPE type_from, type_to;
		{
			if( aIdEA_from[iiea] != aIdEA_to[iiea] ){
				assert(0);
				return false;
			}
			const unsigned int id_ea = aIdEA_from[iiea];
			type_from = field_dis.GetInterpolationType(id_ea,world);
			type_to   = field_str.GetInterpolationType(id_ea,world);
		}

		unsigned int nnoes, ndim;
		if( type_from==TRI11 && type_to==TRI1001 ){
			nnoes = 3; ndim = 2;
		}
		else{
			std::cout << "Error!-->Not Implimented!" << std::endl;
			std::cout << type_from << " " << type_to << std::endl;
			assert(0);
			getchar();
		}

		unsigned int id_ea = aIdEA_to[iiea];
		const CElemAry& ea = world.GetEA(id_ea);
		const CElemAry::CElemSeg& es_c_co = field_dis.GetElemSeg(id_ea,CORNER,false,world);
		const CElemAry::CElemSeg& es_c_va = field_dis.GetElemSeg(id_ea,CORNER,true, world);
		const CElemAry::CElemSeg& es_b_va = field_str.GetElemSeg(id_ea,BUBBLE,true, world);

		Fem::Field::CField::CNodeSegInNodeAry nans_c = field_dis.GetNodeSegInNodeAry(CORNER);
		assert( world.IsIdNA(nans_c.id_na_co) );
		assert( world.IsIdNA(nans_c.id_na_va) );
		const CNodeAry& na_c_co = world.GetNA( nans_c.id_na_co);
		const CNodeAry::CNodeSeg& ns_c_co = na_c_co.GetSeg( nans_c.id_ns_co);
		const CNodeAry& na_c_va = world.GetNA( nans_c.id_na_va);
		const CNodeAry::CNodeSeg& ns_c_va = na_c_va.GetSeg( nans_c.id_ns_va);

		Fem::Field::CField::CNodeSegInNodeAry nans_b = field_str.GetNodeSegInNodeAry(BUBBLE);
		unsigned int id_na_b_va = nans_b.id_na_va;
		unsigned int id_ns_b_va = nans_b.id_ns_va;

		assert( world.IsIdNA(id_na_b_va) );
		CNodeAry& na_b_va = world.GetNA(id_na_b_va);
		CNodeAry::CNodeSeg& ns_b_va = na_b_va.GetSeg(id_ns_b_va);

    const unsigned int nnoes_c = es_c_co.Length();
    assert( nnoes_c < 64 );
		unsigned int noes[64];
    const CEqn_Solid2D& eqn = this->GetEquation(id_ea);
		for(unsigned int ielem=0;ielem<ea.Size();ielem++){
			////////////////
			double coord[16][3];
			// ç¿ïW(coord)Ç∆íl(value)ÇçÏÇÈ
			es_c_co.GetNodes(ielem,noes);
			for(unsigned int inoes=0;inoes<nnoes;inoes++){
				unsigned int ipoi0 = noes[inoes];
				assert( ipoi0 < na_c_co.Size() );
				ns_c_co.GetValue(ipoi0,coord[inoes]);
			}
			////////////////
			double disp[16][3];	// êflì_ïœà 
			es_c_va.GetNodes(ielem,noes);
			for(unsigned int inoes=0;inoes<nnoes;inoes++){
				unsigned int ipoi0 = noes[inoes];
				assert( ipoi0 < na_c_va.Size() );
				for(unsigned int idim=0;idim<ndim;idim++){
					ns_c_va.GetValue(ipoi0,disp[inoes]);
				}
			}
			////////////////
			double dudx[2][2];	// ïœå`å˘îz
			if( type_from == TRI11 ){
				double dldx[3][2];
				double const_term[3];
				TriDlDx(dldx,const_term, coord[0],coord[1],coord[2]);
				for(unsigned int i=0;i<ndim*ndim;i++){ (&dudx[0][0])[i] = 0.0; }
				for(unsigned int knoes=0;knoes<nnoes;knoes++){
					dudx[0][0] += disp[knoes][0]*dldx[knoes][0];
					dudx[0][1] += disp[knoes][0]*dldx[knoes][1];
					dudx[1][0] += disp[knoes][1]*dldx[knoes][0];
					dudx[1][1] += disp[knoes][1]*dldx[knoes][1];
//.........这里部分代码省略.........

bool Ls::CLinearSystem_RigidField2::UpdateValueOfField_NewmarkBeta(
		double gamma, double beta, double dt, 
		unsigned int id_field_val, Fem::Field::CFieldWorld& world, bool IsInitial )
{

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

    if( field_val.GetNodeSegInNodeAry(Fem::Field::CORNER).id_na_va != 0 ){	// 角節点を更新
		unsigned int id_na_val=0, id_ns_a,id_ns_v,id_ns_u;
		{
            const Fem::Field::CField::CNodeSegInNodeAry& nsna = field_val.GetNodeSegInNodeAry(Fem::Field::CORNER);
			id_na_val = nsna.id_na_va;
			id_ns_u = nsna.id_ns_va;
			id_ns_v = nsna.id_ns_ve;
			id_ns_a = nsna.id_ns_ac;
		}
        const int ilss = this->FindIndexArray_Seg(id_field_val,Fem::Field::CORNER,world);
		assert( ilss != -1 );
		const CLinSysSegRF& lss = this->m_aSegRF[ilss];
        const MatVec::CVector_Blk& upd = m_ls.GetVector(-2,ilss);
        Fem::Field::CNodeAry& na_val = world.GetNA(id_na_val);
		// 値を更新する
		if( IsInitial ){
			na_val.AddValueToNodeSegment(id_ns_u ,id_ns_v, dt);
			na_val.AddValueToNodeSegment(id_ns_u ,id_ns_a, 0.5*dt*dt);
		}
		na_val.AddValueToNodeSegment(id_ns_u, upd, beta*dt*dt);
		// 速度を更新する．
		if( IsInitial ){ 
		    na_val.AddValueToNodeSegment(id_ns_v ,id_ns_a, dt); 
		}
		na_val.AddValueToNodeSegment(id_ns_v, upd, gamma*dt);
		// 加速度を更新する
		na_val.AddValueToNodeSegment(id_ns_a, upd, 1.0);
	}
    if( field_val.GetNodeSegInNodeAry(Fem::Field::BUBBLE).id_na_va != 0 ){	// バブル節点を更新
		unsigned int id_na_val=0, id_ns_a,id_ns_v,id_ns_u;
		{
            const Fem::Field::CField::CNodeSegInNodeAry& nsna = field_val.GetNodeSegInNodeAry(Fem::Field::BUBBLE);
			id_na_val = nsna.id_na_va;
			id_ns_u = nsna.id_ns_va;
			id_ns_v = nsna.id_ns_ve;
			id_ns_a = nsna.id_ns_ac;
		}
        const int ilss = this->FindIndexArray_Seg(id_field_val,Fem::Field::BUBBLE,world);
		assert( ilss != -1 );
		const CLinSysSegRF& lss = this->m_aSegRF[ilss];
        const MatVec::CVector_Blk& upd = m_ls.GetVector(-2,ilss);
        Fem::Field::CNodeAry& na_val = world.GetNA(id_na_val);
		// 値を更新する
		if( IsInitial ){
			na_val.AddValueToNodeSegment(id_ns_u ,id_ns_v, dt);
			na_val.AddValueToNodeSegment(id_ns_u ,id_ns_a, 0.5*dt*dt);
		}
		na_val.AddValueToNodeSegment(id_ns_u, upd, beta*dt*dt);
		// 速度を更新する．
		if( IsInitial ){ 
		    na_val.AddValueToNodeSegment(id_ns_v ,id_ns_a, dt); 
		}
		na_val.AddValueToNodeSegment(id_ns_v, upd, gamma*dt);
		// 加速度を更新する
		na_val.AddValueToNodeSegment(id_ns_a, upd, 1.0);
    }
    if( field_val.GetNodeSegInNodeAry(Fem::Field::EDGE).id_na_va != 0 ){	// 角節点を更新
        assert(0);
    }
	return true;
}

bool CDrawerImageBasedFlowVis::Update(const Fem::Field::CFieldWorld& world)
{
	if( !world.IsIdField(this->m_IdFieldVelo) ){ return false; }
	const Fem::Field::CField& fv = world.GetField(this->m_IdFieldVelo);
	////////////////////////////////
	{	// 要素配列の設定
		const std::vector<unsigned int>& aIdEA = fv.GetAryIdEA();
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			assert( world.IsIdEA(id_ea) );
			const unsigned int id_es_v_c = fv.GetIdElemSeg(id_ea,CORNER,true, world);
			assert( id_es_v_c != 0 );
			CIndexArrayElem* pIAE = new CIndexArrayElem(id_ea,id_es_v_c,world);
			pIAE->ilayer = fv.GetLayer(id_ea);
			this->m_apIndexArrayElem.push_back( pIAE );
		}
	}
/*	const std::vector<unsigned int>& aIdEA = fv.GetAryIdEA();
	{
		unsigned int nelem0 = 0;
		for(unsigned int iiea=0;iiea<aIdEA.size();iiea++){
			const unsigned int id_ea = aIdEA[iiea];
			assert( world.IsIdEA(id_ea) );
			const Fem::Field::CElemAry& ea = world.GetEA(id_ea);
			nelem0 += ea.Size();
		}
		if( nelem != nelem0 ){
			nelem = nelem0;
			if( aXYVeloElem != 0 ){ delete[] aXYVeloElem; }
			if( aXYElem     != 0 ){ delete[] aXYElem; }
			aXYVeloElem = new double [nelem*6];
			aXYElem     = new double [nelem*6];
			if( world.IsIdField(m_IdFieldColor) ){
				if( aColorElem != 0 ){ delete[] aColorElem; }
				aColorElem = new double [nelem*3];
			}
		}
	}*/
	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);
	nnode = ns_v.Size();
	if( aVelo ==0 ){ aVelo  = new double [nnode*2]; }
	if( aCoord == 0 ){ aCoord = new double [nnode*2]; }
	for(unsigned int inode=0;inode<nnode;inode++){
		ns_v.GetValue(inode,aVelo+inode*2);
		unsigned int inoc = fv.GetMapVal2Co(inode);
		ns_c.GetValue(inoc,aCoord+inode*2);
	}
	/*
	unsigned int ielem_cur = 0;
	for(unsigned int iiea=0;iiea<aIdEA.size();iiea++)
	{
		const unsigned int id_ea = aIdEA[iiea];
		assert( world.IsIdEA(id_ea) );
		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);
		for(unsigned int ielem=0;ielem<es_v.GetSizeElem();ielem++){
			unsigned int noes[3];
			es_c.GetNodes(ielem,noes);
			double co[3][2];
			ns_c.GetValue(noes[0],co[0]);
			ns_c.GetValue(noes[1],co[1]);
			ns_c.GetValue(noes[2],co[2]);
			////////////////
			es_v.GetNodes(ielem,noes);
			double velo[3][2];
			ns_v.GetValue(noes[0],velo[0]);
			ns_v.GetValue(noes[1],velo[1]);
			ns_v.GetValue(noes[2],velo[2]);
			////////////////
			aXYElem[ielem_cur*6+0] = co[0][0];
			aXYElem[ielem_cur*6+1] = co[0][1];
			aXYElem[ielem_cur*6+2] = co[1][0];
			aXYElem[ielem_cur*6+3] = co[1][1];
			aXYElem[ielem_cur*6+4] = co[2][0];
			aXYElem[ielem_cur*6+5] = co[2][1];
			////////////////
            const double r = velo_scale;
			aXYVeloElem[ielem_cur*6+0] = co[0][0] + r*velo[0][0];
			aXYVeloElem[ielem_cur*6+1] = co[0][1] + r*velo[0][1];
			aXYVeloElem[ielem_cur*6+2] = co[1][0] + r*velo[1][0];
			aXYVeloElem[ielem_cur*6+3] = co[1][1] + r*velo[1][1];
			aXYVeloElem[ielem_cur*6+4] = co[2][0] + r*velo[2][0];
			aXYVeloElem[ielem_cur*6+5] = co[2][1] + r*velo[2][1];
			ielem_cur++;
		}
	}
    assert( ielem_cur == nelem );
	*/
	if( world.IsIdField(m_IdFieldColor) ){
		if( aValColor==0 ){ aValColor = new double [nnode]; }
		const Fem::Field::CField& fc = world.GetField(this->m_IdFieldColor);
		assert( fc.GetNodeSegInNodeAry(CORNER).id_na_va == fv.GetNodeSegInNodeAry(CORNER).id_na_va );
		const Fem::Field::CNodeAry::CNodeSeg& nsv_c = fc.GetNodeSeg(CORNER,true, world,VALUE);
		for(unsigned int inode=0;inode<nnode;inode++){
			nsv_c.GetValue(inode,aValColor+inode);
		}
	}
	else{
		if( aValColor != 0 ){ delete[] aValColor; aValColor = 0; }
//.........这里部分代码省略.........

bool CZLinearSystem::UpdateValueOfField( 
	unsigned int id_field, Fem::Field::CFieldWorld& world, Fem::Field::FIELD_DERIVATION_TYPE fdt )
{
	if( !world.IsIdField(id_field) ) return false;
	const CField& field = world.GetField(id_field);
	{
		unsigned int id_na_val=0, id_ns_val=0;
		{	// CORNERについて更新
			const CField::CNodeSegInNodeAry& nsna_c = field.GetNodeSegInNodeAry(CORNER);
			id_na_val = nsna_c.id_na_va;
			if(      fdt == VALUE )        id_ns_val = nsna_c.id_ns_va;
			else if( fdt == VELOCITY )     id_ns_val = nsna_c.id_ns_ve;
			else if( fdt == ACCELERATION ) id_ns_val = nsna_c.id_ns_ac;
			else{ assert(0); }
		}
		if( id_na_val != 0 ){
			CZVector_Blk* m_Update = this->GetUpdatePtr(id_field,CORNER,world);
			assert( m_Update != 0 );
			CNodeAry& na = world.GetNA(id_na_val);
			assert( na.IsSegID(id_ns_val) );
			na.AddValueToNodeSegment(id_ns_val,*m_Update,1.0);
//			CNodeAry::CNodeSeg& ns = na.GetSeg(id_ns_val);
		}
	}
	{
		unsigned int id_na_val=0, id_ns_val=0;
		{	// EDGEについて更新
			const CField::CNodeSegInNodeAry& nsna_c = field.GetNodeSegInNodeAry(EDGE);
			id_na_val = nsna_c.id_na_va;
			if(      fdt == VALUE )        id_ns_val = nsna_c.id_ns_va;
			else if( fdt == VELOCITY )     id_ns_val = nsna_c.id_ns_ve;
			else if( fdt == ACCELERATION ) id_ns_val = nsna_c.id_ns_ac;
			else{ assert(0); }
		}
		if( id_na_val != 0 ){
			CZVector_Blk* m_Update = this->GetUpdatePtr(id_field,EDGE,world);
			assert( m_Update != 0 );
			CNodeAry& na = world.GetNA(id_na_val);
			assert( na.IsSegID(id_ns_val) );
			na.AddValueToNodeSegment(id_ns_val,*m_Update,1.0);
//			CNodeAry::CNodeSeg& ns = na.GetSeg(id_ns_val);
//			ns.is_updated = true;
		}
	}
	{
		unsigned int id_na_val=0, id_ns_val=0;
		{	// BUBBLEについて更新
			const CField::CNodeSegInNodeAry& nsna_c = field.GetNodeSegInNodeAry(BUBBLE);
			id_na_val = nsna_c.id_na_va;
			if(      fdt == VALUE )        id_ns_val = nsna_c.id_ns_va;
			else if( fdt == VELOCITY )     id_ns_val = nsna_c.id_ns_ve;
			else if( fdt == ACCELERATION ) id_ns_val = nsna_c.id_ns_ac;
			else{ assert(0); }
		}
		if( id_na_val != 0 ){
			CZVector_Blk* m_Update = this->GetUpdatePtr(id_field,BUBBLE,world);
			assert( m_Update != 0 );
			CNodeAry& na = world.GetNA(id_na_val);
			assert( na.IsSegID(id_ns_val) );
			na.AddValueToNodeSegment(id_ns_val,*m_Update,1.0);
//			CNodeAry::CNodeSeg& ns = na.GetSeg(id_ns_val);
//			ns.is_updated = true;
		}
	}
	return true;
}