C++ BFPtr::stiffnessMatrix方法代码示例

//.........这里部分代码省略.........
  bf->addTerm( -u, beta * v->grad());
  bf->addTerm( beta_n_u_minus_sigma_hat, v);

  // ==================== SET INITIAL GUESS ==========================
  mesh->registerSolution(backgroundFlow);

  map<int, Teuchos::RCP<AbstractFunction> > functionMap;
  functionMap[BurgersBilinearForm::U] = Teuchos::rcp(new InitialGuess());
  functionMap[BurgersBilinearForm::SIGMA_1] = Teuchos::rcp(new ZeroFunction());
  functionMap[BurgersBilinearForm::SIGMA_2] = Teuchos::rcp(new ZeroFunction());

  backgroundFlow->projectOntoMesh(functionMap);

  // ==================== END SET INITIAL GUESS ==========================
  // compare stiffness matrices for first linear step:
  int trialOrder = 1;
  pToAdd = 0;
  int testOrder = trialOrder + pToAdd;
  CellTopoPtr quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));
  DofOrderingFactory dofOrderingFactory(bf);
  DofOrderingPtr testOrdering = dofOrderingFactory.testOrdering(testOrder, *quadTopoPtr);
  DofOrderingPtr trialOrdering = dofOrderingFactory.trialOrdering(trialOrder, *quadTopoPtr);

  int numCells = 1;
  // just use testOrdering for both trial and test spaces (we only use to define BasisCache)
  ElementTypePtr elemType  = Teuchos::rcp( new ElementType(trialOrdering, testOrdering, quadTopoPtr) );
  BasisCachePtr basisCache = Teuchos::rcp( new BasisCache(elemType) );
  quadPoints.resize(1,quadPoints.dimension(0),quadPoints.dimension(1));
  basisCache->setPhysicalCellNodes(quadPoints,vector<int>(1),true); // true: do create side cache
  FieldContainer<double> cellSideParities(numCells,quadTopoPtr->getSideCount());
  cellSideParities.initialize(1.0); // not worried here about neighbors actually having opposite parity -- just want the two BF implementations to agree...
  FieldContainer<double> expectedValues(numCells, testOrdering->totalDofs(), trialOrdering->totalDofs() );
  FieldContainer<double> actualValues(numCells, testOrdering->totalDofs(), trialOrdering->totalDofs() );
  oldBurgersBF->stiffnessMatrix(expectedValues, elemType, cellSideParities, basisCache);
  bf->stiffnessMatrix(actualValues, elemType, cellSideParities, basisCache);

  // compare beta's as well
  FieldContainer<double> expectedBeta = oldBurgersBF->getBeta(basisCache);
  Teuchos::Array<int> dim;
  expectedBeta.dimensions(dim);
  FieldContainer<double> actualBeta(dim);
  beta->values(actualBeta,basisCache);

  double tol = 1e-14;
  double maxDiff;
  if (rank == 0)
  {
    if ( ! TestSuite::fcsAgree(expectedBeta,actualBeta,tol,maxDiff) )
    {
      cout << "Test failed: old Burgers beta differs from new; maxDiff " << maxDiff << ".\n";
      cout << "Old beta: \n" << expectedBeta;
      cout << "New beta: \n" << actualBeta;
    }
    else
    {
      cout << "Old and new Burgers beta agree!!\n";
    }

    if ( ! TestSuite::fcsAgree(expectedValues,actualValues,tol,maxDiff) )
    {
      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;
    }

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;
}