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

//---------------------------------------------------------
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;
}

// load {R,S} output nodes
void OutputSampleNodes2D(int sample_N, DVec &R, DVec &S)
{
  int Npts = OutputSampleNpts2D(sample_N);
  R.resize(Npts); S.resize(Npts);
  double denom = (double)(sample_N);
  int sampleNq = sample_N+1;
  for (int sk=0, i=0; i<sampleNq; ++i) {
    for (int j=0; j<sampleNq-i; ++j, ++sk) {
      R[sk] = -1. + (2.*j)/denom;
      S[sk] = -1. + (2.*i)/denom;
    }
  }
}

// compute eigensystem of a real symmetric matrix
//---------------------------------------------------------
void eig_sym(const DMat& A, DVec& ev, DMat& Q, bool bDoEVecs)
//---------------------------------------------------------
{
  if (!A.is_square()) { umERROR("eig_sym(A)", "matrix is not square."); }

  int N = A.num_rows();
  int LDA=N, LDVL=N, LDVR=N, ldwork=10*N, info=0;
  DVec work(ldwork, 0.0, OBJ_temp, "work_TMP");

  Q = A;          // Calculate eigenvectors in Q (optional)
  ev.resize(N);   // Calculate eigenvalues in ev

  char jobV = bDoEVecs ? 'V' : 'N';

  SYEV (jobV,'U', N, Q.data(), LDA, ev.data(), work.data(), ldwork, info);  

  if (info < 0) { 
    umERROR("eig_sym(A, Re,Im)", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    umLOG(1, "eig_sym(A, W): ...\n"
             "\nthe algorithm failed to converge;"
             "\n%d off-diagonal elements of an intermediate"
             "\ntridiagonal form did not converge to zero.\n", info);
  }
}

//==================================================================
bool ParamsFindP(	ParamList &params,
					const SymbolList &globalSymbols,
					DVec<Float3> &out_vectorP,
					int fromIdx )
{
	bool	gotP = false;

	for (size_t i=fromIdx; i < params.size(); i += 2)
	{
		DASSERT( params[i].type == Param::STR );

		const Symbol* pSymbol = globalSymbols.FindSymbol( params[i] );
		if ( pSymbol && pSymbol->IsName( "P" ) )
		{
			DASSTHROW( (i+1) < params.size(), "Invalid number of arguments" );
			
			const FltVec	&fltVec = params[ i+1 ].NumVec();
			
			DASSTHROW( (fltVec.size() % 3) == 0, "Invalid number of arguments" );
			
			out_vectorP.resize( fltVec.size() / 3 );

			for (size_t iv=0, id=0; iv < fltVec.size(); iv += 3)
				out_vectorP[id++] = Float3( &fltVec[ iv ] );

			return true;
		}
	}
	
	return false;
}

//---------------------------------------------------------
void eig(const DMat& A, DVec& Re, DMat& VL, DMat& VR, bool bL, bool bR)
//---------------------------------------------------------
{
  // Compute eigensystem of a real general matrix
  // Currently NOT returning imaginary components

  static DMat B;

  if (!A.is_square()) { umERROR("eig(A)", "matrix is not square."); }

  int N = A.num_rows();
  int LDA=N, LDVL=N, LDVR=N, ldwork=10*N, info=0;

  Re.resize(N);     // store REAL components of eigenvalues in Re
  VL.resize(N,N);   // storage for LEFT eigenvectors
  VR.resize(N,N);   // storage for RIGHT eigenvectors
  DVec Im(N);     // NOT returning imaginary components
  DVec work(ldwork, 0.0);

  // Work on a copy of A
  B = A;

  char jobL = bL ? 'V' : 'N';   // calc LEFT eigenvectors?
  char jobR = bR ? 'V' : 'N';   // calc RIGHT eigenvectors?

  GEEV (jobL,jobR, N, B.data(), LDA, Re.data(), Im.data(), 
        VL.data(), LDVL, VR.data(), LDVR, work.data(), ldwork, info);

  if (info < 0) { 
    umERROR("eig(A, Re,Im)", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    umLOG(1, "eig(A, Re,Im): ...\n"
             "\nThe QR algorithm failed to compute all the"
             "\neigenvalues, and no eigenvectors have been" 
             "\ncomputed;  elements %d+1:N of WR and WI contain"
             "\neigenvalues which have converged.\n", info);
  }

#if (0)
  // Return (Re,Imag) parts of eigenvalues as columns of Ev
  Ev.resize(N,2);
  Ev.set_col(1, Re);
  Ev.set_col(2, Im);
#endif

#ifdef _DEBUG
    //#####################################################
    // Check for imaginary components in eigenvalues
    //#####################################################
    double im_max = Im.max_val_abs();
    if (im_max > 1e-6) {
      umERROR("eig(A)", "imaginary components in eigenvalues.");
    }
    //#####################################################
#endif
}

//==================================================================
bool FileManagerAndroid::GrabFile(const char* pFileName, DVec<U8> &out_data, const char* pMode)
{
    DFUNCTION();
    DLOG("File: %s", pFileName);

#if defined(ANDROID)
	bool prefs = isPrefsMode( pMode );

    if (prefs)
    {
        // Read from the prefs storage

        DLOG("Grabbing file (prefs): %s", pFileName);
        char *contents = ::DoReadPrefs(pFileName);
        if (0 == contents)
        {
            DLOG("Failed to open the file");
            return false;
        }

        // Decode into binary

        size_t binLength = 0;
        void *binData = FromBase64(contents, &binLength);
        free(contents);
        if (binData == 0)
        {
            DPRINT("!! File '%s' was not base64 format.  Rejecting.", pFileName);
            return false;
        }

        // Copy into output array and free original

        // TODO: Could be sped up by pre-calculating length and not
        // copying.

        out_data.resize(binLength);
        memcpy(&out_data[0], binData, binLength);
        free(binData);

        DLOG("copied %d bytes", binLength);
        return true;
    }

    ::FileLoad(pFileName);
    int datasize;
    void *data;
    if (!FileGetData(pFileName, &data, &datasize))
    {
        //DASSTHROW(0, ("failed loading file '%s'", pFileName));
        return false;
    }

    out_data.resize(datasize);

    DVERBOSE("GrabFile: memcpy-ing '%s', %d bytes", pFileName, datasize);
    memcpy(&out_data[0], data, datasize);

    FileRemove(pFileName);
    DVERBOSE("GrabFile: released");
#endif

	return true;
}

//---------------------------------------------------------
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));
  }
}