C++ SolutionPtr::setSolnCoeffsForCellID方法代码示例

bool FunctionTests::testBasisSumFunction()
{
  bool success = true;
  // on a single-element mesh, the BasisSumFunction should be identical to
  // the Solution with those coefficients

  // define a new mesh: more interesting if we're not on the ref cell
  int spaceDim = 2;
  FieldContainer<double> quadPoints(4,2);

  quadPoints(0,0) = 0.0; // x1
  quadPoints(0,1) = 0.0; // y1
  quadPoints(1,0) = 2.0;
  quadPoints(1,1) = 0.0;
  quadPoints(2,0) = 1.0;
  quadPoints(2,1) = 1.0;
  quadPoints(3,0) = 0.0;
  quadPoints(3,1) = 1.0;

  int H1Order = 1, pToAdd = 0;
  int horizontalCells = 1, verticalCells = 1;

  // create a pointer to a new mesh:
  MeshPtr spectralConfusionMesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells,
                                  _confusionBF, H1Order, H1Order+pToAdd);

  BCPtr bc = BC::bc();
  SolutionPtr soln = Teuchos::rcp( new Solution(spectralConfusionMesh, bc) );

  soln->initializeLHSVector();

  int cellID = 0;
  double tol = 1e-16; // overly restrictive, just for now.

  DofOrderingPtr trialSpace = spectralConfusionMesh->getElementType(cellID)->trialOrderPtr;
  set<int> trialIDs = trialSpace->getVarIDs();

  BasisCachePtr volumeCache = BasisCache::basisCacheForCell(spectralConfusionMesh, cellID);

  for (set<int>::iterator trialIt=trialIDs.begin(); trialIt != trialIDs.end(); trialIt++)
  {
    int trialID = *trialIt;
    const vector<int>* sidesForVar = &trialSpace->getSidesForVarID(trialID);
    bool boundaryValued = sidesForVar->size() != 1;
    // note that for volume trialIDs, sideIndex = 0, and numSides = 1…
    for (vector<int>::const_iterator sideIt = sidesForVar->begin(); sideIt != sidesForVar->end(); sideIt++)
    {
      int sideIndex = *sideIt;

      BasisCachePtr sideCache = volumeCache->getSideBasisCache(sideIndex);
      BasisPtr basis = trialSpace->getBasis(trialID, sideIndex);
      int basisCardinality = basis->getCardinality();
      for (int basisOrdinal = 0; basisOrdinal<basisCardinality; basisOrdinal++)
      {
        FieldContainer<double> basisCoefficients(basisCardinality);
        basisCoefficients(basisOrdinal) = 1.0;
        soln->setSolnCoeffsForCellID(basisCoefficients, cellID, trialID, sideIndex);

        VarPtr v = Var::varForTrialID(trialID, spectralConfusionMesh->bilinearForm());
        FunctionPtr solnFxn = Function::solution(v, soln, false);
        FunctionPtr basisSumFxn = Teuchos::rcp( new BasisSumFunction(basis, basisCoefficients, Teuchos::rcp((BasisCache*)NULL), OP_VALUE, boundaryValued) );
        if (!boundaryValued)
        {
          double l2diff = (solnFxn - basisSumFxn)->l2norm(spectralConfusionMesh);
//          cout << "l2diff = " << l2diff << endl;
          if (l2diff > tol)
          {
            success = false;
            cout << "testBasisSumFunction: l2diff of " << l2diff << " exceeds tol of " << tol << endl;
            cout << "l2norm of basisSumFxn: " << basisSumFxn->l2norm(spectralConfusionMesh) << endl;
            cout << "l2norm of solnFxn: " << solnFxn->l2norm(spectralConfusionMesh) << endl;
          }
          l2diff = (solnFxn->dx() - basisSumFxn->dx())->l2norm(spectralConfusionMesh);
          //          cout << "l2diff = " << l2diff << endl;
          if (l2diff > tol)
          {
            success = false;
            cout << "testBasisSumFunction: l2diff of dx() " << l2diff << " exceeds tol of " << tol << endl;
            cout << "l2norm of basisSumFxn->dx(): " << basisSumFxn->dx()->l2norm(spectralConfusionMesh) << endl;
            cout << "l2norm of solnFxn->dx(): " << solnFxn->dx()->l2norm(spectralConfusionMesh) << endl;
          }

          // test that the restriction to a side works
          int numSides = volumeCache->cellTopology()->getSideCount();
          for (int i=0; i<numSides; i++)
          {
            BasisCachePtr mySideCache = volumeCache->getSideBasisCache(i);
            if (! solnFxn->equals(basisSumFxn, mySideCache, tol))
            {
              success = false;
              cout << "testBasisSumFunction: on side 0, l2diff of " << l2diff << " exceeds tol of " << tol << endl;
              reportFunctionValueDifferences(solnFxn, basisSumFxn, mySideCache, tol);
            }
            if (! solnFxn->grad(spaceDim)->equals(basisSumFxn->grad(spaceDim), mySideCache, tol))
            {
              success = false;
              cout << "testBasisSumFunction: on side 0, l2diff of dx() " << l2diff << " exceeds tol of " << tol << endl;
              reportFunctionValueDifferences(solnFxn->grad(spaceDim), basisSumFxn->grad(spaceDim), mySideCache, tol);
            }
          }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // make residual for riesz representation function
  LinearTermPtr residual = Teuchos::rcp(new LinearTerm);// residual
  FunctionPtr parity = Function::sideParity();
  residual->addTerm(-fnhatFxn*v + (beta*uFxn)*v->grad());
  Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(mesh, ip, residual));
  riesz->computeRieszRep();
  map<int,FunctionPtr> err_rep_map;
  err_rep_map[v->ID()] = RieszRep::repFunction(v,riesz);

  ////////////////////   GET BOUNDARY CONDITION DATA    ///////////////////////

  FieldContainer<GlobalIndexType> bcGlobalIndices;
  FieldContainer<double> bcGlobalValues;
  mesh->boundary().bcsToImpose(bcGlobalIndices,bcGlobalValues,*(solution->bc()), NULL);
  set<int> bcInds;
  for (int i=0; i<bcGlobalIndices.dimension(0); i++)
  {
    bcInds.insert(bcGlobalIndices(i));
  }

  ////////////////////   CHECK GALERKIN ORTHOGONALITY   ///////////////////////

  BCPtr nullBC;
  RHSPtr nullRHS;
  IPPtr nullIP;
  SolutionPtr solnPerturbation = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );

  map< int, vector<DofInfo> > infoMap = constructGlobalDofToLocalDofInfoMap(mesh);

  for (map< int, vector<DofInfo> >::iterator mapIt = infoMap.begin();
       mapIt != infoMap.end(); mapIt++)
  {
    int dofIndex = mapIt->first;
    vector< DofInfo > dofInfoVector = mapIt->second; // all the local dofs that map to dofIndex
    // create perturbation in direction du
    solnPerturbation->clear(); // clear all solns
    // set each corresponding local dof to 1.0
    for (vector< DofInfo >::iterator dofInfoIt = dofInfoVector.begin();
         dofInfoIt != dofInfoVector.end(); dofInfoIt++)
    {
      DofInfo info = *dofInfoIt;
      FieldContainer<double> solnCoeffs(info.basisCardinality);
      solnCoeffs(info.basisOrdinal) = 1.0;
      solnPerturbation->setSolnCoeffsForCellID(solnCoeffs, info.cellID, info.trialID, info.sideIndex);
    }
    //    solnPerturbation->setSolnCoeffForGlobalDofIndex(1.0,dofIndex);

    LinearTermPtr b_du =  convectionBF->testFunctional(solnPerturbation);
    FunctionPtr gradient = b_du->evaluate(err_rep_map, TestingUtilities::isFluxOrTraceDof(mesh,dofIndex)); // use boundary part only if flux
    double grad = gradient->integrate(mesh,10);
    if (!TestingUtilities::isFluxOrTraceDof(mesh,dofIndex) && abs(grad)>tol)  // if we're not single-precision zero FOR FIELDS
    {
      //      int cellID = mesh->getGlobalToLocalMap()[dofIndex].first;
      cout << "Failed testGalerkinOrthogonality() for fields with diff " << abs(grad) << " at dof " << dofIndex << "; info:" << endl;
      cout << dofInfoString(infoMap[dofIndex]);
      success = false;
    }
  }
  FieldContainer<double> errorJumps(mesh->numGlobalDofs()); //initialized to zero
  // just test fluxes ON INTERNAL SKELETON here
  set<GlobalIndexType> activeCellIDs = mesh->getActiveCellIDsGlobal();
  for (GlobalIndexType activeCellID : activeCellIDs)
  {
    ElementPtr elem = mesh->getElement(activeCellID);
    for (int sideIndex = 0; sideIndex < 4; sideIndex++)
    {
      ElementTypePtr elemType = elem->elementType();
      vector<int> localDofIndices = elemType->trialOrderPtr->getDofIndices(beta_n_u->ID(), sideIndex);
      for (int i = 0; i<localDofIndices.size(); i++)
      {
        int globalDofIndex = mesh->globalDofIndex(elem->cellID(), localDofIndices[i]);
        vector< DofInfo > dofInfoVector = infoMap[globalDofIndex];

        solnPerturbation->clear();
        TestingUtilities::setSolnCoeffForGlobalDofIndex(solnPerturbation,1.0,globalDofIndex);
        // also add in BCs
        for (int i = 0; i<bcGlobalIndices.dimension(0); i++)
        {
          TestingUtilities::setSolnCoeffForGlobalDofIndex(solnPerturbation,bcGlobalValues(i),bcGlobalIndices(i));
        }

        LinearTermPtr b_du =  convectionBF->testFunctional(solnPerturbation);
        FunctionPtr gradient = b_du->evaluate(err_rep_map, TestingUtilities::isFluxOrTraceDof(mesh,globalDofIndex)); // use boundary part only if flux
        double jump = gradient->integrate(mesh,10);
        errorJumps(globalDofIndex) += jump;
      }
    }
  }
  for (int i = 0; i<mesh->numGlobalDofs(); i++)
  {
    if (abs(errorJumps(i))>tol)
    {
      cout << "Failing Galerkin orthogonality test for fluxes with diff " << errorJumps(i) << " at dof " << i << endl;
      cout << dofInfoString(infoMap[i]);
      success = false;
    }
  }

  return success;
}