C++ ParameterList::set方法代码示例

void
Albany::ResponseFactory::
createResponseFunction(
  const std::string& name,
  Teuchos::ParameterList& responseParams,
  Teuchos::Array< Teuchos::RCP<AbstractResponseFunction> >& responses) const
{
  using std::string;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;
  using Teuchos::Array;

  RCP<const Epetra_Comm> comm = app->getComm();

  if (name == "Solution Average") {
    responses.push_back(rcp(new Albany::SolutionAverageResponseFunction(comm)));
  }

  else if (name == "Solution Two Norm") {
    responses.push_back(rcp(new Albany::SolutionTwoNormResponseFunction(comm)));
  }

  else if (name == "Solution Values") {
    responses.push_back(rcp(new Albany::SolutionValuesResponseFunction(app, responseParams)));
  }

  else if (name == "Solution Max Value") {
    int eq = responseParams.get("Equation", 0);
    int neq = responseParams.get("Num Equations", 1);
    bool inor =  responseParams.get("Interleaved Ordering", true);

    responses.push_back(
      rcp(new Albany::SolutionMaxValueResponseFunction(comm, neq, eq, inor)));
  }

  else if (name == "Solution Two Norm File") {
    responses.push_back(
      rcp(new Albany::SolutionFileResponseFunction<Albany::NormTwo>(comm)));
  }

  else if (name == "Solution Inf Norm File") {
    responses.push_back(
      rcp(new Albany::SolutionFileResponseFunction<Albany::NormInf>(comm)));
  }

  else if (name == "Aggregated") {
    int num_aggregate_responses = responseParams.get<int>("Number");
    Array< RCP<AbstractResponseFunction> > aggregated_responses;
    Array< RCP<ScalarResponseFunction> > scalar_responses;
    for (int i=0; i<num_aggregate_responses; i++) {
      std::string id = Albany::strint("Response",i);
      std::string name = responseParams.get<std::string>(id);
      std::string sublist_name = Albany::strint("ResponseParams",i);
      createResponseFunction(name, responseParams.sublist(sublist_name),
			     aggregated_responses);

    }
    scalar_responses.resize(aggregated_responses.size());
    for (int i=0; i<aggregated_responses.size(); i++) {
      TEUCHOS_TEST_FOR_EXCEPTION(
	aggregated_responses[i]->isScalarResponse() != true, std::logic_error,
	"Response function " << i << " is not a scalar response function." <<
	std::endl <<
	"The aggregated response can only aggregate scalar response " <<
	"functions!");
      scalar_responses[i] =
	Teuchos::rcp_dynamic_cast<ScalarResponseFunction>(
	  aggregated_responses[i]);
    }
    responses.push_back(
      rcp(new Albany::AggregateScalarResponseFunction(comm, scalar_responses)));
  }

  else if (name == "Field Integral" ||
	   name == "Field Value" ||
	   name == "Field Average" ||
	   name == "Surface Velocity Mismatch" ||
	   name == "Center Of Mass" ||
	   name == "Save Field" ||
	   name == "Region Boundary" ||
	   name == "Element Size Field" ||
	   name == "IP to Nodal Field" ||
	   name == "Save Nodal Fields" ||
	   name == "PHAL Field Integral") {
    responseParams.set("Name", name);
    for (int i=0; i<meshSpecs.size(); i++) {
      responses.push_back(
	rcp(new Albany::FieldManagerScalarResponseFunction(
	      app, prob, meshSpecs[i], stateMgr, responseParams)));
    }
  }

  else if (name == "Solution") {
    responses.push_back(
      rcp(new Albany::SolutionResponseFunction(app, responseParams)));
  }

  else if (name == "KL") {
    Array< RCP<AbstractResponseFunction> > base_responses;
//.........这里部分代码省略.........

  TEUCHOS_UNIT_TEST(gather_orientation, gather_constr)
  {

    const std::size_t workset_size = 4;
    const std::string fieldName_q1 = "U";
    const std::string fieldName_qedge1 = "V";

    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::DOFManager<int,int> > dofManager = Teuchos::rcp(new panzer::DOFManager<int,int>(conn_manager,MPI_COMM_WORLD));
    dofManager->addField(fieldName_q1,Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis())));
    dofManager->addField(fieldName_qedge1,Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_qedge1->getIntrepid2Basis())));
    dofManager->setOrientationsRequired(true);
    dofManager->buildGlobalUnknowns();

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

    Teuchos::RCP<PHX::FieldTag> evalField_q1, evalField_qedge1;
    {
       Teuchos::RCP<std::vector<std::string> > dofNames = Teuchos::rcp(new std::vector<std::string>);
       dofNames->push_back(fieldName_q1);

       Teuchos::ParameterList pl;
       pl.set("Indexer Names",dofNames);
       pl.set("DOF Names",dofNames);
       pl.set("Basis",basis_q1);

       Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator  
          = Teuchos::rcp(new panzer::GatherOrientation<panzer::Traits::Residual,panzer::Traits,int,int>(dofManager,pl));

       TEST_EQUALITY(evaluator->evaluatedFields().size(),1);
       evalField_q1 = evaluator->evaluatedFields()[0];

       TEST_EQUALITY(evalField_q1->name(),basis_q1->name()+" Orientation");
       TEST_EQUALITY(evalField_q1->dataLayout().dimension(0),basis_q1->functional->dimension(0));
       TEST_EQUALITY(evalField_q1->dataLayout().dimension(1),basis_q1->functional->dimension(1));

       fm.registerEvaluator<panzer::Traits::Residual>(evaluator);
       fm.requireField<panzer::Traits::Residual>(*evaluator->evaluatedFields()[0]);
    }
    {
       Teuchos::RCP<std::vector<std::string> > dofNames = Teuchos::rcp(new std::vector<std::string>);
       dofNames->push_back(fieldName_qedge1);

       Teuchos::ParameterList pl;
       pl.set("Indexer Names",dofNames);
       pl.set("DOF Names",dofNames);
       pl.set("Basis",basis_qedge1);

       Teuchos::RCP<PHX::Evaluator<panzer::Traits> > evaluator  
          = Teuchos::rcp(new panzer::GatherOrientation<panzer::Traits::Residual,panzer::Traits,int,int>(dofManager,pl));

       TEST_EQUALITY(evaluator->evaluatedFields().size(),1);
       evalField_qedge1 = evaluator->evaluatedFields()[0];

       TEST_EQUALITY(evalField_qedge1->name(),basis_qedge1->name()+" Orientation");
       TEST_EQUALITY(evalField_qedge1->dataLayout().dimension(0),basis_qedge1->functional->dimension(0));
       TEST_EQUALITY(evalField_qedge1->dataLayout().dimension(1),basis_qedge1->functional->dimension(1));

       fm.registerEvaluator<panzer::Traits::Residual>(evaluator);
       fm.requireField<panzer::Traits::Residual>(*evaluator->evaluatedFields()[0]);
    }

    panzer::Traits::SetupData sd;
    fm.postRegistrationSetup(sd);

    // run tests
    /////////////////////////////////////////////////////////////

    panzer::Workset & workset = (*work_sets)[0];
    workset.alpha = 0.0;
    workset.beta = 0.0;
    workset.time = 0.0;
    workset.evaluate_transient_terms = false;

    fm.evaluateFields<panzer::Traits::Residual>(workset);

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

  void MLParameterListInterpreter<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::SetParameterList(const Teuchos::ParameterList & paramList_in) {
    Teuchos::ParameterList paramList = paramList_in;

    RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); // TODO: use internal out (GetOStream())

    //
    // Read top-level of the parameter list
    //

    // hard-coded default values == ML defaults according to the manual
    MUELU_READ_PARAM(paramList, "ML output",                                int,                   0,       verbosityLevel);
    MUELU_READ_PARAM(paramList, "max levels",                               int,                  10,       maxLevels);
    MUELU_READ_PARAM(paramList, "PDE equations",                            int,                   1,       nDofsPerNode);

    MUELU_READ_PARAM(paramList, "coarse: max size",                         int,                 128,       maxCoarseSize);

    MUELU_READ_PARAM(paramList, "aggregation: type",                std::string,         "Uncoupled",       agg_type);
    //MUELU_READ_PARAM(paramList, "aggregation: threshold",                double,                 0.0,       agg_threshold);
    MUELU_READ_PARAM(paramList, "aggregation: damping factor",           double, (double)4/(double)3,       agg_damping);
    //MUELU_READ_PARAM(paramList, "aggregation: smoothing sweeps",            int,                   1,       agg_smoothingsweeps);
    MUELU_READ_PARAM(paramList, "aggregation: nodes per aggregate",         int,                   1,       minPerAgg);
    MUELU_READ_PARAM(paramList, "aggregation: keep Dirichlet bcs",         bool,               false,       bKeepDirichletBcs); // This is a MueLu specific extension that does not exist in ML
    MUELU_READ_PARAM(paramList, "aggregation: max neighbours already aggregated", int,             0,       maxNbrAlreadySelected); // This is a MueLu specific extension that does not exist in ML

    MUELU_READ_PARAM(paramList, "null space: type",                 std::string,   "default vectors",       nullspaceType);
    MUELU_READ_PARAM(paramList, "null space: dimension",                    int,                  -1,       nullspaceDim); // TODO: ML default not in documentation
    MUELU_READ_PARAM(paramList, "null space: vectors",                   double*,               NULL,       nullspaceVec); // TODO: ML default not in documentation

    MUELU_READ_PARAM(paramList, "energy minimization: enable",             bool,               false,       bEnergyMinimization);

    MUELU_READ_PARAM(paramList, "RAP: fix diagonal",                       bool,               false,       bFixDiagonal); // This is a MueLu specific extension that does not exist in ML

    //
    // Move smoothers/aggregation/coarse parameters to sublists
    //

    // ML allows to have level-specific smoothers/aggregation/coarse parameters at the top level of the list or/and defined in sublists:
    // See also: ML Guide section 6.4.1, MueLu::CreateSublists, ML_CreateSublists
    ParameterList paramListWithSubList;
    MueLu::CreateSublists(paramList, paramListWithSubList);
    paramList = paramListWithSubList; // swap

    // std::cout << std::endl << "Parameter list after CreateSublists" << std::endl;
    // std::cout << paramListWithSubList << std::endl;

    //
    // Validate parameter list
    //

    {
      bool validate = paramList.get("ML validate parameter list", true); /* true = default in ML */
      if (validate) {

#if defined(HAVE_MUELU_ML) && defined(HAVE_MUELU_EPETRA)
        // Validate parameter list using ML validator
        int depth = paramList.get("ML validate depth", 5); /* 5 = default in ML */
        TEUCHOS_TEST_FOR_EXCEPTION(! ML_Epetra::ValidateMLPParameters(paramList, depth), Exceptions::RuntimeError,
                                   "ERROR: ML's Teuchos::ParameterList contains incorrect parameter!");
#else
        // If no validator available: issue a warning and set parameter value to false in the output list
        *out << "Warning: MueLu_ENABLE_ML=OFF. The parameter list cannot be validated." << std::endl;
        paramList.set("ML validate parameter list", false);

#endif // HAVE_MUELU_ML
      } // if(validate)
    } // scope


    //
    //
    //



    // Matrix option
    blksize_ = nDofsPerNode;

    // Translate verbosity parameter

    // Translate verbosity parameter
    MsgType eVerbLevel = None;
    if (verbosityLevel == 0) eVerbLevel = None;
    if (verbosityLevel >  0) eVerbLevel = Low;
    if (verbosityLevel >  4) eVerbLevel = Medium;
    if (verbosityLevel >  7) eVerbLevel = High;
    if (verbosityLevel >  9) eVerbLevel = Extreme;
    if (verbosityLevel >  9) eVerbLevel = Test;
    this->verbosity_ = eVerbLevel;


    TEUCHOS_TEST_FOR_EXCEPTION(agg_type != "Uncoupled" && agg_type != "Coupled", Exceptions::RuntimeError, "MueLu::MLParameterListInterpreter::Setup(): parameter \"aggregation: type\": only 'Uncoupled' or 'Coupled' aggregation is supported.");

    // Create MueLu factories
    RCP<CoalesceDropFactory> dropFact = rcp(new CoalesceDropFactory());

    RCP<FactoryBase> CoupledAggFact = Teuchos::null;
    if(agg_type == "Uncoupled") {
      // Uncoupled aggregation
      RCP<UncoupledAggregationFactory> CoupledAggFact2 = rcp(new UncoupledAggregationFactory());
      /*CoupledAggFact2->SetMinNodesPerAggregate(minPerAgg); //TODO should increase if run anything other than 1D
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifndef HAVE_TPETRA_COMPLEX_DOUBLE
#  error "Anasazi: This test requires Scalar = std::complex<double> to be enabled in Tpetra."
#else
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::tuple;
  using std::cout;
  using std::endl;

  typedef std::complex<double>                ST;
  typedef Teuchos::ScalarTraits<ST>          SCT;
  typedef SCT::magnitudeType                  MT;
  typedef Tpetra::MultiVector<ST>             MV;
  typedef MV::global_ordinal_type             GO;
  typedef Tpetra::Operator<ST>                OP;
  typedef Anasazi::MultiVecTraits<ST,MV>     MVT;
  typedef Anasazi::OperatorTraits<ST,MV,OP>  OPT;

  Tpetra::ScopeGuard tpetraScope (&argc, &argv);

  bool success = false;

  const ST ONE = SCT::one ();

  int info = 0;

  RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm ();

  const int MyPID = comm->getRank ();

  bool verbose = false;
  bool debug = false;
  bool insitu = false;
  bool herm = false;
  std::string which("LM");
  std::string filename;
  int nev = 4;
  int blockSize = 4;
  MT tol = 1.0e-6;

  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
  cmdp.setOption("insitu","exsitu",&insitu,"Perform in situ restarting.");
  cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
  cmdp.setOption("herm","nonherm",&herm,"Solve Hermitian or non-Hermitian problem.");
  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix (assumes non-Hermitian unless specified otherwise).");
  cmdp.setOption("nev",&nev,"Number of eigenvalues to compute.");
  cmdp.setOption("blockSize",&blockSize,"Block size for the algorithm.");
  cmdp.setOption("tol",&tol,"Tolerance for convergence.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }
  if (debug) verbose = true;
  if (filename == "") {
    // get default based on herm
    if (herm) {
      filename = "mhd1280b.cua";
    }
    else {
      filename = "mhd1280a.cua";
    }
  }

  if (MyPID == 0) {
    cout << Anasazi::Anasazi_Version() << endl << endl;
  }

  // Get the data from the HB file
  int dim,dim2,nnz;
  int rnnzmax;
  double *dvals;
  int *colptr,*rowind;
  nnz = -1;
  if (MyPID == 0) {
    info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,&colptr,&rowind,&dvals);
    // find maximum NNZ over all rows
    vector<int> rnnz(dim,0);
    for (int *ri=rowind; ri<rowind+nnz; ++ri) {
      ++rnnz[*ri-1];
    }
    rnnzmax = *std::max_element(rnnz.begin(),rnnz.end());
  }
  else {
    // address uninitialized data warnings
    dvals = NULL;
    colptr = NULL;
    rowind = NULL;
  }
  Teuchos::broadcast(*comm,0,&info);
  Teuchos::broadcast(*comm,0,&nnz);
  Teuchos::broadcast(*comm,0,&dim);
  Teuchos::broadcast(*comm,0,&rnnzmax);
  if (info == 0 || nnz < 0) {
    if (MyPID == 0) {
      cout << "Error reading '" << filename << "'" << endl
           << "End Result: TEST FAILED" << endl;
    }
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
#include "MueLu_UseShortNames.hpp"

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

  //
  // MPI initialization
  //

  Teuchos::oblackholestream blackhole;
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, &blackhole);

  bool success = false;
  bool verbose = true;
  try {
    RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();

    //
    // Process command line arguments
    //
    Teuchos::CommandLineProcessor  clp(false);
    Galeri::Xpetra::Parameters<GO> matrixParameters(clp, 81); // manage parameters of the test case
    Xpetra::Parameters             xpetraParameters(clp);     // manage parameters of xpetra

    switch (clp.parse(argc,argv)) {
      case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:        return EXIT_SUCCESS;
      case Teuchos::CommandLineProcessor::PARSE_ERROR:
      case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE;
      case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:          break;
      default:;
    }

    if (comm->getRank() == 0) std::cout << xpetraParameters << matrixParameters;

    //
    // Setup test case (Ax = b)
    //

    // Linear Algebra Library
    Xpetra::UnderlyingLib lib = xpetraParameters.GetLib();

    // Distribution
    RCP<const Map> map = MapFactory::Build(lib, matrixParameters.GetNumGlobalElements(), 0, comm);

    // Matrix
    RCP<Galeri::Xpetra::Problem<Map,CrsMatrixWrap,MultiVector> > Pr =
      Galeri::Xpetra::BuildProblem<SC, LO, GO, Map, CrsMatrixWrap, MultiVector>(matrixParameters.GetMatrixType(), map, matrixParameters.GetParameterList());
    RCP<Matrix> A = Pr->BuildMatrix();

    // User defined nullspace
    RCP<MultiVector> nullSpace = VectorFactory::Build(map,1); nullSpace->putScalar((SC) 1.0);

    // Define B
    RCP<Vector> X = VectorFactory::Build(map,1);
    RCP<Vector> B = VectorFactory::Build(map,1);
    X->setSeed(846930886);
    X->randomize();
    A->apply(*X, *B, Teuchos::NO_TRANS, (SC)1.0, (SC)0.0);

    // X = 0
    X->putScalar((SC) 0.0);

    //
    // Create a multigrid configuration
    //

    // Transfer operators
    RCP<TentativePFactory> TentativePFact = rcp( new TentativePFactory() );
    RCP<SaPFactory>        SaPFact        = rcp( new SaPFactory() );
    RCP<TransPFactory>     RFact          = rcp( new TransPFactory());

    FactoryManager M;
    M.SetFactory("Ptent", TentativePFact);
    M.SetFactory("P",     SaPFact);
    M.SetFactory("R",     RFact);

    M.SetFactory("Smoother", Teuchos::null);      //skips smoother setup
    M.SetFactory("CoarseSolver", Teuchos::null);  //skips coarsest solve setup

    //
    // Multigrid setup phase
    //

    int startLevel = 0;
    int maxLevels = 10;

    std::cout << "=============== Setup transfers only ====================" << std::endl;

    Hierarchy H;
    H.SetDefaultVerbLevel(MueLu::Medium);

    RCP<Level> finestLevel = H.GetLevel();
    finestLevel->Set("A", A);
    finestLevel->Set("Nullspace", nullSpace);

    // Indicate which Hierarchy operators we want to keep
    H.Keep("P", SaPFact.get());  //SaPFact is the generating factory for P.
    H.Keep("R", RFact.get());    //RFact is the generating factory for R.
    H.Keep("Ptent", TentativePFact.get());  //SaPFact is the generating factory for P.
//.........这里部分代码省略.........

bool CompareWithAztecOO(Epetra_LinearProblem& Problem, const std::string what,
                       int Overlap, int ival)
{
  using std::cout;
  using std::endl;

  AztecOO AztecOOSolver(Problem);
  AztecOOSolver.SetAztecOption(AZ_solver,AZ_gmres);
  AztecOOSolver.SetAztecOption(AZ_output,AZ_none);
  AztecOOSolver.SetAztecOption(AZ_overlap,Overlap);
  AztecOOSolver.SetAztecOption(AZ_graph_fill,ival);
  AztecOOSolver.SetAztecOption(AZ_poly_ord, ival);
  AztecOOSolver.SetAztecParam(AZ_drop, 0.0);
  AztecOOSolver.SetAztecParam(AZ_athresh, 0.0);
  AztecOOSolver.SetAztecParam(AZ_rthresh, 0.0);

  Epetra_MultiVector& RHS = *(Problem.GetRHS());
  Epetra_MultiVector& LHS = *(Problem.GetLHS());
  Teuchos::RefCountPtr<Epetra_RowMatrix> A = Teuchos::rcp(Problem.GetMatrix(), false);

  LHS.Random();
  A->Multiply(false,LHS,RHS);

  Teuchos::ParameterList List;
  List.set("fact: level-of-fill", ival);
  List.set("relaxation: sweeps", ival);
  List.set("relaxation: damping factor", 1.0);
  List.set("relaxation: zero starting solution", true);

  //default combine mode is as for AztecOO
  List.set("schwarz: combine mode", Zero);

  Epetra_Time Time(A->Comm());

  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec;

  if (what == "Jacobi") {
    Prec = Teuchos::rcp( new Ifpack_PointRelaxation(&*A) );
    List.set("relaxation: type", "Jacobi");
    AztecOOSolver.SetAztecOption(AZ_precond,AZ_Jacobi);
    AztecOOSolver.SetAztecOption(AZ_reorder,0);
  }
  else if (what == "IC no reord") {
    Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_IC>(&*A,Overlap) );
    AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
    AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_icc);
    AztecOOSolver.SetAztecOption(AZ_reorder,0);
  }
  else if (what == "IC reord") {
    Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_IC>(&*A,Overlap) );
    List.set("schwarz: use reordering", true);
    AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
    AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_icc);
    AztecOOSolver.SetAztecOption(AZ_reorder,1);
  }
  else if (what == "ILU no reord") {
    Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_ILU>(&*A,Overlap) );
    AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
    AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_ilu);
    AztecOOSolver.SetAztecOption(AZ_reorder,0);
  }
  else if (what == "ILU reord") {
    Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_ILU>(&*A,Overlap) );
    List.set("schwarz: use reordering", true);
    AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
    AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_ilu);
    AztecOOSolver.SetAztecOption(AZ_reorder,1);
  }
#ifdef HAVE_IFPACK_AMESOS
  else if (what == "LU") {
    Prec = Teuchos::rcp( new Ifpack_AdditiveSchwarz<Ifpack_Amesos>(&*A,Overlap) );
    List.set("amesos: solver type", "Klu");
    AztecOOSolver.SetAztecOption(AZ_precond,AZ_dom_decomp);
    AztecOOSolver.SetAztecOption(AZ_subdomain_solve,AZ_lu);
  }
#endif
  else {
    cerr << "Option not recognized" << endl;
    exit(EXIT_FAILURE);
  }

  // ==================================== //
  // Solve with AztecOO's preconditioners //
  // ==================================== //

  LHS.PutScalar(0.0);

  Time.ResetStartTime();
  AztecOOSolver.Iterate(150,1e-5);

  if (verbose) {
    cout << endl;
    cout << "==================================================" << endl;
    cout << "Testing `" << what << "', Overlap = "
         << Overlap << ", ival = " << ival << endl;
    cout << endl;
    cout << "[AztecOO] Total time = " << Time.ElapsedTime() << " (s)" << endl;
    cout << "[AztecOO] Residual   = " << AztecOOSolver.TrueResidual() << " (s)" << endl;
    cout << "[AztecOO] Iterations = " << AztecOOSolver.NumIters() << endl;
    cout << endl;
//.........这里部分代码省略.........

// The solve is done in the felix_driver_run function, and the solution is passed back to Glimmer-CISM 
// IK, 12/3/13: time_inc_yr and cur_time_yr are not used here... 
void felix_driver_run(FelixToGlimmer * ftg_ptr, double& cur_time_yr, double time_inc_yr)
{
    //IK, 12/9/13: how come FancyOStream prints an all processors??    
    Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());

    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) {
      std::cout << "In felix_driver_run, cur_time, time_inc = " << cur_time_yr 
                << "   " << time_inc_yr << std::endl;
    }
    
    // ---------------------------------------------
    // get u and v velocity solution from Glimmer-CISM 
    // IK, 11/26/13: need to concatenate these into a single solve for initial condition for Albany/FELIX solve 
    // IK, 3/14/14: moved this step to felix_driver_run from felix_driver init, since we still want to grab and u and v velocities for CISM if the mesh hasn't changed, 
    // in which case only felix_driver_run will be called, not felix_driver_init.   
    // ---------------------------------------------
    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) 
      std::cout << "In felix_driver_run: grabbing pointers to u and v velocities in CISM..." << std::endl; 
    uVel_ptr = ftg_ptr ->getDoubleVar("uvel", "velocity"); 
    vVel_ptr = ftg_ptr ->getDoubleVar("vvel", "velocity"); 

    // ---------------------------------------------
    // Set restart solution to the one passed from CISM
    // IK, 3/14/14: moved this from felix_driver_init to felix_driver_run.  
    // ---------------------------------------------
    
    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) 
      std::cout << "In felix_driver_run: setting initial condition from CISM..." << std::endl;
    //Check what kind of ordering you have in the solution & create solutionField object.
    interleavedOrdering = meshStruct->getInterleavedOrdering();
    Albany::AbstractSTKFieldContainer::VectorFieldType* solutionField;
    if(interleavedOrdering)
      solutionField = Teuchos::rcp_dynamic_cast<Albany::OrdinarySTKFieldContainer<true> >(meshStruct->getFieldContainer())->getSolutionField();
    else
      solutionField = Teuchos::rcp_dynamic_cast<Albany::OrdinarySTKFieldContainer<false> >(meshStruct->getFieldContainer())->getSolutionField();

     //Create vector used to renumber nodes on each processor from the Albany convention (horizontal levels first) to the CISM convention (vertical layers first)
     nNodes2D = (global_ewn + 1)*(global_nsn+1); //number global nodes in the domain in 2D 
     nNodesProc2D = (nsn-2*nhalo+1)*(ewn-2*nhalo+1); //number of nodes on each processor in 2D  
     cismToAlbanyNodeNumberMap.resize(upn*nNodesProc2D);
     for (int j=0; j<nsn-2*nhalo+1;j++) { 
       for (int i=0; i<ewn-2*nhalo+1; i++) {
         for (int k=0; k<upn; k++) { 
           int index = k+upn*i + j*(ewn-2*nhalo+1)*upn; 
           cismToAlbanyNodeNumberMap[index] = k*nNodes2D + global_node_id_owned_map_Ptr[i+j*(ewn-2*nhalo+1)]; 
           //if (mpiComm->MyPID() == 0) 
           //  std::cout << "index: " << index << ", cismToAlbanyNodeNumberMap: " << cismToAlbanyNodeNumberMap[index] << std::endl; 
          }
        }
      }

     //The way it worked out, uVel_ptr and vVel_ptr have more nodes than the nodes in the mesh passed to Albany/CISM for the solve.  In particular, 
     //there is 1 row of halo elements in uVel_ptr and vVel_ptr.  To account for this, we copy uVel_ptr and vVel_ptr into std::vectors, which do not have the halo elements. 
     std::vector<double> uvel_vec(upn*nNodesProc2D); 
     std::vector<double> vvel_vec(upn*nNodesProc2D); 
     int counter1 = 0; 
     int counter2 = 0; 
     int local_nodeID; 
     for (int j=0; j<nsn-1; j++) {
       for (int i=0; i<ewn-1; i++) { 
         for (int k=0; k<upn; k++) {
           if (j >= nhalo-1 & j < nsn-nhalo) {
             if (i >= nhalo-1 & i < ewn-nhalo) {
#ifdef CISM_USE_EPETRA 
               local_nodeID = node_map->LID(cismToAlbanyNodeNumberMap[counter1]); 
#else
               local_nodeID = node_map->getLocalElement(cismToAlbanyNodeNumberMap[counter1]);
#endif
               uvel_vec[counter1] = uVel_ptr[counter2]; 
               vvel_vec[counter1] = vVel_ptr[counter2]; 
               counter1++;
            }
            }
            counter2++; 
         }
        }
     }
     //Loop over all the elements to find which nodes are active.  For the active nodes, copy uvel and vvel from CISM into Albany solution array to 
     //use as initial condition.
     //NOTE: there is some inefficiency here by looping over all the elements.  TO DO? pass only active nodes from Albany-CISM to improve this? 
     double velScale = seconds_per_year*vel_scaling_param;  
     for (int i=0; i<nElementsActive; i++) {
       for (int j=0; j<8; j++) {
        int node_GID =  global_element_conn_active_Ptr[i + nElementsActive*j]; //node_GID is 1-based
#ifdef CISM_USE_EPETRA      
        int node_LID =  node_map->LID(node_GID); //node_LID is 0-based
#else
        int node_LID =  node_map->getLocalElement(node_GID); //node_LID is 0-based
#endif
        stk::mesh::Entity node = meshStruct->bulkData->get_entity(stk::topology::NODE_RANK, node_GID);
        double* sol = stk::mesh::field_data(*solutionField, node);
        //IK, 3/18/14: added division by velScale to convert uvel and vvel from dimensionless to having units of m/year (the Albany units)  
        sol[0] = uvel_vec[node_LID]/velScale;
        sol[1] = vvel_vec[node_LID]/velScale;
      }
    }
    // ---------------------------------------------------------------------------------------------------
    // Solve 
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  // MEMORY_CHECK(true, "Before initializing MPI");

  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
  RCP<const Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
  int rank = comm->getRank();
  int nprocs = comm->getSize();

  MEMORY_CHECK(rank==0, "After initializing MPI");

  if (rank==0)
    cout << "Number of processes: " << nprocs << endl;

  // Default values
  double numGlobalCoords = 1000;
  int numTestCuts = 1;
  int nWeights = 0;
  string timingType("no_timers");
  string debugLevel("basic_status");
  string memoryOn("memoryOn");
  string memoryOff("memoryOff");
  string memoryProcs("0");
  bool doMemory=false;
  int numGlobalParts = nprocs;

  CommandLineProcessor commandLine(false, true);
  commandLine.setOption("size", &numGlobalCoords, 
    "Approximate number of global coordinates.");
  commandLine.setOption("testCuts", &numTestCuts, 
    "Number of test cuts to make when looking for bisector.");
  commandLine.setOption("numParts", &numGlobalParts, 
    "Number of parts (default is one per proc).");
  commandLine.setOption("nWeights", &nWeights, 
    "Number of weights per coordinate, zero implies uniform weights.");

  string balanceCount("balance_object_count");
  string balanceWeight("balance_object_weight");
  string mcnorm1("multicriteria_minimize_total_weight");
  string mcnorm2("multicriteria_balance_total_maximum");
  string mcnorm3("multicriteria_minimize_maximum_weight");

  string objective(balanceWeight);   // default

  string doc(balanceCount); doc.append(": ignore weights\n");

  doc.append(balanceWeight); doc.append(": balance on first weight\n");

  doc.append(mcnorm1);
  doc.append(": given multiple weights, balance their total.\n");

  doc.append(mcnorm3);
  doc.append(": given multiple weights, balance the maximum for each coordinate.\n");

  doc.append(mcnorm2);
  doc.append(": given multiple weights, balance the L2 norm of the weights.\n");

  commandLine.setOption("objective", &objective,  doc.c_str());

  commandLine.setOption("timers", &timingType,
    "no_timers, micro_timers, macro_timers, both_timers, test_timers");

  commandLine.setOption("debug", &debugLevel,
   "no_status, basic_status, detailed_status, verbose_detailed_status");

  commandLine.setOption(memoryOn.c_str(), memoryOff.c_str(), &doMemory,
    "do memory profiling");

  commandLine.setOption("memoryProcs", &memoryProcs,
   "list of processes that output memory usage");

  CommandLineProcessor::EParseCommandLineReturn rc = 
    commandLine.parse(argc, argv);

  if (rc != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL){
    if (rc == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED){
      if (rank==0)
        cout << "PASS" << endl;
      return 1;
    }
    else{
      if (rank==0)
        cout << "FAIL" << endl;
      return 0;
    }
  }

  //MEMORY_CHECK(doMemory && rank==0, "After processing parameters");

  zgno_t globalSize = static_cast<zgno_t>(numGlobalCoords);

  RCP<tMVector_t> coordinates = getMeshCoordinates(comm, globalSize);
  size_t numLocalCoords = coordinates->getLocalLength();

#if 0
  comm->barrier();
  for (int p=0; p < nprocs; p++){
    if (p==rank){
      cout << "Rank " << rank << ", " << numLocalCoords << "coords" << endl;
      const zscalar_t *x = coordinates->getData(0).getRawPtr();
//.........这里部分代码省略.........

  const RCP<const FactoryBase> FactoryManager<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetDefaultFactory(const std::string& varName) const {
    if (defaultFactoryTable_.count(varName)) {
      // The factory for this name was already created (possibly, for previous level, if we reuse factory manager)
      return defaultFactoryTable_.find(varName)->second;

    } else {
      // No factory was created for this name, but we may know which one to create
      if (varName == "A")                               return SetAndReturnDefaultFactory(varName, rcp(new RAPFactory()));
      if (varName == "RAP Pattern")                     return GetFactory("A");
      if (varName == "AP Pattern")                      return GetFactory("A");
      if (varName == "Ptent")                           return SetAndReturnDefaultFactory(varName, rcp(new TentativePFactory()));
      if (varName == "P") {
        // GetFactory("Ptent"): we need to use the same factory instance for both "P" and "Nullspace"
        RCP<Factory> factory = rcp(new SaPFactory());
        factory->SetFactory("P", GetFactory("Ptent"));
        return SetAndReturnDefaultFactory(varName, factory);
      }
      if (varName == "Nullspace") {
        // GetFactory("Ptent"): we need to use the same factory instance for both "P" and "Nullspace"
        RCP<Factory> factory = rcp(new NullspaceFactory());
        factory->SetFactory("Nullspace", GetFactory("Ptent"));
        return SetAndReturnDefaultFactory(varName, factory);
      }

      if (varName == "R")                               return SetAndReturnDefaultFactory(varName, rcp(new TransPFactory()));
#if defined(HAVE_MUELU_ZOLTAN) && defined(HAVE_MPI)
      if (varName == "Partition")                       return SetAndReturnDefaultFactory(varName, rcp(new ZoltanInterface()));
#endif //ifdef HAVE_MPI

      if (varName == "Importer") {
#ifdef HAVE_MPI
                                                        return SetAndReturnDefaultFactory(varName, rcp(new RepartitionFactory()));
#else
                                                        return SetAndReturnDefaultFactory(varName, NoFactory::getRCP());
#endif
      }

      if (varName == "Graph")                           return SetAndReturnDefaultFactory(varName, rcp(new CoalesceDropFactory()));
      if (varName == "UnAmalgamationInfo")              return SetAndReturnDefaultFactory(varName, rcp(new AmalgamationFactory())); //GetFactory("Graph"));
      if (varName == "Aggregates")                      return SetAndReturnDefaultFactory(varName, rcp(new UncoupledAggregationFactory()));
      if (varName == "CoarseMap")                       return SetAndReturnDefaultFactory(varName, rcp(new CoarseMapFactory()));
      if (varName == "DofsPerNode")                     return GetFactory("Graph");
      if (varName == "Filtering")                       return GetFactory("Graph");
      if (varName == "LineDetection_VertLineIds")       return SetAndReturnDefaultFactory(varName, rcp(new LineDetectionFactory()));
      if (varName == "LineDetection_Layers")            return GetFactory("LineDetection_VertLineIds");
      if (varName == "CoarseNumZLayers")                return GetFactory("LineDetection_VertLineIds");

      // Same factory for both Pre and Post Smoother. Factory for key "Smoother" can be set by users.
      if (varName == "PreSmoother")                     return GetFactory("Smoother");
      if (varName == "PostSmoother")                    return GetFactory("Smoother");

      if (varName == "Ppattern") {
        RCP<PatternFactory> PpFact = rcp(new PatternFactory);
        PpFact->SetFactory("P", GetFactory("Ptent"));
        return SetAndReturnDefaultFactory(varName, PpFact);
      }
      if (varName == "Constraint")                      return SetAndReturnDefaultFactory(varName, rcp(new ConstraintFactory()));

      if (varName == "Smoother") {
        Teuchos::ParameterList smootherParamList;
        smootherParamList.set("relaxation: type",           "Symmetric Gauss-Seidel");
        smootherParamList.set("relaxation: sweeps",         Teuchos::OrdinalTraits<LO>::one());
        smootherParamList.set("relaxation: damping factor", Teuchos::ScalarTraits<Scalar>::one());
        return SetAndReturnDefaultFactory(varName, rcp(new SmootherFactory(rcp(new TrilinosSmoother("RELAXATION", smootherParamList)))));
      }
      if (varName == "CoarseSolver")                    return SetAndReturnDefaultFactory(varName, rcp(new SmootherFactory(rcp(new DirectSolver()), Teuchos::null)));

      TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "MueLu::FactoryManager::GetDefaultFactory(): No default factory available for building '" + varName + "'.");
    }
  }

int main(int argc, char *argv[]) {
#include <MueLu_UseShortNames.hpp>

  using Teuchos::RCP; // reference count pointers
  using Teuchos::rcp;
  using Teuchos::TimeMonitor;

  // =========================================================================
  // MPI initialization using Teuchos
  // =========================================================================
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
  RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
  int MyPID   = comm->getRank();
  int NumProc = comm->getSize();

  const Teuchos::RCP<Epetra_Comm> epComm = Teuchos::rcp_const_cast<Epetra_Comm>(Xpetra::toEpetra(comm));

  // =========================================================================
  // Convenient definitions
  // =========================================================================
  //SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();

  // Instead of checking each time for rank, create a rank 0 stream
  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
  Teuchos::FancyOStream& fancyout = *fancy;
  fancyout.setOutputToRootOnly(0);



  // =========================================================================
  // Parameters initialization
  // =========================================================================
  Teuchos::CommandLineProcessor clp(false);
  GO nx = 100, ny = 100;
  clp.setOption("nx",                   &nx,              "mesh size in x direction");
  clp.setOption("ny",                   &ny,              "mesh size in y direction");
  std::string xmlFileName = "xml/s3a.xml"; clp.setOption("xml", &xmlFileName,     "read parameters from a file. Otherwise, this example uses by default 'tutorial1a.xml'");
  int mgridSweeps = 1; clp.setOption("mgridSweeps", &mgridSweeps, "number of multigrid sweeps within Multigrid solver.");
  std::string printTimings = "no";   clp.setOption("timings", &printTimings,     "print timings to screen [yes/no]");
  double tol               = 1e-12;  clp.setOption("tol",                   &tol,              "solver convergence tolerance");

  switch (clp.parse(argc,argv)) {
    case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:        return EXIT_SUCCESS; break;
    case Teuchos::CommandLineProcessor::PARSE_ERROR:
    case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; break;
    case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:                               break;
  }

  // =========================================================================
  // Problem construction
  // =========================================================================
  RCP<TimeMonitor> globalTimeMonitor = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time"))), tm;

  comm->barrier();
  tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));

  Teuchos::ParameterList GaleriList;
  GaleriList.set("nx", nx);
  GaleriList.set("ny", ny);
  GaleriList.set("mx", epComm->NumProc());
  GaleriList.set("my", 1);
  GaleriList.set("lx", 1.0); // length of x-axis
  GaleriList.set("ly", 1.0); // length of y-axis
  GaleriList.set("diff", 1e-5);
  GaleriList.set("conv", 1.0);

  // create map
  Teuchos::RCP<Epetra_Map> epMap = Teuchos::rcp(Galeri::CreateMap("Cartesian2D", *epComm, GaleriList));

  // create coordinates
  Teuchos::RCP<Epetra_MultiVector> epCoord = Teuchos::rcp(Galeri::CreateCartesianCoordinates("2D", epMap.get(), GaleriList));

  // create matrix
  Teuchos::RCP<Epetra_CrsMatrix> epA = Teuchos::rcp(Galeri::CreateCrsMatrix("Recirc2D", epMap.get(), GaleriList));

  // Epetra -> Xpetra
  Teuchos::RCP<CrsMatrix> exA = Teuchos::rcp(new Xpetra::EpetraCrsMatrix(epA));
  Teuchos::RCP<CrsMatrixWrap> exAWrap = Teuchos::rcp(new CrsMatrixWrap(exA));

  RCP<Matrix> A = Teuchos::rcp_dynamic_cast<Matrix>(exAWrap);
  int numPDEs = 1;
  A->SetFixedBlockSize(numPDEs);

  // set rhs and solution vector
  RCP<Epetra_Vector> B = Teuchos::rcp(new Epetra_Vector(*epMap));
  RCP<Epetra_Vector> X = Teuchos::rcp(new Epetra_Vector(*epMap));
  B->PutScalar(1.0);
  X->PutScalar(0.0);

  // Epetra -> Xpetra
  RCP<Vector> xB = Teuchos::rcp(new Xpetra::EpetraVector(B));
  RCP<Vector> xX = Teuchos::rcp(new Xpetra::EpetraVector(X));

  xX->setSeed(100);
  xX->randomize();

  // build null space vector
  RCP<const Map> map = A->getRowMap();
  RCP<MultiVector> nullspace = MultiVectorFactory::Build(map, numPDEs);

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

  void testInitialzation(const Teuchos::RCP<Teuchos::ParameterList>& ipb,
			 std::vector<panzer::BC>& bcs)
  {
    // Physics block
    Teuchos::ParameterList& physics_block = ipb->sublist("test physics");
    {
      Teuchos::ParameterList& p = physics_block.sublist("a");
      p.set("Type","Energy");
      p.set("Prefix","");
      p.set("Model ID","solid");
      p.set("Basis Type","HGrad");
      p.set("Basis Order",2);
    }
    {
      Teuchos::ParameterList& p = physics_block.sublist("b");
      p.set("Type","Energy");
      p.set("Prefix","ION_");
      p.set("Model ID","solid");
      p.set("Basis Type","HGrad");
      p.set("Basis Order",1);
    }

    {
      std::size_t bc_id = 0;
      panzer::BCType neumann = BCT_Dirichlet;
      std::string sideset_id = "left";
      std::string element_block_id = "eblock-0_0";
      std::string dof_name = "TEMPERATURE";
      std::string strategy = "Constant";
      double value = 5.0;
      Teuchos::ParameterList p;
      p.set("Value",value);
      panzer::BC bc(bc_id, neumann, sideset_id, element_block_id, dof_name, 
		    strategy, p);
      bcs.push_back(bc);
    } 
    {
      std::size_t bc_id = 1;
      panzer::BCType neumann = BCT_Dirichlet;
      std::string sideset_id = "right";
      std::string element_block_id = "eblock-1_0";
      std::string dof_name = "TEMPERATURE";
      std::string strategy = "Constant";
      double value = 5.0;
      Teuchos::ParameterList p;
      p.set("Value",value);
      panzer::BC bc(bc_id, neumann, sideset_id, element_block_id, dof_name, 
		    strategy, p);
      bcs.push_back(bc);
    }   
    {
      std::size_t bc_id = 2;
      panzer::BCType neumann = BCT_Dirichlet;
      std::string sideset_id = "top";
      std::string element_block_id = "eblock-1_0";
      std::string dof_name = "TEMPERATURE";
      std::string strategy = "Constant";
      double value = 5.0;
      Teuchos::ParameterList p;
      p.set("Value",value);
      panzer::BC bc(bc_id, neumann, sideset_id, element_block_id, dof_name, 
		    strategy, p);
      bcs.push_back(bc);
    }

  }

  TEUCHOS_UNIT_TEST(scatter_field_evaluators, cell_field)
  {

    const std::size_t workset_size = 5;
    linBasis = buildLinearBasis(workset_size);

    Teuchos::RCP<panzer_stk::STK_Interface> mesh = buildMesh(5,5,false);

    RCP<Epetra_Comm> Comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));

     Teuchos::RCP<shards::CellTopology> topo = 
        Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));
    panzer::CellData cellData(workset_size,topo);
    Teuchos::RCP<panzer::IntegrationRule> intRule = Teuchos::rcp(new panzer::IntegrationRule(1,cellData));

    Teuchos::RCP<Teuchos::ParameterList> ipb = Teuchos::parameterList("Physics Blocks");
    std::vector<panzer::BC> bcs;
    testInitialzation(ipb, bcs);

    Teuchos::RCP<PHX::FieldManager<panzer::Traits> > fm = 
      Teuchos::rcp(new PHX::FieldManager<panzer::Traits>);
    {
       Teuchos::ParameterList pl;
       pl.set("Data Layout",linBasis->functional);
       fm->registerEvaluator<panzer::Traits::Residual>(Teuchos::rcp(new XCoordinate<panzer::Traits::Residual,panzer::Traits>(pl)));
    }

    {
       Teuchos::ParameterList pl;
       pl.set("Nodes",1);
       pl.set("Data Layout",intRule->dl_scalar);
       fm->registerEvaluator<panzer::Traits::Residual>(Teuchos::rcp(new XCoordinate<panzer::Traits::Residual,panzer::Traits>(pl)));
    }

    {
       Teuchos::RCP<std::vector<std::string> > fieldNames
             = Teuchos::rcp(new std::vector<std::string>);
       fieldNames->push_back("x-coord");
    
       Teuchos::ParameterList pl;
       pl.set("Mesh",mesh);
       pl.set("IR",intRule);
       pl.set("Field Names",fieldNames);
       pl.set("Scatter Name", "xcoord-scatter-cell-residual");
       Teuchos::RCP<PHX::Evaluator<panzer::Traits> > eval
             = Teuchos::rcp(new panzer_stk::ScatterCellAvgQuantity<panzer::Traits::Residual,panzer::Traits>(pl));
       fm->registerEvaluator<panzer::Traits::Residual>(eval);
       fm->requireField<panzer::Traits::Residual>(*eval->evaluatedFields()[0]);
    }

    {
       Teuchos::RCP<std::vector<std::string> > fieldNames
             = Teuchos::rcp(new std::vector<std::string>);
       fieldNames->push_back("x-coord");

       std::string scatterName = "xcoord-scatter-residual";
       Teuchos::RCP<PHX::Evaluator<panzer::Traits> > eval
             = Teuchos::rcp(new panzer_stk::ScatterFields<panzer::Traits::Residual,panzer::Traits>(scatterName,mesh,linBasis,*fieldNames));
       fm->registerEvaluator<panzer::Traits::Residual>(eval);
       fm->requireField<panzer::Traits::Residual>(*eval->evaluatedFields()[0]);
    }

    // build physics blocks
    //////////////////////////////////////////////////////////////
    std::vector<Teuchos::RCP<panzer::PhysicsBlock> > physicsBlocks;
    {
      Teuchos::RCP<user_app::MyFactory> eqset_factory = Teuchos::rcp(new user_app::MyFactory);
      user_app::BCFactory bc_factory;
      const int default_integration_order = 1;

      std::map<std::string,std::string> block_ids_to_physics_ids;
      block_ids_to_physics_ids["eblock-0_0"] = "test physics";

      std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;
      block_ids_to_cell_topo["eblock-0_0"] = mesh->getCellTopology("eblock-0_0");
      
      Teuchos::RCP<panzer::GlobalData> gd = panzer::createGlobalData();

      panzer::buildPhysicsBlocks(block_ids_to_physics_ids,
				 block_ids_to_cell_topo,
				 ipb,
				 default_integration_order,
				 workset_size,
				 eqset_factory,
				 gd,
				 false,
				 physicsBlocks);
    }

    // register jacobian size
    //////////////////////////////////////////////////////////////
    std::vector<PHX::index_size_type> derivative_dimensions;
    derivative_dimensions.push_back(9+4);
    fm->setKokkosExtendedDataTypeDimensions<panzer::Traits::Jacobian>(derivative_dimensions);

    // build worksets
    //////////////////////////////////////////////////////////////

    Teuchos::RCP<panzer::PhysicsBlock> physics_block_one = panzer::findPhysicsBlock("eblock-0_0",physicsBlocks);
    Teuchos::RCP<std::vector<panzer::Workset> > volume_worksets = panzer_stk::buildWorksets(*mesh,*physics_block_one);
//.........这里部分代码省略.........

QCAD::ResponseCenterOfMass<EvalT, Traits>::
ResponseCenterOfMass(Teuchos::ParameterList& p,
		     const Teuchos::RCP<Albany::Layouts>& dl) :
  coordVec("Coord Vec", dl->qp_vector),
  weights("Weights", dl->qp_scalar)
{
  // get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  // number of quad points per cell and dimension of space
  Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
  Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
  
  std::vector<PHX::DataLayout::size_type> dims;
  vector_dl->dimensions(dims);
  numQPs  = dims[1];
  numDims = dims[2];

  //! Get material DB from parameters passed down from problem (if given)
  Teuchos::RCP<QCAD::MaterialDatabase> materialDB;
  Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem = 
    p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");
  if(paramsFromProblem != Teuchos::null)
    materialDB = paramsFromProblem->get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
  else materialDB = Teuchos::null;
  
  // User-specified parameters
  fieldName = plist->get<std::string>("Field Name");
  opRegion  = Teuchos::rcp( new QCAD::MeshRegion<EvalT, Traits>("Coord Vec","Weights",*plist,materialDB,dl) );
  
  // setup field
  PHX::MDField<ScalarT> f(fieldName, scalar_dl); field = f;

  // add dependent fields
  this->addDependentField(field);
  this->addDependentField(coordVec);
  this->addDependentField(weights);
  opRegion->addDependentFields(this);

  this->setName(fieldName+" Response Center of Mass" );

  using PHX::MDALayout;

  // Setup scatter evaluator
  p.set("Stand-alone Evaluator", false);
  std::string local_response_name = 
    fieldName + " Local Response Center of Mass";
  std::string global_response_name = 
    fieldName + " Global Response Center of Mass";
  int worksetSize = scalar_dl->dimension(0);
  int responseSize = 4;
  Teuchos::RCP<PHX::DataLayout> local_response_layout =
    Teuchos::rcp(new MDALayout<Cell,Dim>(worksetSize, responseSize));
  Teuchos::RCP<PHX::DataLayout> global_response_layout =
    Teuchos::rcp(new MDALayout<Dim>(responseSize));
  PHX::Tag<ScalarT> local_response_tag(local_response_name, 
				       local_response_layout);
  PHX::Tag<ScalarT> global_response_tag(global_response_name, 
					global_response_layout);
  p.set("Local Response Field Tag", local_response_tag);
  p.set("Global Response Field Tag", global_response_tag);
  PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);
}

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

#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  int MyPID = Comm.MyPID();

  // Number of dimension of the domain
  int space_dim = 2;
  
  // Size of each of the dimensions of the domain
  std::vector<double> brick_dim( space_dim );
  brick_dim[0] = 1.0;
  brick_dim[1] = 1.0;
  
  // Number of elements in each of the dimensions of the domain
  std::vector<int> elements( space_dim );
  elements[0] = 10;
  elements[1] = 10;
  
  // Create problem
  Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
  
  // Get the stiffness and mass matrices
  Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
  Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
	
  //
  // *******************************************************
  // Set up Amesos direct solver for inner iteration
  // *******************************************************
  //

  // Create the shifted system K - sigma * M.
  // For the buckling transformation, this shift must be nonzero.

  double sigma = 1.0;
  Epetra_CrsMatrix Kcopy( *K );

  int addErr = EpetraExt::MatrixMatrix::Add( *M, false, -sigma, Kcopy, 1.0 );
  if (addErr != 0) {
    if (MyPID == 0) {
      std::cout << "EpetraExt::MatrixMatrix::Add returned with error: " << addErr << std::endl;
    }
#ifdef HAVE_MPI
    MPI_Finalize() ;
#endif
    return -1;
  }

  // Create Epetra linear problem class to solve "x = b"
  Epetra_LinearProblem AmesosProblem;
  AmesosProblem.SetOperator(&Kcopy);
  
  // Create Amesos factory and solver for solving "(K - sigma*M)x = b" using a direct factorization
  Amesos amesosFactory;
  Teuchos::RCP<Amesos_BaseSolver> AmesosSolver = 
    Teuchos::rcp( amesosFactory.Create( "Klu", AmesosProblem ) );

  // The AmesosBucklingOp class assumes that the symbolic/numeric factorizations have already
  // been performed on the linear problem.
  AmesosSolver->SymbolicFactorization();
  AmesosSolver->NumericFactorization();

  //
  // ************************************
  // Start the block Arnoldi iteration
  // ************************************
  //
  //  Variables used for the Block Arnoldi Method
  //
  int nev = 10;
  int blockSize = 3;  
  int numBlocks = 3*nev/blockSize;
  int maxRestarts = 5;
  //int step = 5;
  double tol = 1.0e-8;
  std::string which = "LM";
  int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
  //
  // Create parameter list to pass into solver
  //
  Teuchos::ParameterList MyPL;
  MyPL.set( "Verbosity", verbosity );
  MyPL.set( "Which", which );
  MyPL.set( "Block Size", blockSize );
  MyPL.set( "Num Blocks", numBlocks );
  MyPL.set( "Maximum Restarts", maxRestarts );
  MyPL.set( "Convergence Tolerance", tol );
  //MyPL.set( "Step Size", step );
  
  typedef Epetra_MultiVector MV;
  typedef Epetra_Operator OP;
  typedef Anasazi::MultiVecTraits<double, MV> MVT;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    indices[0] = 1; indices[1] = 4; indices[2] = 5; indices[3] = 2;
    matrix->InsertGlobalValues(nodesPerElem, &indices[0], &coefs[0]);
  }
  else if (localProc == 1) {
    //element E2:
    indices[0] = 3; indices[1] = 2; indices[2] = 7; indices[3] = 8;
    matrix->InsertGlobalValues(nodesPerElem, &indices[0], &coefs[0]);
  }
  else { //localProc==2
    //element E3:
    indices[0] = 2; indices[1] = 5; indices[2] = 6; indices[3] = 7;
    matrix->InsertGlobalValues(nodesPerElem, &indices[0], &coefs[0]);
  }

  int err = matrix->GlobalAssemble();
  if (err != 0) {
    std::cout << "err="<<err<<" returned from matrix->GlobalAssemble()"
      << std::endl;
  }

//  std::cout << "matrix: " << std::endl;
//  std::cout << *matrix << std::endl;

  //We'll need a Teuchos::ParameterList object to pass to the
  //Isorropia::Epetra::Partitioner class.
  Teuchos::ParameterList paramlist;

#ifdef HAVE_ISORROPIA_ZOLTAN
  // If Zoltan is available, we'll specify that the Zoltan package be
  // used for the partitioning operation, by creating a parameter
  // sublist named "Zoltan".
  // In the sublist, we'll set parameters that we want sent to Zoltan.

  paramlist.set("PARTITIONING METHOD", "GRAPH");
  paramlist.set("PRINT ZOLTAN METRICS", "2");
  Teuchos::ParameterList& sublist = paramlist.sublist("Zoltan");
  sublist.set("GRAPH_PACKAGE", "PHG");

  //sublist.set("DEBUG_LEVEL", "1"); // Zoltan will print out parameters
  //sublist.set("DEBUG_LEVEL", "5");   // proc 0 will trace Zoltan calls
  //sublist.set("DEBUG_MEMORY", "2");  // Zoltan will trace alloc & free

#else
  // If Zoltan is not available, a simple linear partitioner will be
  // used to partition such that the number of nonzeros is equal (or
  // close to equal) on each processor. No parameter is necessary to
  // specify this.
#endif


  Teuchos::RCP<Isorropia::Epetra::CostDescriber> costs =
    Teuchos::rcp(new Isorropia::Epetra::CostDescriber);

  //Next create a matrix which is a copy of the matrix we just
  //assembled, but we'll replace the values with graph edge weights.

  Teuchos::RCP<Epetra_FECrsMatrix> ge_weights =
    Teuchos::rcp(new Epetra_FECrsMatrix(*matrix));

  Teuchos::RCP<Epetra_CrsMatrix> crs_ge_weights;
  crs_ge_weights = ge_weights;

  //Fill the matrix with a "default" weight of 1.0.
  crs_ge_weights->PutScalar(1.0);

  //Now we'll put a "large" weight on edges that connect nodes