C++ DMat::pCol方法代码示例

//---------------------------------------------------------
void Euler3D::CouetteBC3D
(
    const DVec& xi,
    const DVec& yi,
    const DVec& zi,
    const DVec& nxi,
    const DVec& nyi,
    const DVec& nzi,
    const IVec& tmapI,
    const IVec& tmapO,
    const IVec& tmapW,
    const IVec& tmapC,
    double ti,
    DMat& Qio
)
//---------------------------------------------------------
{
    //#######################################################
    // FIXME: what about top and bottom of 3D annulus?
    //#######################################################


    // gmapB = concat(mapI,mapO,mapW, ...);
    // Check for no boundary faces
    if (gmapB.size() < 1) {
        return;
    }

    // function [Q] = CouetteBC3D(xin, yin, nxin, nyin, mapI, mapO, mapW, mapC, Q, time);
    // Purpose: evaluate solution for Couette flow

    // Couette flow (mach .2 at inner cylinder)
    // gamma = 1.4;

    int Nr = Qio.num_rows();
    DVec rho,rhou,rhov,rhow,Ener;

    // wrap current state data (columns of Qio)
    rho.borrow (Nr, Qio.pCol(1));
    rhou.borrow(Nr, Qio.pCol(2));
    rhov.borrow(Nr, Qio.pCol(3));
    rhow.borrow(Nr, Qio.pCol(4));
    Ener.borrow(Nr, Qio.pCol(5));

    // update boundary nodes of Qio with boundary data
    // pre-calculated in function precalc_bdry_data()

    rho (gmapB) = rhoB;
    rhou(gmapB) = rhouB;
    rhov(gmapB) = rhovB;
    rhow(gmapB) = rhowB;
    Ener(gmapB) = EnerB;
}

//---------------------------------------------------------
void CurvedCNS2D::BoxFlowIC2D
(
  const DVec&   xi,   // [in]
  const DVec&   yi,   // [in]
        double  ti,   // [in]
        DMat&   Qo    // [out]
)
//---------------------------------------------------------
{
  // function Q = BoxFlowIC2D(x, y, time)
  //
  // Purpose: compute plane flow configuration 

  //-------------------------------------
  // adjust parameters
  //-------------------------------------
  this->gamma = 1.4;
  this->gm1   = 0.4;  // (gamma-1)

  //-------------------------------------
  // use wrappers to update Qo in-place
  //-------------------------------------
  DVec rho,rhou,rhov,Ener; int Nr=Np*K;
  rho.borrow (Nr, Qo.pCol(1));
  rhou.borrow(Nr, Qo.pCol(2));
  rhov.borrow(Nr, Qo.pCol(3));
  Ener.borrow(Nr, Qo.pCol(4));

  if (1) 
  {
    this->pref = 12.0;

    rho  =  1.0;
    rhou = -sin(2.0*pi*y);
    rhov =  sin(4.0*pi*x);
    Ener = pref/gm1 + 0.5*(sqr(rhou) + sqr(rhov)).dd(rho);
  } 
  else 
  {
    this->pref = 12.0;

    rho  = 1.0;
    rhou = 0.0;
    rhov = 0.0;
    Ener = pref/gm1 + 0.5 * rho.dm(exp(-4.0*(sqr(cos(pi*x))+sqr(cos(pi*y)))));
  }
}

//---------------------------------------------------------
void NDG2D::OutputSampleXYZ
(
        int sample_N,
        DMat &newX, 
        DMat &newY, 
        DMat &newZ,     // e.g. triangles on a sphere
  const DMat &FData,    // old field data
        DMat &newFData, // new field data
        int zfield      // if>0, use as z-elevation
)
//---------------------------------------------------------
{
  DVec newR, newS, newT;
  DMat newVDM;
  int newNpts = 0;

  // Triangles
  OutputSampleNodes2D(sample_N, newR, newS);
  newNpts = newR.size();
  newVDM = Vandermonde2D(this->N, newR, newS);

  const DMat& oldV = this->V;
  DMat oldtonew(newNpts, this->Np, "OldToNew");
  oldtonew = trans(trans(oldV) | trans(newVDM));

  //-----------------------------------
  // interpolate the field data
  //-----------------------------------
  int Nfields = FData.num_cols();
  newFData.resize(newNpts*this->K, Nfields);
  //DVec scales(Nfields);

  // For each field, use tOldF to wrap field i.
  // Use tNewF to load the interpolated field
  // directly into column i of the output array.
  DMat tOldF, tNewF;
  for (int i=1; i<=Nfields; ++i) {
    tOldF.borrow(this->Np, this->K, (double*)   FData.pCol(i));
    tNewF.borrow(newNpts,  this->K, (double*)newFData.pCol(i));
    tNewF = oldtonew * tOldF;
  //scales(i) = tNewF.max_col_val_abs(i);
  }

  //-----------------------------------
  // interpolate the vertices
  //-----------------------------------
  newX = oldtonew * this->x;
  newY = oldtonew * this->y;

  if (this->bCoord3D) {
    newZ = oldtonew * this->z;
  } 
  else 
  {
    if (zfield>=1 && zfield<=Nfields) {
      // use field data for z-height
      newZ.load(newNpts, K, newFData.pCol(Nfields));
    } else {
      // set z-data to 0.0
      newZ.resize(newNpts, K, true, 0.0);
    }
  }
}