C++ IPPtr::computeInnerProductMatrix方法代码示例

double RieszRep::computeAlternativeNormSqOnCell(IPPtr ip, ElementPtr elem){
  GlobalIndexType cellID = elem->cellID();
  Teuchos::RCP<DofOrdering> testOrdering= elem->elementType()->testOrderPtr;
  bool testVsTest = true;
  Teuchos::RCP<BasisCache> basisCache =   BasisCache::basisCacheForCell(_mesh, cellID, testVsTest,1);

  int numDofs = testOrdering->totalDofs();
  FieldContainer<double> ipMat(1,numDofs,numDofs);
  ip->computeInnerProductMatrix(ipMat,testOrdering,basisCache);

  double sum = 0.0;
  for (int i = 0;i<numDofs;i++){
    for (int j = 0;j<numDofs;j++){
      sum += _rieszRepDofsGlobal[cellID](i)*_rieszRepDofsGlobal[cellID](j)*ipMat(0,i,j);
    }
  }
  
  return sum;
}

//.........这里部分代码省略.........
      cout << "Test failed: old Burgers stiffness differs from new; maxDiff " << maxDiff << ".\n";
      cout << "Old: \n" << expectedValues;
      cout << "New: \n" << actualValues;
      cout << "TrialDofOrdering: \n" << *trialOrdering;
      cout << "TestDofOrdering:\n" << *testOrdering;
    }
    else
    {
      cout << "Old and new Burgers stiffness agree!!\n";
    }
  }

  // define our inner product:
  // Teuchos::RCP<BurgersInnerProduct> ip = Teuchos::rcp( new BurgersInnerProduct( bf, mesh ) );

  // function to scale the squared guy by epsilon/h
  FunctionPtr epsilonOverHScaling = Teuchos::rcp( new EpsilonScaling(epsilon) );
  IPPtr ip = Teuchos::rcp( new IP );
  ip->addTerm(tau);
  ip->addTerm(tau->div());
  ip->addTerm( epsilonOverHScaling * v );
  ip->addTerm( sqrt(sqrt(epsilon)) * v->grad() );
  ip->addTerm( beta * v->grad() );

  // use old IP instead, for now...
  Teuchos::RCP<BurgersInnerProduct> oldIP = Teuchos::rcp( new BurgersInnerProduct( oldBurgersBF, mesh ) );

  expectedValues.resize(numCells, testOrdering->totalDofs(), testOrdering->totalDofs() );
  actualValues.resize  (numCells, testOrdering->totalDofs(), testOrdering->totalDofs() );

  BasisCachePtr ipBasisCache = Teuchos::rcp( new BasisCache(elemType, true) ); // true: test vs. test
  ipBasisCache->setPhysicalCellNodes(quadPoints,vector<int>(1),false); // false: don't create side cache

  oldIP->computeInnerProductMatrix(expectedValues,testOrdering,ipBasisCache);
  ip->computeInnerProductMatrix(actualValues,testOrdering,ipBasisCache);

  tol = 1e-14;
  maxDiff = 0.0;
  if (rank==0)
  {
    if ( ! TestSuite::fcsAgree(expectedValues,actualValues,tol,maxDiff) )
    {
      cout << "Test failed: old inner product differs from new IP; maxDiff " << maxDiff << ".\n";
      cout << "Old: \n" << expectedValues;
      cout << "New IP: \n" << actualValues;
      cout << "testOrdering: \n" << *testOrdering;
    }
    else
    {
      cout << "Old inner product and new IP agree!!\n";
    }
  }

  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  // the RHS as implemented by BurgersProblem divides the first component of beta by 2.0
  // so we follow that.  I've not done the math; just imitating the code...
  Teuchos::RCP<RHSEasy> otherRHS = Teuchos::rcp( new RHSEasy );
  vector<double> e1_div2 = e1;
  e1_div2[0] /= 2.0;
  FunctionPtr rhsBeta = (e1_div2 * beta * e1 + Teuchos::rcp( new ConstantVectorFunction( e2 ) )) * u_prev;
  otherRHS->addTerm( rhsBeta * v->grad() - u_prev * tau->div() );

  Teuchos::RCP<BurgersProblem> problem = Teuchos::rcp( new BurgersProblem(oldBurgersBF) );

  expectedValues.resize(numCells, testOrdering->totalDofs() );
  actualValues.resize  (numCells, testOrdering->totalDofs() );

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