C++ MultiIndex::firstOrderDirection方法代码示例

Expr SpectralPreprocessor::takeDeriv(const Expr& f, const MultiIndex& mi)
{
  TEUCHOS_TEST_FOR_EXCEPT(mi.order() > 1);

  
  const SpectralExpr* se 
    = dynamic_cast<const SpectralExpr*>(f.ptr().get());

  Expr d = new Derivative(mi.firstOrderDirection());
  int n = se->getSpectralBasis().nterms();
  
  if (se)
  {
    Array<Expr> c(n);
    for (int i=0; i<n; i++)
    {
      c[i] = d*se->getCoeff(i);
    }
    return new SpectralExpr(se->getSpectralBasis(), c);
  }
  else
  {
    return d*f;
  }
}

void EdgeLocalizedBasis::evalOnTriangle(const Point& pt, 
  const MultiIndex& deriv,
  Array<double>& result) const
{
  ADReal x = ADReal(pt[0], 0, 2);
	ADReal y = ADReal(pt[1], 1, 2);
	ADReal one(1.0, 2);
	ADReal zero(0.0, 2);

  Array<ADReal> tmp;

  SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
    << y.value());

  result.resize(3);
  tmp.resize(3);

  bool onEdge2 = std::fabs(pt[1]) < 1.0e-14;
  bool onEdge0 = std::fabs(1.0-pt[0]-pt[1]) < 1.0e-14;
  bool onEdge1 = std::fabs(pt[0]) < 1.0e-14;
  
  TEUCHOS_TEST_FOR_EXCEPTION(!(onEdge0 || onEdge1 || onEdge2),
    std::runtime_error,
    "EdgeLocalizedBasis should not be evaluated at points not on edges");
  
  TEUCHOS_TEST_FOR_EXCEPTION((onEdge0 && onEdge1) || (onEdge1 && onEdge2)
    || (onEdge2 && onEdge0), std::runtime_error,
    "Ambiguous edge in EdgeLocalizedBasis::evalOnTriangle()");

  if (onEdge0)
  {
    tmp[0] = one;
    tmp[1] = zero;
    tmp[2] = zero;
  }
  if (onEdge1)
  {
    tmp[0] = zero;
    tmp[1] = one;
    tmp[2] = zero;
  }
  if (onEdge2)
  {
    tmp[0] = zero;
    tmp[1] = zero;
    tmp[2] = one;
  }


	for (int i=0; i<tmp.length(); i++)
  {
    SUNDANCE_OUT(this->verb() > 3,
      "tmp[" << i << "]=" << tmp[i].value() 
      << " grad=" << tmp[i].gradient());
    if (deriv.order()==0) result[i] = tmp[i].value();
    else 
      result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
  }
  
}

void CubicHermite::evalOnTriangle(const Point& pt, 
  const MultiIndex& deriv,
  Array<double>& result) const

{
  result.resize(10);
  ADReal x = ADReal(pt[0], 0, 2);
  ADReal y = ADReal(pt[1], 1, 2);
  ADReal one(1.0, 2);
  
  Array<ADReal> tmp(10);

  SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
    << y.value());
 
  tmp[0] = 1 - 3*x*x + 2 * x*x*x - 13*x*y + 13*x*x*y - 3*y*y + 13 *x*y*y + 2 *y*y*y;
  tmp[1] = x - 2 *x*x + x*x*x - 3*x*y + 3*x*x*y + 2*x*y*y;
  tmp[2] = y - 3 *x *y + 2* x*x* y - 2* y*y + 3* x*y*y + y*y*y;
  tmp[3] = 3 * x*x - 2* x*x*x - 7* x* y + 7* x*x *y + 7* x*y*y;
  tmp[4] = -x*x + x*x*x + 2*x *y - 2* x*x* y - 2* x* y*y;
  tmp[5] = -x* y + 2* x*x* y + x* y*y;
  tmp[6] = -7* x* y + 7* x*x*y + 3* y*y + 7* x* y*y - 2* y*y*y;
  tmp[7] = -x *y + x*x* y + 2* x* y*y;
  tmp[8] = 2 *x *y - 2* x*x* y - y*y - 2* x* y*y + y*y*y;
  tmp[9] = 27* x *y - 27* x*x* y - 27* x* y*y;

  for (int i=0; i<tmp.length(); i++)
  {
    if (deriv.order()==0) result[i] = tmp[i].value();
    else 
      result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
  }
}

void Bernstein::evalOnTet(const Point& pt, 
  const MultiIndex& deriv,
  Array<double>& result) const
{
  ADReal x = ADReal(pt[0], 0, 3);
  ADReal y = ADReal(pt[1], 1, 3);
  ADReal z = ADReal(pt[2], 2, 3);
  ADReal one(1.0, 3);
  
  Array<ADReal> tmp(result.length());

  if(order_==0)
  {
    tmp.resize(1);
    result.resize(1);
    tmp[0] = one;
  }
  else
  {
  }

  for (int i=0; i<tmp.length(); i++)
  {
    if (deriv.order()==0) result[i] = tmp[i].value();
    else 
      result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
  }
}

Set<MultipleDeriv> Xx(const MultiIndex& x)
{
  Set<MultipleDeriv> rtn;

  TEUCHOS_TEST_FOR_EXCEPTION(x.order() < 0 || x.order() > 1, std::logic_error,
    "invalid multiindex " << x << " in this context");

  MultipleDeriv xmd = makeMultiDeriv(coordDeriv(x.firstOrderDirection()));
  rtn.put(xmd);
  return rtn;
}

void EdgeLocalizedBasis::evalOnLine(const Point& pt, 
													const MultiIndex& deriv,
													Array<double>& result) const
{
	ADReal one(1.0, 1);
	result.resize(1);
	Array<ADReal> tmp(result.length());

  tmp[0] = one;

	for (int i=0; i<tmp.length(); i++)
		{
			if (deriv.order()==0) result[i] = tmp[i].value();
			else result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
		}
}

double TestEvalMediator::evalDummyBasis(int m, const MultiIndex& mi) const
{
  TEUCHOS_TEST_FOR_EXCEPTION(mi.order() > 1, std::runtime_error, 
                     "TestEvalMediator::evalDummyBasis found multiindex "
                     "order > 1. The bad multiindex was " << mi.toString());

  ADReal result = fields_[m].basis().evaluate(ADField::evalPoint());
  SUNDANCE_MSG3(verb(), "basis.value() " << result.value());
  SUNDANCE_MSG3(verb(), "basis.gradient() " << result.gradient());

  if (mi.order()==0)
    {
      return result.value();
    }
  else 
    {
      return result.gradient()[mi.firstOrderDirection()];
    }
}

void Legendre::evalOnQuad(const Point& pt,
  const MultiIndex& deriv,
  Array<double>& result) const

{
  result.resize( 4 + 4*nrDOF_edge_ + nrDOF_face_);
  ADReal x = ADReal(pt[0], 0, 2);
  ADReal y = ADReal(pt[1], 1, 2);
  ADReal one(1.0, 2);
  
  Array<ADReal> refAllx(7);
  Array<ADReal> refAlly(7);

  refAllx[0] = 1-x;
  refAllx[1] = x;
  refAllx[2] = 2.44948974278318 * ( (2*x-1)*(2*x-1) - 1 ) / 4;
  refAllx[3] = 3.16227766016838 * ( (2*x-1)*(2*x-1) - 1 ) * (2*x-1) / 4;
  refAllx[4] = 3.74165738677394 * ( 5*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 6*(2*x-1)*(2*x-1) + 1) / 16;
  refAllx[5] = 4.24264068711929 * (2*x-1) * (7*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 10*(2*x-1)*(2*x-1) + 3) / 16;
  refAllx[6] = 4.69041575982343 * (21*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) -
		                          35*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) + 15*(2*x-1)*(2*x-1) - 1) / 32;

  refAlly[0] = 1-y;
  refAlly[1] = y;
  refAlly[2] = 2.44948974278318 * ( (2*y-1)*(2*y-1) - 1 ) / 4;
  refAlly[3] = 3.16227766016838 * ( (2*y-1)*(2*y-1) - 1 ) * (2*y-1) / 4;
  refAlly[4] = 3.74165738677394 * ( 5*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) - 6*(2*y-1)*(2*y-1) + 1) / 16;
  refAlly[5] = 4.24264068711929 * (2*y-1) * (7*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) - 10*(2*y-1)*(2*y-1) + 3) / 16;
  refAlly[6] = 4.69041575982343 * (21*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) -
		                          35*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) + 15*(2*y-1)*(2*y-1) - 1) / 32;

  SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
    << y.value());

  int sideIndex[4][2] = { {0,0} , {1,0} , {0,1} , {1,1}};
  int edgeIndex[4]    = { 0 , 1 , 1 , 0};
  int edgeMultI[4]    = { 0 , 0 , 1 , 1};
  int offs = 0;
  Array<ADReal> tmp(4 + 4*nrDOF_edge_ + nrDOF_face_);

  // loop over vertexes
  for (int i=0; i < 4 ; i++, offs++){
	  tmp[offs] = refAllx[sideIndex[i][0]] * refAlly[sideIndex[i][1]];
  }

  // loop over edges
  for (int i=0; i < 4 ; i++){
	  // loop over each DOF on the edge
	  if (edgeIndex[i] == 0){
		  for (int j = 0 ; j < nrDOF_edge_ ; j++ , offs++){
			  tmp[offs] = refAllx[2+j] * refAlly[edgeMultI[i]];
		  }
	  }
	  else
	  {
		  for (int j = 0 ; j < nrDOF_edge_ ; j++ , offs++){
			  tmp[offs] = refAllx[edgeMultI[i]] * refAlly[2+j];
		  }
	  }
  }

  // loop over all internal DOFs
  if ( nrDOF_face_ > 0 ){
	  // loop for each hierarchical layer
	  for (int hierarch = 0 ; hierarch < order_ - 3 ; hierarch++)
	  {
		  //SUNDANCE_OUT( true , "Legendre::evalOnQuad hierarch:" << hierarch );
		  // for each layer add the basis function
		  for (int i=0 ; i < hierarch+1 ; i++ , offs++)
		  {
			  //SUNDANCE_OUT( true , "Legendre::evalOnQuad offs:" << offs << " 2+i:" << 2+i << " , 2+(hierarch-1-i):" << 2+(hierarch-i));
			  tmp[offs] = refAllx[2+i] * refAlly[2+(hierarch-i)];
		  }
	  }
  }

  // compute the results
  for (int i=0; i<result.length(); i++)
  {
    if (deriv.order()==0) result[i] = tmp[i].value();
    else 
      result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
  }
  //SUNDANCE_OUT( true , "Legendre::evalOnQuad result.length():" << result.length() );
}

void Bernstein::evalOnTriangle(const Point& pt, 
  const MultiIndex& deriv,
  Array<double>& result) const

{
  ADReal x = ADReal(pt[0], 0, 2);
  ADReal y = ADReal(pt[1], 1, 2);
  ADReal one(1.0, 2);
  
  Array<ADReal> tmp;

  SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
    << y.value());
 
  if(order_==0) {
    result.resize(1);
    tmp.resize(1);
    tmp[0] = one;
  }
  else {
    int N = (order()+1)*(order()+2)/2;
    result.resize(N);
    tmp.resize(N);
    // these are the barycentric coordinates
    ADReal b1 = 1.0 - x - y;
    ADReal b2 = x;
    ADReal b3 = y;

    // will hold \binom{n}{\alpha_1}
    int bfcur = 0;

    for (int alpha1=order();alpha1>=0;alpha1--) 
    {
      for (int alpha2 = order()-alpha1;alpha2>=0;alpha2--)
      {
        int alpha3 = order() - alpha1 - alpha2;
        tmp[bfcur] = one;
        for (int i=0;i<alpha1;i++)
	      {
          tmp[bfcur] *= b1;
	      }
        for (int i=0;i<alpha2;i++) 
	      {
          tmp[bfcur] *= b2;
	      }
        for (int i=0;i<alpha3;i++)
	      {
          tmp[bfcur] *= b3;
	      }
        for (int i=2;i<=order();i++)
	      {
          tmp[bfcur] *= (double) i;
	      }
        for (int i=2;i<=alpha1;i++) 
	      {
          tmp[bfcur] /= (double) i;
	      }
        for (int i=2;i<=alpha2;i++) 
	      {
          tmp[bfcur] /= (double) i;
	      }
        for (int i=2;i<=alpha3;i++) 
	      {
          tmp[bfcur] /= (double) i;
	      }
        bfcur++;
      }
    }
  }

  for (int i=0; i<tmp.length(); i++)
  {
    if (deriv.order()==0) result[i] = tmp[i].value();
    else 
      result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
  }
}

DiffOp::DiffOp(const MultiIndex& op, 
  const RCP<ScalarExpr>& arg)
  : UnaryExpr(arg), mi_(op), myCoordDeriv_(coordDeriv(op.firstOrderDirection())), requiredFunctions_(),
    ignoreFuncTerms_(false)
{}