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

void UnitNormalFunction::values(Intrepid::FieldContainer<double> &values, BasisCachePtr basisCache)
{
    this->CHECK_VALUES_RANK(values);
    int numCells = values.dimension(0);
    int numPoints = values.dimension(1);
    int spaceDim = basisCache->getSpaceDim();
    if (_comp == -1)
    {
        // check the the "D" dimension of values is correct:
        if (_spaceTime)
        {
            TEUCHOS_TEST_FOR_EXCEPTION(values.dimension(2) != spaceDim+1, std::invalid_argument, "For space-time normals, values.dimension(2) should be spaceDim + 1.");
        }
        else
        {
            TEUCHOS_TEST_FOR_EXCEPTION(values.dimension(2) != spaceDim, std::invalid_argument, "For spatial normals, values.dimension(2) should be spaceDim.");
        }
    }
    const Intrepid::FieldContainer<double> *sideNormals = _spaceTime ? &(basisCache->getSideNormalsSpaceTime()) : &(basisCache->getSideNormals());

    int comp = _comp;
    if (comp == -2)
    {
        // want to select the temporal component, t()
        comp = spaceDim;
    }
    for (int cellIndex=0; cellIndex<numCells; cellIndex++)
    {
        for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
        {
            if (comp == -1)
            {
                for (int d=0; d<spaceDim; d++)
                {
                    double nd = (*sideNormals)(cellIndex,ptIndex,d);
                    values(cellIndex,ptIndex,d) = nd;
                }
                if (_spaceTime)
                {
                    double nd = (*sideNormals)(cellIndex,ptIndex,spaceDim);
                    values(cellIndex,ptIndex,spaceDim) = nd;
                }
            }
            else
            {
                double ni = (*sideNormals)(cellIndex,ptIndex,comp);
                values(cellIndex,ptIndex) = ni;
            }
        }
    }
}

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 values(FieldContainer<double> &values, BasisCachePtr basisCache)
  {
    // sets values(_cellIndex,P,D)
    TEUCHOS_TEST_FOR_EXCEPTION(_cellIndex == -1, std::invalid_argument, "must call setCellIndex before calling values!");

//    cout << "_basisCoefficients:\n" << _basisCoefficients;

    BasisCachePtr spaceTimeBasisCache;
    if (basisCache->cellTopologyIsSpaceTime())
    {
      // then we require that the basisCache provided be a space-time basis cache
      SpaceTimeBasisCache* spaceTimeCache = dynamic_cast<SpaceTimeBasisCache*>(basisCache.get());
      TEUCHOS_TEST_FOR_EXCEPTION(!spaceTimeCache, std::invalid_argument, "space-time requires a SpaceTimeBasisCache");
      spaceTimeBasisCache = basisCache;
      basisCache = spaceTimeCache->getSpatialBasisCache();
    }

    int numDofs = _basis->getCardinality();
    int spaceDim = basisCache->getSpaceDim();

    bool basisIsVolumeBasis = (spaceDim == _basis->domainTopology()->getDimension());
    bool useCubPointsSideRefCell = basisIsVolumeBasis && basisCache->isSideCache();

    int numPoints = values.dimension(1);

    // check if we're taking a temporal derivative
    int component;
    Intrepid::EOperator relatedOp = BasisEvaluation::relatedOperator(_op, _basis->functionSpace(), spaceDim, component);
    if ((relatedOp == Intrepid::OPERATOR_GRAD) && (component==spaceDim)) {
      // then we are taking the temporal part of the Jacobian of the reference to curvilinear-reference space
      // based on our assumptions that curvilinearity is just in the spatial direction (and is orthogonally extruded in the
      // temporal direction), this is always the identity.
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        for (int d=0; d<values.dimension(2); d++)
        {
          if (d < spaceDim)
            values(_cellIndex,ptIndex,d) = 0.0;
          else
            values(_cellIndex,ptIndex,d) = 1.0;
        }
      }
      return;
    }
    constFCPtr transformedValues = basisCache->getTransformedValues(_basis, _op, useCubPointsSideRefCell);

    // transformedValues has dimensions (C,F,P,[D,D])
    // therefore, the rank of the sum is transformedValues->rank() - 3
    int rank = transformedValues->rank() - 3;
    TEUCHOS_TEST_FOR_EXCEPTION(rank != values.rank()-2, std::invalid_argument, "values rank is incorrect.");


    int spaceTimeSideOrdinal = (spaceTimeBasisCache != Teuchos::null) ? spaceTimeBasisCache->getSideIndex() : -1;
    // I'm pretty sure much of this treatment of the time dimension could be simplified by taking advantage of SpaceTimeBasisCache::getTemporalBasisCache()...
    double t0 = -1, t1 = -1;
    if ((spaceTimeSideOrdinal != -1) && (!spaceTimeBasisCache->cellTopology()->sideIsSpatial(spaceTimeSideOrdinal)))
    {
      unsigned sideTime0 = spaceTimeBasisCache->cellTopology()->getTemporalSideOrdinal(0);
      unsigned sideTime1 = spaceTimeBasisCache->cellTopology()->getTemporalSideOrdinal(1);
      // get first node of each of the time-orthogonal sides, and use that to determine t0 and t1:
      unsigned spaceTimeNodeTime0 = spaceTimeBasisCache->cellTopology()->getNodeMap(spaceDim, sideTime0, 0);
      unsigned spaceTimeNodeTime1 = spaceTimeBasisCache->cellTopology()->getNodeMap(spaceDim, sideTime1, 0);
      t0 = spaceTimeBasisCache->getPhysicalCellNodes()(_cellIndex,spaceTimeNodeTime0,spaceDim);
      t1 = spaceTimeBasisCache->getPhysicalCellNodes()(_cellIndex,spaceTimeNodeTime1,spaceDim);
    }

    // initialize the values we're responsible for setting
    if (_op == OP_VALUE)
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
      {
        for (int d=0; d<values.dimension(2); d++)
        {
          if (d < spaceDim)
            values(_cellIndex,ptIndex,d) = 0.0;
          else if ((spaceTimeBasisCache != Teuchos::null) && (spaceTimeSideOrdinal == -1))
            values(_cellIndex,ptIndex,spaceDim) = spaceTimeBasisCache->getPhysicalCubaturePoints()(_cellIndex,ptIndex,spaceDim);
          else if ((spaceTimeBasisCache != Teuchos::null) && (spaceTimeSideOrdinal != -1))
          {
            if (spaceTimeBasisCache->cellTopology()->sideIsSpatial(spaceTimeSideOrdinal))
            {
              values(_cellIndex,ptIndex,spaceDim) = spaceTimeBasisCache->getPhysicalCubaturePoints()(_cellIndex,ptIndex,spaceDim-1);
            }
            else
            {
              double temporalPoint;
              unsigned temporalNode = spaceTimeBasisCache->cellTopology()->getTemporalComponentSideOrdinal(spaceTimeSideOrdinal);
              if (temporalNode==0)
                temporalPoint = t0;
              else
                temporalPoint = t1;
              values(_cellIndex,ptIndex,spaceDim) = temporalPoint;
            }
          }
        }
      }
    }
    else if ((_op == OP_DX) || (_op == OP_DY) || (_op == OP_DZ))
    {
      for (int ptIndex=0; ptIndex<numPoints; ptIndex++)
//.........这里部分代码省略.........