C++ DVec::size方法代码示例

// DPOSV uses Cholesky factorization A=U^T*U, A=L*L^T 
// to compute the solution to a real system of linear 
// equations A*X=B, where A is a square, (N,N) symmetric 
// positive definite matrix and X and B are (N,NRHS).
//---------------------------------------------------------
void umSOLVE_CH(const DMat& mat, const DVec& b, DVec& x)
//---------------------------------------------------------
{
  // check args
  assert(mat.is_square());            // symmetric
  assert(b.size() >= mat.num_rows()); // is b consistent?
  assert(b.size() <= x.size());       // can x store solution?
  
  DMat A(mat);    // work with copy of input
  x = b;          // allocate solution vector

  int rows=A.num_rows(), LDA=A.num_rows(), cols=A.num_cols();
  int  LDB=b.size(), NRHS=1, info=0;
  if (rows<1) {umWARNING("umSOLVE_CH()", "system is empty"); return;}

  // Solve the system.
  POSV('U', rows, NRHS, A.data(), LDA, x.data(), LDB, info);

  if (info < 0) { 
    x = 0.0;
    umERROR("umSOLVE_CH(A,b, x)", 
            "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    x = 0.0;
    umERROR("umSOLVE_CH(A,b, x)", 
            "\nINFO = %d.  The leading minor of order %d of A"
            "\nis not positive definite, so the factorization" 
            "\ncould not be completed. No solution computed.", 
              info, info);
  }
}

// access //////////////////////////////////////////////////////////////////////
void TabFunction::setData(const DVec &x, const DVec &y)
{
    if (x.size() != y.size())
    {
        LATAN_ERROR(Size, "tabulated function x/y data size mismatch");
    }
    FOR_VEC(x, i)
    {
        value_[x(i)] = y(i);
    }

//---------------------------------------------------------
void xyztorst
(
  const DVec& X,  // [in]
  const DVec& Y,  // [in]
  const DVec& Z,  // [in]
        DVec& r,  // [out]
        DVec& s,  // [out]
        DVec& t   // [out]
)
//---------------------------------------------------------
{
  // function [r,s,t] = xyztorst(x, y, z)
  // Purpose : Transfer from (x,y,z) in equilateral tetrahedron
  //           to (r,s,t) coordinates in standard tetrahedron

  double sqrt3=sqrt(3.0), sqrt6=sqrt(6.0); int Nc=X.size();
  DVec v1(3),v2(3),v3(3),v4(3);
  DMat tmat1(3,Nc), A(3,3), rhs; 
  
  v1(1)=(-1.0);  v1(2)=(-1.0/sqrt3);  v1(3)=(-1.0/sqrt6);
  v2(1)=( 1.0);  v2(2)=(-1.0/sqrt3);  v2(3)=(-1.0/sqrt6);
  v3(1)=( 0.0);  v3(2)=( 2.0/sqrt3);  v3(3)=(-1.0/sqrt6);
  v4(1)=( 0.0);  v4(2)=( 0.0      );  v4(3)=( 3.0/sqrt6);

  // back out right tet nodes
  tmat1.set_row(1,X); tmat1.set_row(2,Y); tmat1.set_row(3,Z);
  rhs = tmat1 - 0.5*outer(v2+v3+v4-v1, ones(Nc));
  A.set_col(1,0.5*(v2-v1)); A.set_col(2,0.5*(v3-v1)); A.set_col(3,0.5*(v4-v1));

  DMat RST = A|rhs;

  r=RST.get_row(1); s=RST.get_row(2); t=RST.get_row(3);
}

//==================================================================
void MemFile::InitExclusiveOwenership( DVec<U8> &fromData )
{
	mDataSize = fromData.size();
	mOwnData.get_ownership( fromData );

	mpData		= &mOwnData[0];
	mReadPos	= 0;
	mIsReadOnly	= true;
}

BasicStats::BasicStats(DVec &raw) {
  this->n = raw.size();
  double fcount = raw.size();
  DVec data(raw.begin(), raw.end());
  std::sort(data.begin(), data.end());
  this->min = data.front();
  this->max = data.back();
  this->median = data.at(data.size() / 2);

  double s = 0;
  for (decimal_t v : data)
    s += v;
  this->mean = s / fcount;

  s = 0;
  for (decimal_t v : data)
    s += ::pow(v - this->mean, 2);
  this->stdev = ::sqrt(s / fcount);
}

//---------------------------------------------------------
void CurvedINS2D::INScylinderBC2D
(
  const DVec&   xin,    // [in]
  const DVec&   yin,    // [in]
  const DVec&   nxi,    // [in]
  const DVec&   nyi,    // [in]
  const IVec&   MAPI,   // [in]
  const IVec&   MAPO,   // [in]
  const IVec&   MAPW,   // [in]
  const IVec&   MAPC,   // [in]
        double  ti,     // [in]
        double  nu,     // [in]
        DVec&   BCUX,   // [out]
        DVec&   BCUY,   // [out]
        DVec&   BCPR,   // [out]
        DVec&   BCDUNDT // [out]
)
//---------------------------------------------------------
{
  // function [bcUx, bcUy, bcPR, bcdUndt] = INScylinderBC2D(x, y, nx, ny, mapI, mapO, mapW, mapC, time, nu)
  // Purpose: evaluate boundary conditions for channel bounded cylinder flow with walls at y=+/- .15

  // TEST CASE: from 
  // V. John "Reference values for drag and lift of a two-dimensional time dependent flow around a cylinder", 
  // Int. J. Numer. Meth. Fluids 44, 777 - 788, 2004

  DVec yI("yI");  int len = xin.size();
  BCUX.resize(len); BCUY.resize(len);     // resize result arrays
  BCPR.resize(len); BCDUNDT.resize(len);  // and set to zero

  // inflow
#ifdef _MSC_VER
  yI = yin(MAPI);  yI += 0.20;
#else
  int Ni=MAPI.size(); yI.resize(Ni);
  for(int n=1;n<=Ni;++n) {yI(n)=yin(MAPI(n))+0.2;}
#endif

  BCUX(MAPI)    =  SQ(1.0/0.41)*6.0 * yI.dm(0.41 - yI);
  BCUY(MAPI)    =  0.0;
  BCDUNDT(MAPI) = -SQ(1.0/0.41)*6.0 * yI.dm(0.41 - yI);

  // wall
  BCUX(MAPW) = 0.0;
  BCUY(MAPW) = 0.0;

  // cylinder
  BCUX(MAPC) = 0.0;
  BCUY(MAPC) = 0.0;

  // outflow
  BCUX(MAPO)    = 0.0;
  BCUY(MAPO)    = 0.0;
  BCDUNDT(MAPO) = 0.0;
}

//---------------------------------------------------------
bool isInf(const DVec& V)
//---------------------------------------------------------
{
  for (int i=1; i<=V.size(); ++i) {

    double Vi = V(i);
    if (isinf(Vi)) {
      return true;    // v has a non-finite element
    }
  }
  return false;       // all elements are finite 
}

//---------------------------------------------------------
void VertexAngles
(
  const DVec& x1, const DVec& x2, const DVec& x3,
  const DVec& y1, const DVec& y2, const DVec& y3,
        DVec& a1,       DVec& a2,       DVec& a3
)
//---------------------------------------------------------
{

  //------------------------------------------
  // Expand definitions from ElmTools
  //------------------------------------------
  // a1 = acos ( -a23() / sqrt(a22()*a33()) );
  // a2 = acos ( -a13() / sqrt(a11()*a33()) );
  // a3 = acos ( -a12() / sqrt(a11()*a22()) );

  DVec g2x=(y3-y1), g2y=(x1-x3), g3x=(y1-y2), g3y=(x2-x1);
  
  DVec det = g3y*g2x - g3x*g2y;  
  DVec d   = 1.0/det;
  
  g2x *= d;  g2y *= d;  g3x *= d;  g3y *= d;

  DVec g1x =  - g2x - g3x;
  DVec g1y =  - g2y - g3y;

  a1 = acos( -(g2x*g3x + g2y*g3y) / sqrt((sqr(g2x)+sqr(g2y)) * (sqr(g3x)+sqr(g3y)) ));
  a2 = acos( -(g1x*g3x + g1y*g3y) / sqrt((sqr(g1x)+sqr(g1y)) * (sqr(g3x)+sqr(g3y)) ));
  a3 = acos( -(g1x*g2x + g1y*g2y) / sqrt((sqr(g1x)+sqr(g1y)) * (sqr(g2x)+sqr(g2y)) ));

#if (0)
  // check that the angles in each element sum to 180
  int Ni = x1.size(); double sum=0.0;
  umMSG(1, "\nChecking sum of angles in %d element\n", Ni);
  for (int i=1; i<=Ni; ++i) {
    sum = fabs(a1(i)) + fabs(a2(i)) + fabs(a3(i));
    if ( fabs(sum-M_PI) > 1e-15) {
      umMSG(1, "element %4d: %12.5e\n", i, fabs(sum-M_PI)); 
    }
  }
#endif
}

//---------------------------------------------------------
void umPOLISH(DVec& V, double eps)
//---------------------------------------------------------
{
  // round elements close to certain values

  int N = V.size();
  double *p = V.data();

  for (int i=0; i<N; ++i) 
  {
    if (fabs(p[i]) < eps) 
    {
      p[i] = 0.0;
    }
    else
    {
      if (p[i] > 0.0) 
      {
        // check for proximity to certain positive values
        if      (fabs (p[i] - 0.10) < eps) { p[i] = 0.10; }
        else if (fabs (p[i] - 0.20) < eps) { p[i] = 0.20; }
        else if (fabs (p[i] - 0.25) < eps) { p[i] = 0.25; }
        else if (fabs (p[i] - 0.50) < eps) { p[i] = 0.50; }
        else if (fabs (p[i] - 0.75) < eps) { p[i] = 0.75; }
        else if (fabs (p[i] - 0.80) < eps) { p[i] = 0.80; }
        else if (fabs (p[i] - 0.90) < eps) { p[i] = 0.90; }
        else if (fabs (p[i] - 1.00) < eps) { p[i] = 1.00; }
        else if (fabs (p[i] - 2.00) < eps) { p[i] = 2.00; }
        else if (fabs (p[i] - 4.00) < eps) { p[i] = 4.00; }
        else if (fabs (p[i] - 4.50) < eps) { p[i] = 4.50; }
        else if (fabs (p[i] - 5.00) < eps) { p[i] = 5.00; }

        else if (fabs (p[i] - M_PI  ) < eps) { p[i] = M_PI  ; }
        else if (fabs (p[i] - M_PI_2) < eps) { p[i] = M_PI_2; }
        else if (fabs (p[i] - M_PI_4) < eps) { p[i] = M_PI_4; }
        else if (fabs (p[i] - M_E   ) < eps) { p[i] = M_E   ; }
      }
      else
      {
        // check for proximity to certain negative values
        if      (fabs (p[i] + 0.10) < eps) { p[i] = -0.10; }
        else if (fabs (p[i] + 0.20) < eps) { p[i] = -0.20; }
        else if (fabs (p[i] + 0.25) < eps) { p[i] = -0.25; }
        else if (fabs (p[i] + 0.50) < eps) { p[i] = -0.50; }
        else if (fabs (p[i] + 0.75) < eps) { p[i] = -0.75; }
        else if (fabs (p[i] + 0.80) < eps) { p[i] = -0.80; }
        else if (fabs (p[i] + 0.90) < eps) { p[i] = -0.90; }
        else if (fabs (p[i] + 1.00) < eps) { p[i] = -1.00; }
        else if (fabs (p[i] + 2.00) < eps) { p[i] = -2.00; }
        else if (fabs (p[i] + 4.00) < eps) { p[i] = -4.00; }
        else if (fabs (p[i] + 4.50) < eps) { p[i] = -4.50; }
        else if (fabs (p[i] + 5.00) < eps) { p[i] = -5.00; }

        else if (fabs (p[i] + M_PI  ) < eps) { p[i] = -M_PI  ; }
        else if (fabs (p[i] + M_PI_2) < eps) { p[i] = -M_PI_2; }
        else if (fabs (p[i] + M_PI_4) < eps) { p[i] = -M_PI_4; }
        else if (fabs (p[i] + M_E   ) < eps) { p[i] = -M_E   ; }
      }
    }
  }
}

int OCCEdge::createNURBS(OCCVertex *start, OCCVertex *end, std::vector<OCCStruct3d> points,
                          DVec knots, DVec weights, IVec mult)
{
    try {
        Standard_Boolean periodic = false;
        
        int vertices = 0;
        if (start != NULL && end != NULL) {
            vertices = 2;
            periodic = true;
        }
        
        int nbControlPoints = points.size() + vertices;
        TColgp_Array1OfPnt  ctrlPoints(1, nbControlPoints);
        
        TColStd_Array1OfReal _knots(1, knots.size());
        TColStd_Array1OfReal _weights(1, weights.size());
        TColStd_Array1OfInteger  _mult(1, mult.size());
        
        for (unsigned i = 0; i < knots.size(); i++) {
            _knots.SetValue(i+1, knots[i]);
        }
        
        for (unsigned i = 0; i < weights.size(); i++) {
            _weights.SetValue(i+1, weights[i]);
        }
        
        int totKnots = 0;
        for (unsigned i = 0; i < mult.size(); i++) {
            _mult.SetValue(i+1, mult[i]);   
            totKnots += mult[i];
        }

        const int degree = totKnots - nbControlPoints - 1;

        int index = 1;
        
        if (!periodic) {
            ctrlPoints.SetValue(index++, gp_Pnt(start->X(), start->Y(), start->Z()));
        }
        
        for (unsigned i = 0; i < points.size(); i++) {
            gp_Pnt aP(points[i].x,points[i].y,points[i].z);
            ctrlPoints.SetValue(index++, aP);
        }
        
        if (!periodic) {
            ctrlPoints.SetValue(index++, gp_Pnt(end->X(), end->Y(), end->Z()));
        }
        
        Handle(Geom_BSplineCurve) NURBS = new Geom_BSplineCurve
        (ctrlPoints, _weights, _knots, _mult, degree, periodic);
        
        if (!periodic) {
            this->setShape(BRepBuilderAPI_MakeEdge(NURBS, start->vertex, end->vertex));
        } else {
            this->setShape(BRepBuilderAPI_MakeEdge(NURBS));
        }
    } catch(Standard_Failure &err) {
        Handle_Standard_Failure e = Standard_Failure::Caught();
        const Standard_CString msg = e->GetMessageString();
        if (msg != NULL && strlen(msg) > 1) {
            setErrorMessage(msg);
        } else {
            setErrorMessage("Failed to create nurbs");
        }
        return 1;
    }
    return 0;
}

// function call ///////////////////////////////////////////////////////////////
double DoubleModel::operator()(const DVec &arg, const DVec &par) const
{
    checkSize(arg.size(), par.size());
    
    return (*this)(arg.data(), par.data());
}

//---------------------------------------------------------
void CurvedINS2D::KovasznayBC2D
(
  const DVec&   xin,    // [in]
  const DVec&   yin,    // [in]
  const DVec&   nxi,    // [in]
  const DVec&   nyi,    // [in]
  const IVec&   MAPI,   // [in]
  const IVec&   MAPO,   // [in]
  const IVec&   MAPW,   // [in]
  const IVec&   MAPC,   // [in]
        double  ti,     // [in]
        double  nu,     // [in]
        DVec&   BCUX,   // [out]
        DVec&   BCUY,   // [out]
        DVec&   BCPR,   // [out]
        DVec&   BCDUNDT // [out]
)
//---------------------------------------------------------
{
  // function [bcUx, bcUy, bcPR, bcdUndt] = KovasznayBC2D(x, y, nx, ny, MAPI, MAPO, MAPW, MAPC, time, nu)
  // Purpose: evaluate boundary conditions for Kovasznay flow 

  static DVec xI("xI"), yI("yI"), xO("xO"), yO("yO");

  int len = xin.size();
  BCUX.resize(len); BCUY.resize(len);     // resize result arrays
  BCPR.resize(len); BCDUNDT.resize(len);  // and set to zero

  double lam = (0.5/nu) - sqrt( (0.25/SQ(nu)) + 4.0*SQ(pi));

  // inflow
#ifdef _MSC_VER
  xI = xin(MAPI);  yI = yin(MAPI);
#else
  int Ni=MAPI.size(), n=0; xI.resize(Ni); yI.resize(Ni);
  for (n=1;n<=Ni;++n) {xI(n)=xin(MAPI(n));}
  for (n=1;n<=Ni;++n) {yI(n)=yin(MAPI(n));}
#endif

  DVec elamXI = exp(lam*xI), twopiyI = 2.0*pi*yI;

  BCUX(MAPI) =          1.0 - elamXI.dm(cos(twopiyI));
  BCUY(MAPI) = (0.5*lam/pi) * elamXI.dm(sin(twopiyI));


  // outflow
#ifdef _MSC_VER
  xO = xin(MAPO); yO = yin(MAPO);
#else
  int No=MAPO.size(); xO.resize(No); yO.resize(No);
  for (n=1;n<=No;++n) {xO(n)=xin(MAPO(n));}
  for (n=1;n<=No;++n) {yO(n)=yin(MAPO(n));}
#endif

  DVec elamXO = exp(lam*xO), twopiyO = 2.0*pi*yO;

  BCPR   (MAPO) =  0.5*(1.0-exp(2.0*lam*xO));
//BCDUNDT(MAPO) = -lam*elamXO.dm(cos(twopiyO));

  if (0) {
    // THIS WORKED
    BCUX(MAPO) =         1.0  - elamXO.dm(cos(twopiyO));
    BCUY(MAPO) = (0.5*lam/pi) * elamXO.dm(sin(twopiyO));
  } else {
    // Neumann data for each velocity
    BCUX(MAPO) = -lam*               elamXO.dm(cos(twopiyO));
    BCUY(MAPO) =  lam*(0.5*lam/pi) * elamXO.dm(sin(twopiyO));
  }
}

//==================================================================
void MakeTree( TokNode *pRoot, DVec<Token> &tokens, u_int &out_blockCnt )
{
	//TokNode	*pParent = pRoot;
	//DVec<TokNode*>	pParentsMemory;
	//pParentsMemory.push_back( pParent );
	//pPrev = pParent = pParent->AddNewChild( &tokens[i] );
	//pParentsMemory.pop_back();

	out_blockCnt = 0;

	TokNode			*pNode = pRoot;
	DVec<TokenID>	bracketsMemory;

	for (size_t i=0; i < tokens.size(); ++i)
	{
		switch ( tokens[i].id )
		{
		case RSLC::T_OP_LFT_BRACKET		:
		case RSLC::T_OP_LFT_CRL_BRACKET	:
		case RSLC::T_OP_LFT_SQ_BRACKET	:
			bracketsMemory.push_back( tokens[i].id );

			if ( tokens[i].id == RSLC::T_OP_LFT_SQ_BRACKET )
			{
				pNode->AddNewChild( &tokens[i] );

				Token	*pNewBrkTok =
					DNEW Token(
							"(",
							RSLC::T_OP_LFT_BRACKET,
							RSLC::T_TYPE_OPERATOR,
							&tokens[i] );

				pNode = pNode->AddNewChild( pNewBrkTok );

				//pNode = pNewParent;
			}
			else
			// square bracket is indents itself
			//if ( tokens[i].id == RSLC::T_OP_LFT_SQ_BRACKET )
			//{		
			//	if NOT( pNode->mpChilds.size() )
			//		throw Exception( "Misplaced square bracket ?", &tokens[i] );

			//	pNode = pNode->mpChilds[pNode->mpChilds.size()-1]->AddNewChild( &tokens[i] );
			//}
			//else
			{
				pNode = pNode->AddNewChild( &tokens[i] );
			}
			break;

		case RSLC::T_OP_RGT_BRACKET		:
		case RSLC::T_OP_RGT_CRL_BRACKET	:
		case RSLC::T_OP_RGT_SQ_BRACKET	:
			if (!pNode->mpParent ||
				!bracketsMemory.size() ||
				!doBracketsMatch( bracketsMemory.back(), tokens[i].id ) )
			{
				throw Exception( "Mismatched brackets ?", &tokens[i] );
			}

			pNode = pNode->mpParent;
			// square bracket is indented itself
			//if ( tokens[i].id == RSLC::T_OP_RGT_SQ_BRACKET )
			//	pNode = pNode->mpParent;

			bracketsMemory.pop_back();
			break;

		default:
			pNode->AddNewChild( &tokens[i] );
			break;
		}
	}

	defineBlockTypeAndID( pRoot, out_blockCnt );
}