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

Teuchos::RCP<const Teuchos::ParameterList>
IfpackPreconditionerFactory::getValidParameters() const
{
  using Teuchos::rcp;
  using Teuchos::rcp_implicit_cast;
  typedef Teuchos::ParameterEntryValidator PEV;
  static Teuchos::RCP<Teuchos::ParameterList> validParamList;
  if(validParamList.get()==NULL) {
    validParamList = Teuchos::rcp(new Teuchos::ParameterList(Ifpack_name));
    {
      // Create the validator for the preconditioner type!
      Teuchos::Array<std::string>
        precTypeNames;
      precTypeNames.insert(
        precTypeNames.begin(),
        &Ifpack::precTypeNames[0],
        &Ifpack::precTypeNames[0] + Ifpack::numPrecTypes
        );
      Teuchos::Array<Ifpack::EPrecType>
        precTypeValues;
      precTypeValues.insert(
        precTypeValues.begin(),
        &Ifpack::precTypeValues[0],
        &Ifpack::precTypeValues[0] + Ifpack::numPrecTypes
        );
      precTypeValidator = rcp(
        new Teuchos::StringToIntegralParameterEntryValidator<Ifpack::EPrecType>(
          precTypeNames,precTypeValues,PrecType_name
          )
        );
    }
    validParamList->set(
      PrecType_name, PrecTypeName_default,
      "Type of Ifpack preconditioner to use.",
      rcp_implicit_cast<const PEV>(precTypeValidator)
      );
    validParamList->set(
      Overlap_name, Overlap_default,
      "Number of rows/columns overlapped between subdomains in different"
      "\nprocesses in the additive Schwarz-type domain-decomposition preconditioners."
      );
    validParamList->sublist(
      IfpackSettings_name, false,
      "Preconditioner settings that are passed onto the Ifpack preconditioners themselves."
      ).setParameters(Ifpack_GetValidParameters());
    // Note that in the above setParameterList(...) function that we actually
    // validate down into the first level of this sublist.  Really the
    // Ifpack_Preconditioner objects themselves should do validation but we do
    // it ourselves taking the return from the Ifpack_GetValidParameters()
    // function as gospel!
    Teuchos::setupVerboseObjectSublist(&*validParamList);
  }
  return validParamList;
}

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

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

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

/*!
 * \brief Find the set of entities a point neighbors.
 */
void CoarseLocalSearch::search( const Teuchos::ArrayView<const double>& point,
				const Teuchos::ParameterList& parameters,
				Teuchos::Array<Entity>& neighbors ) const
{
    // Find the leaf of nearest neighbors.
    int num_neighbors = 100;
    if ( parameters.isParameter("Coarse Local Search kNN") )
    {
	num_neighbors = parameters.get<int>("Coarse Local Search kNN");
    }
    num_neighbors = 
	std::min( num_neighbors, Teuchos::as<int>(d_entity_map.size()) );
    Teuchos::Array<unsigned> local_neighbors = 
	d_tree->nnSearch( point, num_neighbors );

    // Extract the neighbors.
    neighbors.resize( local_neighbors.size() );
    Teuchos::Array<unsigned>::const_iterator local_it;
    Teuchos::Array<Entity>::iterator entity_it;
    for ( local_it = local_neighbors.begin(),
	 entity_it = neighbors.begin();
	  local_it != local_neighbors.end();
	  ++local_it, ++entity_it )
    {
	DTK_CHECK( d_entity_map.count(*local_it) );
	*entity_it = d_entity_map.find(*local_it)->second;
    }
}

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

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

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

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

//---------------------------------------------------------------------------//
// Given a range entity id, get the ids of the domain entities that it mapped to.
void ParallelSearch::getDomainEntitiesFromRange( 
    const EntityId range_id,
    Teuchos::Array<EntityId>& domain_ids ) const
{
    DTK_REQUIRE( !d_empty_range );
    auto range_pair = d_range_to_domain_map.equal_range( range_id );
    domain_ids.resize( std::distance(range_pair.first,range_pair.second) );
    auto domain_it = domain_ids.begin();
    for ( auto range_it = range_pair.first;
	  range_it != range_pair.second;
	  ++range_it, ++domain_it )
    {
	*domain_it = range_it->second;
    }
}

//---------------------------------------------------------------------------//
// Given a domain entity id, get the ids of the range entities that mapped to it.
void ParallelSearch::getRangeEntitiesFromDomain( 
    const EntityId domain_id,
    Teuchos::Array<EntityId>& range_ids ) const
{
    DTK_REQUIRE( !d_empty_domain );
    auto domain_pair = d_domain_to_range_map.equal_range( domain_id );
    range_ids.resize( std::distance(domain_pair.first,domain_pair.second) );
    auto range_it = range_ids.begin();
    for ( auto domain_it = domain_pair.first;
	  domain_it != domain_pair.second;
	  ++domain_it, ++range_it )
    {
	*range_it = domain_it->second;
    }
}

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

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