C++ BasisPtr::getCardinality方法代码示例

bool checkVertexOrdinalsQuad(BasisPtr basis, vector<int> &vertexOrdinals)
{
  // check that the given indices are exactly the vertex basis functions:
  // a) these are (1,0) or (0,1) at the corresponding vertices
  // b) others are (0,0) at the vertices

  int numVertices = 4;

  FieldContainer<double> refCellPoints(numVertices,2); // vertices, in order
  refCellPoints(0,0) = -1;
  refCellPoints(0,1) = -1;
  refCellPoints(1,0) =  1;
  refCellPoints(1,1) = -1;
  refCellPoints(2,0) =  1;
  refCellPoints(2,1) =  1;
  refCellPoints(3,0) = -1;
  refCellPoints(3,1) =  1;

  // assume scalar basis for now -- we'll throw an exception if not...
  FieldContainer<double> values(basis->getCardinality(), numVertices); // F, P
  basis->getValues(values, refCellPoints, OPERATOR_VALUE);

  double tol = 1e-14;
  for (int vertexIndex=0; vertexIndex<numVertices; vertexIndex++)
  {
    int vertexOrdinal = vertexOrdinals[vertexIndex];
    for (int fieldIndex=0; fieldIndex<basis->getCardinality(); fieldIndex++)
    {
      double value = values(fieldIndex,vertexIndex);
      if (fieldIndex==vertexOrdinal)
      {
        // expect non-zero
        if (value < tol)
        {
          return false;
        }
      }
      else
      {
        // expect zero
        if (value > tol)
        {
          return false;
        }
      }
    }
  }
  return true;
}

void DofOrdering::addEntry(int varID, BasisPtr basis, int basisRank, int sideIndex) {
  // test to see if we already have one matching this.  (If so, that's an error.)
  pair<int,int> key = make_pair(varID, sideIndex);
  if ( indices.find(key) != indices.end() ) {
    TEUCHOS_TEST_FOR_EXCEPTION( true,
                       std::invalid_argument,
                       "Already have an entry in DofOrdering for this varID, sideIndex pair.");
  } else {
    indices[key] = vector<int>(basis->getCardinality());
  }
  
  vector<int>* dofIndices = &(indices[key]);
  
  for (vector<int>::iterator dofEntryIt = dofIndices->begin(); dofEntryIt != dofIndices->end(); dofEntryIt++) {
    *dofEntryIt = _nextIndex;
    _nextIndex++;
  }
  
  varIDs.insert(varID);
  
  numSidesForVarID[varID]++;
  //cout << "numSidesForVarID[" << varID << "]" << numSidesForVarID[varID] << endl;
  pair<int, int> basisKey = make_pair(varID,sideIndex);
//  _nextIndex += basis->getCardinality();
  bases[basisKey] = basis;
  basisRanks[varID] = basisRank;
}

void DofOrdering::copyLikeCoefficients( FieldContainer<double> &newValues, Teuchos::RCP<DofOrdering> oldDofOrdering,
                          const FieldContainer<double> &oldValues ) {
  // copy the coefficients for the bases that agree between the two DofOrderings
  // requires that "like" bases are actually pointers to the same memory location
  TEUCHOS_TEST_FOR_EXCEPTION( newValues.rank() != 1, std::invalid_argument, "newValues.rank() != 1");
  TEUCHOS_TEST_FOR_EXCEPTION( newValues.size() != totalDofs(), std::invalid_argument, "newValues.size() != totalDofs()");
  TEUCHOS_TEST_FOR_EXCEPTION( oldValues.rank() != 1, std::invalid_argument, "oldValues.rank() != 1");
  TEUCHOS_TEST_FOR_EXCEPTION( oldValues.size() != oldDofOrdering->totalDofs(), std::invalid_argument, "oldValues.size() != oldDofOrdering->totalDofs()");
  
  newValues.initialize(0.0);
  
  for (set<int>::iterator varIDIt = varIDs.begin(); varIDIt != varIDs.end(); varIDIt++) {
    int varID = *varIDIt;
    int numSides = getNumSidesForVarID(varID);
    if ( numSides == oldDofOrdering->getNumSidesForVarID(varID) ) {
      for (int sideIndex=0; sideIndex < numSides; sideIndex++) {
        BasisPtr basis = getBasis(varID,sideIndex);
        if (basis.get() == oldDofOrdering->getBasis(varID,sideIndex).get() ) {
          // bases alike: copy coefficients
          int cardinality = basis->getCardinality();
          for (int dofOrdinal=0; dofOrdinal < cardinality; dofOrdinal++) {
            int dofIndex = getDofIndex(varID,dofOrdinal,sideIndex);
            newValues(dofIndex) = oldValues( oldDofOrdering->getDofIndex(varID,dofOrdinal,sideIndex) );
          }
        }
      }
    }
  }
}

int DofOrdering::getTotalBasisCardinality() { // total number of *distinct* basis functions
  int totalBasisCardinality = 0;
  set< ::Camellia::Basis<>* > basesCounted;
  for (map< pair<int, int>, BasisPtr >::iterator basisIt = bases.begin(); basisIt != bases.end(); basisIt++) {
    BasisPtr basis = basisIt->second;
    if (basesCounted.find(basis.get())==basesCounted.end()) {
      basesCounted.insert(basis.get());
      totalBasisCardinality += basis->getCardinality();
    }
  }
  return totalBasisCardinality;
}

void DofOrdering::rebuildIndex() {
  if (dofIdentifications.size() == 0) {
    // nothing to do
    return;
  }
  set<int>::iterator varIDIterator;
  int numIdentificationsProcessed = 0;
  for (varIDIterator = varIDs.begin(); varIDIterator != varIDs.end(); varIDIterator++) {
    int varID = *varIDIterator;
    //cout << "rebuildIndex: varID=" << varID << endl;
    for (int sideIndex=0; sideIndex<numSidesForVarID[varID]; sideIndex++) {
      BasisPtr basis = getBasis(varID,sideIndex);
      int cellTopoSideIndex = (numSidesForVarID[varID]==1) ? -1 : sideIndex;
      if ( _cellTopologyForSide.find(cellTopoSideIndex) == _cellTopologyForSide.end() ) {
        CellTopoPtr cellTopoPtr = basis->domainTopology();
        _cellTopologyForSide[cellTopoSideIndex] = cellTopoPtr;
      }
      for (int dofOrdinal=0; dofOrdinal < basis->getCardinality(); dofOrdinal++) {
        pair<int, pair<int,int> > key = make_pair(varID, make_pair(sideIndex,dofOrdinal));
        pair<int, int> indexKey = make_pair(key.first,key.second.first); // key into indices container
        if ( dofIdentifications.find(key) != dofIdentifications.end() ) {
          int earlierSideIndex  = dofIdentifications[key].first;
          int earlierDofOrdinal = dofIdentifications[key].second;
          pair<int,int> earlierIndexKey = make_pair(varID,earlierSideIndex);
          if (indices[indexKey][dofOrdinal] != indices[earlierIndexKey][earlierDofOrdinal]) {
            indices[indexKey][dofOrdinal] = indices[earlierIndexKey][earlierDofOrdinal];
//            cout << "processed identification for varID " << varID << ": (" << sideIndex << "," << dofOrdinal << ")";
//            cout << " --> " << "(" << earlierSideIndex << "," << earlierDofOrdinal << ")" << endl;
            numIdentificationsProcessed++;
          }
        } else {
          // modify the index according to the number of dofs we've consolidated
//          cout << "Reducing indices for key (varID=" << indexKey.first << ", sideIndex " << indexKey.second << ") for dofOrdinal " << dofOrdinal << " from ";
//          cout << indices[indexKey][dofOrdinal] << " to ";
          indices[indexKey][dofOrdinal] -= numIdentificationsProcessed;
//          cout << indices[indexKey][dofOrdinal] << "\n";
        }
      }
    }
  }
  _nextIndex -= numIdentificationsProcessed;
  //cout << "index rebuilt; _nextIndex = " << _nextIndex << "; numIdentificationsProcessed: " << numIdentificationsProcessed << endl;
  _indexNeedsToBeRebuilt = false;
}

int DofOrdering::getBasisCardinality(int varID, int sideIndex) {
  BasisPtr basis = getBasis(varID,sideIndex);
  return basis->getCardinality();
//  return getBasis(varID,sideIndex)->getCardinality();
}

bool VectorizedBasisTestSuite::testVectorizedBasis()
{
  bool success = true;

  string myName = "testVectorizedBasis";

  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );

  int polyOrder = 3, numPoints = 5, spaceDim = 2;

  BasisPtr hgradBasis = BasisFactory::basisFactory()->getBasis(polyOrder, quad_4.getKey(), Camellia::FUNCTION_SPACE_HGRAD);

  // first test: make a single-component vector basis.  This should agree in every entry with the basis itself, but its field container will have one higher rank...
  VectorizedBasis<> oneComp(hgradBasis, 1);

  FieldContainer<double> linePoints(numPoints, spaceDim);
  for (int i=0; i<numPoints; i++)
  {
    for (int j=0; j<spaceDim; j++)
    {
      linePoints(i,j) = ((double)(i + j)) / (numPoints + spaceDim);
    }
  }

  FieldContainer<double> compValues(hgradBasis->getCardinality(),numPoints);
  hgradBasis->getValues(compValues, linePoints, Intrepid::OPERATOR_VALUE);

  FieldContainer<double> values(hgradBasis->getCardinality(),linePoints.dimension(0),1); // one component
  oneComp.getValues(values, linePoints, Intrepid::OPERATOR_VALUE);

  for (int i=0; i<compValues.size(); i++)
  {
    double diff = abs(values[i]-compValues[i]);
    if (diff != 0.0)
    {
      success = false;
      cout << myName << ": one-component vector basis doesn't produce same values as component basis." << endl;
      cout << "difference: " << diff << " in enumerated value " << i << endl;
      cout << "values:\n" << values;
      cout << "compValues:\n" << compValues;
      return success;
    }
  }

  vector< BasisPtr > twoComps;
  twoComps.push_back( Teuchos::rcp( new VectorizedBasis<>(hgradBasis, 2) ) );
  twoComps.push_back( BasisFactory::basisFactory()->getBasis( polyOrder,
                      quad_4.getKey(), Camellia::FUNCTION_SPACE_VECTOR_HGRAD) );


  vector< BasisPtr >::iterator twoCompIt;
  for (twoCompIt = twoComps.begin(); twoCompIt != twoComps.end(); twoCompIt++)
  {
    BasisPtr twoComp = *twoCompIt;

    int componentCardinality = hgradBasis->getCardinality();

    if (twoComp->getCardinality() != 2 * hgradBasis->getCardinality() )
    {
      success = false;
      cout << myName << ": two-component vector basis cardinality != one-component cardinality * 2." << endl;
      cout << "twoComp->getCardinality(): " << twoComp->getCardinality() << endl;
      cout << "oneComp->getCardinality(): " << oneComp.getCardinality() << endl;
    }

    values.resize(twoComp->getCardinality(),linePoints.dimension(0),2); // two components
    twoComp->getValues(values, linePoints, Intrepid::OPERATOR_VALUE);
    for (int basisIndex=0; basisIndex<twoComp->getCardinality(); basisIndex++)
    {
      for (int k=0; k<numPoints; k++)
      {
        double xValueExpected = (basisIndex < componentCardinality) ? compValues(basisIndex,k) : 0;
        double xValueActual = values(basisIndex,k,0);
        double yValueExpected = (basisIndex >= componentCardinality) ? compValues(basisIndex - componentCardinality,k) : 0;
        double yValueActual = values(basisIndex,k,1);
        if ( ( abs(xValueActual - xValueExpected) != 0) || ( abs(yValueActual - yValueExpected) != 0) )
        {
          success = false;
          cout << myName << ": expected differs from actual\n";
          cout << "component\n" << compValues;
          cout << "vector values:\n" << values;
          return success;
        }
      }
    }

    // test the mapping from oneComp dofOrdinal to twoComp:
    VectorizedBasis<>* twoCompAsVectorBasis = (VectorizedBasis<>  *) twoComp.get();

    for (int compDofOrdinal=0; compDofOrdinal<oneComp.getCardinality(); compDofOrdinal++)
    {
      int dofOrdinal_0 = twoCompAsVectorBasis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 0);
      int dofOrdinal_1 = twoCompAsVectorBasis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 1);
      // we expect the lists to be stacked (this is implicit in the test above)
      // dofOrdinal_0 we expect to be == compDofOrdinal
      // dofOrdinal_1 we expect to be == compDofOrdinal + oneComp.getCardinality()
      if (dofOrdinal_0 != compDofOrdinal)
      {
        success = false;
        cout << "getDofOrdinalFromComponentDofOrdinal() not returning expected value in first component.\n";
//.........这里部分代码省略.........

void Projector::projectFunctionOntoBasis(FieldContainer<double> &basisCoefficients, FunctionPtr fxn, 
                                         BasisPtr basis, BasisCachePtr basisCache, IPPtr ip, VarPtr v,
                                         set<int> fieldIndicesToSkip) {
  CellTopoPtr cellTopo = basis->domainTopology();
  DofOrderingPtr dofOrderPtr = Teuchos::rcp(new DofOrdering());
  
  if (! fxn.get()) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "fxn cannot be null!");
  }
  
  int cardinality = basis->getCardinality();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numDofs = cardinality - fieldIndicesToSkip.size();
  if (numDofs==0) {
    // we're skipping all the fields, so just initialize basisCoefficients to 0 and return
    basisCoefficients.resize(numCells,cardinality);
    basisCoefficients.initialize(0);
    return;
  }
  
  FieldContainer<double> gramMatrix(numCells,cardinality,cardinality);
  FieldContainer<double> ipVector(numCells,cardinality);

  // fake a DofOrdering
  DofOrderingPtr dofOrdering = Teuchos::rcp( new DofOrdering );
  if (! basisCache->isSideCache()) {
    dofOrdering->addEntry(v->ID(), basis, v->rank());
  } else {
    dofOrdering->addEntry(v->ID(), basis, v->rank(), basisCache->getSideIndex());
  }
  
  ip->computeInnerProductMatrix(gramMatrix, dofOrdering, basisCache);
  ip->computeInnerProductVector(ipVector, v, fxn, dofOrdering, basisCache);
  
//  cout << "physical points for projection:\n" << basisCache->getPhysicalCubaturePoints();
//  cout << "gramMatrix:\n" << gramMatrix;
//  cout << "ipVector:\n" << ipVector;
  
  map<int,int> oldToNewIndices;
  if (fieldIndicesToSkip.size() > 0) {
    // the code to do with fieldIndicesToSkip might not be terribly efficient...
    // (but it's not likely to be called too frequently)
    int i_indices_skipped = 0;
    for (int i=0; i<cardinality; i++) {
      int new_index;
      if (fieldIndicesToSkip.find(i) != fieldIndicesToSkip.end()) {
        i_indices_skipped++;
        new_index = -1;
      } else {
        new_index = i - i_indices_skipped;
      }
      oldToNewIndices[i] = new_index;
    }
    
    FieldContainer<double> gramMatrixFiltered(numCells,numDofs,numDofs);
    FieldContainer<double> ipVectorFiltered(numCells,numDofs);
    // now filter out the values that we're to skip
    
    for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
      for (int i=0; i<cardinality; i++) {
        int i_filtered = oldToNewIndices[i];
        if (i_filtered == -1) {
          continue;
        }
        ipVectorFiltered(cellIndex,i_filtered) = ipVector(cellIndex,i);
        
        for (int j=0; j<cardinality; j++) {
          int j_filtered = oldToNewIndices[j];
          if (j_filtered == -1) {
            continue;
          }
          gramMatrixFiltered(cellIndex,i_filtered,j_filtered) = gramMatrix(cellIndex,i,j);
        }
      }
    }
//    cout << "gramMatrixFiltered:\n" << gramMatrixFiltered;
//    cout << "ipVectorFiltered:\n" << ipVectorFiltered;
    gramMatrix = gramMatrixFiltered;
    ipVector = ipVectorFiltered;
  }
  
  for (int cellIndex=0; cellIndex<numCells; cellIndex++){
    
    // TODO: rewrite to take advantage of SerialDenseWrapper...
    Epetra_SerialDenseSolver solver;
    
    Epetra_SerialDenseMatrix A(Copy,
                               &gramMatrix(cellIndex,0,0),
                               gramMatrix.dimension(2), 
                               gramMatrix.dimension(2),  
                               gramMatrix.dimension(1)); // stride -- fc stores in row-major order (a.o.t. SDM)
    
    Epetra_SerialDenseVector b(Copy,
                               &ipVector(cellIndex,0),
                               ipVector.dimension(1));
    
    Epetra_SerialDenseVector x(gramMatrix.dimension(1));
    
    solver.SetMatrix(A);
    int info = solver.SetVectors(x,b);
//.........这里部分代码省略.........

bool testBasisClassifications(BasisPtr basis) {
  bool success = true;
  
  CellTopoPtr cellTopo = basis->domainTopology();
  
  int numVertices = cellTopo->getVertexCount();
  int numEdges = cellTopo->getEdgeCount();
  
  int degree = basis->getDegree();
  
  // TODO: finish this
  vector<int> vertexOrdinals;
  for (int vertexIndex=0; vertexIndex < numVertices; vertexIndex++) {
    set<int> dofOrdinals = basis->dofOrdinalsForVertex(vertexIndex);
    if (dofOrdinals.size() == 0) TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "No dofOrdinal for vertex...");
    if (dofOrdinals.size() > 1) TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "More than one dofOrdinal per vertex...");
    vertexOrdinals.push_back(*(dofOrdinals.begin()));
  }
//  
//  if (! checkVertexOrdinalsQuad(basis, vertexOrdinals) ) {
//    success = false;
//    cout << "vertex dof ordinals don't match expected\n";
//    cout << "ordinals: ";
//    for (int vertexIndex=0; vertexIndex < numVertices; vertexIndex++) {
//      cout << vertexOrdinals[vertexIndex] << " ";
//    }
//    cout << endl;
//  }
  
  // get the points in reference space for each vertex
  FieldContainer<double> points;
  if (numVertices == 2) { // line
    points.resize(2,1);
    points(0,0) = -1;
    points(1,0) = 1;
  } else if (numVertices == 3) { // triangle
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "triangles not yet supported");
  } else if (numVertices == 4) { // quad
    points.resize(4,2);
    points(0,0) = -1;
    points(0,1) = -1;
    points(1,0) =  1;
    points(1,1) = -1;
    points(2,0) =  1;
    points(2,1) =  1;
    points(3,0) = -1;
    points(3,1) =  1;
  } else {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unsupported topology");
  }
  
  FieldContainer<double> vertexValues;
  if (basis->rangeRank() > 0) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "rank > 0 bases not yet supported");
  } else {
    vertexValues.resize(basis->getCardinality(),numVertices);
  }
  
  basis->getValues(vertexValues, points, Intrepid::OPERATOR_VALUE);
  
  // test that the points are correctly classified
  for (int fieldIndex=0; fieldIndex<basis->getCardinality(); fieldIndex++) {
    for (int ptIndex=0; ptIndex<numVertices; ptIndex++) {
      int dofOrdinalForPoint = vertexOrdinals[ptIndex];
      bool expectZero = (dofOrdinalForPoint != fieldIndex);
      if (expectZero) {
        if (vertexValues(fieldIndex,ptIndex) != 0) {
          success = false;
          cout << "Expected 0 for fieldIndex " << fieldIndex << " and ptIndex " << ptIndex;
          cout << ", but got " << vertexValues(fieldIndex,ptIndex) << endl;
        }
      } else {
        if (vertexValues(fieldIndex,ptIndex) == 0) {
          cout << "Expected nonzero for fieldIndex " << fieldIndex << " and ptIndex " << ptIndex << endl;
          success = false;
        }
      }
    }
  }
  
  if (!success) {
    cout << "Failed testBasisClassifications; vertexValues:\n" << vertexValues;
  }
  
  return success;
}

void Boundary::bcsToImpose( map<  GlobalIndexType, Scalar > &globalDofIndicesAndValues, TBC<Scalar> &bc,
                            GlobalIndexType cellID, DofInterpreter* dofInterpreter)
{
  // this is where we actually compute the BCs; the other bcsToImpose variants call this one.
  CellPtr cell = _mesh->getTopology()->getCell(cellID);

  // define a couple of important inner products:
  TIPPtr<Scalar> ipL2 = Teuchos::rcp( new TIP<Scalar> );
  TIPPtr<Scalar> ipH1 = Teuchos::rcp( new TIP<Scalar> );
  VarFactoryPtr varFactory = VarFactory::varFactory();
  VarPtr trace = varFactory->traceVar("trace");
  VarPtr flux = varFactory->traceVar("flux");
  ipL2->addTerm(flux);
  ipH1->addTerm(trace);
  ipH1->addTerm(trace->grad());
  ElementTypePtr elemType = _mesh->getElementType(cellID);
  DofOrderingPtr trialOrderingPtr = elemType->trialOrderPtr;
  vector< int > trialIDs = _mesh->bilinearForm()->trialIDs();

  vector<unsigned> boundarySides = cell->boundarySides();
  if (boundarySides.size() > 0)
  {
    BasisCachePtr basisCache = BasisCache::basisCacheForCell(_mesh, cellID);
    for (vector< int >::iterator trialIt = trialIDs.begin(); trialIt != trialIDs.end(); trialIt++)
    {
      int trialID = *(trialIt);
      if ( bc.bcsImposed(trialID) )
      {
//        // DEBUGGING: keep track of which sides we impose BCs on:
//        set<unsigned> bcImposedSides;
//
        // Determine global dof indices and values, in one pass per side
        for (int i=0; i<boundarySides.size(); i++)
        {
          unsigned sideOrdinal = boundarySides[i];
          // TODO: here, we need to treat the volume basis case.
          /*
           To do this:
           1. (Determine which dofs in the basis have support on the side.)
           2. (Probably should resize dirichletValues to match number of dofs with support on the side.)
           3. (Within coefficientsForBC, and the projection method it calls, when it's a side cache, check whether the basis being projected has a higher dimension.  If so, do the same determination regarding the support of basis on the side as #1.)
           4. DofInterpreter::interpretLocalBasisCoefficients() needs to handle the case that trialID has volume support, and in this case interpret the provided data appropriately.
           */
          
          BasisPtr basis;
          int numDofsSide;
          if (trialOrderingPtr->getSidesForVarID(trialID).size() == 1)
          {
            // volume basis
            basis = trialOrderingPtr->getBasis(trialID);
            // get the dof ordinals for the side (interpreted as a "continuous" basis)
            numDofsSide = basis->dofOrdinalsForSide(sideOrdinal).size();
          }
          else if (! trialOrderingPtr->hasBasisEntry(trialID, sideOrdinal))
          {
            continue;
          }
          else
          {
            basis = trialOrderingPtr->getBasis(trialID,sideOrdinal);
            numDofsSide = basis->getCardinality();
          }
          
          GlobalIndexType numCells = 1;
          if (numCells > 0)
          {
            FieldContainer<double> dirichletValues(numCells,numDofsSide);
            // project bc function onto side basis:
            BCPtr bcPtr = Teuchos::rcp(&bc, false);
            Teuchos::RCP<BCFunction<double>> bcFunction = BCFunction<double>::bcFunction(bcPtr, trialID);
            bcPtr->coefficientsForBC(dirichletValues, bcFunction, basis, basisCache->getSideBasisCache(sideOrdinal));
            dirichletValues.resize(numDofsSide);
            if (bcFunction->imposeOnCell(0))
            {
              FieldContainer<double> globalData;
              FieldContainer<GlobalIndexType> globalDofIndices;

              dofInterpreter->interpretLocalBasisCoefficients(cellID, trialID, sideOrdinal, dirichletValues, globalData, globalDofIndices);
              for (int globalDofOrdinal=0; globalDofOrdinal<globalDofIndices.size(); globalDofOrdinal++)
              {
                GlobalIndexType globalDofIndex = globalDofIndices(globalDofOrdinal);
                Scalar value = globalData(globalDofOrdinal);
                
                // sanity check: if this has been previously set, do the two values roughly agree?
                if (globalDofIndicesAndValues.find(globalDofIndex) != globalDofIndicesAndValues.end())
                {
                  double tol = 1e-10;
                  Scalar prevValue = globalDofIndicesAndValues[globalDofIndex];
                  double absDiff = abs(prevValue - value);
                  if (absDiff > tol)
                  {
                    double relativeDiff = absDiff / max(abs(prevValue),abs(value));
                    int rank = _mesh->Comm()->MyPID();
                    if (relativeDiff > tol)
                    {
                      cout << "WARNING: in Boundary::bcsToImpose(), inconsistent values for BC: " << prevValue << " and ";
                      cout << value << " prescribed for global dof index " << globalDofIndex;
                      cout << " on rank " << rank << endl;
                    }
                  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                      cout << neighborCellID << ", " << neighborLocalDofIndex << ") -- ";
                      cout << globalDofIndex << " != " << neighborsGlobalDofIndex << "\n";
                      success = false;
                    }
                  }
                }
              }
              else if (neighbor->getNeighborCellID(ancestralSideIndexInNeighbor) != cellID)
              {
                // elem is small, neighbor big
                // first, find my leaf index in neighbor:
                int ancestorCellID = neighbor->getNeighborCellID(ancestralSideIndexInNeighbor);
                Teuchos::RCP<Element> ancestor = mesh->getElement(ancestorCellID);
                int ancestorSideIndex = neighbor->getSideIndexInNeighbor(ancestralSideIndexInNeighbor);
                vector< pair<int,int> > descendantsForSide = ancestor->getDescendantsForSide(ancestorSideIndex);
                int descendantIndex = 0;
                int leafIndexInNeighbor = -1;
                for (vector< pair<int,int> >::iterator entryIt = descendantsForSide.begin();
                     entryIt != descendantsForSide.end();  entryIt++, descendantIndex++)
                {
                  if (entryIt->first == cellID)
                  {
                    leafIndexInNeighbor = GDAMaximumRule2D::neighborChildPermutation(descendantIndex, descendantsForSide.size());
                    break;
                  }
                }
                if (leafIndexInNeighbor == -1)
                {
                  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Could not determine leafIndexInNeigbor.");
                }
                // check whether the basis is the right size:
                MultiBasis<>* neighborMultiBasis = (MultiBasis<>*) neighbor->elementType()->trialOrderPtr->getBasis(trialID,ancestralSideIndexInNeighbor).get();
                BasisPtr neighborLeafBasis = neighborMultiBasis->getLeafBasis(leafIndexInNeighbor);
                if (numBasisDofs != neighborLeafBasis->getCardinality())
                {
                  success = false;
                  cout << "FAILURE: cellID " << cellID << "'s basis for trialID " << trialID;
                  cout << " along sideIndex " << sideIndex << " has cardinality " << numBasisDofs;
                  cout << ", but neighbor leaf basis along that side (cellID " << neighbor->cellID();
                  cout << ", sideIndex " << ancestralSideIndexInNeighbor;
                  cout << ", leaf node " << leafIndexInNeighbor;
                  cout << ") has cardinality " << neighborLeafBasis->getCardinality() << endl;
                }
                else
                {
                  // cardinalities match, check that global dofs line up
                  for (int dofOrdinal = 0; dofOrdinal < numBasisDofs; dofOrdinal++)
                  {
                    int permutedDofOrdinal = GDAMaximumRule2D::neighborDofPermutation(dofOrdinal, numBasisDofs);
                    int neighborDofOrdinal = neighborMultiBasis->relativeToAbsoluteDofOrdinal(permutedDofOrdinal,
                                             leafIndexInNeighbor);
                    int neighborLocalDofIndex = neighbor->elementType()->trialOrderPtr->getDofIndex(trialID, neighborDofOrdinal,ancestralSideIndexInNeighbor);
                    int neighborGlobalDofIndex = mesh->globalDofIndex(neighbor->cellID(), neighborLocalDofIndex);
                    int myLocalDofIndex = elem->elementType()->trialOrderPtr->getDofIndex(trialID, dofOrdinal, sideIndex);
                    int myGlobalDofIndex = mesh->globalDofIndex(cellID, myLocalDofIndex);
                    if (neighborGlobalDofIndex != myGlobalDofIndex)
                    {
                      success = false;
                      cout << "FAILURE: checkDofConnectivities--(cellID, localDofIndex) : (" << cellID << ", ";
                      cout << myLocalDofIndex << ") != (";
                      cout << neighbor->cellID() << ", " << neighborLocalDofIndex << ") -- ";
                      cout << myGlobalDofIndex << " != " << neighborGlobalDofIndex << "\n";
                    }
                  }
                }
              }

bool MeshTestUtility::checkLocalGlobalConsistency(MeshPtr mesh, double tol)
{
  bool success = true;

  set<GlobalIndexType> cellIDs = mesh->getActiveCellIDs();

  GlobalDofAssignmentPtr gda = mesh->globalDofAssignment();

  // TODO: make this only check the locally-owned cells (right now does the whole mesh on every rank)

  int numGlobalDofs = gda->globalDofCount();
  FieldContainer<double> globalCoefficients(numGlobalDofs);
  for (int i=0; i<numGlobalDofs; i++)
  {
    globalCoefficients(i) = 2*i + 1; // arbitrary cofficients
  }
  FieldContainer<double> globalCoefficientsExpected = globalCoefficients;
  FieldContainer<double> globalCoefficientsActual(numGlobalDofs);
  FieldContainer<double> localCoefficients;

  Epetra_SerialComm Comm;
  Epetra_BlockMap map(numGlobalDofs, 1, 0, Comm);
  Epetra_Vector globalCoefficientsVector(map);
  for (int i=0; i<numGlobalDofs; i++)
  {
    globalCoefficientsVector[i] = globalCoefficients(i);
  }

  cellIDs = mesh->getActiveCellIDs();
  for (set<GlobalIndexType>::iterator cellIDIt = cellIDs.begin(); cellIDIt != cellIDs.end(); cellIDIt++)
  {
    GlobalIndexType cellID = *cellIDIt;
    gda->interpretGlobalCoefficients(cellID, localCoefficients, globalCoefficientsVector);
    FieldContainer<GlobalIndexType> globalDofIndices;
    FieldContainer<double> globalCoefficientsForCell;

    DofOrderingPtr trialOrder = mesh->getElementType(cellID)->trialOrderPtr;
    set<int> varIDs = trialOrder->getVarIDs();

    for (set<int>::iterator varIt=varIDs.begin(); varIt != varIDs.end(); varIt++)
    {
      int varID = *varIt;
      const vector<int>* sidesForVar = &trialOrder->getSidesForVarID(varID);
      for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++)
      {
        int sideOrdinal = *sideIt;

        BasisPtr basis;
        if (sidesForVar->size() == 1)
        {
          basis = trialOrder->getBasis(varID);
        }
        else
        {
          basis = trialOrder->getBasis(varID,sideOrdinal);
        }

        FieldContainer<double> basisCoefficients(basis->getCardinality());

        for (int dofOrdinal=0; dofOrdinal<basis->getCardinality(); dofOrdinal++)
        {
          int localDofIndex;
          if (sidesForVar->size() == 1)
          {
            localDofIndex = trialOrder->getDofIndex(varID, dofOrdinal);
          }
          else
          {
            localDofIndex = trialOrder->getDofIndex(varID, dofOrdinal, sideOrdinal);
          }

          basisCoefficients(dofOrdinal) = localCoefficients(localDofIndex);
        }

        gda->interpretLocalBasisCoefficients(cellID, varID, sideOrdinal,
                                             basisCoefficients, globalCoefficientsForCell, globalDofIndices);

        //        if ( (cellID==2) && (sideOrdinal==1) && (varID==0) ) {
        //          cout << "DEBUGGING: (cellID==2) && (sideOrdinal==1).\n";
        //          cout << "globalDofIndices:\n" << globalDofIndices;
        //          cout << "globalCoefficientsForCell:\n" << globalCoefficientsForCell;
        //          cout << "basisCoefficients:\n" << basisCoefficients;
        //        }

        for (int dofOrdinal=0; dofOrdinal < globalDofIndices.size(); dofOrdinal++)
        {
          GlobalIndexType dofIndex = globalDofIndices(dofOrdinal);
          globalCoefficientsActual(dofIndex) = globalCoefficientsForCell(dofOrdinal);

          double diff = abs(globalCoefficientsForCell(dofOrdinal) - globalCoefficientsExpected(dofIndex)) / globalCoefficientsExpected(dofIndex);
          if (diff > tol)
          {
            cout << "In mapping for cell " << cellID << " and var " << varID << " on side " << sideOrdinal << ", ";
            cout << "expected coefficient for global dof index " << dofIndex << " to be " << globalCoefficientsExpected(dofIndex);
            cout << ", but was " << globalCoefficientsForCell(dofOrdinal);
            cout << " (relative diff = " << diff << "; tol = " << tol << ")\n";
            success = false;
          }
        }
      }
//.........这里部分代码省略.........

FCPtr BasisEvaluation::getValues(BasisPtr basis, Camellia::EOperator op,
                                 const FieldContainer<double> &refPoints)
{
  int numPoints = refPoints.dimension(0);
  int spaceDim = refPoints.dimension(1);  // points dimensions are (numPoints, spaceDim)
  int basisCardinality = basis->getCardinality();
  int componentOfInterest = -1;
  Camellia::EFunctionSpace fs = basis->functionSpace();
  Intrepid::EOperator relatedOp = relatedOperator(op, fs, spaceDim, componentOfInterest);

  int opCardinality = spaceDim; // for now, we assume basis values are in the same spaceDim as points (e.g. vector 1D has just 1 component)
  
  bool relatedOpIsDkOperator = (relatedOp >= OPERATOR_D1) && (relatedOp <= OPERATOR_D10);
  if (relatedOpIsDkOperator)
  {
    opCardinality = Intrepid::getDkCardinality(relatedOp, spaceDim);
  }
  
  if ((Camellia::EOperator)relatedOp != op)
  {
    // we can assume relatedResults has dimensions (numPoints,basisCardinality,spaceDimOut)
    FCPtr relatedResults = Teuchos::rcp(new FieldContainer<double>(basisCardinality,numPoints,opCardinality));
    basis->getValues(*(relatedResults.get()), refPoints, relatedOp);
    FCPtr result = getComponentOfInterest(relatedResults,op,fs,componentOfInterest);
    if ( result.get() == 0 )
    {
      result = relatedResults;
    }
    return result;
  }
  // if we get here, we should have a standard Intrepid operator, in which case we should
  // be able to: size a FieldContainer appropriately, and then call basis->getValues

  // But let's do just check that we have a standard Intrepid operator
  if ( (op >= Camellia::OP_X) || (op <  Camellia::OP_VALUE) )
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Unknown operator.");
  }

  // result dimensions should be either (numPoints,basisCardinality) or (numPoints,basisCardinality,spaceDimOut);
  Teuchos::Array<int> dimensions;
  dimensions.push_back(basisCardinality);
  dimensions.push_back(numPoints);
  int basisRank = basis->rangeRank();
  if (((basisRank == 0) && (op ==  Camellia::OP_VALUE)))
  {
    // scalar; leave as is
  }
  else if ((basisRank == 0) && relatedOpIsDkOperator)
  {
    dimensions.push_back(opCardinality);
  }
  else if ( ( ( basisRank == 1) && (op ==  Camellia::OP_VALUE) ) )
  {
    dimensions.push_back(opCardinality);
  }
  else if (
    ( ( basisRank == 0) && (op == Camellia::OP_GRAD) )
    ||
    ( ( basisRank == 0) && (op == Camellia::OP_CURL) ) )
  {
    dimensions.push_back(spaceDim);
  }
  else if ( (basis->rangeRank() == 1) && (op == Camellia::OP_GRAD) )
  {
    // grad of vector: a tensor
    dimensions.push_back(spaceDim);
    dimensions.push_back(opCardinality);
  }
  FCPtr result = Teuchos::rcp(new FieldContainer<double>(dimensions));
  basis->getValues(*(result.get()), refPoints, (Intrepid::EOperator)op);
  return result;
}

//.........这里部分代码省略.........
        break;
      case Camellia::FUNCTION_SPACE_VECTOR_HVOL:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD_DISC: // shouldn't take the transform so late
      default:
        TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
        break;
    }
      break;
    case(Intrepid::OPERATOR_DIV):
      switch(fs)
    {
      case Camellia::FUNCTION_SPACE_HDIV:
      case Camellia::FUNCTION_SPACE_HDIV_DISC:
      case Camellia::FUNCTION_SPACE_HDIV_FREE:
        fst::HDIVtransformDIV<double>(*transformedValues,basisCache->getJacobianDet(),*referenceValues);
        break;
      case Camellia::FUNCTION_SPACE_VECTOR_HVOL:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD:
      case Camellia::FUNCTION_SPACE_VECTOR_HGRAD_DISC: // shouldn't take the transform so late
      default:
        TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
        break;
    }
      break;
    case(Intrepid::OPERATOR_D2):
      switch(fs)
    {
      case Camellia::FUNCTION_SPACE_HVOL:
      case Camellia::FUNCTION_SPACE_HGRAD:
      case Camellia::FUNCTION_SPACE_HGRAD_DISC:
      {
        // first term in the sum is:
        //      J^{-1} * OPD2 * J^{-T}
        // where J^{-1} is the inverse Jacabian, and OPD2 is the matrix of reference-space values formed from OP_D2
        
        const FieldContainer<double> *J_inv = &basisCache->getJacobianInv();
        transformedValues->initialize(0.0);
        int basisCardinality = basis->getCardinality();
        
        Teuchos::Array<int> multiplicities(spaceDim);
        Teuchos::Array<int> multiplicitiesTransformed(spaceDim);
        int dkCardinality = Intrepid::getDkCardinality(Intrepid::OPERATOR_D2, spaceDim);
        int numPoints = referenceValues->dimension(1);
        
        for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
        {
          for (int fieldOrdinal=0; fieldOrdinal<basisCardinality; fieldOrdinal++)
          {
            for (int pointOrdinal=0; pointOrdinal<numPoints; pointOrdinal++)
            {
              for (int dkOrdinal=0; dkOrdinal<dkCardinality; dkOrdinal++) // ref values dkOrdinal
              {
                double refValue = (*referenceValues)(fieldOrdinal,pointOrdinal,dkOrdinal);
                Intrepid::getDkMultiplicities(multiplicities, dkOrdinal, Intrepid::OPERATOR_D2, spaceDim);
                int k,l; // ref space coordinate directions for refValue (columns in J_inv)
                getD2_ij_FromMultiplicities(k,l,multiplicities);
                
                for (int dkOrdinalTransformed=0; dkOrdinalTransformed<dkCardinality; dkOrdinalTransformed++) // physical values dkOrdinal
                {
                  Intrepid::getDkMultiplicities(multiplicitiesTransformed, dkOrdinalTransformed, Intrepid::OPERATOR_D2, spaceDim);
                  int i,j; // physical space coordinate directions for refValue (rows in J_inv)
                  getD2_ij_FromMultiplicities(i,j,multiplicitiesTransformed);
                  double J_inv_ik = (*J_inv)(cellOrdinal,pointOrdinal,i,k);
                  double J_inv_jl = (*J_inv)(cellOrdinal,pointOrdinal,j,l);
                  (*transformedValues)(cellOrdinal,fieldOrdinal,pointOrdinal,dkOrdinalTransformed) += refValue * J_inv_ik * J_inv_jl;
                }
              }
            }
          }
        }
        
        // do we need the second term in the sum?  (So far, we don't support this)
        if (!basisCache->neglectHessian())
        {
          // then we need to do include the gradient of the basis times the (inverse) of the Hessian
          // this seems complicated, and we don't need it for computing the Laplacian in physical space
          // (at least not unless we have curvilinear geometry) so we don't support it for now...
          TEUCHOS_TEST_FOR_EXCEPTION(!basisCache->neglectHessian(), std::invalid_argument, "Support for the Hessian of the reference-to-physical mapping not yet implemented.");
        }
      }
        break;
      default:
        TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
        break;
    }
      break;

    default:
      TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "unhandled transformation");
      break;
  }
  
  FCPtr result = getComponentOfInterest(transformedValues,op,fs,componentOfInterest);
  if ( result.get() == 0 )
  {
    result = transformedValues;
  }
  return result;
}

void ParametricSurface::basisWeightsForEdgeInterpolant(FieldContainer<double> &edgeInterpolationCoefficients,
    VectorBasisPtr basis,
    MeshPtr mesh, int cellID)
{
  vector< ParametricCurvePtr > curves = mesh->parametricEdgesForCell(cellID);
  Teuchos::RCP<TransfiniteInterpolatingSurface> exactSurface = Teuchos::rcp( new TransfiniteInterpolatingSurface(curves) );
  exactSurface->setNeglectVertices(false);

  int basisDegree = basis->getDegree();
  shards::CellTopology line_2(shards::getCellTopologyData<shards::Line<2> >() );
  BasisPtr basis1D = BasisFactory::basisFactory()->getBasis(basisDegree, line_2.getKey(),
                     Camellia::FUNCTION_SPACE_HGRAD);

  BasisPtr compBasis = basis->getComponentBasis();
  int numComponents = basis->getNumComponents();
  if (numComponents != 2)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Only 2D surfaces supported right now");
  }

  edgeInterpolationCoefficients.resize(basis->getCardinality());

  set<int> edgeNodeFieldIndices = BasisFactory::basisFactory()->sideFieldIndices(basis,true); // true: include vertex dofs

  FieldContainer<double> dofCoords(compBasis->getCardinality(),2);
  IntrepidBasisWrapper< double, Intrepid::FieldContainer<double> >* intrepidBasisWrapper = dynamic_cast< IntrepidBasisWrapper< double, Intrepid::FieldContainer<double> >* >(compBasis.get());
  if (!intrepidBasisWrapper)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "compBasis does not appear to be an instance of IntrepidBasisWrapper");
  }
  Basis_HGRAD_QUAD_Cn_FEM<double, Intrepid::FieldContainer<double> >* intrepidBasis = dynamic_cast< Basis_HGRAD_QUAD_Cn_FEM<double, Intrepid::FieldContainer<double> >* >(intrepidBasisWrapper->intrepidBasis().get());
  if (!intrepidBasis)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "IntrepidBasisWrapper does not appear to wrap Basis_HGRAD_QUAD_Cn_FEM.");
  }
  intrepidBasis->getDofCoords(dofCoords);

  int edgeDim = 1;
  int vertexDim = 0;

  // set vertex dofs:
  for (int vertexIndex=0; vertexIndex<curves.size(); vertexIndex++)
  {
    double x = exactSurface->vertices()[vertexIndex].first;
    double y = exactSurface->vertices()[vertexIndex].second;
    int compDofOrdinal = compBasis->getDofOrdinal(vertexDim, vertexIndex, 0);
    int basisDofOrdinal_x = basis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 0);
    int basisDofOrdinal_y = basis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, 1);
    edgeInterpolationCoefficients[basisDofOrdinal_x] = x;
    edgeInterpolationCoefficients[basisDofOrdinal_y] = y;
  }

  for (int edgeIndex=0; edgeIndex<curves.size(); edgeIndex++)
  {
    bool edgeDofsFlipped = edgeIndex >= 2; // because Intrepid's ordering of dofs on the quad is not CCW but tensor-product, we need to flip for the opposite edges
    // (what makes things worse is that the vertex/edge numbering *is* CCW)
    if (curves.size() != 4)
    {
      cout << "WARNING: have not worked out the rule for flipping or not flipping edge dofs for anything but quads.\n";
    }
    double edgeLength = curves[edgeIndex]->linearLength();

    //    cout << "edgeIndex " << edgeIndex << endl;
    for (int comp=0; comp<numComponents; comp++)
    {
      FieldContainer<double> basisCoefficients_comp;
      bool useH1ForEdgeInterpolant = true; // an experiment
      curves[edgeIndex]->projectionBasedInterpolant(basisCoefficients_comp, basis1D, comp, edgeLength, useH1ForEdgeInterpolant);
      //      cout << "for edge " << edgeIndex << " and comp " << comp << ", projection-based interpolant dofs:\n";
      //      cout << basisCoefficients_comp;
      ////      cout << "basis dof coords:\n" << dofCoords;
      //      int basisDofOrdinal = basis->getDofOrdinalFromComponentDofOrdinal(v0_dofOrdinal_comp, comp);
      //      edgeInterpolationCoefficients[basisDofOrdinal] = basisCoefficients_comp[v0_dofOrdinal_1D];

      if (compBasis->getDegree() >= 2)   // then there are some "middle" nodes on the edge
      {
        // get the first dofOrdinal for the edge, so we can check the number of edge basis functions
        int firstEdgeDofOrdinal = compBasis->getDofOrdinal(edgeDim, edgeIndex, 0);

        //        cout << "first edge dofOrdinal: " << firstEdgeDofOrdinal << endl;

        int numEdgeDofs = compBasis->getDofTag(firstEdgeDofOrdinal)[3];
        if (numEdgeDofs != basis1D->getCardinality() - 2)
        {
          TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "numEdgeDofs does not match 1D basis cardinality");
        }
        for (int edgeDofOrdinal=0; edgeDofOrdinal<numEdgeDofs; edgeDofOrdinal++)
        {
          // determine the index into basisCoefficients_comp:
          int edgeDofOrdinalIn1DBasis = edgeDofsFlipped ? numEdgeDofs - 1 - edgeDofOrdinal : edgeDofOrdinal;
          int dofOrdinal1D = basis1D->getDofOrdinal(edgeDim, 0, edgeDofOrdinalIn1DBasis);
          // determine the ordinal of the edge dof in the component basis:
          int compDofOrdinal = compBasis->getDofOrdinal(edgeDim, edgeIndex, edgeDofOrdinal);
          // now, determine its ordinal in the vector basis
          int basisDofOrdinal = basis->getDofOrdinalFromComponentDofOrdinal(compDofOrdinal, comp);

          //          cout << "edge dof ordinal " << edgeDofOrdinal << " has basis weight " << basisCoefficients_comp[dofOrdinal1D] << " for component " << comp << endl;
          //          cout << "node on cell is at (" << dofCoords(compDofOrdinal,0) << ", " << dofCoords(compDofOrdinal,1) << ")\n";
          //          cout << "mapping to basisDofOrdinal " << basisDofOrdinal << endl;

//.........这里部分代码省略.........