C++ DofOrderingPtr::addEntry方法代码示例

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 FunctionTests::testValuesDottedWithTensor()
{
  bool success = true;

  vector< FunctionPtr > vectorFxns;

  double xValue = 3, yValue = 4;
  FunctionPtr simpleVector = Function::vectorize(Function::constant(xValue), Function::constant(yValue));
  vectorFxns.push_back(simpleVector);
  FunctionPtr x = Function::xn(1);
  FunctionPtr y = Function::yn(1);
  vectorFxns.push_back( Function::vectorize(x*x, x*y) );

  VGPStokesFormulation vgpStokes = VGPStokesFormulation(1.0);
  BFPtr bf = vgpStokes.bf();

  int h1Order = 1;
  MeshPtr mesh = MeshFactory::quadMesh(bf, h1Order);

  int cellID=0; // the only cell
  BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, cellID);

  for (int i=0; i<vectorFxns.size(); i++)
  {
    FunctionPtr vectorFxn_i = vectorFxns[i];
    for (int j=0; j<vectorFxns.size(); j++)
    {
      FunctionPtr vectorFxn_j = vectorFxns[j];
      FunctionPtr dotProduct = vectorFxn_i * vectorFxn_j;
      FunctionPtr expectedDotProduct = vectorFxn_i->x() * vectorFxn_j->x() + vectorFxn_i->y() * vectorFxn_j->y();
      if (! expectedDotProduct->equals(dotProduct, basisCache))
      {
        cout << "testValuesDottedWithTensor() failed: expected dot product does not match dotProduct.\n";
        success = false;
        double tol = 1e-14;
        reportFunctionValueDifferences(dotProduct, expectedDotProduct, basisCache, tol);
      }
    }
  }

  // now, let's try the same thing, but for a LinearTerm dot product
  VarFactoryPtr vf = VarFactory::varFactory();
  VarPtr v = vf->testVar("v", HGRAD);

  DofOrderingPtr dofOrdering = Teuchos::rcp( new DofOrdering(CellTopology::quad()) );
  shards::CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
  BasisPtr basis = BasisFactory::basisFactory()->getBasis(h1Order, quad_4.getKey(), Camellia::FUNCTION_SPACE_HGRAD);
  dofOrdering->addEntry(v->ID(), basis, v->rank());

  int numCells = 1;
  int numFields = basis->getCardinality();

  for (int i=0; i<vectorFxns.size(); i++)
  {
    FunctionPtr f_i = vectorFxns[i];
    LinearTermPtr lt_i = f_i * v;
    LinearTermPtr lt_i_x = f_i->x() * v;
    LinearTermPtr lt_i_y = f_i->y() * v;
    for (int j=0; j<vectorFxns.size(); j++)
    {
      FunctionPtr f_j = vectorFxns[j];
      LinearTermPtr lt_j = f_j * v;
      LinearTermPtr lt_j_x = f_j->x() * v;
      LinearTermPtr lt_j_y = f_j->y() * v;
      FieldContainer<double> values(numCells,numFields,numFields);
      lt_i->integrate(values, dofOrdering, lt_j, dofOrdering, basisCache);
      FieldContainer<double> values_expected(numCells,numFields,numFields);
      lt_i_x->integrate(values_expected,dofOrdering,lt_j_x,dofOrdering,basisCache);
      lt_i_y->integrate(values_expected,dofOrdering,lt_j_y,dofOrdering,basisCache);
      double tol = 1e-14;
      double maxDiff = 0;
      if (!fcsAgree(values, values_expected, tol, maxDiff))
      {
        cout << "FunctionTests::testValuesDottedWithTensor: ";
        cout << "dot product and sum of the products of scalar components differ by maxDiff " << maxDiff;
        cout << " in LinearTerm::integrate().\n";
        success = false;
      }
    }
  }

//  // finally, let's try the same sort of thing, but now with a vector-valued basis
//  BasisPtr vectorBasisTemp = BasisFactory::basisFactory()->getBasis(h1Order, quad_4.getKey(), Camellia::FUNCTION_SPACE_VECTOR_HGRAD);
//  VectorBasisPtr vectorBasis = Teuchos::rcp( (VectorizedBasis<double, FieldContainer<double> > *)vectorBasisTemp.get(),false);
//
//  BasisPtr compBasis = vectorBasis->getComponentBasis();
//
//  // create a new v, and a new dofOrdering
//  VarPtr v_vector = vf->testVar("v_vector", VECTOR_HGRAD);
//  dofOrdering = Teuchos::rcp( new DofOrdering );
//  dofOrdering->addEntry(v_vector->ID(), vectorBasis, v_vector->rank());
//
//  DofOrderingPtr dofOrderingComp = Teuchos::rcp( new DofOrdering );
//  dofOrderingComp->addEntry(v->ID(), compBasis, v->rank());
//

  return success;
}

bool LobattoBasisTests::testSimpleStiffnessMatrix() {
  bool success = true;
  
  int rank = Teuchos::GlobalMPISession::getRank();
  
  VarFactory varFactory;
  VarPtr u = varFactory.fieldVar("u");
  VarPtr un = varFactory.fluxVar("un_hat");
  VarPtr v = varFactory.testVar("v", HGRAD);
  
  BFPtr bf = Teuchos::rcp( new BF(varFactory) );
  vector<double> beta;
  beta.push_back(1);
  beta.push_back(1);
  bf->addTerm(beta * u, v->grad());
  bf->addTerm(un, v);
  
  DofOrderingPtr trialSpace = Teuchos::rcp( new DofOrdering );
  DofOrderingPtr testSpace = Teuchos::rcp( new DofOrdering );
  
  const int numSides = 4;
  const int spaceDim = 2;
  
  int fieldOrder = 3;
  int testOrder = fieldOrder+2;
  BasisPtr fieldBasis = Camellia::intrepidQuadHGRAD(fieldOrder);
  BasisPtr fluxBasis = Camellia::intrepidLineHGRAD(fieldOrder);
  trialSpace->addEntry(u->ID(), fieldBasis, fieldBasis->rangeRank());
  for (int i=0; i<numSides; i++) {
    trialSpace->addEntry(un->ID(), fluxBasis, fluxBasis->rangeRank(), i);
  }
  
  BasisPtr testBasis = Camellia::lobattoQuadHGRAD(testOrder+1,false); // +1 because it lives in HGRAD
  testSpace->addEntry(v->ID(), testBasis, testBasis->rangeRank());
  
  int numTrialDofs = trialSpace->totalDofs();
  int numTestDofs = testSpace->totalDofs();
  int numCells = 1;
  
  FieldContainer<double> cellNodes(numCells,numSides,spaceDim);
  cellNodes(0,0,0) = 0;
  cellNodes(0,0,1) = 0;
  cellNodes(0,1,0) = 1;
  cellNodes(0,1,1) = 0;
  cellNodes(0,2,0) = 1;
  cellNodes(0,2,1) = 1;
  cellNodes(0,3,0) = 0;
  cellNodes(0,3,1) = 1;
  
  FieldContainer<double> stiffness(numCells,numTestDofs,numTrialDofs);
  
  FieldContainer<double> cellSideParities(numCells,numSides);
  cellSideParities.initialize(1.0);
  
  Teuchos::RCP<shards::CellTopology> quad_4 = Teuchos::rcp( new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ) );
  Teuchos::RCP<ElementType> elemType = Teuchos::rcp( new ElementType(trialSpace, testSpace, quad_4));
  
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType) );
  vector<GlobalIndexType> cellIDs;
  cellIDs.push_back(0);
  basisCache->setPhysicalCellNodes(cellNodes, cellIDs, true);
  bf->stiffnessMatrix(stiffness, elemType, cellSideParities, basisCache);

  // TODO: finish this test
  
//  cout << stiffness;
  if (rank==0)
    cout << "Warning: testSimpleStiffnessMatrix() unfinished.\n";
  
  return success;
}

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

//    DofOrderingPtr testOrderTrace = mesh->getDofOrderingFactory().getRelabeledDofOrdering(traceOrder, traceTestMap);
//    testOrderTrace = mesh->getDofOrderingFactory().setBasisDegree(testOrderTrace, testDegreeTrace);

        int numTestInteriorDofs = testOrderInterior->totalDofs();
        int numTestTraceDofsIncludingInterior = testOrderTrace->totalDofs();
        int numTestTraceDofs = testTraceIndex;
        int numTestDofs = numTestTraceDofs + numTestInteriorDofs;

        timer.ResetStartTime();

        bool printTimings = true;

        if (printTimings) {
            cout << "numCells: " << numCells << endl;
            cout << "numTestDofs: " << numTestDofs << endl;
        }

        FieldContainer<double> rhsVectorTest(numCells,testOrderInterior->totalDofs()); // rhsVector is zero for the "trace" test dofs
        {
            // project the load f onto the space of interior test dofs.
            LinearTermPtr f = rhs->linearTerm();
            set<int> testIDs = f->varIDs();
            for (int eqn=0; eqn<numEquations; eqn++) {
                VarPtr v = _virtualTerms.getFieldTestVars()[eqn];

                if (testIDs.find(v->ID()) != testIDs.end()) {
                    BasisPtr testInteriorBasis = testOrderInterior->getBasis(v->ID());
                    FieldContainer<double> fValues(numCells,testInteriorBasis->getCardinality());
//          DofOrderingPtr oneVarOrderingTest = Teuchos::rcp(new DofOrdering(testInteriorBasis->domainTopology()));
                    DofOrderingPtr oneVarOrderingTest = Teuchos::rcp(new DofOrdering);
                    oneVarOrderingTest->addEntry(v->ID(), testInteriorBasis, testInteriorBasis->rangeRank());

                    LinearTermPtr f_v = Teuchos::rcp( new LinearTerm );
                    typedef pair< FunctionPtr, VarPtr > LinearSummand;
                    vector<LinearSummand> summands = f->summands();
                    for (int i=0; i<summands.size(); i++) {
                        FunctionPtr f = summands[i].first;
                        if (v->ID() == summands[i].second->ID()) {
                            f_v->addTerm(f * v);
                            f_v->integrate(fValues, oneVarOrderingTest, basisCache);
                        }
                    }

                    LinearTermPtr v_lt = 1.0 * v;
                    FieldContainer<double> l2(numCells,testInteriorBasis->getCardinality(),testInteriorBasis->getCardinality());
                    v_lt->integrate(l2,oneVarOrderingTest,v_lt,oneVarOrderingTest,basisCache,basisCache->isSideCache());

                    Teuchos::Array<int> testTestDim(2), testOneDim(2);
                    testTestDim[0] = testInteriorBasis->getCardinality();
                    testTestDim[1] = testInteriorBasis->getCardinality();
                    testOneDim[0] = testInteriorBasis->getCardinality();
                    testOneDim[1] = 1;
                    FieldContainer<double> projection(testOneDim);
                    for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++) {
                        FieldContainer<double> l2cell(testTestDim,&l2(cellOrdinal,0,0));
                        FieldContainer<double> f_cell(testOneDim,&fValues(cellOrdinal,0));

                        SerialDenseWrapper::solveSystemUsingQR(projection, l2cell, f_cell);

                        // rows in projection correspond to Ae_i, columns to the e_j.  I.e. projection coefficients for e_i are found in the ith row
                        for (int basisOrdinal_j=0; basisOrdinal_j<projection.dimension(0); basisOrdinal_j++) {
                            int testIndex = testOrderInterior->getDofIndex(v->ID(), basisOrdinal_j);
                            rhsVectorTest(cellOrdinal,testIndex) = projection(basisOrdinal_j,0);