C++ Array::end方法代码示例

Adapter *
Adapter::initialize(const Teuchos::RCP<Teuchos::ParameterList>& catalystParams)
{
  // Validate parameters against list for this specific class
  catalystParams->validateParameters(*getValidAdapterParameters(),0);

  if (Adapter::instance)
    delete Adapter::instance;
  Adapter::instance = new Adapter();

  // Register our Grid class with Catalyst so that it can be used in a pipeline.
  if (vtkClientServerInterpreterInitializer *intInit =
      vtkClientServerInterpreterInitializer::GetInitializer()) {
    if (vtkClientServerInterpreter *interp = intInit->GetGlobalInterpreter()) {
      interp->AddNewInstanceFunction("Grid", Private::MakeGrid);
    }
  }

  // Load pipeline file
  Teuchos::Array<std::string> files =
      catalystParams->get<Teuchos::Array<std::string> >("Pipeline Files");
  typedef Teuchos::Array<std::string>::const_iterator FileIterT;
  for (FileIterT it = files.begin(), itEnd = files.end(); it != itEnd; ++it)
    Adapter::instance->addPythonScriptPipeline(*it);

  return Adapter::instance;
}

Teuchos::RCP<Epetra_Map>
Albany::SolutionResponseFunction::
buildCulledMap(const Epetra_Map& x_map,
	       const Teuchos::Array<int>& keepDOF) const
{
  int numKeepDOF = std::accumulate(keepDOF.begin(), keepDOF.end(), 0);
  int Neqns = keepDOF.size();
  int N = x_map.NumMyElements(); // x_map is map for solution vector

  TEUCHOS_ASSERT( !(N % Neqns) ); // Assume that all the equations for
                                  // a given node are on the assigned
                                  // processor. I.e. need to ensure
                                  // that N is exactly Neqns-divisible

  int nnodes = N / Neqns;          // number of fem nodes
  int N_new = nnodes * numKeepDOF; // length of local x_new

  int *gids = x_map.MyGlobalElements(); // Fill local x_map into gids array
  Teuchos::Array<int> gids_new(N_new);
  int idx = 0;
  for ( int inode = 0; inode < N/Neqns ; ++inode) // For every node
    for ( int ieqn = 0; ieqn < Neqns; ++ieqn )  // Check every dof on the node
      if ( keepDOF[ieqn] == 1 )  // then want to keep this dof
	gids_new[idx++] = gids[(inode*Neqns)+ieqn];
  // end cull

  Teuchos::RCP<Epetra_Map> x_new_map =
    Teuchos::rcp( new Epetra_Map( -1, N_new, &gids_new[0], 0, x_map.Comm() ) );

  return x_new_map;
}

Teuchos::Array<GO>
Albany::NodeGIDsSolutionCullingStrategy::
selectedGIDsT(Teuchos::RCP<const Tpetra_Map> sourceMapT) const
{
  Teuchos::Array<GO> result;
  {
    Teuchos::Array<GO> mySelectedGIDs;

    // Subract 1 to convert exodus GIDs to our GIDs
    for (int i=0; i<nodeGIDs_.size(); i++)
      if (sourceMapT->isNodeGlobalElement(nodeGIDs_[i] -1) ) mySelectedGIDs.push_back(nodeGIDs_[i] - 1);

    Teuchos::RCP<const Teuchos::Comm<int> >commT = sourceMapT->getComm(); 
    {
      GO selectedGIDCount;
      {
        GO mySelectedGIDCount = mySelectedGIDs.size();
        Teuchos::reduceAll<LO, GO>(*commT, Teuchos::REDUCE_SUM, 1, &mySelectedGIDCount, &selectedGIDCount); 
      }
      result.resize(selectedGIDCount);
    }

    const int ierr = Tpetra::GatherAllV(
        commT,
        mySelectedGIDs.getRawPtr(), mySelectedGIDs.size(),
        result.getRawPtr(), result.size());
    TEUCHOS_ASSERT(ierr == 0);
  }

  std::sort(result.begin(), result.end());

  return result;
}

Teuchos::Array<EpetraGlobalIndex> getMyLIDs(
    const Epetra_BlockMap &map,
    const Teuchos::ArrayView<const EpetraGlobalIndex> &selectedGIDs)
{
  Teuchos::Array<EpetraGlobalIndex> sortedMyGIDs(map.MyGlobalElements(), map.MyGlobalElements() + map.NumMyElements());
  std::sort(sortedMyGIDs.begin(), sortedMyGIDs.end());

  Teuchos::Array<EpetraGlobalIndex> sortedSelectedGIDs(selectedGIDs);
  std::sort(sortedSelectedGIDs.begin(), sortedSelectedGIDs.end());

  Teuchos::Array<EpetraGlobalIndex> mySelectedGIDs;
  std::set_intersection(sortedMyGIDs.begin(), sortedMyGIDs.end(),
                        sortedSelectedGIDs.begin(), sortedSelectedGIDs.end(),
                        std::back_inserter(mySelectedGIDs));

  Teuchos::Array<EpetraGlobalIndex> result;
  result.reserve(mySelectedGIDs.size());

  std::transform(
      mySelectedGIDs.begin(), mySelectedGIDs.end(),
      std::back_inserter(result),
      std::bind1st(std::mem_fun_ref(static_cast<int(Epetra_BlockMap::*)(EpetraGlobalIndex) const>(&Epetra_BlockMap::LID)), map));

  return result;
}

void epetraFromThyra(
  const Teuchos::RCP<const Epetra_Comm> &comm,
  const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > &thyraResponses,
  const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > &thyraSensitivities,
  Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > &responses,
  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > &sensitivities)
{
  responses.clear();
  responses.reserve(thyraResponses.size());
  typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > ThyraResponseArray;
  for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
      it_end = thyraResponses.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    responses.push_back(epetraVectorFromThyra(comm, *it));
  }

  sensitivities.clear();
  sensitivities.reserve(thyraSensitivities.size());
  typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > ThyraSensitivityArray;
  for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
      it_end = thyraSensitivities.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    ThyraSensitivityArray::const_reference sens_thyra = *it;
    Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > sens;
    sens.reserve(sens_thyra.size());
    for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
        jt_end = sens_thyra.end();
        jt != jt_end;
        ++jt) {
        sens.push_back(epetraMultiVectorFromThyra(comm, *jt));
    }
    sensitivities.push_back(sens);
  }
}

size_type computeSize(const Teuchos::Array< DIM_TYPE > & dimensions)
{
  // In the MDArray<T>(const MDArrayView<T> &) constructor, I try to
  // pass the MDArrayView dimensions to computeSize(), but they come
  // in as ArrayView<const T> (for reasons I can't determine) and
  // cause all sorts of const-correctness problems.  So I copy them
  // into a new Array<T> and pass its view to the main computeSize()
  // function.  Fortunately, the array of dimensions is small.
  Teuchos::Array< DIM_TYPE > nonConstDims(0);
  nonConstDims.insert(nonConstDims.begin(),
                      dimensions.begin(),
                      dimensions.end());
  return computeSize(nonConstDims());
}

BoundingBox GeometryManager<Geometry,GlobalOrdinal>::localBoundingBox() const
{
    double global_x_min = Teuchos::ScalarTraits<double>::rmax();
    double global_y_min = Teuchos::ScalarTraits<double>::rmax();
    double global_z_min = Teuchos::ScalarTraits<double>::rmax();
    double global_x_max = -Teuchos::ScalarTraits<double>::rmax();
    double global_y_max = -Teuchos::ScalarTraits<double>::rmax();
    double global_z_max = -Teuchos::ScalarTraits<double>::rmax();

    // Get the local bounding boxes compute the local bounding box.
    BoundingBox local_box;
    Teuchos::Tuple<double,6> box_bounds;
    Teuchos::Array<BoundingBox> boxes = boundingBoxes();
    Teuchos::Array<BoundingBox>::const_iterator box_iterator;
    DTK_CHECK( !boxes.empty() );
    for ( box_iterator = boxes.begin();
	  box_iterator != boxes.end();
	  ++box_iterator )
    {
	box_bounds = box_iterator->getBounds();

	if ( box_bounds[0] < global_x_min )
	{
	    global_x_min = box_bounds[0];
	}
	if ( box_bounds[1] < global_y_min )
	{
	    global_y_min = box_bounds[1];
	}
	if ( box_bounds[2] < global_z_min )
	{
	    global_z_min = box_bounds[2];
	}
	if ( box_bounds[3] > global_x_max )
	{
	    global_x_max = box_bounds[3];
	}
	if ( box_bounds[4] > global_y_max )
	{
	    global_y_max = box_bounds[4];
	}
	if ( box_bounds[5] > global_z_max )
	{
	    global_z_max = box_bounds[5];
	}
    }

    return BoundingBox( global_x_min, global_y_min, global_z_min,
			global_x_max, global_y_max, global_z_max );
}

void tpetraFromThyra(
  const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<ST> > > &thyraResponses,
  const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<ST> > > > &thyraSensitivities,
  Teuchos::Array<Teuchos::RCP<const Tpetra_Vector> > &responses,
  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > > &sensitivities)
{
  responses.clear();
  responses.reserve(thyraResponses.size());
  typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<ST> > > ThyraResponseArray;
  for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
      it_end = thyraResponses.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    responses.push_back(Teuchos::nonnull(*it) ? ConverterT::getConstTpetraVector(*it) : Teuchos::null);
  }

  sensitivities.clear();
  sensitivities.reserve(thyraSensitivities.size());
  typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<ST> > > > ThyraSensitivityArray;
  for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
      it_end = thyraSensitivities.end(),
      it = it_begin;
      it != it_end;
      ++it) {
  ThyraSensitivityArray::const_reference sens_thyra = *it;
    Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > sens;
    sens.reserve(sens_thyra.size());
    for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
        jt_end = sens_thyra.end();
        jt != jt_end;
        ++jt) {
        sens.push_back(Teuchos::nonnull(*jt) ? ConverterT::getConstTpetraMultiVector(*jt) : Teuchos::null);
    }
    sensitivities.push_back(sens);
  }
}

  bool IsInRangeList(const Integral val, const Teuchos::Array<Integral> &valList, bool sorted=true)
{
  if (allValuesAreInRangeList(valList))
    return true;
  else if (noValuesAreInRangeList(valList))
    return false;

  if (sorted){
    typename Teuchos::Array<Integral>::const_iterator flag = valList.end();
    --flag;
    if (std::binary_search(valList.begin(), flag, val))
      return true;
    else
      return false;
  }
  else{
    for (typename Teuchos::Array<Integral>::size_type i=0; i < valList.size()-1; i++){
      if (valList[i] == val)
        return true;
    }
    return false;
  }
}

//.........这里部分代码省略.........
    const size_t nVtx = model->getLocalNumVertices();
    // RCM constructs invPerm, not perm
    ArrayRCP<lno_t> invPerm = solution->getPermutationRCP(true);
  
    // Check if there are actually edges to reorder.
    // If there are not, then just use the natural ordering.
    if (numEdges == 0) {
      for (size_t i = 0; i < nVtx; ++i) {
        invPerm[i] = i;
      }
      return 0;
    }
  
    // Set the label of each vertex to invalid.
    Tpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
    for (size_t i = 0; i < nVtx; ++i) {
      invPerm[i] = INVALID;
    }
  
    // Loop over all connected components.
    // Do BFS within each component.
    lno_t root;
    std::queue<lno_t> Q;
    size_t count = 0; // CM label, reversed later
    size_t next = 0;  // next unmarked vertex
    Teuchos::Array<std::pair<lno_t, size_t> >  children; // children and their degrees
  
    while (count < nVtx) {
  
      // Find suitable root vertex for this component.
      // First find an unmarked vertex, use to find root in next component.
      while ((next < nVtx) && (static_cast<Tpetra::global_size_t>(invPerm[next]) != INVALID)) next++;

      // Select root method. Pseudoperipheral usually gives the best
      // ordering, but the user may choose a faster method.
      std::string root_method = pl->get("root_method", "pseudoperipheral");
      if (root_method == string("first"))
        root = next;
      else if (root_method == string("smallest_degree"))
        root = findSmallestDegree(next, nVtx, edgeIds, offsets);
      else if (root_method == string("pseudoperipheral"))
        root = findPseudoPeripheral(next, nVtx, edgeIds, offsets);
      else {
        // This should never happen if pl was validated.
      }

      // Label connected component starting at root
      Q.push(root);
      //cout << "Debug: invPerm[" << root << "] = " << count << endl;
      invPerm[root] = count++;
  
      while (Q.size()){
        // Get a vertex from the queue
        lno_t v = Q.front();
        Q.pop();
        //cout << "Debug: v= " << v << ", offsets[v] = " << offsets[v] << endl;
  
        // Add unmarked children to list of pairs, to be added to queue.
        children.resize(0);
        for (lno_t ptr = offsets[v]; ptr < offsets[v+1]; ++ptr){
          lno_t child = edgeIds[ptr];
          if (static_cast<Tpetra::global_size_t>(invPerm[child]) == INVALID){
            // Not visited yet; add child to list of pairs.
            std::pair<lno_t,size_t> newchild;
            newchild.first = child;
            newchild.second = offsets[child+1] - offsets[child];
            children.push_back(newchild); 
          }
        }
        // Sort children by increasing degree
        // TODO: If edge weights, sort children by decreasing weight,
        SortPairs<lno_t,size_t> zort;
        zort.sort(children);

        typename Teuchos::Array<std::pair<lno_t,size_t> >::iterator it = children.begin ();
        for ( ; it != children.end(); ++it){
          // Push children on the queue in sorted order.
          lno_t child = it->first;
          invPerm[child] = count++; // Label as we push on Q
          Q.push(child);
          //cout << "Debug: invPerm[" << child << "] = " << count << endl;
        }
      }
    }
  
    // Reverse labels for RCM
    bool reverse = true; // TODO: Make parameter
    if (reverse) {
      lno_t temp;
      for (size_t i=0; i < nVtx/2; ++i) {
        // Swap (invPerm[i], invPerm[nVtx-i])
        temp = invPerm[i];
        invPerm[i] = invPerm[nVtx-1-i];
        invPerm[nVtx-1-i] = temp;
      }
    }
  
    solution->setHaveInverse(true);
    return ierr;
  }

ATO::StiffnessObjectiveBase<EvalT, Traits>::
StiffnessObjectiveBase(Teuchos::ParameterList& p,
		  const Teuchos::RCP<Albany::Layouts>& dl,
                  const Albany::MeshSpecsStruct* meshSpecs) :
  qp_weights ("Weights", dl->qp_scalar     ),
  BF         ("BF",      dl->node_qp_scalar)

{
  Teuchos::ParameterList* responseParams = p.get<Teuchos::ParameterList*>("Parameter List");
  elementBlockName = meshSpecs->ebName;

  m_excludeBlock = false;
  if(responseParams->isType<Teuchos::Array<std::string>>("Blocks")){
    Teuchos::Array<std::string> 
      blocks = responseParams->get<Teuchos::Array<std::string>>("Blocks");
    if(find(blocks.begin(),blocks.end(),elementBlockName) == blocks.end()){
      m_excludeBlock = true;

      stiffness_objective_tag =
        Teuchos::rcp(new PHX::Tag<ScalarT>(className, dl->dummy));
      this->addEvaluatedField(*stiffness_objective_tag);
      return;
    }
  }

  ATO::PenaltyModelFactory<ScalarT> penaltyFactory;
  penaltyModel = penaltyFactory.create(p, dl, elementBlockName);

  Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem =
    p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");

  topologies = paramsFromProblem->get<Teuchos::RCP<TopologyArray> >("Topologies");

  int nTopos = topologies->size();
  FName = responseParams->get<std::string>("Response Name");
  std::string dFdpBaseName = responseParams->get<std::string>("Response Derivative Name");
  dFdpNames.resize(nTopos);
  for(int itopo=0; itopo<nTopos; itopo++){
    dFdpNames[itopo] = Albany::strint(dFdpBaseName,itopo);
  }

  this->pStateMgr = p.get< Albany::StateManager* >("State Manager Ptr");
  this->pStateMgr->registerStateVariable(FName, dl->workset_scalar, dl->dummy, 
                                         "all", "scalar", 0.0, false, false);
  for(int itopo=0; itopo<nTopos; itopo++){
    this->pStateMgr->registerStateVariable(dFdpNames[itopo], 
                                           dl->node_scalar, dl->dummy, 
                                           "all", "scalar", 0.0, false, false);
  }

  this->addDependentField(qp_weights);
  this->addDependentField(BF);

  Teuchos::Array< PHX::MDField<ScalarT> > depFields;
  penaltyModel->getDependentFields(depFields);

  int nFields = depFields.size();
  for(int ifield=0; ifield<nFields; ifield++)
    this->addDependentField(depFields[ifield].fieldTag());

  stiffness_objective_tag =
    Teuchos::rcp(new PHX::Tag<ScalarT>(className, dl->dummy));
  this->addEvaluatedField(*stiffness_objective_tag);

}

 //! Iterator pointing to end of array
 const_iterator end() const { 
   return const_iterator(indices.end(), values.end()); 
 }

int main(int argc, char *argv[])
{
  Teuchos::oblackholestream blackhole;
  Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);

  using Teuchos::TimeMonitor;
  using TpetraExamples::FDStencil;
  using TpetraExamples::ScaleKernel;

  // 
  // Get the default communicator and node
  //
  auto &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
  auto comm = platform.getComm();
  auto node = platform.getNode();
  const int myImageID = comm->getRank();
  const int numImages = comm->getSize();

  //
  // Get example parameters from command-line processor
  //  
  bool verbose = (myImageID==0);
  int numGlobal_user = 100*comm->getSize();
  int numTimeTrials = 3;
  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("global-size",&numGlobal_user,"Global test size.");
  cmdp.setOption("num-time-trials",&numTimeTrials,"Number of trials in timing loops.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }

  // 
  // Say hello, print some communicator info
  //
  if (verbose) {
    std::cout << "\n" << Tpetra::version() << std::endl;
    std::cout << "Comm info: " << *comm;
    std::cout << "Node type: " << Teuchos::typeName(*node) << std::endl;
    std::cout << std::endl;
  }

  // 
  // Create a simple map for domain and range
  // 
  Tpetra::global_size_t numGlobalRows = numGlobal_user;
  auto map = Tpetra::createUniformContigMapWithNode<int,int>(numGlobalRows, comm, node);
  // const size_t numLocalRows = map->getNodeNumElements();
  auto x = Tpetra::createVector<double>(map),
       y = Tpetra::createVector<double>(map);

  // Create a simple diagonal operator using lambda function
  auto fTwoOp = Tpetra::RTI::binaryOp<double>( [](double /*y*/, double x) { return 2.0 * x; } , map );
  // y = 3*fTwoOp*x + 2*y = 3*2*1 + 2*1 = 8
  x->putScalar(1.0);
  y->putScalar(1.0);
  fTwoOp->apply( *x, *y, Teuchos::NO_TRANS, 3.0, 2.0 );
  // check that y == eights
  double norm = y->norm1();
  if (verbose) {
    std::cout << "Tpetra::RTI::binaryOp" << std::endl
              << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
              << ", expected value is " << numGlobalRows * 8.0 << std::endl;
  }

  // Create the same diagonal operator using a Kokkos kernel
  auto kTwoOp = Tpetra::RTI::kernelOp<double>( ScaleKernel<double>(2.0), map );
  // y = 0.5*kTwop*x + 0.75*y = 0.5*2*1 + 0.75*8 = 7
  kTwoOp->apply( *x, *y, Teuchos::NO_TRANS, 0.5, 0.75 );
  // check that y == sevens
  norm = y->norm1();
  if (verbose) {
    std::cout << "Tpetra::RTI::kernelOp" << std::endl
              << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
              << ", expected value is " << numGlobalRows * 7.0 << std::endl;
  }

  //
  // Create a finite-difference stencil using a Kokkos kernel and non-trivial maps
  //
  decltype(map) colmap;
  if (numImages > 1) {
    Teuchos::Array<int>           colElements;
    Teuchos::ArrayView<const int> rowElements = map->getNodeElementList();
    // This isn't the most efficient Map/Import layout, but it makes for a very straight-forward kernel
    if (myImageID != 0) colElements.push_back( map->getMinGlobalIndex() - 1 );
    colElements.insert(colElements.end(), rowElements.begin(), rowElements.end());
    if (myImageID != numImages-1) colElements.push_back( map->getMaxGlobalIndex() + 1 );
    colmap = Tpetra::createNonContigMapWithNode<int,int>(colElements(), comm, node);
  }
  else {
    colmap = map;
  }
  auto importer = createImport(map,colmap);
  // Finite-difference kernel = tridiag(-1, 2, -1)
  FDStencil<double> kern(myImageID, numImages, map->getNodeNumElements(), -1.0, 2.0, -1.0);
  auto FDStencilOp = Tpetra::RTI::kernelOp<double>( kern, map, map, importer );
  // x = ones(), FD(x) = [1 zeros() 1]
  auto timeFDStencil = TimeMonitor::getNewTimer("FD RTI Stencil");
  {
//.........这里部分代码省略.........

/*
 * \brief This constructor will pull the mesh data DTK needs out of Moab,
 * partition it for the example, and build a DataTransferKit::MeshContainer
 * object from the local data in the partition. You can directly write the
 * traits interface yourself, but this is probably the easiest way to get
 * started (although potentially inefficient).
 */
MoabMesh::MoabMesh( const RCP_Comm& comm,
		    const std::string& filename,
		    const moab::EntityType& block_topology,
		    const int partitioning_type )
    : d_comm( comm )
{
    // Compute the node dimension.
    int node_dim = 0;
    if ( block_topology == moab::MBTRI )
    {
	node_dim = 2;
    }
    else if ( block_topology == moab::MBQUAD )
    {
	node_dim = 2;
    }
    else if ( block_topology == moab::MBTET )
    {
	node_dim = 3;
    }
    else if ( block_topology == moab::MBHEX )
    {
	node_dim = 3;
    }
    else if ( block_topology == moab::MBPYRAMID )
    {
	node_dim = 3;
    }
    else
    {
	node_dim = 0;
    }

    // Create a moab instance.
    moab::ErrorCode error;
    d_moab = Teuchos::rcp( new moab::Core() );

    std::cout<<"Filename: "<<filename<<std::endl;

    // Load the mesh.
    d_moab->load_mesh( &filename[0] );
    moab::EntityHandle root_set = d_moab->get_root_set();

    // Extract the elements with this block's topology.
    std::vector<moab::EntityHandle> global_elements;
    error = d_moab->get_entities_by_type(
	root_set, block_topology, global_elements );
    assert( error == moab::MB_SUCCESS );

    std::cout<<"Global elements: "<<global_elements.size()<<std::endl;

    // Partition the mesh.
    int comm_size = d_comm->getSize();
    int comm_rank = d_comm->getRank();

    // Get the number of nodes in an element.
    std::vector<moab::EntityHandle> elem_vertices;
    error = d_moab->get_adjacencies( &global_elements[0],
				     1,
				     0,
				     false,
				     elem_vertices );
    assert( error == moab::MB_SUCCESS );
    int nodes_per_element = elem_vertices.size();

    // Get the global element coordinates.
    std::vector<double> global_coords;
    error = d_moab->get_vertex_coordinates( global_coords );
    assert( error == moab::MB_SUCCESS );

    // Get the global max and min values for the coordinates. This problem is
    // symmetric.
    double min = *(std::min_element( global_coords.begin(),
				     global_coords.end() ) );
    double max = *(std::max_element( global_coords.begin(),
				     global_coords.end() ) );
    double width = max - min;

    Teuchos::Array<moab::EntityHandle> elements;
    elem_vertices.resize( nodes_per_element );
    std::vector<double> elem_coords( 3*nodes_per_element );
    std::vector<moab::EntityHandle>::const_iterator global_elem_iterator;
    for ( global_elem_iterator = global_elements.begin();
	  global_elem_iterator != global_elements.end();
	  ++global_elem_iterator )
    {
	// Get the individual element vertices.
	error = d_moab->get_adjacencies( &*global_elem_iterator,
					 1,
					 0,
					 false,
					 elem_vertices );
	assert( error == moab::MB_SUCCESS );
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  Teuchos::oblackholestream blackhole;
  Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);

  // 
  // Get the default communicator and node
  //
  auto &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
  auto comm = platform.getComm();
  auto node = platform.getNode();
  const int myImageID = comm->getRank();
  const int numImages = comm->getSize();
  const bool verbose = (myImageID==0);

  // 
  // Say hello, print some communicator info
  //
  if (verbose) {
    std::cout << "\n" << Tpetra::version() << std::endl;
    std::cout << "Comm info: " << *comm;
    std::cout << "Node type: " << Teuchos::typeName(*node) << std::endl;
    std::cout << std::endl;
  }

  // 
  // Create a simple map for domain and range
  // 
  Tpetra::global_size_t numGlobalRows = 1000*numImages;
  auto map = Tpetra::createUniformContigMapWithNode<int,int>(numGlobalRows, comm, node);
  auto x = Tpetra::createVector<double>(map),
       y = Tpetra::createVector<double>(map);

  // Create a simple diagonal operator using lambda function
  auto fTwoOp = Tpetra::RTI::binaryOp<double>( [](double /*y*/, double x) { return 2.0 * x; } , map );
  // y = 3*fTwoOp*x + 2*y = 3*2*1 + 2*1 = 8
  x->putScalar(1.0);
  y->putScalar(1.0);
  fTwoOp->apply( *x, *y, Teuchos::NO_TRANS, 3.0, 2.0 );
  // check that y == eights
  double norm = y->norm1();
  if (verbose) {
    std::cout << "Tpetra::RTI::binaryOp" << std::endl
              << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
              << ", expected value is " << numGlobalRows * 8.0 << std::endl;
  }

  //
  // Create a finite-difference stencil using a Kokkos kernel and non-trivial maps
  //
  decltype(map) colmap;
  if (numImages > 1) {
    Teuchos::Array<int>           colElements;
    Teuchos::ArrayView<const int> rowElements = map->getNodeElementList();
    // This isn't the most efficient Map/Import layout, but it makes for a very straight-forward kernel
    if (myImageID != 0) colElements.push_back( map->getMinGlobalIndex() - 1 );
    colElements.insert(colElements.end(), rowElements.begin(), rowElements.end());
    if (myImageID != numImages-1) colElements.push_back( map->getMaxGlobalIndex() + 1 );
    colmap = Tpetra::createNonContigMapWithNode<int,int>(colElements(), comm, node);
  }
  else {
    colmap = map;
  }
  auto importer = createImport(map,colmap);
  // Finite-difference kernel = tridiag(-1, 2, -1)
  FDStencil<double> kern(myImageID, numImages, map->getNodeNumElements(), -1.0, 2.0, -1.0);
  auto FDStencilOp = Tpetra::RTI::kernelOp<double>( kern, map, map, importer );
  // x = ones(), FD(x) = [1 zeros() 1]
  FDStencilOp->apply( *x, *y );
  norm = y->norm1();
  if (verbose) {
    std::cout << std::endl
              << "TpetraExamples::FDStencil" << std::endl
              << "norm(y) result == " << std::setprecision(2) << std::scientific << norm 
              << ", expected value is " << 2.0 << std::endl;
  }

  std::cout << "\nEnd Result: TEST PASSED" << std::endl;
  return 0;
}