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

int main(int narg, char **arg)
{
  Teuchos::GlobalMPISession session(&narg, &arg);
  Teuchos::RCP<const Teuchos::Comm<int> > tcomm = 
               Teuchos::DefaultComm<int>::getComm();
  Teuchos::RCP<const Epetra_Comm> ecomm = Xpetra::toEpetra(tcomm);

  const int nGlobRows = 50;
  const Epetra_Map emap(nGlobRows, 0, *ecomm);
  Teuchos::RCP<const Epetra_BlockMap> ebmap = Teuchos::rcpFromRef(emap);

  typedef Xpetra::EpetraMapT<int, Tpetra::Map<>::node_type> xemap_t;
  Teuchos::RCP<const xemap_t> xmap(new xemap_t(ebmap));
  
  const Teuchos::RCP<const Teuchos::Comm<int> > &xcomm = xmap->getComm();

  std::cout << "Teuchos:  Hello from " 
            << tcomm->getRank() << " of " 
            << tcomm->getSize() << std::endl;
  std::cout << "Epetra:   Hello from " 
            << ecomm->MyPID() << " of " 
            << ecomm->NumProc() << std::endl;
  std::cout << "Xpetra:   Hello from " 
            << xcomm->getRank() << " of " 
            << xcomm->getSize() << std::endl;
}

void Iterative_Inverse_Operator::operator () (const Epetra_MultiVector &b, Epetra_MultiVector &x)
{
  // Initialize the solution to zero
  x.PutScalar( 0.0 );

  // Reset the solver, problem, and status test for next solve (HKT)
  pProb->setProblem( Teuchos::rcp(&x, false), Teuchos::rcp(&b, false) );

  timer.start();
  Belos::ReturnType ret = pBelos->solve();
  timer.stop();

  int pid = pComm->MyPID();

  if (pid == 0 && print) {
    if (ret == Belos::Converged)
    {
      std::cout << std::endl << "pid[" << pid << "] Block GMRES converged" << std::endl;
      std::cout << "Solution time: " << timer.totalElapsedTime() << std::endl;

    }
    else
      std::cout << std::endl << "pid[" << pid << "] Block GMRES did not converge" << std::endl;
  }
}

Real random(const Teuchos::RCP<Epetra_Comm> &comm) {
  Real val = 0.0;
  if ( comm->MyPID()==0 ) {
    val = (Real)rand()/(Real)RAND_MAX;
  }
  comm->Broadcast(&val,1,0);
  return val;
}

TEUCHOS_UNIT_TEST(tSGEpetraLinearObjFactory, basic)
{
   // build global (or serial communicator)
   #ifdef HAVE_MPI
      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
   #else
      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
   #endif

   using Teuchos::RCP;
   using Teuchos::rcp;
   using Teuchos::rcp_dynamic_cast;

   int myRank = eComm->MyPID();
   int numProc = eComm->NumProc();

   // panzer::pauseToAttach();

   RCP<Stokhos::OrthogPolyExpansion<int,double> > sgExpansion = buildExpansion(3,5);
   RCP<panzer::UniqueGlobalIndexer<int,int> > indexer 
         = rcp(new unit_test::UniqueGlobalIndexer<int>(myRank,numProc));
 
   // setup factory
   RCP<panzer::EpetraLinearObjFactory<panzer::Traits,int> > epetraFactory
         = rcp(new panzer::EpetraLinearObjFactory<panzer::Traits,int>(eComm.getConst(),indexer));
   RCP<panzer::SGEpetraLinearObjFactory<panzer::Traits,int> > la_factory
         = rcp(new panzer::SGEpetraLinearObjFactory<panzer::Traits,int>(epetraFactory,sgExpansion,Teuchos::null));
   
   RCP<panzer::LinearObjContainer> ghostedContainer = la_factory->buildGhostedLinearObjContainer(); 
   RCP<panzer::LinearObjContainer> container = la_factory->buildLinearObjContainer(); 

   TEST_NOTHROW(rcp_dynamic_cast<panzer::SGEpetraLinearObjContainer>(ghostedContainer,true));
   TEST_NOTHROW(rcp_dynamic_cast<panzer::SGEpetraLinearObjContainer>(container,true));

   RCP<panzer::SGEpetraLinearObjContainer> sgGhostedContainer 
         = rcp_dynamic_cast<panzer::SGEpetraLinearObjContainer>(ghostedContainer,true);
   RCP<panzer::SGEpetraLinearObjContainer> sgContainer 
         = rcp_dynamic_cast<panzer::SGEpetraLinearObjContainer>(container,true);

   // make sure we have correct number of sub containers
   TEST_EQUALITY(sgExpansion->size(),sgGhostedContainer->end()-sgGhostedContainer->begin());
   TEST_EQUALITY(sgExpansion->size(),sgContainer->end()-sgContainer->begin());

   // make sure all "sub-containers" are epetra containers
   panzer::SGEpetraLinearObjContainer::const_iterator itr;
   for(itr=sgGhostedContainer->begin();itr!=sgGhostedContainer->end();++itr)
      TEST_NOTHROW(rcp_dynamic_cast<EpetraLinearObjContainer>(*itr));
   for(itr=sgContainer->begin();itr!=sgContainer->end();++itr)
      TEST_NOTHROW(rcp_dynamic_cast<EpetraLinearObjContainer>(*itr));

   la_factory->ghostToGlobalContainer(*ghostedContainer,*container,0);
   la_factory->globalToGhostContainer(*container,*ghostedContainer,0);
}

TEUCHOS_UNIT_TEST(tCloneLOF, blocked_epetra)
{

   // build global (or serial communicator)
   #ifdef HAVE_MPI
      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
   #else
      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
   #endif

   Teuchos::RCP<const Teuchos::MpiComm<int> > tComm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));


   int myRank = eComm->MyPID();
   int numProc = eComm->NumProc();

   RCP<ConnManager<int,int> > connManager = rcp(new unit_test::ConnManager(myRank,numProc));

   RCP<const FieldPattern> patternC1
         = buildFieldPattern<Intrepid2::Basis_HGRAD_HEX_C1_FEM<double,FieldContainer> >();

   RCP<panzer::BlockedDOFManager<int,int> > indexer = rcp(new panzer::BlockedDOFManager<int,int>());
   indexer->setConnManager(connManager,MPI_COMM_WORLD);
   indexer->addField("U",patternC1);
   indexer->addField("V",patternC1);

   std::vector<std::vector<std::string> > fieldOrder(2);
   fieldOrder[0].push_back("U");
   fieldOrder[1].push_back("V");
   indexer->setFieldOrder(fieldOrder);
   indexer->buildGlobalUnknowns();

   RCP<panzer::BlockedDOFManager<int,int> > control_indexer = rcp(new panzer::BlockedDOFManager<int,int>());
   control_indexer->setConnManager(connManager,MPI_COMM_WORLD);
   control_indexer->addField("Z",patternC1);
   fieldOrder[0][0] = "Z";
   control_indexer->buildGlobalUnknowns();
 
   // setup factory
   RCP<BlockedEpetraLinearObjFactory<Traits,int> > ep_lof
         = Teuchos::rcp(new BlockedEpetraLinearObjFactory<Traits,int>(tComm,indexer));
   
   // this is the member we are testing!
   RCP<const LinearObjFactory<Traits> > control_lof = cloneWithNewDomain(*ep_lof,control_indexer);
   RCP<const BlockedEpetraLinearObjFactory<Traits,int> > ep_control_lof 
       = rcp_dynamic_cast<const BlockedEpetraLinearObjFactory<Traits,int> >(control_lof);

/*
   TEST_ASSERT(ep_control_lof->getMap()->SameAs(*ep_lof->getMap()));
   TEST_EQUALITY(ep_control_lof->getColMap()->NumMyElements(),Teuchos::as<int>(control_owned.size()));
*/
}

void
Albany::IossSTKMeshStruct::readSerialMesh(const Teuchos::RCP<const Epetra_Comm>& comm,
                                          std::vector<std::string>& entity_rank_names){

#ifdef ALBANY_ZOLTAN // rebalance needs Zoltan

  MPI_Group group_world;
  MPI_Group peZero;
  MPI_Comm peZeroComm;

  // Read a single exodus mesh on Proc 0 then rebalance it across the machine

  MPI_Comm theComm = Albany::getMpiCommFromEpetraComm(*comm);

  int process_rank[1]; // the reader process

  process_rank[0] = 0;
  int my_rank = comm->MyPID();

  //get the group under theComm
  MPI_Comm_group(theComm, &group_world);
  // create the new group. This group includes only processor zero - that is the only processor that reads the file
  MPI_Group_incl(group_world, 1, process_rank, &peZero);
  // create the new communicator - it just contains processor zero
  MPI_Comm_create(theComm, peZero, &peZeroComm);

  // Note that peZeroComm == MPI_COMM_NULL on all processors but processor 0

  if(my_rank == 0){

    *out << "Albany_IOSS: Loading serial STKMesh from Exodus file  " 
         << params->get<std::string>("Exodus Input File Name") << std::endl;

  }

  /* 
   * This checks the existence of the file, checks to see if we can open it, builds a handle to the region
   * and puts it in mesh_data (in_region), and reads the metaData into metaData.
   */

  stk::io::create_input_mesh("exodusii",
//  create_input_mesh("exodusii",
                             params->get<std::string>("Exodus Input File Name"), 
                             peZeroComm, 
                             *metaData, *mesh_data,
                             entity_rank_names); 

  // Here, all PEs have read the metaData from the input file, and have a pointer to in_region in mesh_data

#endif

}

Albany::AsciiSTKMesh2D::AsciiSTKMesh2D(
		const Teuchos::RCP<Teuchos::ParameterList>& params,
		const Teuchos::RCP<const Epetra_Comm>& comm) :
		GenericSTKMeshStruct(params, Teuchos::null, 2), out(
				Teuchos::VerboseObjectBase::getDefaultOStream()), periodic(false), sh(0) {
	std::string fname = params->get("Ascii Input Mesh File Name",
			"greenland.msh");

	std::string shape;
	if (comm->MyPID() == 0) {
		std::ifstream ifile;

		NumElemNodes = 0;
		ifile.open(fname.c_str());
		if (ifile.is_open()) {
		  ifile >> shape >> NumElemNodes;
		  if(shape == "Triangle") {
		    TEUCHOS_TEST_FOR_EXCEPTION(NumElemNodes != 3, Teuchos::Exceptions::InvalidParameter,
		              std::endl << "Error in AsciiSTKMesh2D: Triangles must be linear. Number of nodes per element in file " << fname << " is: " << NumElemNodes << ". Should be 3!" << std::endl);
		  }
		  else if(shape == "Quadrilateral") {
        TEUCHOS_TEST_FOR_EXCEPTION(NumElemNodes != 4, Teuchos::Exceptions::InvalidParameter,
                  std::endl << "Error in AsciiSTKMesh2D: Quadrilaterals must be bilinear. Number of nodes per element in file " << fname << " is: "  << " is: " << NumElemNodes << ". Should be 4!" << std::endl);
      }
		  else
		    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
		                      std::endl << "Error in AsciiSTKMesh2D: Only Triangle or Quadrilateral grids can be imported. Shape in in file " << fname << " is: " << shape << ". Should be Triangle or Quadrialteral" << std::endl);


			ifile >> NumNodes >> NumEles >> NumBdEdges;
			//read in nodes coordinates
			xyz = new double[NumNodes][3];
			for (int i = 0; i < NumNodes; i++) {
				ifile >> xyz[i][0] >> xyz[i][1] >> xyz[i][2];
			}

			//read in element connectivity
			eles = new int[NumEles][4];
			int temp;

			if(shape == "Triangle")
        for (int i = 0; i < NumEles; i++) {
          ifile >> eles[i][0] >> eles[i][1] >> eles[i][2] >> temp;
			  }
			else
			  for (int i = 0; i < NumEles; i++) {

Albany::AsciiSTKMesh2D::AsciiSTKMesh2D(
		const Teuchos::RCP<Teuchos::ParameterList>& params,
		const Teuchos::RCP<const Epetra_Comm>& comm) :
		GenericSTKMeshStruct(params, Teuchos::null, 2), out(
				Teuchos::VerboseObjectBase::getDefaultOStream()), periodic(false), sh(0) {
	std::string fname = params->get("Ascii Input Mesh File Name",
			"greenland.msh");

	if (comm->MyPID() == 0) {
		std::ifstream ifile;
		ifile.open(fname.c_str());
		if (ifile.is_open()) {
			ifile >> NumNodes >> NumEles >> NumBdEdges;
			//read in nodes coordinates
			xyz = new double[NumNodes][3];
			for (int i = 0; i < NumNodes; i++) {
				ifile >> xyz[i][0] >> xyz[i][1] >> xyz[i][2];
				*out << "i: " << i << ", x: " << xyz[i][0] << ", y: " << xyz[i][1]
						<< ", z: " << xyz[i][2] << std::endl;
			}
			//read in element connectivity
			eles = new int[NumEles][4];
			int temp;
			for (int i = 0; i < NumEles; i++) {
				ifile >> eles[i][0] >> eles[i][1] >> eles[i][2] >> temp;
				*out << "elm" << i << ": " << eles[i][0] << " " << eles[i][1] << " "
						<< eles[i][2] << std::endl;
			}
			//read in boundary edges connectivity
			be = new int[NumBdEdges][2];
			for (int i = 0; i < NumBdEdges; i++) {
				ifile >> be[i][0] >> be[i][1] >> temp;
				*out << "edge #:" << i << " " << be[i][0] << " " << be[i][1]
						<< std::endl;
			}
			ifile.close();
		} else {

void
Albany::IossSTKMeshStruct::setFieldAndBulkData(
                                               const Teuchos::RCP<const Epetra_Comm>& comm,
                                               const Teuchos::RCP<Teuchos::ParameterList>& params,
                                               const unsigned int neq_,
                                               const AbstractFieldContainer::FieldContainerRequirements& req,
                                               const Teuchos::RCP<Albany::StateInfoStruct>& sis,
                                               const unsigned int worksetSize)
{
  this->SetupFieldData(comm, neq_, req, sis, worksetSize);

  *out << "IOSS-STK: number of node sets = " << nsPartVec.size() << std::endl;
  *out << "IOSS-STK: number of side sets = " << ssPartVec.size() << std::endl;

  metaData->commit();

  // Restart index to read solution from exodus file.
  int index = params->get("Restart Index",-1); // Default to no restart
  double res_time = params->get<double>("Restart Time",-1.0); // Default to no restart
  Ioss::Region *region = mesh_data->m_input_region;

  /*
   * The following code block reads a single mesh on PE 0, then distributes the mesh across
   * the other processors. stk_rebalance is used, which requires Zoltan
   *
   * This code is only compiled if ALBANY_MPI and ALBANY_ZOLTAN are true
   */

#ifdef ALBANY_ZOLTAN // rebalance needs Zoltan

  if(useSerialMesh){

    bulkData->modification_begin();

    if(comm->MyPID() == 0){ // read in the mesh on PE 0

      stk::io::process_mesh_bulk_data(region, *bulkData);

      // Read solution from exodus file.
      if (index >= 0) { // User has specified a time step to restart at
        *out << "Restart Index set, reading solution index : " << index << std::endl;
        stk::io::input_mesh_fields(region, *bulkData, index);
        m_restartDataTime = region->get_state_time(index);
        m_hasRestartSolution = true;
      }
      else if (res_time >= 0) { // User has specified a time to restart at
        *out << "Restart solution time set, reading solution time : " << res_time << std::endl;
        stk::io::input_mesh_fields(region, *bulkData, res_time);
        m_restartDataTime = res_time;
        m_hasRestartSolution = true;
      }
      else {

        *out << "Neither restart index or time are set. Not reading solution data from exodus file"<< std::endl;

      }
    }

    bulkData->modification_end();

  } // End UseSerialMesh - reading mesh on PE 0

  else 
#endif

    /*
     * The following code block reads a single mesh when Albany is compiled serially, or a
     * Nemspread fileset if ALBANY_MPI is true.
     *
     */

  { // running in Serial or Parallel read from Nemspread files

    stk::io::populate_bulk_data(*bulkData, *mesh_data);

    if (!usePamgen)  {

      // Read solution from exodus file.
      if (index >= 0) { // User has specified a time step to restart at
        *out << "Restart Index set, reading solution index : " << index << std::endl;
        stk::io::process_input_request(*mesh_data, *bulkData, index);
        m_restartDataTime = region->get_state_time(index);
        m_hasRestartSolution = true;
      }
      else if (res_time >= 0) { // User has specified a time to restart at
        *out << "Restart solution time set, reading solution time : " << res_time << std::endl;
        stk::io::process_input_request(*mesh_data, *bulkData, res_time);
        m_restartDataTime = res_time;
        m_hasRestartSolution = true;
      }
      else {
        *out << "Restart Index not set. Not reading solution from exodus (" 
             << index << ")"<< std::endl;

      }
    }

    bulkData->modification_end();

  } // End Parallel Read - or running in serial
//.........这里部分代码省略.........

Albany::AsciiSTKMeshStruct::AsciiSTKMeshStruct(
                                             const Teuchos::RCP<Teuchos::ParameterList>& params, 
                                             const Teuchos::RCP<const Epetra_Comm>& comm) :
  GenericSTKMeshStruct(params,Teuchos::null,3),
  out(Teuchos::VerboseObjectBase::getDefaultOStream()),
  periodic(false)
{
   int numProc = comm->NumProc(); //total number of processors
   contigIDs = params->get("Contiguous IDs", true); 
   std::cout << "Number of processors: " << numProc << std::endl; 
   //names of files giving the mesh
   char meshfilename[100]; 
   char shfilename[100];
   char confilename[100];
   char bffilename[100];
   char geIDsfilename[100];
   char gnIDsfilename[100];
   char bfIDsfilename[100];
   char flwafilename[100]; //flow factor file
   char tempfilename[100]; //temperature file
   char betafilename[100]; //basal friction coefficient file
   if ((numProc == 1) & (contigIDs == true)) { //serial run with contiguous global IDs
     std::cout << "Ascii mesh has contiguous IDs; no bfIDs, geIDs, gnIDs files required." << std::endl;
     sprintf(meshfilename, "%s", "xyz");
     sprintf(shfilename, "%s", "sh");
     sprintf(confilename, "%s", "eles");
     sprintf(bffilename, "%s", "bf");
     sprintf(flwafilename, "%s", "flwa");
     sprintf(tempfilename, "%s", "temp");
     sprintf(betafilename, "%s", "beta");
   }
   else { //parallel run or serial run with non-contiguous global IDs - proc # is appended to file name to indicate what processor the mesh piece is on 
     if ((numProc == 1) & (contigIDs == false))
        std::cout << "1 processor run with non-contiguous IDs; bfIDs0, geIDs0, gnIDs0 files required." << std::endl;
     int suffix = comm->MyPID(); //current processor number 
     sprintf(meshfilename, "%s%i", "xyz", suffix);
     sprintf(shfilename, "%s%i", "sh", suffix);
     sprintf(confilename, "%s%i", "eles", suffix);
     sprintf(bffilename, "%s%i", "bf", suffix);
     sprintf(geIDsfilename, "%s%i", "geIDs", suffix);
     sprintf(gnIDsfilename, "%s%i", "gnIDs", suffix);
     sprintf(bfIDsfilename, "%s%i", "bfIDs", suffix);
     sprintf(flwafilename, "%s%i", "flwa", suffix);
     sprintf(tempfilename, "%s%i", "temp", suffix);
     sprintf(betafilename, "%s%i", "beta", suffix);
   }

    //read in coordinates of mesh -- right now hard coded for 3D
    //assumes mesh file is called "xyz" and its first row is the number of nodes  
    FILE *meshfile = fopen(meshfilename,"r");
    if (meshfile == NULL) { //check if coordinates file exists
      *out << "Error in AsciiSTKMeshStruct: coordinates file " << meshfilename <<" not found!"<< std::endl;
      TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
          std::endl << "Error in AsciiSTKMeshStruct: coordinates file " << meshfilename << " not found!"<< std::endl);
    }
    double temp; 
    fseek(meshfile, 0, SEEK_SET); 
    fscanf(meshfile, "%lf", &temp); 
    NumNodes = int(temp); 
    std::cout << "numNodes: " << NumNodes << std::endl;  
    xyz = new double[NumNodes][3]; 
    char buffer[100];
    fgets(buffer, 100, meshfile); 
    for (int i=0; i<NumNodes; i++){
      fgets(buffer, 100, meshfile); 
      sscanf(buffer, "%lf %lf %lf", &xyz[i][0], &xyz[i][1], &xyz[i][2]); 
      //*out << "i: " << i << ", x: " << xyz[i][0] << ", y: " << xyz[i][1] << ", z: " << xyz[i][2] << std::endl; 
     }
    //read in surface height data from mesh 
    //assumes surface height file is called "sh" and its first row is the number of nodes  
    FILE *shfile = fopen(shfilename,"r");
    have_sh = false;
    if (shfile != NULL) have_sh = true;
    if (have_sh) {
      fseek(shfile, 0, SEEK_SET); 
      fscanf(shfile, "%lf", &temp); 
      int NumNodesSh = int(temp);
      std::cout << "NumNodesSh: " << NumNodesSh<< std::endl; 
      if (NumNodesSh != NumNodes) { 
           *out << "Error in AsciiSTKMeshStruct: sh file must have same number nodes as xyz file!  numNodes in xyz = " << NumNodes <<", numNodes in sh = "<< NumNodesSh  << std::endl;
          TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
            std::endl << "Error in AsciiSTKMeshStruct: sh file must have same number nodes as xyz file!  numNodes in xyz = " << NumNodes << ", numNodes in sh = "<< NumNodesSh << std::endl);
      }
      sh = new double[NumNodes]; 
      fgets(buffer, 100, shfile); 
      for (int i=0; i<NumNodes; i++){
        fgets(buffer, 100, shfile); 
        sscanf(buffer, "%lf", &sh[i]); 
        //*out << "i: " << i << ", sh: " << sh[i] << std::endl; 
       }
     }
     //read in connectivity file -- right now hard coded for 3D hexes
     //assumes mesh file is called "eles" and its first row is the number of elements  
     FILE *confile = fopen(confilename,"r"); 
     if (confile == NULL) { //check if element connectivity file exists
      *out << "Error in AsciiSTKMeshStruct: element connectivity file " << confilename <<" not found!"<< std::endl;
      TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
          std::endl << "Error in AsciiSTKMeshStruct: element connectivity file " << confilename << " not found!"<< std::endl);
     }
     fseek(confile, 0, SEEK_SET); 
//.........这里部分代码省略.........

void
createEpetraProblem (const Teuchos::RCP<const Epetra_Comm>& epetraComm,
		     const std::string& filename,
		     Teuchos::RCP<Epetra_Map>& rowMap,
		     Teuchos::RCP<Epetra_CrsMatrix>& A,
		     Teuchos::RCP<Epetra_MultiVector>& B,
		     Teuchos::RCP<Epetra_MultiVector>& X,
		     int& numRHS)
{
  using Teuchos::inOutArg;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::set_extra_data;

  const int MyPID = epetraComm->MyPID();

  int n_nonzeros, N_update;
  int *bindx = NULL, *update = NULL, *col_inds = NULL;
  double *val = NULL, *row_vals = NULL;
  double *xguess = NULL, *b = NULL, *xexact = NULL;

  //
  // Set up the problem to be solved
  //
  int NumGlobalElements;  // total # of rows in matrix

  try {
    // Read in matrix from HB file
    Trilinos_Util_read_hb (const_cast<char *> (filename.c_str()), MyPID, 
			   &NumGlobalElements, &n_nonzeros, &val, &bindx, 
			   &xguess, &b, &xexact);

    // Distribute data among processors
    Trilinos_Util_distrib_msr_matrix (*epetraComm, &NumGlobalElements, 
				      &n_nonzeros, &N_update, &update, &val, 
				      &bindx, &xguess, &b, &xexact);
    //
    // Construct the matrix
    //
    int NumMyElements = N_update; // # local rows of matrix on processor
    //
    // Create an int array NumNz that is used to build the Petra Matrix.
    // NumNz[i] is the number of OFF-DIAGONAL terms for the i-th global
    // equation on this processor.
    //
    std::vector<int> NumNz (NumMyElements);
    for (int i = 0; i < NumMyElements; ++i) {
      NumNz[i] = bindx[i+1] - bindx[i] + 1;
    }
    rowMap = rcp (new Epetra_Map (NumGlobalElements, NumMyElements, 
				  update, 0, *epetraComm));
    //set_extra_data (epetraComm, "Map::Comm", inOutArg (rowMap));

    // Create an Epetra sparse matrix.
    if (NumMyElements == 0) {
      A = rcp (new Epetra_CrsMatrix (Copy, *rowMap, static_cast<int*>(NULL)));
    } else {
      A = rcp (new Epetra_CrsMatrix (Copy, *rowMap, &NumNz[0]));
    }
    //set_extra_data (rowMap, "Operator::Map", ptr (A));

    // Add rows to the sparse matrix one at a time.
    int NumEntries;
    for (int i = 0; i < NumMyElements; ++i) {
      row_vals = val + bindx[i];
      col_inds = bindx + bindx[i];
      NumEntries = bindx[i+1] - bindx[i];
      int info = A->InsertGlobalValues (update[i], NumEntries, row_vals, col_inds);
      TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				  "Failed to insert global value into A." );
      info = A->InsertGlobalValues (update[i], 1, val+i, update+i);
      TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				  "Failed to insert global value into A." );
    }
    // Finish initializing the sparse matrix.
    int info = A->FillComplete();
    TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				"FillComplete() failed on the sparse matrix A.");
    info = A->OptimizeStorage();
    TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				"OptimizeStorage() failed on the sparse matrix A.");
    A->SetTracebackMode (1); // Shut down Epetra warning tracebacks
    //
    // Construct the right-hand side and solution multivectors.
    //
    if (false && b != NULL) {
      B = rcp (new Epetra_MultiVector (::Copy, *rowMap, b, NumMyElements, 1));
      numRHS = 1;
    } else {
      B = rcp (new Epetra_MultiVector (*rowMap, numRHS));
      B->Random ();
    }
    X = rcp (new Epetra_MultiVector (*rowMap, numRHS));
    X->PutScalar (0.0);

    //set_extra_data (rowMap, "X::Map", Teuchos::ptr (X));
    //set_extra_data (rowMap, "B::Map", Teuchos::ptr (B));
  } 
  catch (std::exception& e) {
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

  int status=0; // 0 = pass, failures are incremented
  bool success = true;
  Teuchos::GlobalMPISession mpiSession(&argc,&argv);

  using Teuchos::RCP;
  using Teuchos::rcp;

  RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
  
  //***********************************************************
  // Command-line argument for input file
  //***********************************************************

  Albany::CmdLineArgs cmd;
  cmd.parse_cmdline(argc, argv, *out);
  std::string xmlfilename_coupled = cmd.xml_filename;


  try {

    RCP<Teuchos::Time> totalTime = 
      Teuchos::TimeMonitor::getNewTimer("AlbanySG: ***Total Time***");
    RCP<Teuchos::Time> setupTime = 
      Teuchos::TimeMonitor::getNewTimer("AlbanySG: Setup Time");
    Teuchos::TimeMonitor totalTimer(*totalTime); //start timer

    //***********************************************************
    // Set up coupled solver first to setup comm's
    //***********************************************************
    Teuchos::RCP<Epetra_Comm> globalComm = 
      Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
    // Connect vtune for performance profiling
    if (cmd.vtune) {
      Albany::connect_vtune(globalComm->MyPID());
    }
    Albany::SolverFactory coupled_slvrfctry(xmlfilename_coupled, 
					    Albany::createTeuchosCommFromEpetraComm(globalComm));
    Teuchos::ParameterList& coupledParams = coupled_slvrfctry.getParameters();
    Teuchos::ParameterList& coupledSystemParams = 
      coupledParams.sublist("Coupled System");
    Teuchos::Array<std::string> model_filenames =
      coupledSystemParams.get<Teuchos::Array<std::string> >("Model XML Files");
    int num_models = model_filenames.size();
    Teuchos::Array< RCP<Albany::Application> > apps(num_models);
    Teuchos::Array< RCP<EpetraExt::ModelEvaluator> > models(num_models);
    Teuchos::Array< RCP<Teuchos::ParameterList> > piroParams(num_models);
    Teuchos::RCP< Teuchos::ParameterList> coupledPiroParams = 
      Teuchos::rcp(&(coupledParams.sublist("Piro")),false);
    Teuchos::RCP<Piro::Epetra::StokhosSolver> coupledSolver =
      Teuchos::rcp(new Piro::Epetra::StokhosSolver(coupledPiroParams, 
						   globalComm));
    Teuchos::RCP<const Epetra_Comm> app_comm = coupledSolver->getSpatialComm();

    // Set up each model
    Teuchos::Array< Teuchos::RCP<NOX::Epetra::Observer> > observers(num_models);
    for (int m=0; m<num_models; m++) {
      Albany::SolverFactory slvrfctry(
	model_filenames[m], 
	Albany::createTeuchosCommFromEpetraComm(app_comm));
      models[m] = slvrfctry.createAlbanyAppAndModel(apps[m], app_comm);
      Teuchos::ParameterList& appParams = slvrfctry.getParameters();
      piroParams[m] = Teuchos::rcp(&(appParams.sublist("Piro")),false);
      observers[m] = Teuchos::rcp(new Albany_NOXObserver(apps[m]));
    }

    // Setup network model
    std::string network_name = 
      coupledSystemParams.get("Network Model", "Param To Response");
    RCP<Piro::Epetra::AbstractNetworkModel> network_model;
    if (network_name == "Param To Response")
      network_model = rcp(new Piro::Epetra::ParamToResponseNetworkModel);
    else if (network_name == "Reactor Network")
      network_model = rcp(new Albany::ReactorNetworkModel(1));
    else
      TEUCHOS_TEST_FOR_EXCEPTION(
	true, std::logic_error, "Invalid network model name " << network_name);
    RCP<EpetraExt::ModelEvaluator> coupledModel =
      rcp(new Piro::Epetra::NECoupledModelEvaluator(models, piroParams,
						    network_model,
						    coupledPiroParams, 
						    globalComm,
						    observers));
    coupledSolver->setup(coupledModel);

    // Solve coupled system
    EpetraExt::ModelEvaluator::InArgs inArgs = coupledSolver->createInArgs();
    EpetraExt::ModelEvaluator::OutArgs outArgs = coupledSolver->createOutArgs();
    for (int i=0; i<inArgs.Np(); i++)
      if (inArgs.supports(EpetraExt::ModelEvaluator::IN_ARG_p_sg, i))
	inArgs.set_p_sg(i, coupledSolver->get_p_sg_init(i));
    for (int i=0; i<outArgs.Ng(); i++) 
      if (outArgs.supports(EpetraExt::ModelEvaluator::OUT_ARG_g_sg, i)) {
	RCP<Stokhos::EpetraVectorOrthogPoly> g_sg = 
	  coupledSolver->create_g_sg(i);
	outArgs.set_g_sg(i, g_sg);
      }
    coupledSolver->evalModel(inArgs, outArgs);

//.........这里部分代码省略.........

void
Albany::MpasSTKMeshStruct::constructMesh(
                                               const Teuchos::RCP<const Epetra_Comm>& comm,
                                               const Teuchos::RCP<Teuchos::ParameterList>& params,
                                               const unsigned int neq_,
                                               const Albany::AbstractFieldContainer::FieldContainerRequirements& req,
                                               const Teuchos::RCP<Albany::StateInfoStruct>& sis,
                                               const std::vector<int>& indexToVertexID, const std::vector<double>& verticesCoords, const std::vector<bool>& isVertexBoundary, int nGlobalVertices,
                                               const std::vector<int>& verticesOnTria,
                                               const std::vector<bool>& isBoundaryEdge, const std::vector<int>& trianglesOnEdge, const std::vector<int>& trianglesPositionsOnEdge,
                                               const std::vector<int>& verticesOnEdge,
                                               const std::vector<int>& indexToEdgeID, int nGlobalEdges,
                                               const unsigned int worksetSize)
{
  this->SetupFieldData(comm, neq_, req, sis, worksetSize);

  metaData->commit();

  bulkData->modification_begin(); // Begin modifying the mesh

  stk::mesh::PartVector nodePartVec;
  stk::mesh::PartVector singlePartVec(1);
  stk::mesh::PartVector emptyPartVec;
  std::cout << "elem_map # elments: " << elem_map->NumMyElements() << std::endl;
  unsigned int ebNo = 0; //element block #??? 
  int sideID = 0;

  AbstractSTKFieldContainer::IntScalarFieldType* proc_rank_field = fieldContainer->getProcRankField();
  AbstractSTKFieldContainer::VectorFieldType* coordinates_field = fieldContainer->getCoordinatesField();

  for (int i=0; i<indexToVertexID.size(); i++)
  {
	  stk::mesh::Entity node = bulkData->declare_entity(stk::topology::NODE_RANK, indexToVertexID[i]+1, nodePartVec);

	  double* coord;
	  coord = stk::mesh::field_data(*coordinates_field, node);
	  coord[0] = verticesCoords[3*i];   coord[1] = verticesCoords[3*i+1]; coord[2] = verticesCoords[3*i+2];
  }

  for (int i=0; i<elem_map->NumMyElements(); i++)
  {

     singlePartVec[0] = partVec[ebNo];
     stk::mesh::Entity elem  = bulkData->declare_entity(stk::topology::ELEMENT_RANK, elem_map->GID(i)+1, singlePartVec);

     for(int j=0; j<3; j++)
     {
    	 stk::mesh::Entity node = bulkData->get_entity(stk::topology::NODE_RANK, indexToVertexID[verticesOnTria[3*i+j]]+1);
    	 bulkData->declare_relation(elem, node, j);
     }
    
     int* p_rank = (int*)stk::mesh::field_data(*proc_rank_field, elem);
     p_rank[0] = comm->MyPID();
  }

  for (int i=0; i<indexToEdgeID.size(); i++) {

	 if(isBoundaryEdge[i])
	 {

		 singlePartVec[0] = ssPartVec["lateralside"];
		 stk::mesh::Entity side = bulkData->declare_entity(metaData->side_rank(), indexToEdgeID[i]+1, singlePartVec);
		 stk::mesh::Entity elem  = bulkData->get_entity(stk::topology::ELEMENT_RANK,  elem_map->GID(trianglesOnEdge[2*i])+1);
		 bulkData->declare_relation(elem, side,  trianglesPositionsOnEdge[2*i] /*local side id*/);
		 for(int j=0; j<2; j++)
		 {
			 stk::mesh::Entity node = bulkData->get_entity(stk::topology::NODE_RANK, indexToVertexID[verticesOnEdge[2*i+j]]+1);
			 bulkData->declare_relation(side, node, j);
		 }
	 }
  }

  bulkData->modification_end();
}

int main(int argc, char *argv[]) {

// Initialize MPI
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  int MyPID;

  try {

    // Create a communicator for Epetra objects
    Teuchos::RCP<const Epetra_Comm> globalComm;
#ifdef HAVE_MPI
    globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
#else
    globalComm = Teuchos::rcp(new Epetra_SerialComm);
#endif
    MyPID = globalComm->MyPID();
    
    // Create Stochastic Galerkin basis and expansion
    int p = 5;
    Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(1); 
    bases[0] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p));
    Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = 
      Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
    int sz = basis->size();
    Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
    Cijk = basis->computeTripleProductTensor();
    Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion = 
      Teuchos::rcp(new Stokhos::AlgebraicOrthogPolyExpansion<int,double>(basis,
									 Cijk));
    if (MyPID == 0)
      std::cout << "Stochastic Galerkin expansion size = " << sz << std::endl;

    // Create stochastic parallel distribution
    int num_spatial_procs = -1;
    Teuchos::ParameterList parallelParams;
    parallelParams.set("Number of Spatial Processors", num_spatial_procs);
    Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
      Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
					     parallelParams));
    Teuchos::RCP<const EpetraExt::MultiComm> sg_comm = 
      sg_parallel_data->getMultiComm();
    Teuchos::RCP<const Epetra_Comm> app_comm = 
      sg_parallel_data->getSpatialComm();
    
    // Create application model evaluator
    Teuchos::RCP<EpetraExt::ModelEvaluator> model = 
      Teuchos::rcp(new SimpleME(app_comm));
    
    // Setup stochastic Galerkin algorithmic parameters
    Teuchos::RCP<Teuchos::ParameterList> sgParams = 
      Teuchos::rcp(new Teuchos::ParameterList);
    sgParams->set("Jacobian Method", "Matrix Free");
    sgParams->set("Mean Preconditioner Type", "Ifpack");

    // Create stochastic Galerkin model evaluator
    Teuchos::RCP<Stokhos::SGModelEvaluator> sg_model =
      Teuchos::rcp(new Stokhos::SGModelEvaluator(model, basis, Teuchos::null,
    						 expansion, sg_parallel_data, 
						 sgParams));

    // Stochastic Galerkin initial guess
    // Set the mean to the deterministic initial guess, higher-order terms
    // to zero
    Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_init_sg = 
      sg_model->create_x_sg();
    x_init_sg->init(0.0);
    (*x_init_sg)[0] = *(model->get_x_init());
    sg_model->set_x_sg_init(*x_init_sg);

    // Stochastic Galerkin parameters
    // Linear expansion with the mean given by the deterministic initial
    // parameter values, linear terms equal to 1, and higher order terms
    // equal to zero.
    Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> p_init_sg =
      sg_model->create_p_sg(0);
    p_init_sg->init(0.0);
    (*p_init_sg)[0] = *(model->get_p_init(0));
    for (int i=0; i<model->get_p_map(0)->NumMyElements(); i++)
      (*p_init_sg)[i+1][i] = 1.0;
    sg_model->set_p_sg_init(0, *p_init_sg);
    std::cout << "Stochatic Galerkin parameter expansion = " << std::endl
	      << *p_init_sg << std::endl;

    // Set up NOX parameters
    Teuchos::RCP<Teuchos::ParameterList> noxParams =
      Teuchos::rcp(new Teuchos::ParameterList);

    // Set the nonlinear solver method
    noxParams->set("Nonlinear Solver", "Line Search Based");

    // Set the printing parameters in the "Printing" sublist
    Teuchos::ParameterList& printParams = noxParams->sublist("Printing");
    printParams.set("MyPID", MyPID); 
    printParams.set("Output Precision", 3);
    printParams.set("Output Processor", 0);
    printParams.set("Output Information", 
                    NOX::Utils::OuterIteration + 
//.........这里部分代码省略.........

  TEUCHOS_UNIT_TEST(block_assembly, scatter_dirichlet_residual)
  {
    PHX::KokkosDeviceSession session;

#ifdef HAVE_MPI
    Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
#else
    Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
#endif

    int myRank = eComm->MyPID();

    const std::size_t workset_size = 4;
    const std::string fieldName1_q1 = "U";
    const std::string fieldName2_q1 = "V";
    const std::string fieldName_qedge1 = "B";

    Teuchos::RCP<panzer_stk_classic::STK_Interface> mesh = buildMesh(2,2);

    // build input physics block
    Teuchos::RCP<panzer::PureBasis> basis_q1 = buildBasis(workset_size,"Q1");
    Teuchos::RCP<panzer::PureBasis> basis_qedge1 = buildBasis(workset_size,"QEdge1");

    Teuchos::RCP<Teuchos::ParameterList> ipb = Teuchos::parameterList();
    testInitialization(ipb);

    const int default_int_order = 1;
    std::string eBlockID = "eblock-0_0";    
    Teuchos::RCP<user_app::MyFactory> eqset_factory = Teuchos::rcp(new user_app::MyFactory);
    panzer::CellData cellData(workset_size,mesh->getCellTopology("eblock-0_0"));
    Teuchos::RCP<panzer::GlobalData> gd = panzer::createGlobalData();
    Teuchos::RCP<panzer::PhysicsBlock> physicsBlock = 
      Teuchos::rcp(new PhysicsBlock(ipb,eBlockID,default_int_order,cellData,eqset_factory,gd,false));

    Teuchos::RCP<std::vector<panzer::Workset> > work_sets = panzer_stk_classic::buildWorksets(*mesh,*physicsBlock); 
    TEST_EQUALITY(work_sets->size(),1);

    // build connection manager and field manager
    const Teuchos::RCP<panzer::ConnManager<int,int> > conn_manager = Teuchos::rcp(new panzer_stk_classic::STKConnManager<int>(mesh));
    RCP<panzer::BlockedDOFManager<int,int> > dofManager = Teuchos::rcp(new panzer::BlockedDOFManager<int,int>(conn_manager,MPI_COMM_WORLD));

    dofManager->addField(fieldName1_q1,Teuchos::rcp(new panzer::IntrepidFieldPattern(basis_q1->getIntrepidBasis())));
    dofManager->addField(fieldName2_q1,Teuchos::rcp(new panzer::IntrepidFieldPattern(basis_q1->getIntrepidBasis())));
    dofManager->addField(fieldName_qedge1,Teuchos::rcp(new panzer::IntrepidFieldPattern(basis_qedge1->getIntrepidBasis())));

    std::vector<std::vector<std::string> > fieldOrder(3);
    fieldOrder[0].push_back(fieldName1_q1);
    fieldOrder[1].push_back(fieldName_qedge1);
    fieldOrder[2].push_back(fieldName2_q1);
    dofManager->setFieldOrder(fieldOrder);

    // dofManager->setOrientationsRequired(true);
    dofManager->buildGlobalUnknowns();

    // setup linear object factory
    /////////////////////////////////////////////////////////////

    Teuchos::RCP<BlockedEpetraLinearObjFactory<panzer::Traits,int> > be_lof 
       = Teuchos::rcp(new BlockedEpetraLinearObjFactory<panzer::Traits,int>(eComm.getConst(),dofManager));
    Teuchos::RCP<LinearObjFactory<panzer::Traits> > lof = be_lof;
    Teuchos::RCP<LinearObjContainer> dd_loc = be_lof->buildGhostedLinearObjContainer();
    Teuchos::RCP<LinearObjContainer> loc = be_lof->buildGhostedLinearObjContainer();
    be_lof->initializeGhostedContainer(LinearObjContainer::F,*dd_loc);
    dd_loc->initialize();

    be_lof->initializeGhostedContainer(LinearObjContainer::X | LinearObjContainer::F,*loc);
    loc->initialize();

    Teuchos::RCP<BlockedEpetraLinearObjContainer> b_dd_loc = Teuchos::rcp_dynamic_cast<BlockedEpetraLinearObjContainer>(dd_loc);
    Teuchos::RCP<BlockedEpetraLinearObjContainer> b_loc = Teuchos::rcp_dynamic_cast<BlockedEpetraLinearObjContainer>(loc);
    Teuchos::RCP<Thyra::ProductVectorBase<double> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(b_loc->get_x());
    Thyra::assign(p_vec->getNonconstVectorBlock(0).ptr(),123.0+myRank);
    Thyra::assign(p_vec->getNonconstVectorBlock(1).ptr(),456.0+myRank);
    Thyra::assign(p_vec->getNonconstVectorBlock(2).ptr(),789.0+myRank);

    // setup field manager, add evaluator under test
    /////////////////////////////////////////////////////////////
 
    PHX::FieldManager<panzer::Traits> fm;

    std::string resName = "";
    Teuchos::RCP<std::map<std::string,std::string> > names_map =
       Teuchos::rcp(new std::map<std::string,std::string>);
    names_map->insert(std::make_pair(fieldName1_q1,resName+fieldName1_q1));
    names_map->insert(std::make_pair(fieldName2_q1,resName+fieldName2_q1));
    names_map->insert(std::make_pair(fieldName_qedge1,resName+fieldName_qedge1));

    // evaluators under test
    {
       using Teuchos::RCP;
       using Teuchos::rcp;
       RCP<std::vector<std::string> > names = rcp(new std::vector<std::string>);
       names->push_back(resName+fieldName1_q1);
       names->push_back(resName+fieldName2_q1);

       Teuchos::ParameterList pl; 
       pl.set("Scatter Name", "ScatterQ1");
       pl.set("Basis", basis_q1);
       pl.set("Dependent Names", names);
       pl.set("Dependent Map", names_map);
//.........这里部分代码省略.........