C++ BasisCachePtr::getPhysicalCubaturePoints方法代码示例

 void values(FieldContainer<double> &values, BasisCachePtr sliceBasisCache) {
   vector<GlobalIndexType> sliceCellIDs = sliceBasisCache->cellIDs();
   
   Teuchos::Array<int> dim;
   values.dimensions(dim);
   dim[0] = 1; // one cell
   Teuchos::Array<int> offset(dim.size());
   
   for (int cellOrdinal = 0; cellOrdinal < sliceCellIDs.size(); cellOrdinal++) {
     offset[0] = cellOrdinal;
     int enumeration = values.getEnumeration(offset);
     FieldContainer<double>valuesForCell(dim,&values[enumeration]);
     GlobalIndexType sliceCellID = sliceCellIDs[cellOrdinal];
     int numPoints = sliceBasisCache->getPhysicalCubaturePoints().dimension(1);
     int spaceDim = sliceBasisCache->getPhysicalCubaturePoints().dimension(2);
     FieldContainer<double> spaceTimePhysicalPoints(1,numPoints,spaceDim+1);
     for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++) {
       for (int d=0; d<spaceDim; d++) {
         spaceTimePhysicalPoints(0,ptOrdinal,d) = sliceBasisCache->getPhysicalCubaturePoints()(cellOrdinal,ptOrdinal,d);
       }
       spaceTimePhysicalPoints(0,ptOrdinal,spaceDim) = _t;
     }
     
     GlobalIndexType cellID = _cellIDMap[sliceCellID];
     BasisCachePtr spaceTimeBasisCache = BasisCache::basisCacheForCell(_spaceTimeMesh, cellID);
     
     FieldContainer<double> spaceTimeRefPoints(1,numPoints,spaceDim+1);
     CamelliaCellTools::mapToReferenceFrame(spaceTimeRefPoints, spaceTimePhysicalPoints, _spaceTimeMesh, cellID);
     spaceTimeRefPoints.resize(numPoints,spaceDim+1);
     spaceTimeBasisCache->setRefCellPoints(spaceTimeRefPoints);
     _spaceTimeFunction->values(valuesForCell, spaceTimeBasisCache);
   }
 }

void FunctionTests::checkFunctionsAgree(FunctionPtr f1, FunctionPtr f2, BasisCachePtr basisCache) {
  ASSERT_EQ(f1->rank(), f2->rank())
    << "f1->rank() " << f1->rank() << " != f2->rank() " << f2->rank() << endl;

  int rank = f1->rank();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numPoints = basisCache->getPhysicalCubaturePoints().dimension(1);
  int spaceDim = basisCache->getPhysicalCubaturePoints().dimension(2);
  Teuchos::Array<int> dim;
  dim.append(numCells);
  dim.append(numPoints);
  for (int i=0; i<rank; i++) {
    dim.append(spaceDim);
  }
  FieldContainer<double> f1Values(dim);
  FieldContainer<double> f2Values(dim);
  f1->values(f1Values,basisCache);
  f2->values(f2Values,basisCache);
  
  double tol = 1e-14;
  double maxDiff;
  EXPECT_TRUE(fcsAgree(f1Values,f2Values,tol,maxDiff))
    << "Test failed: f1 and f2 disagree; maxDiff " << maxDiff << ".\n"
    << "f1Values: \n" << f1Values
    << "f2Values: \n" << f2Values;
}

    void values(FieldContainer<double> &values, BasisCachePtr basisCache) {
      int numCells = values.dimension(0);
      int numPoints = values.dimension(1);

      const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
      for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
        double xCenter = 0;
        double yCenter = 0;
        int nPts = 0;
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
          double x = (*points)(cellIndex,ptIndex,0);
          double y = (*points)(cellIndex,ptIndex,1);
          xCenter += x;
          yCenter += y;
          nPts++;
        }
        xCenter /= nPts;
        yCenter /= nPts;
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
          double x = (*points)(cellIndex,ptIndex,0);
          double y = (*points)(cellIndex,ptIndex,1);
          if (abs(y) <= halfWidth && abs(yCenter) < halfWidth)
            values(cellIndex, ptIndex) = 1.0;
          else
            values(cellIndex, ptIndex) = 0.0;
        }
      }
    }

  // bool matchesPoint(double x, double y) {
  //   return true;
  // }
  bool matchesPoints(FieldContainer<bool> &pointsMatch, BasisCachePtr basisCache)
  {
    const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
    const FieldContainer<double> *normals = &(basisCache->getSideNormals());
    int numCells = (*points).dimension(0);
    int numPoints = (*points).dimension(1);

    FieldContainer<double> beta_pts(numCells,numPoints,2);
    _beta->values(beta_pts,basisCache);

    double tol=1e-14;
    bool somePointMatches = false;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        double n1 = (*normals)(cellIndex,ptIndex,0);
        double n2 = (*normals)(cellIndex,ptIndex,1);
        double beta_n = beta_pts(cellIndex,ptIndex,0)*n1 + beta_pts(cellIndex,ptIndex,1)*n2 ;
        // cout << "beta = (" << beta_pts(cellIndex,ptIndex,0)<<", "<<
        //   beta_pts(cellIndex,ptIndex,1)<<"), n = ("<<n1<<", "<<n2<<")"<<endl;
        pointsMatch(cellIndex,ptIndex) = false;
        if (beta_n < 0)
        {
          pointsMatch(cellIndex,ptIndex) = true;
          somePointMatches = true;
        }
      }
    }
    return somePointMatches;
  }

  bool matchesPoints(FieldContainer<bool> &pointsMatch, BasisCachePtr basisCache)
  {
    const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
    const FieldContainer<double> *normals = &(basisCache->getSideNormals());
    int numCells = (*points).dimension(0);
    int numPoints = (*points).dimension(1);

    double tol = 1e-3;
    bool somePointMatches = false;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        double x = (*points)(cellIndex,ptIndex,0);
        double y = (*points)(cellIndex,ptIndex,1);
        double n1 = (*normals)(cellIndex,ptIndex,0);
        double n2 = (*normals)(cellIndex,ptIndex,1);
        double beta_n = _beta[0]*n1 + _beta[1]*n2 ;
        pointsMatch(cellIndex,ptIndex) = false;
        if (abs((x-.5)*(x-.5)+y*y) < 0.75+tol && beta_n < 0)
        {
          pointsMatch(cellIndex,ptIndex) = true;
          somePointMatches = true;
        }
      }
    }
    return somePointMatches;
  }

void PreviousSolutionFunction::values(FieldContainer<double> &values, BasisCachePtr basisCache) {
  int rank = Teuchos::GlobalMPISession::getRank();
  if (_overrideMeshCheck) {
    _solnExpression->evaluate(values, _soln, basisCache);
    return;
  }
  if (!basisCache.get()) cout << "basisCache is nil!\n";
  if (!_soln.get()) cout << "_soln is nil!\n";
  // values are stored in (C,P,D) order
  if (basisCache->mesh().get() == _soln->mesh().get()) {
    _solnExpression->evaluate(values, _soln, basisCache);
  } else {
    static bool warningIssued = false;
    if (!warningIssued) {
      if (rank==0)
        cout << "NOTE: In PreviousSolutionFunction, basisCache's mesh doesn't match solution's.  If this is not what you intended, it would be a good idea to make sure that the mesh is passed in on BasisCache construction; the evaluation will be a lot slower without it...\n";
      warningIssued = true;
    }
    // get the physicalPoints, and make a basisCache for each...
    FieldContainer<double> physicalPoints = basisCache->getPhysicalCubaturePoints();
    FieldContainer<double> value(1,1); // assumes scalar-valued solution function.
    int numCells = physicalPoints.dimension(0);
    int numPoints = physicalPoints.dimension(1);
    int spaceDim = physicalPoints.dimension(2);
    physicalPoints.resize(numCells*numPoints,spaceDim);
    vector< ElementPtr > elements = _soln->mesh()->elementsForPoints(physicalPoints, false); // false: don't make elements null just because they're off-rank.
    FieldContainer<double> point(1,spaceDim);
    FieldContainer<double> refPoint(1,spaceDim);
    int combinedIndex = 0;
    vector<GlobalIndexType> cellID;
    cellID.push_back(-1);
    BasisCachePtr basisCacheOnePoint;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++, combinedIndex++) {
        if (elements[combinedIndex].get()==NULL) continue; // no element found for point; skip it…
        ElementTypePtr elemType = elements[combinedIndex]->elementType();
        for (int d=0; d<spaceDim; d++) {
          point(0,d) = physicalPoints(combinedIndex,d);
        }
        if (elements[combinedIndex]->cellID() != cellID[0]) {
          cellID[0] = elements[combinedIndex]->cellID();
          basisCacheOnePoint = Teuchos::rcp( new BasisCache(elemType, _soln->mesh()) );
          basisCacheOnePoint->setPhysicalCellNodes(_soln->mesh()->physicalCellNodesForCell(cellID[0]),cellID,false); // false: don't createSideCacheToo
        }
        // compute the refPoint:
        typedef CellTools<double>  CellTools;
        int whichCell = 0;
        CellTools::mapToReferenceFrame(refPoint,point,_soln->mesh()->physicalCellNodesForCell(cellID[0]),
                                       *(elemType->cellTopoPtr),whichCell);
        basisCacheOnePoint->setRefCellPoints(refPoint);
        //          cout << "refCellPoints:\n " << refPoint;
        //          cout << "physicalCubaturePoints:\n " << basisCacheOnePoint->getPhysicalCubaturePoints();
        _solnExpression->evaluate(value, _soln, basisCacheOnePoint);
        //          cout << "value at point (" << point(0,0) << ", " << point(0,1) << ") = " << value(0,0) << endl;
        values(cellIndex,ptIndex) = value(0,0);
      }
    }
  }
}

void MeshTransformationFunction::values(FieldContainer<double> &values, BasisCachePtr basisCache)
{
  CHECK_VALUES_RANK(values);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
  // identity map is the right thing most of the time
  // we'll do something different only where necessary
  int spaceDim = values.dimension(2);
  TEUCHOS_TEST_FOR_EXCEPTION(basisCache->cellTopology()->getDimension() != spaceDim, std::invalid_argument, "cellTopology dimension does not match the shape of the values container");
  if (_op == OP_VALUE)
  {
    values = basisCache->getPhysicalCubaturePoints(); // identity
  }
  else if (_op == OP_DX)
  {
    // identity map is 1 in all the x slots, 0 in all others
    int mod_value = 0; // the x slots are the mod spaceDim = 0 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
  else if (_op == OP_DY)
  {
    // identity map is 1 in all the y slots, 0 in all others
    int mod_value = 1; // the y slots are the mod spaceDim = 1 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
  else if (_op == OP_DZ)
  {
    // identity map is 1 in all the z slots, 0 in all others
    int mod_value = 2; // the z slots are the mod spaceDim = 2 slots;
    for (int i=0; i<values.size(); i++)
    {
      values[i] = (i%spaceDim == mod_value) ? 1.0 : 0.0;
    }
  }
//  if (_op == OP_DX) {
//    cout << "values before cellTransformation:\n" << values;
//  }
  for (int cellIndex=0; cellIndex < cellIDs.size(); cellIndex++)
  {
    GlobalIndexType cellID = cellIDs[cellIndex];
    if (_cellTransforms.find(cellID) == _cellTransforms.end()) continue;
    TFunctionPtr<double> cellTransformation = _cellTransforms[cellID];
    ((CellTransformationFunction*)cellTransformation.get())->setCellIndex(cellIndex);
    cellTransformation->values(values, basisCache);
  }
//  if (_op == OP_DX) {
//    cout << "values after cellTransformation:\n" << values;
//  }
}

    void values(FieldContainer<double> &values, BasisCachePtr basisCache) {
      int numCells = values.dimension(0);
      int numPoints = values.dimension(1);

      const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
      for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
          double x = (*points)(cellIndex,ptIndex,0);
          double y = (*points)(cellIndex,ptIndex,1);
          values(cellIndex, ptIndex) = 0;
        }
      }
    }

 void values(FieldContainer<double> &values, BasisCachePtr basisCache)
 {
   vector<GlobalIndexType> cellIDs = basisCache->cellIDs();
   int numPoints = values.dimension(1);
   FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
   for (int i = 0; i<cellIDs.size(); i++)
   {
     for (int j = 0; j<numPoints; j++)
     {
       values(i,j) = cellIDs[i];
     }
   }
 }

bool FunctionTests::functionsAgree(FunctionPtr f1, FunctionPtr f2, BasisCachePtr basisCache)
{
  if (f2->rank() != f1->rank() )
  {
    cout << "f1->rank() " << f1->rank() << " != f2->rank() " << f2->rank() << endl;
    return false;
  }
  int rank = f1->rank();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numPoints = basisCache->getPhysicalCubaturePoints().dimension(1);
  int spaceDim = basisCache->getSpaceDim();
  Teuchos::Array<int> dim;
  dim.append(numCells);
  dim.append(numPoints);
  for (int i=0; i<rank; i++)
  {
    dim.append(spaceDim);
  }
  FieldContainer<double> f1Values(dim);
  FieldContainer<double> f2Values(dim);
  f1->values(f1Values,basisCache);
  f2->values(f2Values,basisCache);

  double tol = 1e-14;
  double maxDiff;
  bool functionsAgree = TestSuite::fcsAgree(f1Values,f2Values,tol,maxDiff);
  if ( ! functionsAgree )
  {
    functionsAgree = false;
    cout << "Test failed: f1 and f2 disagree; maxDiff " << maxDiff << ".\n";
    cout << "f1Values: \n" << f1Values;
    cout << "f2Values: \n" << f2Values;
  }
  else
  {
//    cout << "f1 and f2 agree!" << endl;
  }
  return functionsAgree;
}

void ExactSolution::solutionValues(FieldContainer<double> &values, int trialID, BasisCachePtr basisCache) {
  if (_exactFunctions.find(trialID) != _exactFunctions.end() ) {
    _exactFunctions[trialID]->values(values,basisCache);
    return;
  }
  
  // TODO: change ExactSolution::solutionValues (below) to take a *const* points FieldContainer, to avoid this copy:
  FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
  if (basisCache->getSideIndex() >= 0) {
    FieldContainer<double> unitNormals = basisCache->getSideNormals();
    this->solutionValues(values,trialID,points,unitNormals);
  } else {
    this->solutionValues(values,trialID,points);
  }
}

  void values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    FieldContainer<double> points  = basisCache->getPhysicalCubaturePoints();
    FieldContainer<double> normals = basisCache->getSideNormals();
    int numCells = points.dimension(0);
    int numPoints = points.dimension(1);

    FieldContainer<double> Tv(numCells,numPoints);
    FieldContainer<double> u1v(numCells,numPoints);;
    _u1->values(u1v,basisCache);
    _T->values(Tv,basisCache);

    bool isSubsonic = false;
    double min_y = YTOP;
    double max_y = 0.0;
    values.initialize(0.0);
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        double x = points(cellIndex,ptIndex,0);
        double y = points(cellIndex,ptIndex,1);

        double T = Tv(cellIndex,ptIndex);
        double un = u1v(cellIndex,ptIndex); // WARNING: ASSUMES NORMAL AT OUTFLOW = (1,0)
        double c = sqrt(_gamma * (_gamma-1.0) * _cv * T);

        bool outflowMatch = ((abs(x-2.0) < _tol) && (y > 0.0) && (y < YTOP));
        bool subsonicMatch = (un < c) && (un > 0.0);
        if (subsonicMatch && outflowMatch)
        {
          values(cellIndex,ptIndex) = 1.0;
          isSubsonic = true;
          min_y = min(y,min_y);
          max_y = max(y,max_y);
          //	  cout << "y = " << y << endl;
        }
      }
    }
    if (isSubsonic)
    {
      //      cout << "subsonic in interval y =(" << min_y << "," << max_y << ")" << endl;
    }
  }

 void values(FieldContainer<double> &values, BasisCachePtr basisCache) {
   int numCells = values.dimension(0);
   int numPoints = values.dimension(1);
   int spaceDim = values.dimension(2);
   
   const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
   for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
     for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
       for (int d = 0; d < spaceDim; d++) {
         double x = (*points)(cellIndex,ptIndex,0);
         double y = (*points)(cellIndex,ptIndex,1);
         if ((x+.5)*(x+.5)+y*y < .25*.25)
           values(cellIndex,ptIndex) = 1.0;
         else
           values(cellIndex,ptIndex) = 0.0;
       }
     }
   }
 }

bool FunctionTests::testVectorFunctionValuesOrdering()
{
  bool success = true;

  FunctionPtr x = Function::xn(1);
  FunctionPtr x_vector = Function::vectorize(x, Function::zero());

  BasisCachePtr basisCache = BasisCache::parametricQuadCache(10);

  FieldContainer<double> points = basisCache->getPhysicalCubaturePoints();
  int numCells = points.dimension(0);
  int numPoints = points.dimension(1);
  int spaceDim = points.dimension(2);
  FieldContainer<double> values(numCells,numPoints,spaceDim);

  x_vector->values(values, basisCache);

//  cout << "(x,0) function values:\n" << values;

  double tol = 1e-14;
  for (int cellIndex=0; cellIndex<numCells; cellIndex++)
  {
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
    {
      double xValueExpected = points(cellIndex,ptIndex,0);
      double yValueExpected = 0;
      double xValue = values(cellIndex,ptIndex,0);
      double yValue = values(cellIndex,ptIndex,1);
      double xErr = abs(xValue-xValueExpected);
      double yErr = abs(yValue-yValueExpected);
      if ( (xErr > tol) || (yErr > tol) )
      {
        success = false;
        cout << "testVectorFunctionValuesOrdering(): vectorized function values incorrect (presumably out of order).\n";
        cout << "x: " << xValueExpected << " ≠ " << xValue << endl;
        cout << "y: " << yValueExpected << " ≠ " << yValue << endl;
      }
    }
  }

  return success;
}

  void values(FieldContainer<double> &values, BasisCachePtr basisCache) {
    int numCells = values.dimension(0);
    int numPoints = values.dimension(1);

    FieldContainer<double> beta_pts(numCells,numPoints,2);
    _beta->values(beta_pts,basisCache);

    const FieldContainer<double> *points = &(basisCache->getPhysicalCubaturePoints());
    double tol=1e-14;
    for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
        double x = (*points)(cellIndex,ptIndex,0);
        double y = (*points)(cellIndex,ptIndex,1);
	double b1 =  beta_pts(cellIndex,ptIndex,0);
	double b2 =  beta_pts(cellIndex,ptIndex,1);
	double beta_norm =b1*b1 + b2*b2;
	values(cellIndex,ptIndex) = sqrt(beta_norm);
      }
    }
  }