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

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);
//.........这里部分代码省略.........

void Boundary::bcsToImpose(FieldContainer<GlobalIndexType> &globalIndices,
                           FieldContainer<Scalar> &globalValues, TBC<Scalar> &bc,
                           DofInterpreter* dofInterpreter)
{
  set< GlobalIndexType > rankLocalCells = _mesh->cellIDsInPartition();
  map< GlobalIndexType, double> bcGlobalIndicesAndValues;

  for (GlobalIndexType cellID : rankLocalCells)
  {
    bcsToImpose(bcGlobalIndicesAndValues, bc, cellID, dofInterpreter);
  }
  singletonBCsToImpose(bcGlobalIndicesAndValues, bc, dofInterpreter);
  
  // ****** New, tag-based BC imposition follows ******
  map< GlobalIndexType, double> bcTagGlobalIndicesAndValues;
  
  map< int, vector<pair<VarPtr, TFunctionPtr<Scalar>>>> tagBCs = bc.getDirichletTagBCs(); // keys are tags
  
  MeshTopology* meshTopo = dynamic_cast<MeshTopology*>(_mesh->getTopology().get());
  
  TEUCHOS_TEST_FOR_EXCEPTION(!meshTopo, std::invalid_argument, "pure MeshTopologyViews are not yet supported by new tag-based BC imposition");
  
  for (auto tagBC : tagBCs)
  {
    int tagID = tagBC.first;
    
    vector<EntitySetPtr> entitySets = meshTopo->getEntitySetsForTagID(DIRICHLET_SET_TAG_NAME, tagID);
    for (EntitySetPtr entitySet : entitySets)
    {
      // get rank-local cells that match the entity set:
      set<IndexType> matchingCellIDs = entitySet->cellIDsThatMatch(_mesh->getTopology(), rankLocalCells);
      for (IndexType cellID : matchingCellIDs)
      {
        ElementTypePtr elemType = _mesh->getElementType(cellID);
        BasisCachePtr basisCache = BasisCache::basisCacheForCell(_mesh, cellID);
        
        for (auto varFunctionPair : tagBC.second)
        {
          VarPtr var = varFunctionPair.first;
          FunctionPtr f = varFunctionPair.second;
          
          vector<int> sideOrdinals = elemType->trialOrderPtr->getSidesForVarID(var->ID());
          
          for (int sideOrdinal : sideOrdinals)
          {
            BasisPtr basis = elemType->trialOrderPtr->getBasis(var->ID(), sideOrdinal);
            bool isVolume = basis->domainTopology()->getDimension() == _mesh->getDimension();
            for (int d=0; d<_mesh->getDimension(); d++)
            {
              vector<unsigned> matchingSubcells;
              if (isVolume)
                matchingSubcells = entitySet->subcellOrdinals(_mesh->getTopology(), cellID, d);
              else
              {
                CellTopoPtr cellTopo = elemType->cellTopoPtr;
                int sideDim = cellTopo->getDimension() - 1;
                vector<unsigned> matchingSubcellsOnSide = entitySet->subcellOrdinalsOnSide(_mesh->getTopology(), cellID, sideOrdinal, d);
                for (unsigned sideSubcellOrdinal : matchingSubcellsOnSide)
                {
                  unsigned cellSubcellOrdinal = CamelliaCellTools::subcellOrdinalMap(cellTopo, sideDim, sideOrdinal, d, sideSubcellOrdinal);
                  matchingSubcells.push_back(cellSubcellOrdinal);
                }
              }
              
              if (matchingSubcells.size() == 0) continue; // nothing to impose
              
              /*
               What follows - projecting the function onto the basis on the whole domain - is more expensive than necessary,
               in the general case: we can do the projection on just the matching subcells, and if we had a way of taking the
               restriction of a basis to a subcell of the domain, then we could avoid computing the whole basis as well.
               
               But for now, this should work, and it's simple to implement.
               */
              BasisCachePtr basisCacheForImposition = isVolume ? basisCache : basisCache->getSideBasisCache(sideOrdinal);
              int numCells = 1;
              FieldContainer<double> basisCoefficients(numCells,basis->getCardinality());
              Projector<double>::projectFunctionOntoBasisInterpolating(basisCoefficients, f, basis, basisCacheForImposition);
              basisCoefficients.resize(basis->getCardinality());
              
              set<GlobalIndexType> matchingGlobalIndices;
              for (unsigned matchingSubcell : matchingSubcells)
              {
                set<GlobalIndexType> subcellGlobalIndices = dofInterpreter->globalDofIndicesForVarOnSubcell(var->ID(),cellID,d,matchingSubcell);
                matchingGlobalIndices.insert(subcellGlobalIndices.begin(),subcellGlobalIndices.end());
              }
              
              FieldContainer<double> globalData;
              FieldContainer<GlobalIndexType> globalDofIndices;
//              dofInterpreter->interpretLocalBasisCoefficients(cellID, var->ID(), sideOrdinal, basisCoefficientsToImpose, globalData, globalDofIndices);
              dofInterpreter->interpretLocalBasisCoefficients(cellID, var->ID(), sideOrdinal, basisCoefficients, globalData, globalDofIndices);
              for (int globalDofOrdinal=0; globalDofOrdinal<globalDofIndices.size(); globalDofOrdinal++)
              {
                GlobalIndexType globalDofIndex = globalDofIndices(globalDofOrdinal);
                if (matchingGlobalIndices.find(globalDofIndex) != matchingGlobalIndices.end())
                  bcTagGlobalIndicesAndValues[globalDofIndex] = globalData(globalDofOrdinal);
              }
            }
          }
        }
      }
//.........这里部分代码省略.........

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 TRieszRep<Scalar>::computeRepresentationValues(FieldContainer<Scalar> &values, int testID, Camellia::EOperator op, BasisCachePtr basisCache)
{

  if (_repsNotComputed)
  {
    cout << "Computing riesz rep dofs" << endl;
    computeRieszRep();
  }

  int spaceDim = _mesh->getTopology()->getDimension();
  int numCells = values.dimension(0);
  int numPoints = values.dimension(1);
  vector<GlobalIndexType> cellIDs = basisCache->cellIDs();

  // all elems coming in should be of same type
  ElementTypePtr elemTypePtr = _mesh->getElementType(cellIDs[0]);
  DofOrderingPtr testOrderingPtr = elemTypePtr->testOrderPtr;
  int numTestDofsForVarID = testOrderingPtr->getBasisCardinality(testID, VOLUME_INTERIOR_SIDE_ORDINAL);
  BasisPtr testBasis = testOrderingPtr->getBasis(testID);

  bool testBasisIsVolumeBasis = (spaceDim == testBasis->domainTopology()->getDimension());
  bool useCubPointsSideRefCell = testBasisIsVolumeBasis && basisCache->isSideCache();

  Teuchos::RCP< const FieldContainer<double> > transformedBasisValues = basisCache->getTransformedValues(testBasis,op,useCubPointsSideRefCell);

  int rank = values.rank() - 2; // if values are shaped as (C,P), scalar...
  if (rank > 1)
  {
    cout << "ranks greater than 1 not presently supported...\n";
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "ranks greater than 1 not presently supported...");
  }

//  Camellia::print("cellIDs",cellIDs);

  values.initialize(0.0);
  for (int cellIndex = 0; cellIndex<numCells; cellIndex++)
  {
    int cellID = cellIDs[cellIndex];
    for (int j = 0; j<numTestDofsForVarID; j++)
    {
      int dofIndex = testOrderingPtr->getDofIndex(testID, j);
      for (int i = 0; i<numPoints; i++)
      {
        if (rank==0)
        {
          double basisValue = (*transformedBasisValues)(cellIndex,j,i);
          values(cellIndex,i) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
        }
        else
        {
          for (int d = 0; d<spaceDim; d++)
          {
            double basisValue = (*transformedBasisValues)(cellIndex,j,i,d);
            values(cellIndex,i,d) += basisValue*_rieszRepDofsGlobal[cellID](dofIndex);
          }
        }
      }
    }
  }
//  TestSuite::serializeOutput("rep values", values);
}