C++ Tpetra_Vector::get1dView方法代码示例

void Albany::GenericSTKFieldContainer<Interleaved>::saveVectorHelperT(const Tpetra_Vector& solnT,
    ScalarFieldType* solution_field,
    const Teuchos::RCP<const Tpetra_Map>& node_mapT,
    const stk::mesh::Bucket& bucket, int offset) {

  // Fill the result vector
  // Create a multidimensional array view of the
  // solution field data for this bucket of nodes.
  // The array is two dimensional ( Cartesian X NumberNodes )
  // and indexed by ( 0..2 , 0..NumberNodes-1 )

  BucketArray<ScalarFieldType> solution_array(*solution_field, bucket);

  const int num_nodes_in_bucket = solution_array.dimension(0);

  stk::mesh::BulkData& mesh = solution_field->get_mesh();

  //get const (read-only) view of solnT
  Teuchos::ArrayRCP<const ST> solnT_constView = solnT.get1dView();

  for(std::size_t i = 0; i < num_nodes_in_bucket; i++)  {

    const GO node_gid = mesh.identifier(bucket[i]) - 1;

    if(node_mapT->getLocalElement(node_gid) != Teuchos::OrdinalTraits<LO>::invalid()){
      int node_lid = node_mapT->getLocalElement(node_gid);
      solution_array(i) = solnT_constView[getDOF(node_lid, offset)];
    }
  }
}

void
Albany::SolutionMaxValueResponseFunction::
evaluateGradientT(const double current_time,
		 const Tpetra_Vector* xdotT,
		 const Tpetra_Vector* xdotdotT,
		 const Tpetra_Vector& xT,
		 const Teuchos::Array<ParamVec>& p,
		 ParamVec* deriv_p,
		 Tpetra_Vector* gT,
		 Tpetra_MultiVector* dg_dxT,
		 Tpetra_MultiVector* dg_dxdotT,
		 Tpetra_MultiVector* dg_dxdotdotT,
		 Tpetra_MultiVector* dg_dpT)
{
  int global_index;
  double mxv;
  computeMaxValueT(xT, mxv, global_index);
  
  // Evaluate response g
  if (gT != NULL) {
    Teuchos::ArrayRCP<ST> gT_nonconstView = gT->get1dViewNonConst();
    gT_nonconstView[0] = mxv;
  }

  Teuchos::ArrayRCP<const ST> xT_constView = xT.get1dView();
  
  // Evaluate dg/dx
  if (dg_dxT != NULL) {
    Teuchos::ArrayRCP<ST> dg_dxT_nonconstView;
    int im = -1;
    for (int i=0; i<xT.getMap()->getNodeNumElements(); i++) {
       dg_dxT_nonconstView = dg_dxT->getDataNonConst(0); 
       if (xT_constView[i] == mxv) { dg_dxT_nonconstView[i] = 1.0; im = i; }
       else                          dg_dxT_nonconstView[i] = 0.0;
    }

  }

  // Evaluate dg/dxdot
  if (dg_dxdotT != NULL)
    dg_dxdotT->putScalar(0.0);
  if (dg_dxdotdotT != NULL)
    dg_dxdotdotT->putScalar(0.0);

  // Evaluate dg/dp
  if (dg_dpT != NULL)
    dg_dpT->putScalar(0.0);

}

void
Albany::SolutionMaxValueResponseFunction::
computeMaxValueT(const Tpetra_Vector& xT, double& global_max, int& global_index)
{
  //The following is needed b/c Epetra_MaxDouble comes from Trilinos Epetra package.
  double Tpetra_MaxDouble = 1.0E+100; 
  double my_max = -Tpetra_MaxDouble;
  int my_index = -1, index;
  
  Teuchos::ArrayRCP<const ST> xT_constView = xT.get1dView();
  
  // Loop over nodes to find max value for equation eq
  int num_my_nodes = xT.getLocalLength() / neq;
  for (int node=0; node<num_my_nodes; node++) {
    if (interleavedOrdering)  index = node*neq+eq;
    else                      index = node + eq*num_my_nodes;
    if (xT_constView[index] > my_max) {
      my_max = xT_constView[index];
      my_index = index;
    }
  }

  // Check remainder (AGS: NOT SURE HOW THIS CODE GETS CALLED?)
  if (num_my_nodes*neq+eq < xT.getLocalLength()) {
    if (interleavedOrdering)  index = num_my_nodes*neq+eq;
    else                      index = num_my_nodes + eq*num_my_nodes;
    if (xT_constView[index] > my_max) {
      my_max = xT_constView[index];
      my_index = index;
    }
  }

  Teuchos::RCP<const Teuchos::Comm<int> > commT = xT.getMap()->getComm(); 
  // Get max value across all proc's
  Teuchos::reduceAll(*commT, Teuchos::REDUCE_MAX, my_max, Teuchos::ptr(&global_max)); 

  // Compute min of all global indices equal to max value
  if (my_max == global_max)
    my_index = xT.getMap()->getGlobalElement(my_index);
  else
    my_index = xT.getGlobalLength();
  Teuchos::reduceAll(*commT, Teuchos::REDUCE_MIN, my_index, Teuchos::ptr(&global_index)); 
}

void
twoD_diffusion_problem<Scalar,MeshScalar,BasisScalar,LocalOrdinal,GlobalOrdinal,
		       Node>::
computeA(const FuncT& func, const Tpetra_Vector& p, Tpetra_CrsMatrix& jac)
{
  using Teuchos::ArrayView;
  using Teuchos::arrayView;

  jac.resumeFill();
  jac.setAllToScalar(0.0);

  Teuchos::ArrayRCP<const Scalar> p_view = p.get1dView();
  Teuchos::Array<Scalar> rv(p_view());
  size_t NumMyElements = x_map->getNodeNumElements();
  ArrayView<const GlobalOrdinal> MyGlobalElements = 
    x_map->getNodeElementList ();
  MeshScalar h2 = h*h;
  Scalar val;

  for(size_t i=0 ; i<NumMyElements; ++i ) {
      
    // Center
    GlobalOrdinal global_idx = MyGlobalElements[i];
    if (mesh[global_idx].boundary) {
      val = 1.0;
      jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1), 
			      arrayView(&val,1));
    }
    else {
      Scalar a_down = 
	-func(mesh[global_idx].x, mesh[global_idx].y-h/2.0, rv)/h2;
      Scalar a_left = 
	-func(mesh[global_idx].x-h/2.0, mesh[global_idx].y, rv)/h2;
      Scalar a_right = 
	-func(mesh[global_idx].x+h/2.0, mesh[global_idx].y, rv)/h2;
      Scalar a_up = 
	-func(mesh[global_idx].x, mesh[global_idx].y+h/2.0, rv)/h2;
      
      // Center
      val = -(a_down + a_left + a_right + a_up);
      jac.replaceGlobalValues(global_idx, arrayView(&global_idx,1), 
			      arrayView(&val,1));

      // Down
      if (!(eliminate_bcs && mesh[mesh[global_idx].down].boundary))
	jac.replaceGlobalValues(global_idx, 
				arrayView(&mesh[global_idx].down,1), 
				arrayView(&a_down,1));
      
      // Left
      if (!(eliminate_bcs && mesh[mesh[global_idx].left].boundary))
	jac.replaceGlobalValues(global_idx, 
				arrayView(&mesh[global_idx].left,1), 
				arrayView(&a_left,1));
      
      // Right
      if (!(eliminate_bcs && mesh[mesh[global_idx].right].boundary))
	jac.replaceGlobalValues(global_idx, 
				arrayView(&mesh[global_idx].right,1), 
				arrayView(&a_right,1));

      // Up
      if (!(eliminate_bcs && mesh[mesh[global_idx].up].boundary))
	jac.replaceGlobalValues(global_idx, 
				arrayView(&mesh[global_idx].up,1), 
				arrayView(&a_up,1));
    }
  }
  jac.fillComplete();
}