C++ RCP::Map方法代码示例

NOX::Epetra::Vector::
Vector(const Teuchos::RCP<Epetra_Vector>& source,
       NOX::Epetra::Vector::MemoryType memoryType,
       NOX::CopyType type,
       Teuchos::RCP<NOX::Epetra::VectorSpace> vs)
{
  if (Teuchos::is_null(vs))
    vectorSpace = Teuchos::rcp(new NOX::Epetra::VectorSpaceL2);
  else
    vectorSpace = vs;

  if (memoryType == NOX::Epetra::Vector::CreateView)
    epetraVec = source;
  else {

    switch (type) {

    case DeepCopy:        // default behavior

      epetraVec = Teuchos::rcp(new Epetra_Vector(*source));
      break;

    case ShapeCopy:

      epetraVec = Teuchos::rcp(new Epetra_Vector(source->Map()));
      break;
    }

  }
}

/// building the stencil
void PoissonSolver::StencilGeometry(Teuchos::RCP<Epetra_Vector> & RHS,
                                    Teuchos::RCP<Epetra_CrsMatrix> & A)
{
    const Epetra_BlockMap & MyMap = RHS->Map();
    int NumMyElements = MyMap.NumMyElements();
    int* MyGlobalElements = MyMap.MyGlobalElements();
    double * rhsvalues = RHS->Values();

    std::vector<double> Values(5);
    std::vector<int> Indices(5);

    const auto & inside = _bend.getInsideMask();

    for (int lid = 0; lid < NumMyElements; ++ lid) {
        size_t NumEntries = 0;

        const size_t & gid = MyGlobalElements[lid];

        cutoffStencil(Indices,
                      Values,
                      rhsvalues[lid],
                      NumEntries,
                      inside,
                      gid);
        A->InsertGlobalValues(gid, NumEntries, &Values[0], &Indices[0]);
    }

    A->FillComplete();
    A->OptimizeStorage();
}

void PoissonSolver::plotPotential(BinaryVtkFile & vtkFile,
                                  const Teuchos::RCP<Epetra_Vector> & LHS,
                                  const VField_Edge_t & EFD,
                                  const int & component)
{
    const Epetra_BlockMap & Map = LHS->Map();
    const int * MyGlobalElements = Map.MyGlobalElements();
    const int NumMyElements = Map.NumMyElements();

    boost::shared_ptr<SField_Vert_t> scalarField = Utils::getScalarVertField(EFD);
    NDIndex<DIM> elem;
    char fieldName[50];
    ST * values = LHS->Values();

    sprintf(fieldName, "potential_%d", component);

    for (int lid = 0; lid < NumMyElements; ++ lid) {
        const size_t idx = MyGlobalElements[lid] % _Nx;
        const size_t idy = MyGlobalElements[lid] / _Nx;
        elem[0] = Index(idx, idx);
        elem[1] = Index(idy, idy);

        scalarField->localElement(elem) = values[lid];
    }

    vtkFile.addScalarField(*scalarField, fieldName);
}

 void ISVDUDV::Initialize( const Teuchos::RCP< Teuchos::ParameterList >& params,
                           const Teuchos::RCP< const Epetra_MultiVector >& ss,
                           const Teuchos::RCP< RBGen::FileIOHandler< Epetra_Operator > >& fileio) 
 {
   workU_ = Teuchos::rcp( new Epetra_MultiVector(ss->Map(),maxBasisSize_,false) );
   Epetra_LocalMap lclmap(ss->NumVectors(),0,ss->Comm());
   workV_ = Teuchos::rcp( new Epetra_MultiVector(lclmap,maxBasisSize_,false) );
 }

void panzer::ScatterDirichletResidual_Epetra<panzer::Traits::Residual, Traits,LO,GO>::
evaluateFields(typename Traits::EvalData workset)
{ 
   std::vector<GO> GIDs;
   std::vector<int> LIDs;
 
   // for convenience pull out some objects from workset
   std::string blockId = workset.block_id;
   const std::vector<std::size_t> & localCellIds = workset.cell_local_ids;

   Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
   Teuchos::RCP<Epetra_Vector> r = epetraContainer->get_f(); 

   // NOTE: A reordering of these loops will likely improve performance
   //       The "getGIDFieldOffsets may be expensive.  However the
   //       "getElementGIDs" can be cheaper. However the lookup for LIDs
   //       may be more expensive!

   // scatter operation for each cell in workset
   for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
      std::size_t cellLocalId = localCellIds[worksetCellIndex];

      globalIndexer_->getElementGIDs(cellLocalId,GIDs); 

      // caculate the local IDs for this element
      LIDs.resize(GIDs.size());
      for(std::size_t i=0;i<GIDs.size();i++)
         LIDs[i] = r->Map().LID(GIDs[i]);

      // loop over each field to be scattered
      for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
         int fieldNum = fieldIds_[fieldIndex];
   
         // this call "should" get the right ordering accordint to the Intrepid basis
         const std::pair<std::vector<int>,std::vector<int> > & indicePair 
               = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
         const std::vector<int> & elmtOffset = indicePair.first;
         const std::vector<int> & basisIdMap = indicePair.second;
   
         // loop over basis functions
         for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
            int offset = elmtOffset[basis];
            int lid = LIDs[offset];
            if(lid<0) // not on this processor!
               continue;

            int basisId = basisIdMap[basis];
            (*r)[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);

            // record that you set a dirichlet condition
            if(dirichletCounter_!=Teuchos::null)
              (*dirichletCounter_)[lid] = 1.0;
         }
      }
   }
}

void PeridigmNS::BoundaryCondition::evaluateParser(const int & localNodeID, double & currentValue, double & previousValue, const double & timeCurrent, const double & timePrevious){
  Teuchos::RCP<Epetra_Vector> x = peridigm->getX();
  const Epetra_BlockMap& threeDimensionalMap = x->Map();
  TEUCHOS_TEST_FOR_EXCEPT_MSG(threeDimensionalMap.ElementSize() != 3, "**** setVectorValues() must be called with map having element size = 3.\n");

  bool success(true);
  // set the coordinates and set the return value to 0.0
  if(success)
    rtcFunction->varValueFill(0, (*x)[localNodeID*3]);
  if(success)
    success = rtcFunction->varValueFill(1, (*x)[localNodeID*3 + 1]);
  if(success)
    success = rtcFunction->varValueFill(2, (*x)[localNodeID*3 + 2]);
  if(success)
    success = rtcFunction->varValueFill(4, 0.0);
  // evaluate at previous time
  if(success)
    success = rtcFunction->varValueFill(3, timePrevious);
  if(success)
    success = rtcFunction->execute();
  if(success)
    previousValue = rtcFunction->getValueOfVar("value");
  // evaluate at current time
  if(success)
    success = rtcFunction->varValueFill(3, timeCurrent);
  if(success)
    success = rtcFunction->execute();
  if(success)
    currentValue = rtcFunction->getValueOfVar("value");
  if(!success){
    string msg = "\n**** Error in BoundaryCondition::evaluateParser().\n";
    msg += "**** " + rtcFunction->getErrors() + "\n";
    TEUCHOS_TEST_FOR_EXCEPT_MSG(!success, msg);
  }

  // if this is any other boundary condition besides prescribed displacement
  // get the previous value from evaluating the string function as above
  // if it is a perscribed displacement bc, check to see if the increment is
  // zero. This could happen if the user specifies a constant prescribed displacement.
  // If so, the increment should be the current prescribed value
  // minus the existing field value instead of the parser evaluation
  if(bcType==PRESCRIBED_DISPLACEMENT && currentValue - previousValue == 0.0)
  {
    Teuchos::RCP<Epetra_Vector> previousDisplacement = peridigm->getU();
    previousValue = (*previousDisplacement)[localNodeID*3 + coord];
  }

  // TODO: we should revisit how prescribed boundary conditions interact with initial conditions
  // in the case that the prescribed boundary condition at time zero doesn't match the inital condition
  // who wins?
}

Library::
Library(Teuchos::RCP<const Epetra_MultiVector> input_coords, int itype)
  : input_type_(itype),
    numPartSizes(0),
    partGIDs(NULL),
    partSizes(NULL),
    input_graph_(0),
    input_matrix_(0),
    input_coords_(input_coords),
    costs_(0),
    weights_(0)
{
  input_map_ = Teuchos::rcp(&(input_coords->Map()), false);
}

Teuchos::RCP<const Epetra_Vector>
EpetraExt::createInverseModelScalingVector(
  Teuchos::RCP<const Epetra_Vector> const& scalingVector
  )
{
  Teuchos::RCP<Epetra_Vector> invScalingVector =
    Teuchos::rcpWithEmbeddedObj(
      new Epetra_Vector(scalingVector->Map()),
      scalingVector
      );
  invScalingVector->Reciprocal(*scalingVector);
  return invScalingVector;
  // Above, we embedd the forward scaling vector.  This is done in order to
  // achieve the exact same numerics as before this refactoring and to improve
  // runtime speed and accruacy.
}

void panzer::ScatterInitialCondition_Epetra<panzer::Traits::Tangent, Traits,LO,GO>::
evaluateFields(typename Traits::EvalData workset)
{ 
   TEUCHOS_ASSERT(false);

   std::vector<GO> GIDs;
   std::vector<int> LIDs;
 
   // for convenience pull out some objects from workset
   std::string blockId = workset.block_id;
   const std::vector<std::size_t> & localCellIds = workset.cell_local_ids;

   Teuchos::RCP<Epetra_Vector> x = epetraContainer_->get_x(); 

   // NOTE: A reordering of these loops will likely improve performance
   //       The "getGIDFieldOffsets may be expensive.  However the
   //       "getElementGIDs" can be cheaper. However the lookup for LIDs
   //       may be more expensive!

   // scatter operation for each cell in workset
   for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
      std::size_t cellLocalId = localCellIds[worksetCellIndex];

      globalIndexer_->getElementGIDs(cellLocalId,GIDs); 

      // caculate the local IDs for this element
      LIDs.resize(GIDs.size());
      for(std::size_t i=0;i<GIDs.size();i++)
         LIDs[i] = x->Map().LID(GIDs[i]);

      // loop over each field to be scattered
      for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
         int fieldNum = fieldIds_[fieldIndex];
         const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
   
         // loop over basis functions
         for(std::size_t basis=0;basis<elmtOffset.size();basis++) {

            int offset = elmtOffset[basis];
            int lid = LIDs[offset];
	    if (lid != -1)
	      (*x)[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
         }
      }
   }
}

static int run_test(Teuchos::RCP<const Epetra_MultiVector> &coords,
		    const Teuchos::ParameterList &params)
{
  Teuchos::RCP<Isorropia::Epetra::Partitioner> part =
    Teuchos::rcp(new Isorropia::Epetra::Partitioner(coords, params));

  // both Zoltan and Isorropia partitioners will use unit weights,
  // so we create unit weights to compute balances

  Epetra_MultiVector tmp(coords->Map(), 1);
  tmp.PutScalar(1.0);

  Teuchos::RCP<const Epetra_MultiVector> unitweights =
    Teuchos::rcp(new const Epetra_MultiVector(tmp));

  return check_test(coords, unitweights, part);
}

void SeparableScatterScalarResponse<PHAL::AlbanyTraits::DistParamDeriv, Traits>::
postEvaluate(typename Traits::PostEvalData workset)
{
#if defined(ALBANY_EPETRA)
  // Here we scatter the *global* response and its derivatives
  Teuchos::RCP<Epetra_Vector> g = workset.g;
  Teuchos::RCP<Epetra_MultiVector> dgdp = workset.dgdp;
  Teuchos::RCP<Epetra_MultiVector> overlapped_dgdp = workset.overlapped_dgdp;
  if (g != Teuchos::null)
     for (std::size_t res = 0; res < this->global_response.size(); res++) {
       (*g)[res] = this->global_response[res].val();
   }
  if (dgdp != Teuchos::null) {
    Epetra_Export exporter(overlapped_dgdp->Map(), dgdp->Map());
    dgdp->Export(*overlapped_dgdp, exporter, Add);
  }
#endif
}

void PoissonSolver::initializeLHS(Teuchos::RCP<Epetra_Vector> & LHS) const
{
    const Epetra_BlockMap & Map = LHS->Map();
    const int * MyGlobalElements = Map.MyGlobalElements();
    const int NumMyElements = Map.NumMyElements();

    ST * values = LHS->Values();

    // size_t length = _bend.getLengthStraightSection() + 1;
    // size_t width = _bend.getWidthStraightSection() + 1;
    Vector_t position;
    _bunch.get_rmean(position);
    _bend.getPositionInCells(position);

    int lastX = (_bend.getGlobalDomain())[0].last();
    int lastY = (_bend.getGlobalDomain())[1].last();

    // NDIndex<DIM> initDom(Index(0,length), Index(lastY - width, lastY));

    NDIndex<DIM> elem;
    for (int lid = 0; lid < NumMyElements; ++ lid) {
        const size_t idx = MyGlobalElements[lid] % _Nx;
        const size_t idy = MyGlobalElements[lid] / _Nx;
        elem[0] = Index(idx, idx);
        elem[1] = Index(idy, idy);

        Vector_t contr;
        if (idx > position(0)) {
            contr(0) = (idx - position(0)) / (lastX - position(0));
        } else {
            contr(0) = (position(0) - idx) / position(0);
        }
        if (idy > position(1)) {
            contr(1) = (idy - position(1)) / (lastY - position(1));
        } else {
            contr(1) = (position(1) - idy) / position(1);
        }

        const double val = sqrt(dot(contr, contr));
        dbg << std::min(10.0, -log(val)) << std::endl;
        values[lid] = -std::max(0.0, std::min(10.0,-log(val)));
    }
}

Teuchos::RCP<const Epetra_MultiVector>
ReducedSpaceFactory::getProjector(const Teuchos::RCP<Teuchos::ParameterList> &params)
{
  const Teuchos::RCP<const Epetra_MultiVector> basis = this->getBasis(params);

  const Teuchos::RCP<const Epetra_Map> basisMap = mapDowncast(basis->Map());
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(basisMap));

  const Teuchos::RCP<const Epetra_Operator> collocationOperator =
    this->getSamplingOperator(params, *basisMap);

  Teuchos::RCP<const Epetra_MultiVector> result = basis;
  if (Teuchos::nonnull(collocationOperator)) {
    const Teuchos::RCP<Epetra_MultiVector> dualBasis(
        new Epetra_MultiVector(collocationOperator->OperatorRangeMap(), basis->NumVectors(), false));
    dualize(*basis, *collocationOperator, *dualBasis);
    result = dualBasis;
  }
  return result;
}

void PoissonSolver::shiftLHS(SField_t & rho,
                             Teuchos::RCP<Epetra_Vector> & LHS,
                             double tau) const
{
    const Epetra_BlockMap & Map = LHS->Map();
    const int * MyGlobalElements = Map.MyGlobalElements();
    const int NumMyElements = Map.NumMyElements();

    NDIndex<DIM> elem;
    NDIndex<DIM> ldom = rho.getLayout().getLocalNDIndex();
    Index II = ldom[0], JJ = ldom[1];
    ST * values = LHS->Values();

    for (int lid = 0; lid < NumMyElements; ++ lid) {
        const size_t idx = MyGlobalElements[lid] % _Nx;
        const size_t idy = MyGlobalElements[lid] / _Nx;
        elem[0] = Index(idx, idx);
        elem[1] = Index(idy, idy);

        rho.localElement(elem) = values[lid];
    }
    rho.fillGuardCells();

    if (tau > 0) {
        if (tau > 1) tau = 1;
        rho[II][JJ] += tau * rho[II+1][JJ] - tau * rho[II][JJ];
    } else {
        tau = - tau;
        if (tau > 1) tau = 1;
        rho[II][JJ] += tau * rho[II-1][JJ] - tau * rho[II][JJ];
    }

    for (int lid = 0; lid < NumMyElements; ++ lid) {
        const int idx = MyGlobalElements[lid] % _Nx;
        const int idy = MyGlobalElements[lid] / _Nx;
        elem[0] = Index(idx, idx);
        elem[1] = Index(idy, idy);

        values[lid] = rho.localElement(elem);
    }
}

void PoissonSolver::add(const SField_t & rho,
                        Teuchos::RCP<Epetra_Vector> & LHS) const
{
    const Epetra_BlockMap & Map = LHS->Map();
    const int * MyGlobalElements = Map.MyGlobalElements();
    const int NumMyElements = Map.NumMyElements();

    NDIndex<DIM> elem;
    NDIndex<DIM> ldom = rho.getLayout().getLocalNDIndex();
    Index II = ldom[0], JJ = ldom[1];
    ST * values = LHS->Values();

    for (int lid = 0; lid < NumMyElements; ++ lid) {
        const size_t idx = MyGlobalElements[lid] % _Nx;
        const size_t idy = MyGlobalElements[lid] / _Nx;
        elem[0] = Index(idx, idx);
        elem[1] = Index(idy, idy);

        values[lid] += rho.localElement(elem);
    }
}