C++ Array::pop_back方法代码示例

FCPtr BasisEvaluation::getTransformedValuesWithBasisValues(BasisPtr basis, Camellia::EOperator op,
                                                           constFCPtr referenceValues, int numCells,
                                                           BasisCache* basisCache)
{
  typedef FunctionSpaceTools fst;
  //  int numCells = cellJacobian.dimension(0);
  
  int spaceDim = basisCache->getSpaceDim(); // changed 3-21-16
  
  int componentOfInterest;
  Camellia::EFunctionSpace fs = basis->functionSpace();
  Intrepid::EOperator relatedOp = relatedOperator(op,fs,spaceDim, componentOfInterest);
  Teuchos::Array<int> dimensions;
  referenceValues->dimensions(dimensions);
  dimensions.insert(dimensions.begin(), numCells);
  Teuchos::RCP<FieldContainer<double> > transformedValues = Teuchos::rcp(new FieldContainer<double>(dimensions));
  bool vectorizedBasis = functionSpaceIsVectorized(fs);
  if (vectorizedBasis && (op ==  Camellia::OP_VALUE))
  {
    TEUCHOS_TEST_FOR_EXCEPTION( vectorizedBasis && (op ==  Camellia::OP_VALUE),
                               std::invalid_argument, "Vector HGRAD/HVOL with OP_VALUE not supported by getTransformedValuesWithBasisValues.  Please use getTransformedVectorValuesWithComponentBasisValues instead.");
  }
  switch(relatedOp)
  {
    case(Intrepid::OPERATOR_VALUE):
      switch(fs)
    {
      case Camellia::FUNCTION_SPACE_REAL_SCALAR:
        //          cout << "Reference values for FUNCTION_SPACE_REAL_SCALAR: " << *referenceValues;
      case Camellia::FUNCTION_SPACE_HGRAD:
      case Camellia::FUNCTION_SPACE_HGRAD_DISC:
        fst::HGRADtransformVALUE<double>(*transformedValues,*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HCURL:
      case Camellia::FUNCTION_SPACE_HCURL_DISC:
        fst::HCURLtransformVALUE<double>(*transformedValues,basisCache->getJacobianInv(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HDIV:
      case Camellia::FUNCTION_SPACE_HDIV_DISC:
      case Camellia::FUNCTION_SPACE_HDIV_FREE:
        fst::HDIVtransformVALUE<double>(*transformedValues,basisCache->getJacobian(),basisCache->getJacobianDet(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_HVOL:
      case Camellia::FUNCTION_SPACE_HVOL_SPACE_HGRAD_TIME:
      case Camellia::FUNCTION_SPACE_HGRAD_SPACE_HVOL_TIME:
        //        {
        //          static bool haveWarned = false;
        //          if (!haveWarned) {
        //            cout << "WARNING: for the moment, switching to the standard HVOLtransformVALUE method.\n";
        //            haveWarned = true;
        //          }
        //        }
        //          fst::HVOLtransformVALUE<double>(*transformedValues, cellJacobianDet, *referenceValues);
        // for the moment, use the fact that we know the HVOL basis is always an HGRAD basis:
        // (I think using the below amounts to solving for the HVOL variables scaled by Jacobian)
        fst::HGRADtransformVALUE<double>(*transformedValues,*referenceValues);
        break;
      default:
        TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
        break;
    }
      break;
    case(Intrepid::OPERATOR_GRAD):
    case(Intrepid::OPERATOR_D1):
      switch(fs)
    {
      case Camellia::FUNCTION_SPACE_HVOL:
      case Camellia::FUNCTION_SPACE_HGRAD:
      case Camellia::FUNCTION_SPACE_HGRAD_DISC:
        fst::HGRADtransformGRAD<double>(*transformedValues,basisCache->getJacobianInv(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_VECTOR_HVOL:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD_DISC:
        // referenceValues has dimensions (F,P,D1,D2).  D1 is our component dimension, and D2 is the one that came from the gradient.
        // HGRADtransformGRAD expects (F,P,D) for input, and (C,F,P,D) for output.
        // If we split referenceValues into (F,P,D1=0,D2) and (F,P,D1=1,D2), then we can transform each of those, and then interleave the results…
      {
        // block off so we can create new stuff inside the switch case
        Teuchos::Array<int> dimensions;
        referenceValues->dimensions(dimensions);
        int numFields = dimensions[0];
        int numPoints = dimensions[1];
        int D1 = dimensions[dimensions.size()-2];
        int D2 = dimensions[dimensions.size()-1];
        dimensions[dimensions.size()-2] = D2; // put D2 in the D1 spot
        dimensions.pop_back(); // get rid of original D2
        FieldContainer<double> refValuesSlice(dimensions);
        dimensions.insert(dimensions.begin(),numCells);
        FieldContainer<double> transformedValuesSlice(dimensions);
        
        //          int numEntriesPerSlice = refValuesSlice.size();
        //          int numEntriesPerTransformedSlice = transformedValuesSlice.size();
        
        for (int compIndex1=0; compIndex1<D1; compIndex1++)
        {
          // could speed the following along by doing the enumeration arithmetic in place...
          for (int fieldIndex=0; fieldIndex<numFields; fieldIndex++)
          {
            for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
//.........这里部分代码省略.........

FCPtr BasisEvaluation::getComponentOfInterest(constFCPtr values, Camellia::EOperator op, Camellia::EFunctionSpace fs, int componentOfInterest)
{
  FCPtr result;
  bool vectorizedBasis = functionSpaceIsVectorized(fs);
  if (   vectorizedBasis
         && ( (op == Camellia::OP_CURL) || (op == Camellia::OP_DIV)) )
  {
    TEUCHOS_TEST_FOR_EXCEPTION(values->rank() != 5, std::invalid_argument, "rank of values must be 5 for VECTOR_HGRAD_GRAD");
    int numCells  = values->dimension(0);
    int numFields = values->dimension(1);
    int numPoints = values->dimension(2);
    result = Teuchos::rcp(new FieldContainer<double>(numCells,numFields,numPoints));
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      for (int field=0; field<numFields; field++)
      {
        for (int point=0; point<numPoints; point++)
        {
          if (op==Camellia::OP_DIV)
          {
            (*result)(cellIndex,field,point) = (*values)(cellIndex,field,point,0,0) + (*values)(cellIndex,field,point,1,1);
          }
          else if (op==Camellia::OP_CURL)
          {
            // curl of 2D vector: d(Fx)/dy - d(Fy)/dx
            (*result)(cellIndex,field,point) = (*values)(cellIndex,field,point,1,0) - (*values)(cellIndex,field,point,0,1);
          }
        }
      }
    }
    return result;
  }
  else if ((componentOfInterest < 0) && (op != Camellia::OP_LAPLACIAN))     // then just return values
  {
    // the copy is a bit unfortunate, but can't be avoided unless we change a bunch of constFCPtrs to FCPtrs (or vice versa)
    // in the API...
//    cout << "values:\n" << *values;
    return Teuchos::rcp( new FieldContainer<double>(*values));
  }
  Teuchos::Array<int> dimensions;
  values->dimensions(dimensions);
  int opCardinality = dimensions[dimensions.size()-1];
  dimensions.pop_back(); // get rid of last, spatial dimension
  result = Teuchos::rcp(new FieldContainer<double>(dimensions));
  
  if (op == Camellia::OP_LAPLACIAN)
  {
    // back out spaceDim from opCardinality
    int spaceDim = -1;
    if (opCardinality == 1) spaceDim = 1;
    if (opCardinality == 3) spaceDim = 2;
    if (opCardinality == 6) spaceDim = 3;
//    if (opCardinality == 10) spaceDim = 4;
    TEUCHOS_TEST_FOR_EXCEPTION(spaceDim == -1, std::invalid_argument, "unhandled opCardinality and/or spaceDim");
    vector<int> dkEnumeration(spaceDim);
    if (spaceDim == 1)
    {
      dkEnumeration[0] = Intrepid::getDkEnumeration(2);
    }
    else if (spaceDim == 2)
    {
      dkEnumeration[0] = Intrepid::getDkEnumeration(2,0);
      dkEnumeration[1] = Intrepid::getDkEnumeration(0,2);
    }
    else if (spaceDim == 3)
    {
      dkEnumeration[0] = Intrepid::getDkEnumeration(2,0,0);
      dkEnumeration[1] = Intrepid::getDkEnumeration(0,2,0);
      dkEnumeration[2] = Intrepid::getDkEnumeration(0,0,2);
    }
    for (int comp : dkEnumeration)
    {
      int size = result->size();
      int enumeratedLocation;
      if (values->rank() == 3) // (F,P,D)
      {
        enumeratedLocation = values->getEnumeration(0,0,comp);
      }
      else if (values->rank() == 4) // (C,F,P,D)
      {
        enumeratedLocation = values->getEnumeration(0,0,0,comp);
      }
      else
      {
        // TODO: consider computing the enumerated location in a rank-independent way.
        TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unsupported values container rank.");
      }
      for (int i=0; i<size; i++)
      {
        (*result)[i] += (*values)[enumeratedLocation];
        enumeratedLocation += opCardinality;
      }
    }
    return result;
  }
  
  TEUCHOS_TEST_FOR_EXCEPTION(componentOfInterest >= opCardinality, std::invalid_argument, "componentOfInterest is out of bounds!");
//  int numPoints = dimensions[0];
//  int basisCardinality = dimensions[1];
  int size = result->size();
//.........这里部分代码省略.........