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

int RestrictedMultiVectorWrapper::restrict_comm(Teuchos::RCP<Epetra_MultiVector> input_mv){
 input_mv_=input_mv;

  /* Pull the Input Matrix Info */
  const Epetra_MpiComm *InComm = dynamic_cast<const Epetra_MpiComm*>(& input_mv_->Comm());
  const Epetra_BlockMap *InMap = dynamic_cast<const Epetra_BlockMap*>(& input_mv_->Map());

  if(!InComm || !InMap) return -1;
  
  if(!subcomm_is_set){
    /* Build the Split Communicators, If Needed */
    int color;
    if(InMap->NumMyElements()) color=1;
    else color=MPI_UNDEFINED;
    MPI_Comm_split(InComm->Comm(),color,InComm->MyPID(),&MPI_SubComm_);
  }
  else{
    /* Sanity check user-provided subcomm - drop an error if the MPISubComm
       does not include a processor with data. */
    if (input_mv->MyLength() && MPI_SubComm_ == MPI_COMM_NULL)
      return(-2);
  }

  /* Mark active processors */
  if(MPI_SubComm_ == MPI_COMM_NULL) proc_is_active=false;
  else proc_is_active=true;

  
  int Nrows=InMap->NumGlobalElements();
      
  if(proc_is_active){      
    RestrictedComm_=new Epetra_MpiComm(MPI_SubComm_);
    
    /* Build the Restricted Maps */
    ResMap_ = new Epetra_BlockMap(Nrows,InMap->NumMyElements(),InMap->MyGlobalElements(),
                                  InMap->ElementSizeList(),InMap->IndexBase(),*RestrictedComm_);

    /* Allocate the restricted matrix*/
    double *A; int LDA;
    input_mv_->ExtractView(&A,&LDA);
    restricted_mv_ = Teuchos::rcp(new Epetra_MultiVector(View,*ResMap_,A,LDA,input_mv_->NumVectors()));
  }
}/*end restrict_comm*/

void
Albany::ModelEvaluator::evalModel(const InArgs& inArgs,
                                 const OutArgs& outArgs) const
{
  Teuchos::TimeMonitor Timer(*timer); //start timer
  //
  // Get the input arguments
  //
  Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
  Teuchos::RCP<const Epetra_Vector> x_dot;
  Teuchos::RCP<const Epetra_Vector> x_dotdot;
  double alpha     = 0.0;
  double omega     = 0.0;
  double beta      = 1.0;
  double curr_time = 0.0;
  x_dot = inArgs.get_x_dot();
  x_dotdot = inArgs.get_x_dotdot();
  if (x_dot != Teuchos::null || x_dotdot != Teuchos::null) {
    alpha = inArgs.get_alpha();
    omega = inArgs.get_omega();
    beta = inArgs.get_beta();
    curr_time  = inArgs.get_t();
  }
  for (int i=0; i<num_param_vecs; i++) {
    Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i);
    if (p != Teuchos::null) {
      for (unsigned int j=0; j<sacado_param_vec[i].size(); j++)
        sacado_param_vec[i][j].baseValue = (*p)[j];
    }
  }
  for (int i=0; i<num_dist_param_vecs; i++) {
    Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_param_vecs);
    if (p != Teuchos::null) {
      *(distParamLib->get(dist_param_names[i])->vector()) = *p;
    }
  }

  //
  // Get the output arguments
  //
  EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out = outArgs.get_f();
  Teuchos::RCP<Epetra_Operator> W_out = outArgs.get_W();

  // Cast W to a CrsMatrix, throw an exception if this fails
  Teuchos::RCP<Epetra_CrsMatrix> W_out_crs;

  if (W_out != Teuchos::null)
    W_out_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out, true);

int test_var = 0;
if(test_var != 0){
std::cout << "The current solution length is: " << x->MyLength() << std::endl;
x->Print(std::cout);

}

  // Get preconditioner operator, if requested
  Teuchos::RCP<Epetra_Operator> WPrec_out;
  if (outArgs.supports(OUT_ARG_WPrec)) WPrec_out = outArgs.get_WPrec();

  //
  // Compute the functions
  //
  bool f_already_computed = false;

  // W matrix
  if (W_out != Teuchos::null) {
    app->computeGlobalJacobian(alpha, beta, omega, curr_time, x_dot.get(), x_dotdot.get(),*x,
                               sacado_param_vec, f_out.get(), *W_out_crs);
    f_already_computed=true;
if(test_var != 0){
//std::cout << "The current rhs length is: " << f_out->MyLength() << std::endl;
//f_out->Print(std::cout);
std::cout << "The current Jacobian length is: " << W_out_crs->NumGlobalRows() << std::endl;
W_out_crs->Print(std::cout);
}
  }

  if (WPrec_out != Teuchos::null) {
    app->computeGlobalJacobian(alpha, beta, omega, curr_time, x_dot.get(), x_dotdot.get(), *x,
                               sacado_param_vec, f_out.get(), *Extra_W_crs);
    f_already_computed=true;
if(test_var != 0){
//std::cout << "The current rhs length is: " << f_out->MyLength() << std::endl;
//f_out->Print(std::cout);
std::cout << "The current preconditioner length is: " << Extra_W_crs->NumGlobalRows() << std::endl;
Extra_W_crs->Print(std::cout);
}

    app->computeGlobalPreconditioner(Extra_W_crs, WPrec_out);
  }

  // scalar df/dp
  for (int i=0; i<num_param_vecs; i++) {
    Teuchos::RCP<Epetra_MultiVector> dfdp_out =
      outArgs.get_DfDp(i).getMultiVector();
    if (dfdp_out != Teuchos::null) {
      Teuchos::Array<int> p_indexes =
        outArgs.get_DfDp(i).getDerivativeMultiVector().getParamIndexes();
      Teuchos::RCP<ParamVec> p_vec;
//.........这里部分代码省略.........

void
Albany::ModelEvaluator::evalModel(const InArgs& inArgs,
                                 const OutArgs& outArgs) const
{
  Teuchos::TimeMonitor Timer(*timer); //start timer
  //
  // Get the input arguments
  //
  Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
  Teuchos::RCP<const Epetra_Vector> x_dot;
  Teuchos::RCP<const Epetra_Vector> x_dotdot;

  //create comm and node objects for Epetra -> Tpetra conversions
  Teuchos::RCP<const Teuchos::Comm<int> > commT = app->getComm();
  Teuchos::RCP<Epetra_Comm> comm = Albany::createEpetraCommFromTeuchosComm(commT);
  //Create Tpetra copy of x, call it xT
  Teuchos::RCP<const Tpetra_Vector> xT;
  if (x != Teuchos::null)
    xT  = Petra::EpetraVector_To_TpetraVectorConst(*x, commT);

  double alpha     = 0.0;
  double omega     = 0.0;
  double beta      = 1.0;
  double curr_time = 0.0;

  if(num_time_deriv > 0)
    x_dot = inArgs.get_x_dot();
  if(num_time_deriv > 1)
    x_dotdot = inArgs.get_x_dotdot();

  //Declare and create Tpetra copy of x_dot, call it x_dotT
  Teuchos::RCP<const Tpetra_Vector> x_dotT;
  if (Teuchos::nonnull(x_dot))
    x_dotT = Petra::EpetraVector_To_TpetraVectorConst(*x_dot, commT);

  //Declare and create Tpetra copy of x_dotdot, call it x_dotdotT
  Teuchos::RCP<const Tpetra_Vector> x_dotdotT;
  if (Teuchos::nonnull(x_dotdot))
    x_dotdotT = Petra::EpetraVector_To_TpetraVectorConst(*x_dotdot, commT);

  if (Teuchos::nonnull(x_dot)){
    alpha = inArgs.get_alpha();
    beta = inArgs.get_beta();
    curr_time  = inArgs.get_t();
  }
  if (Teuchos::nonnull(x_dotdot)) {
    omega = inArgs.get_omega();
  }

  for (int i=0; i<num_param_vecs; i++) {
    Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i);
    if (p != Teuchos::null) {
      for (unsigned int j=0; j<sacado_param_vec[i].size(); j++) {
        sacado_param_vec[i][j].baseValue = (*p)[j];
      }
    }
  }
  for (int i=0; i<num_dist_param_vecs; i++) {
    Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_param_vecs);
    //create Tpetra copy of p
    Teuchos::RCP<const Tpetra_Vector> pT;
    if (p != Teuchos::null) {
      pT = Petra::EpetraVector_To_TpetraVectorConst(*p, commT);
      //*(distParamLib->get(dist_param_names[i])->vector()) = *p;
      *(distParamLib->get(dist_param_names[i])->vector()) = *pT;
    }
  }

  //
  // Get the output arguments
  //
  EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out = outArgs.get_f();
  Teuchos::RCP<Epetra_Operator> W_out = outArgs.get_W();

  // Cast W to a CrsMatrix, throw an exception if this fails
  Teuchos::RCP<Epetra_CrsMatrix> W_out_crs;
#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
  //IK, 7/15/14: adding object to hold mass matrix to be written to matrix market file
  Teuchos::RCP<Epetra_CrsMatrix> Mass;
  //IK, 7/15/14: needed for writing mass matrix out to matrix market file
  EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> ftmp = outArgs.get_f();
#endif

  if (W_out != Teuchos::null) {
    W_out_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out, true);
#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
    //IK, 7/15/14: adding object to hold mass matrix to be written to matrix market file
    Mass = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out, true);
#endif
  }


int test_var = 0;
if(test_var != 0){
std::cout << "The current solution length is: " << x->MyLength() << std::endl;
x->Print(std::cout);

}

  // Get preconditioner operator, if requested
//.........这里部分代码省略.........

void InitialConditions(const Teuchos::RCP<Epetra_Vector>& soln,
                       const Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >& wsElNodeEqID,
                       const Teuchos::ArrayRCP<std::string>& wsEBNames,
                       const Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > > coords,
                       const int neq, const int numDim,
                       Teuchos::ParameterList& icParams, const bool hasRestartSolution) {
  // Called twice, with x and xdot. Different param lists are sent in.
  icParams.validateParameters(*AAdapt::getValidInitialConditionParameters(wsEBNames), 0);

  // Default function is Constant, unless a Restart solution vector
  // was used, in which case the Init COnd defaults to Restart.
  std::string name;

  if(!hasRestartSolution) name = icParams.get("Function", "Constant");

  else                     name = icParams.get("Function", "Restart");

  if(name == "Restart") return;

  // Handle element block specific constant data

  if(name == "EBPerturb" || name == "EBPerturbGaussian" || name == "EBConstant") {

    bool perturb_values = false;

    Teuchos::Array<double> defaultData(neq);
    Teuchos::Array<double> perturb_mag;

    // Only perturb if the user has told us by how much to perturb
    if(name != "EBConstant" && icParams.isParameter("Perturb IC")) {

      perturb_values = true;

      perturb_mag = icParams.get("Perturb IC", defaultData);

    }

    /* The element block-based IC specification here is currently a hack. It assumes the initial value is constant
     * within each element across the element block (or optionally perturbed somewhat element by element). The
     * proper way to do this would be to project the element integration point values to the nodes using the basis
     * functions and a consistent mass matrix.
     *
     * The current implementation uses a single integration point per element - this integration point value for this
     * element within the element block is specified in the input file (and optionally perturbed). An approximation
     * of the load vector is obtained by accumulating the resulting (possibly perturbed) value into the nodes. Then,
     * a lumped version of the mass matrix is inverted and used to solve for the approximate nodal point initial
     * conditions.
     */

    // Use an Epetra_Vector to hold the lumped mass matrix (has entries only on the diagonal). Zero-ed out.

    Epetra_Vector lumpedMM(soln->Map(), true);

    // Make sure soln is zeroed - we are accumulating into it

    for(int i = 0; i < soln->MyLength(); i++)

      (*soln)[i] = 0;

    // Loop over all worksets, elements, all local nodes: compute soln as a function of coord and wsEBName


    Teuchos::RCP<AAdapt::AnalyticFunction> initFunc;

    for(int ws = 0; ws < wsElNodeEqID.size(); ws++) { // loop over worksets

      Teuchos::Array<double> data = icParams.get(wsEBNames[ws], defaultData);
      // Call factory method from library of initial condition functions

      if(perturb_values) {

        if(name == "EBPerturb")

          initFunc = Teuchos::rcp(new AAdapt::ConstantFunctionPerturbed(neq, numDim, ws, data, perturb_mag));

        else // name == EBGaussianPerturb

          initFunc = Teuchos::rcp(new
                                  AAdapt::ConstantFunctionGaussianPerturbed(neq, numDim, ws, data, perturb_mag));
      }

      else

        initFunc = Teuchos::rcp(new AAdapt::ConstantFunction(neq, numDim, data));

      std::vector<double> X(neq);
      std::vector<double> x(neq);

      for(int el = 0; el < wsElNodeEqID[ws].size(); el++) { // loop over elements in workset

        for(int i = 0; i < neq; i++)
          X[i] = 0;

        for(int ln = 0; ln < wsElNodeEqID[ws][el].size(); ln++) // loop over node local to the element
          for(int i = 0; i < neq; i++)
            X[i] += coords[ws][el][ln][i]; // nodal coords

        for(int i = 0; i < neq; i++)
          X[i] /= (double)neq;

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

int 
RestrictedMultiVectorWrapper::
restrict_comm (Teuchos::RCP<Epetra_MultiVector> input_mv)
{
  using Teuchos::rcp;
  input_mv_ = input_mv;

  // Extract the input MV's communicator and Map.
  const Epetra_MpiComm *InComm = 
    dynamic_cast<const Epetra_MpiComm*> (& (input_mv_->Comm ()));
  const Epetra_BlockMap *InMap = 
    dynamic_cast<const Epetra_BlockMap*> (& (input_mv_->Map ()));

  if (! InComm || ! InMap) {
    return -1; // At least one dynamic cast failed.
  }
  
  if (! subcomm_is_set) {
    /* Build the Split Communicators, If Needed */
    int color;
    if (InMap->NumMyElements()) {
      color = 1;
    }
    else {
      color = MPI_UNDEFINED;
    }
    const int err = 
      MPI_Comm_split (InComm->Comm(), color, InComm->MyPID(), &MPI_SubComm_);
    if (err != MPI_SUCCESS) {
      return -2;
    }
  }
  else {
    /* Sanity check user-provided subcomm - drop an error if the MPISubComm
       does not include a processor with data. */
    if (input_mv->MyLength() && MPI_SubComm_ == MPI_COMM_NULL) {
      return -2;
    }
  }

  /* Mark active processors */
  if (MPI_SubComm_ == MPI_COMM_NULL) {
    proc_is_active = false;
  }
  else {
    proc_is_active = true;
  }

  if (proc_is_active) {

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  if(InMap->GlobalIndicesInt()) {
    int Nrows = InMap->NumGlobalElements ();
    RestrictedComm_ = new Epetra_MpiComm (MPI_SubComm_);
    
    // Build the restricted Maps
    ResMap_ = new Epetra_BlockMap (Nrows, InMap->NumMyElements(), 
				   InMap->MyGlobalElements(),
				   InMap->ElementSizeList(),
				   InMap->IndexBase(), *RestrictedComm_);
  }
  else
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  if(InMap->GlobalIndicesLongLong()) {
    long long Nrows = InMap->NumGlobalElements64 ();
    RestrictedComm_ = new Epetra_MpiComm (MPI_SubComm_);
    
    // Build the restricted Maps
    ResMap_ = new Epetra_BlockMap (Nrows, InMap->NumMyElements(), 
				   InMap->MyGlobalElements64(),
				   InMap->ElementSizeList(),
				   InMap->IndexBase64(), *RestrictedComm_);
  }
  else
#endif
    throw "EpetraExt::RestrictedMultiVectorWrapper::restrict_comm ERROR: GlobalIndices type unknown";

    // Allocate the restricted matrix
    double *A; 
    int LDA;
    input_mv_->ExtractView (&A,&LDA);
    restricted_mv_ = rcp (new Epetra_MultiVector (View, *ResMap_, A, LDA, 
						  input_mv_->NumVectors ()));
  }
  return 0; // Success!
}/*end restrict_comm*/

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

    ICFieldDescriptor densDesc;
    densDesc.fieldName  = "DENSITY";
    densDesc.basisName  = "Const";
    densDesc.basisOrder = 0;

    ICFieldDescriptor condDesc;
    condDesc.fieldName  = "CONDUCTIVITY";
    condDesc.basisName  = "HGrad";
    condDesc.basisOrder = 1;

    RCP<panzer::PureBasis> const_basis = rcp(new panzer::PureBasis(densDesc.basisName,densDesc.basisOrder,cellData));
    RCP<panzer::PureBasis> hgrad_basis = rcp(new panzer::PureBasis(condDesc.basisName,condDesc.basisOrder,cellData));

    RCP<const panzer::FieldPattern> constFP = rcp(new panzer::Intrepid2FieldPattern(const_basis->getIntrepid2Basis()));
    RCP<const panzer::FieldPattern> hgradFP = rcp(new panzer::Intrepid2FieldPattern(hgrad_basis->getIntrepid2Basis()));
 
    // setup DOF manager
    /////////////////////////////////////////////
    RCP<panzer::ConnManager<int,int> > conn_manager 
           = Teuchos::rcp(new panzer_stk_classic::STKConnManager<int>(mesh));

    RCP<panzer::DOFManager<int,int> > dofManager
        = rcp(new panzer::DOFManager<int,int>(conn_manager,MPI_COMM_WORLD));
    dofManager->addField(densDesc.fieldName, constFP);
    dofManager->addField(condDesc.fieldName, hgradFP);
    dofManager->buildGlobalUnknowns();

    Teuchos::RCP<panzer::EpetraLinearObjFactory<panzer::Traits,int> > elof 
          = Teuchos::rcp(new panzer::EpetraLinearObjFactory<panzer::Traits,int>(comm.getConst(),dofManager));

    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof = elof;

    // setup worksets
    /////////////////////////////////////////////
    
    std::map<std::string,panzer::WorksetNeeds> needs;
    needs["eblock-0_0"].cellData = cellData;
    needs["eblock-0_0"].int_rules.push_back(int_rule);
    needs["eblock-0_0"].bases          = { const_basis,    hgrad_basis};
    needs["eblock-0_0"].rep_field_name = {   densDesc.fieldName, condDesc.fieldName};

    needs["eblock-1_0"].cellData = cellData;
    needs["eblock-1_0"].int_rules.push_back(int_rule);
    needs["eblock-1_0"].bases          = { const_basis,    hgrad_basis};
    needs["eblock-1_0"].rep_field_name = {   densDesc.fieldName, condDesc.fieldName};

    Teuchos::RCP<panzer_stk_classic::WorksetFactory> wkstFactory 
        = Teuchos::rcp(new panzer_stk_classic::WorksetFactory(mesh)); // build STK workset factory
    Teuchos::RCP<panzer::WorksetContainer> wkstContainer              // attach it to a workset container (uses lazy evaluation)
        = Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,needs));

    // setup field manager builder
    /////////////////////////////////////////////
      
    // Add in the application specific closure model factory
    panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory; 
    user_app::STKModelFactory_TemplateBuilder cm_builder;
    cm_factory.buildObjects(cm_builder);

    Teuchos::ParameterList user_data("User Data");

    Teuchos::ParameterList ic_closure_models("Initial Conditions");
    ic_closure_models.sublist("eblock-0_0").sublist(densDesc.fieldName).set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-0_0").sublist(condDesc.fieldName).set<double>("Value",9.0);
    ic_closure_models.sublist("eblock-1_0").sublist(densDesc.fieldName).set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-1_0").sublist(condDesc.fieldName).set<double>("Value",9.0);    

    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");
    block_ids_to_cell_topo["eblock-1_0"] = mesh->getCellTopology("eblock-1_0");

    std::map<std::string,std::vector<ICFieldDescriptor> > block_ids_to_fields;
    block_ids_to_fields["eblock-0_0"] = {densDesc,condDesc};
    block_ids_to_fields["eblock-1_0"] = {densDesc,condDesc};
    
    int workset_size = 4;

    Teuchos::RCP<panzer::LinearObjContainer> loc  = lof->buildLinearObjContainer();
    lof->initializeContainer(panzer::LinearObjContainer::X,*loc);
    Teuchos::RCP<panzer::EpetraLinearObjContainer> eloc = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
    Teuchos::RCP<Thyra::VectorBase<double> > vec = eloc->get_x_th();

    // this is the Function under test
    panzer::setupControlInitialCondition(block_ids_to_cell_topo,
                                         block_ids_to_fields,
                                         *wkstContainer,
                                         *lof,cm_factory,ic_closure_models,user_data,
                                         workset_size,
                                         0.0, // t0
                                         vec); 

    Teuchos::RCP<Epetra_Vector> x = eloc->get_x();
    out << x->GlobalLength() << " " << x->MyLength() << std::endl;
    for (int i=0; i < x->MyLength(); ++i) {
      double v = (*x)[i];
      TEST_ASSERT(v==3.0 || v==9.0); 
    }

  }

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

  // Loop to solve the same eigenproblem numtrial times (different initial vectors)
  int numfailed = 0;
  int iter = 0;
  double solvetime = 0;

  // Set random seed to have consistent initial vectors between
  // experiments.  Different seed in each loop iteration.
  ivec->SetSeed(2*(MyPID) +1); // Odd seed

  // Set up initial vectors
  // Using random values as the initial guess.
  if (initvec == "random"){
    MVT::MvRandom(*ivec);
  } 
  else if (initvec == "zero"){
    // All zero initial vector should be essentially the same,
    // but appears slightly worse in practice.
    ivec->PutScalar(0.);
  }
  else if (initvec == "unit"){
    // Orthogonal unit initial vectors.
    ivec->PutScalar(0.);
    for (int i = 0; i < blockSize; i++)
      ivec->ReplaceGlobalValue(i,i,1.);
  }
  else if (initvec == "random2"){
    // Partially random but orthogonal (0,1) initial vectors.
    // ivec(i,*) is zero in all but one column (for each i)
    // Inefficient implementation but this is only done once...
    double rowmax;
    int col;
    ivec->Random();
    for (int i = 0; i < ivec->MyLength(); i++){
      rowmax = -1;
      col = -1;
      for (int j = 0; j < blockSize; j++){
        // Make ivec(i,j) = 1 for largest random value in row i
        if ((*ivec)[j][i] > rowmax){
          rowmax = (*ivec)[j][i];
          col = j;
        }
        ivec->ReplaceMyValue(i,j,0.);
      }
      ivec->ReplaceMyValue(i,col,1.);
    }
  }
  else
    cout << "ERROR: Unknown value for initial vectors." << endl;

  if (verbose && (ivec->GlobalLength64() < TINYMATRIX)) 
    ivec->Print(std::cout);
  
  // Inform the eigenproblem that you are finished passing it information
  
  bool boolret = MyProblem->setProblem();
  if (boolret != true) {
    if (verbose && MyPID == 0) {
      cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." 
           << endl;
    }
    FINALIZE;
    return -1;
  }
 
  Teuchos::RCP<Anasazi::SolverManager<double, MV, OP> > MySolverMgr;

//.........这里部分代码省略.........
    TEST_ASSERT(myGlobalElements[3] == 3);
    TEST_ASSERT(myGlobalElements[4] == 0);
    TEST_ASSERT(myGlobalElements[5] == 2);
    TEST_ASSERT(myGlobalElements[6] == 4);
    TEST_ASSERT(myGlobalElements[7] == 6);
  }

  // check the bond map
  // the horizon was chosen such that each point should have three neighbors
  // note that if the NeighborhoodType parameter is not set to Spherical, this will fail
  Teuchos::RCP<const Epetra_BlockMap> bondMap = discretization->getGlobalBondMap();
  TEST_ASSERT(bondMap->NumGlobalElements() == 8);
  TEST_ASSERT(bondMap->NumMyElements() == 4);
  TEST_ASSERT(bondMap->IndexBase() == 0);
  TEST_ASSERT(bondMap->UniqueGIDs() == true);
  myGlobalElements = bondMap->MyGlobalElements();
  if(rank == 0){
    TEST_ASSERT(myGlobalElements[0] == 0);
    TEST_ASSERT(myGlobalElements[1] == 2);
    TEST_ASSERT(myGlobalElements[2] == 4);
    TEST_ASSERT(myGlobalElements[3] == 6);
  }
  if(rank == 1){
    TEST_ASSERT(myGlobalElements[0] == 5);
    TEST_ASSERT(myGlobalElements[1] == 7);
    TEST_ASSERT(myGlobalElements[2] == 1);
    TEST_ASSERT(myGlobalElements[3] == 3);
  }
  TEST_ASSERT(discretization->getNumBonds() == 4*3);

  // check the initial positions
  // all three coordinates are contained in a single vector
  Teuchos::RCP<Epetra_Vector> initialX = discretization->getInitialX();
  TEST_ASSERT(initialX->MyLength() == 4*3);
  TEST_ASSERT(initialX->GlobalLength() == 8*3);
  if(rank == 0){
    TEST_FLOATING_EQUALITY((*initialX)[0],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[1],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[2],  0.25, 1.0e-16);
 
    TEST_FLOATING_EQUALITY((*initialX)[3],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[4],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[5],  0.25, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[6],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[7],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[8],  0.75, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[9],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[10], 0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[11], 0.75, 1.0e-16);
  }
  if(rank == 1){
    TEST_FLOATING_EQUALITY((*initialX)[0],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[1],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[2],  0.75, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[3],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[4],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[5],  0.75, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[6],  0.75, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[7],  0.25, 1.0e-16);
    TEST_FLOATING_EQUALITY((*initialX)[8],  0.25, 1.0e-16);

    TEST_FLOATING_EQUALITY((*initialX)[9],  0.75, 1.0e-16);

//.........这里部分代码省略.........
      
       panzer::buildPhysicsBlocks(block_ids_to_physics_ids,
                                  block_ids_to_cell_topo,
				  ipb,
				  default_integration_order,
				  workset_size,
                                  eqset_factory,
				  gd,
		    	          false,
                                  physics_blocks);
    }

    // setup worksets
    /////////////////////////////////////////////

    Teuchos::RCP<panzer_stk::WorksetFactory> wkstFactory 
       = Teuchos::rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory
    Teuchos::RCP<panzer::WorksetContainer> wkstContainer     // attach it to a workset container (uses lazy evaluation)
       = Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,physics_blocks,workset_size));

    // get vector of element blocks
    std::vector<std::string> elementBlocks;
    mesh->getElementBlockNames(elementBlocks);
 
    // build volume worksets from container
    std::map<panzer::BC,Teuchos::RCP<std::map<unsigned,panzer::Workset> >,panzer::LessBC> bc_worksets;
    panzer::getSideWorksetsFromContainer(*wkstContainer,bcs,bc_worksets);

    // setup DOF manager
    /////////////////////////////////////////////
    const Teuchos::RCP<panzer::ConnManager<int,int> > conn_manager 
           = Teuchos::rcp(new panzer_stk::STKConnManager<int>(mesh));

    Teuchos::RCP<const panzer::UniqueGlobalIndexerFactory<int,int,int,int> > indexerFactory
          = Teuchos::rcp(new panzer::DOFManagerFactory<int,int>);
    const Teuchos::RCP<panzer::UniqueGlobalIndexer<int,int> > dofManager 
          = indexerFactory->buildUniqueGlobalIndexer(Teuchos::opaqueWrapper(MPI_COMM_WORLD),physics_blocks,conn_manager);

    // and linear object factory
    Teuchos::RCP<const Teuchos::MpiComm<int> > tComm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));

    Teuchos::RCP<panzer::EpetraLinearObjFactory<panzer::Traits,int> > elof 
          = Teuchos::rcp(new panzer::EpetraLinearObjFactory<panzer::Traits,int>(tComm.getConst(),dofManager));

    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof = elof;

    // setup field manager builder
    /////////////////////////////////////////////
      
    // Add in the application specific closure model factory
    panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory; 
    user_app::STKModelFactory_TemplateBuilder cm_builder;
    cm_factory.buildObjects(cm_builder);

    Teuchos::ParameterList closure_models("Closure Models");
    closure_models.sublist("solid").sublist("SOURCE_TEMPERATURE").set<double>("Value",1.0);
    closure_models.sublist("solid").sublist("SOURCE_ELECTRON_TEMPERATURE").set<double>("Value",1.0);
    closure_models.sublist("ion solid").sublist("SOURCE_ION_TEMPERATURE").set<double>("Value",1.0);

    Teuchos::ParameterList user_data("User Data");
    user_data.sublist("Panzer Data").set("Mesh", mesh);
    user_data.sublist("Panzer Data").set("DOF Manager", dofManager);
    user_data.sublist("Panzer Data").set("Linear Object Factory", lof);

    fmb.setWorksetContainer(wkstContainer);
    fmb.setupVolumeFieldManagers(physics_blocks,cm_factory,closure_models,*elof,user_data);
    fmb.setupBCFieldManagers(bcs,physics_blocks,*eqset_factory,cm_factory,bc_factory,closure_models,*elof,user_data);

    Teuchos::ParameterList ic_closure_models("Initial Conditions");
    ic_closure_models.sublist("eblock-0_0").sublist("TEMPERATURE").set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-0_0").sublist("ELECTRON_TEMPERATURE").set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-1_0").sublist("TEMPERATURE").set<double>("Value",3.0);
    ic_closure_models.sublist("eblock-1_0").sublist("ION_TEMPERATURE").set<double>("Value",3.0);    

    std::map<std::string, Teuchos::RCP< PHX::FieldManager<panzer::Traits> > > phx_ic_field_managers;
    panzer::setupInitialConditionFieldManagers(*wkstContainer,
					       physics_blocks,
					       cm_factory,
					       ic_closure_models,
					       *elof,
					       user_data,
					       true,
					       "initial_condition_test",
					       phx_ic_field_managers);

    
    Teuchos::RCP<panzer::LinearObjContainer> loc  = elof->buildLinearObjContainer();
    elof->initializeContainer(panzer::EpetraLinearObjContainer::X,*loc);

    Teuchos::RCP<panzer::EpetraLinearObjContainer> eloc = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
    eloc->get_x()->PutScalar(0.0);

    panzer::evaluateInitialCondition(*wkstContainer, phx_ic_field_managers, loc, *elof, 0.0);

    Teuchos::RCP<Epetra_Vector> x = eloc->get_x();
    out << x->GlobalLength() << " " << x->MyLength() << std::endl;
    for (int i=0; i < x->MyLength(); ++i)
      TEST_FLOATING_EQUALITY((*x)[i], 3.0, 1.0e-10);

  }