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

Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > create_tridiag_graph(LocalOrdinal num_rows_per_proc)
{
  Teuchos::RCP<const Teuchos::Comm<int> > comm = getDefaultComm();

  const global_size_t INVALID = Teuchos::OrdinalTraits<global_size_t>::invalid();
  const LocalOrdinal indexBase = 0;

  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowmap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(INVALID, num_rows_per_proc, indexBase, comm));

  Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > crsgraph = Teuchos::rcp(new Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>(rowmap, 0));

  Teuchos::Array<GlobalOrdinal> cols;
 for(LocalOrdinal l_row = 0; (size_t) l_row<rowmap->getNodeNumElements(); l_row++) {
    GlobalOrdinal g_row = rowmap->getGlobalElement(l_row);
    if (g_row == rowmap->getMinAllGlobalIndex() ||
        g_row == rowmap->getMaxAllGlobalIndex()) cols.resize(2);
    else cols.resize(3);

    size_t coloffset = 0;
    if (g_row > rowmap->getMinAllGlobalIndex()) cols[coloffset++] = g_row-1;
    cols[coloffset++] = g_row;
    if (g_row < rowmap->getMaxAllGlobalIndex()) cols[coloffset++] = g_row+1;

    crsgraph->insertGlobalIndices(g_row, cols());
  }

  crsgraph->fillComplete();

  return crsgraph;
}

void
Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
getQuadPoints(ordinal_type quad_order,
	      Teuchos::Array<value_type>& quad_points,
	      Teuchos::Array<value_type>& quad_weights,
	      Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
{
  // This assumes Gaussian quadrature
  ordinal_type num_points = 
    static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
  const Pecos::RealArray& gp = pecosPoly->collocation_points(num_points);
  const Pecos::RealArray& gw = pecosPoly->type1_collocation_weights(num_points); 
  quad_points.resize(num_points);
  quad_weights.resize(num_points);
  for (ordinal_type i=0; i<num_points; i++) {
    quad_points[i] = gp[i];
    quad_weights[i] = gw[i];
  }
  
  // Evalute basis at gauss points
  quad_values.resize(num_points);
  for (ordinal_type i=0; i<num_points; i++) {
    quad_values[i].resize(p+1);
    evaluateBases(quad_points[i], quad_values[i]);
  }
}

 MeanDeviationFromTarget( Teuchos::ParameterList &parlist )
   : RiskMeasure<Real>(), firstReset_(true) {
   Teuchos::ParameterList &list
     = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation From Target");
   // Get data from parameter list
   Teuchos::Array<Real> target
     = Teuchos::getArrayFromStringParameter<double>(list,"Targets");
   Teuchos::Array<Real> order
     = Teuchos::getArrayFromStringParameter<double>(list,"Orders");
   Teuchos::Array<Real> coeff
     = Teuchos::getArrayFromStringParameter<double>(list,"Coefficients");
   // Check inputs
   NumMoments_ = order.size();
   target_.clear(); order_.clear(); coeff_.clear();
   if ( NumMoments_ != static_cast<uint>(target.size()) ) {
     target.resize(NumMoments_,0.0);
   }
   if ( NumMoments_ != static_cast<uint>(coeff.size()) ) {
     coeff.resize(NumMoments_,1.0);
   }
   for ( uint i = 0; i < NumMoments_; i++ ) {
     target_.push_back(target[i]);
     order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
     coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
   }
   initialize();
   // Build (approximate) positive function
   if ( list.get("Deviation Type","Upper") == "Upper" ) {
     positiveFunction_ = Teuchos::rcp(new PlusFunction<Real>(list));
   }
   else {
     positiveFunction_ = Teuchos::rcp(new AbsoluteValue<Real>(list));
   }
 }

Teuchos::Array<bool> createResponseTable(
    int count,
    const std::string selectionType,
    int index,
    const Teuchos::ArrayView<const int> &list)
{
  Teuchos::Array<bool> result;

  if (count > 0) {
    if (selectionType == "All") {
      result.resize(count, true);
    } else if (selectionType == "Last") {
      result = createResponseTableFromIndex(count - 1, count);
    } else if (selectionType == "AllButLast") {
      result.reserve(count);
      result.resize(count - 1, true);
      result.push_back(false);
    } else if (selectionType == "Index") {
      result = createResponseTableFromIndex(index, count);
    } else if (selectionType == "List") {
      result.resize(count, false);
      for (Teuchos::ArrayView<const int>::const_iterator it = list.begin(), it_end = list.end(); it != it_end; ++it) {
        result.at(*it) = true;
      }
    } else {
      TEUCHOS_TEST_FOR_EXCEPT(false);
    }
  }

  return result;
}

void BCManager::modifyPrimalSystem(
    std::string const& val, int offset, std::string const& set)
{
  // should we fill in BC info?
  bool fillRes = (res != Teuchos::null);
  bool fillJac = (jac != Teuchos::null);
  if ((!fillRes) && (!fillJac)) return;

  // get views of the solution and residual vectors
  ArrayRCP<const ST> x_const_view;
  ArrayRCP<ST> f_nonconst_view;
  if (fillRes) {
    x_const_view = sol->get1dView();
    f_nonconst_view = res->get1dViewNonConst();
  }

  // set up some data for replacing jacobian values
  Teuchos::Array<LO> index(1);
  Teuchos::Array<ST> value(1);
  value[0] = 1.0;
  size_t numEntries;
  Teuchos::Array<ST> matrixEntries;
  Teuchos::Array<LO> matrixIndices;

  // loop over all of the nodes in this node set
  std::vector<Albany::GOALNode> nodes = ns[set];
  for (int i=0; i < nodes.size(); ++i)
  {
    Albany::GOALNode node = nodes[i];
    int lunk = disc->getDOF(node.lid, offset);

    // if the node is higher order, we set the value of the DBC to be 0
    // if it is not, we parse the expression from the input file to get the value
    // note: this assumes that bcs are either constant or linear in space
    // anything else would require a linear solve to find coefficients v
    double v = 0.0;
    if (!node.higherOrder)
      v = parseExpression(val, node.coord[0], node.coord[1], node.coord[2], t);

    // modify the residual if necessary
    if (fillRes)
      f_nonconst_view[lunk] = x_const_view[lunk] - v;

    // modify the jacobian if necessary
    if (fillJac)
    {
      index[0] = lunk;
      numEntries = jac->getNumEntriesInLocalRow(lunk);
      matrixEntries.resize(numEntries);
      matrixIndices.resize(numEntries);
      jac->getLocalRowCopy(lunk, matrixIndices(), matrixEntries(), numEntries);
      for (int i=0; i < numEntries; ++i)
        matrixEntries[i] = 0.0;
      jac->replaceLocalValues(lunk, matrixIndices(), matrixEntries());
      jac->replaceLocalValues(lunk, index(), value());
    }
  }
}

void DirichletField<PHAL::AlbanyTraits::Jacobian, Traits>::
evaluateFields(typename Traits::EvalData dirichletWorkset) {

#if defined(ALBANY_EPETRA)
  const Albany::NodalDOFManager& fieldDofManager = dirichletWorkset.disc->getDOFManager(this->field_name);
  Teuchos::RCP<const Epetra_Map> fiedNodeMap = dirichletWorkset.disc->getNodeMap(this->field_name);
  const std::vector<GO>& nsNodesGIDs = dirichletWorkset.disc->getNodeSetGIDs().find(this->nodeSetID)->second;
#endif
  Teuchos::RCP<Tpetra_Vector> pvecT = dirichletWorkset.distParamLib->get(this->field_name)->vector();
  Teuchos::ArrayRCP<const ST> pT = pvecT->get1dView();
  Teuchos::RCP<Tpetra_Vector> fT = dirichletWorkset.fT;
  Teuchos::RCP<const Tpetra_Vector> xT = dirichletWorkset.xT;
  Teuchos::ArrayRCP<const ST> xT_constView = xT->get1dView();
  Teuchos::RCP<Tpetra_CrsMatrix> jacT = dirichletWorkset.JacT;

  const RealType j_coeff = dirichletWorkset.j_coeff;
  const std::vector<std::vector<int> >& nsNodes = dirichletWorkset.nodeSets->find(this->nodeSetID)->second;

  bool fillResid = (fT != Teuchos::null);
  Teuchos::ArrayRCP<ST> fT_nonconstView;                                         
  if (fillResid) fT_nonconstView = fT->get1dViewNonConst();

  Teuchos::Array<LO> index(1);
  Teuchos::Array<ST> value(1); 
  size_t numEntriesT;  
  value[0] = j_coeff; 
  Teuchos::Array<ST> matrixEntriesT; 
  Teuchos::Array<LO> matrixIndicesT; 

  for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
      int lunk = nsNodes[inode][this->offset];
      index[0] = lunk; 
      numEntriesT = jacT->getNumEntriesInLocalRow(lunk);
      matrixEntriesT.resize(numEntriesT); 
      matrixIndicesT.resize(numEntriesT); 

      jacT->getLocalRowCopy(lunk, matrixIndicesT(), matrixEntriesT(), numEntriesT); 
      for (int i=0; i<numEntriesT; i++) matrixEntriesT[i]=0;
      jacT->replaceLocalValues(lunk, matrixIndicesT(), matrixEntriesT()); 

      jacT->replaceLocalValues(lunk, index(), value()); 
      
      if (fillResid) {
#if defined(ALBANY_EPETRA)
        GO node_gid = nsNodesGIDs[inode];
        int lfield = fieldDofManager.getLocalDOF(fiedNodeMap->LID(node_gid),this->offset);
        fT_nonconstView[lunk] = xT_constView[lunk] - pT[lfield];
#else
        fT_nonconstView[lunk] = xT_constView[lunk] - pT[lunk];
#endif
      }
  }
}

void EquilibriumConcentrationBC<PHAL::AlbanyTraits::Jacobian, Traits>::
evaluateFields(typename Traits::EvalData dirichletWorkset)
{
  Teuchos::RCP<Tpetra_Vector> fT = dirichletWorkset.fT;
  Teuchos::RCP<const Tpetra_Vector> xT = dirichletWorkset.xT;
  Teuchos::ArrayRCP<const ST> xT_constView = xT->get1dView();
  Teuchos::RCP<Tpetra_CrsMatrix> jacT = dirichletWorkset.JacT;

  const RealType j_coeff = dirichletWorkset.j_coeff;
  const std::vector<std::vector<int>>& nsNodes =
    dirichletWorkset.nodeSets->find(this->nodeSetID)->second;

  bool fillResid = (fT != Teuchos::null);
  Teuchos::ArrayRCP<ST> fT_nonconstView;
  if (fillResid) fT_nonconstView = fT->get1dViewNonConst();

  int cunk, punk;
  ScalarT Cval;
  ScalarT pressure;

  Teuchos::Array<LO> index(1);
  Teuchos::Array<ST> value(1);
  size_t numEntriesT;
  value[0] = j_coeff;
  Teuchos::Array<ST> matrixEntriesT;
  Teuchos::Array<LO> matrixIndicesT;


  for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
  {
    cunk = nsNodes[inode][this->coffset_];
    punk = nsNodes[inode][this->poffset_];
    pressure = xT_constView[punk];
    this->computeBCs(pressure, Cval);
    
    // replace jac values for the C dof 
    numEntriesT = jacT->getNumEntriesInLocalRow(cunk);
    matrixEntriesT.resize(numEntriesT);
    matrixIndicesT.resize(numEntriesT);
    jacT->getLocalRowCopy(cunk, matrixIndicesT(), matrixEntriesT(), numEntriesT);
    for (int i=0; i<numEntriesT; i++) matrixEntriesT[i]=0;
    jacT->replaceLocalValues(cunk, matrixIndicesT(), matrixEntriesT());

    index[0] = cunk;
    jacT->replaceLocalValues(cunk, index(), value());

    if (fillResid)
    {
      fT_nonconstView[cunk] = xT_constView[cunk] - Cval.val();
    } 
  }
}

/* Return the kth row  */
template <class Scalar> inline
void MultiVectorOperator<Scalar>
::getRow(const int& k, 
  Teuchos::Array<int>& indices, 
  Teuchos::Array<Scalar>& values) const
{
  int low = this->range()->baseGlobalNaturalIndex();
  indices.resize(cols_.size());
  values.resize(cols_.size());
  for (int j=0; j<cols_.size(); j++)
  {
    indices[j] = j;
    values[j] = cols_[j][k-low];
  }
}

void
Stokhos::GaussPattersonLegendreBasis<ordinal_type,value_type>::
getQuadPoints(ordinal_type quad_order,
	      Teuchos::Array<value_type>& quad_points,
	      Teuchos::Array<value_type>& quad_weights,
	      Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
{
#ifdef HAVE_STOKHOS_DAKOTA
  // Gauss-Patterson points have the following structure 
  // (cf. http://people.sc.fsu.edu/~jburkardt/f_src/patterson_rule/patterson_rule.html):
  //  Level l  Num points n   Precision p
  //  -----------------------------------
  //    0        1                1
  //    1        3                5
  //    2        7                11
  //    3        15               23
  //    4        31               47
  //    5        63               95
  //    6        127              191
  //    7        255              383
  // Thus for l > 0, n = 2^{l+1}-1 and p = 3*2^l-1.  So for a given quadrature
  // order p, we find the smallest l s.t. 3*s^l-1 >= p and then compute the
  // number of points n from the above.  In this case, l = ceil(log2((p+1)/3))
  ordinal_type num_points;
  if (quad_order <= ordinal_type(1))
    num_points = 1;
  else {
    ordinal_type l = std::ceil(std::log((quad_order+1.0)/3.0)/std::log(2.0));
    num_points = (1 << (l+1)) - 1; // std::pow(2,l+1)-1;
  }
  
  quad_points.resize(num_points);
  quad_weights.resize(num_points);
  quad_values.resize(num_points);

  webbur::patterson_lookup(num_points, &quad_points[0], &quad_weights[0]);

  for (ordinal_type i=0; i<num_points; i++) {
    quad_weights[i] *= 0.5;  // scale to unit measure
    quad_values[i].resize(this->p+1);
    this->evaluateBases(quad_points[i], quad_values[i]);
  }

#else
  TEUCHOS_TEST_FOR_EXCEPTION(
    true, std::logic_error, "Clenshaw-Curtis requires TriKota to be enabled!");
#endif
}

void
getRemotePIDs (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
               Teuchos::Array<int>& RemotePIDs)
{
  const Tpetra::Distributor& D = Importer.getDistributor();

  // Get the importer's data
  Teuchos::ArrayView<const LocalOrdinal> RemoteLIDs  = Importer.getRemoteLIDs();

  // Get the distributor's data
  size_t NumReceives                           = D.getNumReceives();
  Teuchos::ArrayView<const int> ProcsFrom      = D.getProcsFrom();
  Teuchos::ArrayView<const size_t> LengthsFrom = D.getLengthsFrom();

  // Resize the outgoing data structure
  RemotePIDs.resize(Importer.getNumRemoteIDs());

  // Now, for each remote ID, record who actually owns it.  This loop
  // follows the operation order in the MpiDistributor so it ought to
  // duplicate that effect.
  size_t i,j,k;
  for (i = 0, j = 0; i < NumReceives; ++i) {
    const int pid = ProcsFrom[i];
    for (k = 0; k < LengthsFrom[i]; ++k) {
      RemotePIDs[j] = pid;
      j++;
    }
  }
}

void SharedDomainMap<Mesh,CoordinateField>::computePointOrdinals(
    const RCP_CoordFieldManager& target_coord_manager,
    Teuchos::Array<GlobalOrdinal>& target_ordinals )
{
    // Set an existence value for the target coords.
    bool target_exists = true;
    if ( target_coord_manager.is_null() ) target_exists = false;
    int comm_rank = d_comm->getRank();
    GlobalOrdinal local_size = 0;

    if ( target_exists )
    {
	int point_dim = CFT::dim( *target_coord_manager->field() );
	local_size = std::distance( 
	    CFT::begin( *target_coord_manager->field() ),
	    CFT::end( *target_coord_manager->field() ) ) / point_dim;
    }

    GlobalOrdinal global_max;
    Teuchos::reduceAll<int,GlobalOrdinal>( *d_comm,
					   Teuchos::REDUCE_MAX,
					   1,
					   &local_size,
					   &global_max );

    target_ordinals.resize( local_size );
    for ( GlobalOrdinal n = 0; n < local_size; ++n )
    {
	target_ordinals[n] = comm_rank*global_max + n;
    }
}

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

/*!
 * \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 FELIX::Hydrology::constructNeumannEvaluators (const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
{
    // Note: we only enter this function if sidesets are defined in the mesh file
    // i.e. meshSpecs.ssNames.size() > 0

    Albany::BCUtils<Albany::NeumannTraits> nbcUtils;

    // Check to make sure that Neumann BCs are given in the input file
    if (!nbcUtils.haveBCSpecified(this->params))
    {
        return;
    }

    // Construct BC evaluators for all side sets and names
    // Note that the string index sets up the equation offset, so ordering is important

    std::vector<std::string> neumannNames(neq + 1);
    Teuchos::Array<Teuchos::Array<int> > offsets;
    offsets.resize(1);

    neumannNames[0] = "Hydraulic Potential";
    neumannNames[1] = "all";
    offsets[0].resize(1);
    offsets[0][0] = 0;

    // Construct BC evaluators for all possible names of conditions
    std::vector<std::string> condNames(1);
    condNames[0] = "neumann";

    nfm.resize(1); // FELIX problem only has one element block

    nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, neumannNames, dof_names, false, 0,
                                            condNames, offsets, dl,
                                            this->params, this->paramLib);
}

//---------------------------------------------------------------------------//
// Get the parameteric center of an entity.
void MoabEntityLocalMap::parametricCenter(
    const Entity& entity,
    Teuchos::Array<double>& center ) const
{
    moab::EntityType moab_type = d_moab_mesh->get_moab()->type_from_handle(
	MoabHelpers::extractEntity(entity) );

    switch( moab_type )
    {
	case moab::MBTRI:
	    center.assign( 2, 1.0 / 3.0 );
	    break;
	case moab::MBQUAD:
	    center.assign( 2, 0.0 );
	    break;
	case moab::MBTET:
	    center.assign( 3, 1.0 / 6.0 );
	    break;
	case moab::MBHEX:
	    center.assign( 3, 0.0 );
	    break;
	default:
	    center.resize( 0 );
	    break;
    }
}