void get_1d_copy_helper<MV,S>::
    do_get (const Teuchos::Ptr<const MV>& mv,
            const Teuchos::ArrayView<S>& vals,
            const size_t ldx,
            EDistribution distribution,
            typename MV::global_ordinal_t indexBase)
    {
      typedef typename MV::local_ordinal_t lo_t;
      typedef typename MV::global_ordinal_t go_t;
      typedef typename MV::global_size_t gs_t;
      typedef typename MV::node_t node_t;

      TEUCHOS_TEST_FOR_EXCEPTION(
        mv.getRawPtr () == NULL, std::invalid_argument,
        "Amesos2::get_1d_copy_helper::do_get(5 args): mv is null.");

      Teuchos::RCP<const Tpetra::Map<lo_t,go_t,node_t> > map
        = Amesos2::Util::getDistributionMap<lo_t,go_t,gs_t,node_t> (distribution,
                                                                    mv->getGlobalLength (),
                                                                    mv->getComm (),
                                                                    indexBase,
                                                                    mv->getMap());

      do_get (mv, vals, ldx, Teuchos::ptrInArg (*map), distribution);
    }

    void diff_type_get_copy<MV,S>::
    apply (const Teuchos::Ptr<const MV>& mv,
           const Teuchos::ArrayView<S>& v,
           const size_t& ldx,
           Teuchos::Ptr<const Tpetra::Map<typename MV::local_ordinal_t, typename MV::global_ordinal_t, typename MV::node_t> > distribution_map,
           EDistribution distribution )
    {
      typedef typename MV::scalar_t mv_scalar_t;
      typedef typename Teuchos::Array<mv_scalar_t>::size_type size_type;

      TEUCHOS_TEST_FOR_EXCEPTION(
        mv.getRawPtr () == NULL, std::invalid_argument,
        "Amesos2::diff_type_get_copy::apply: mv is null.");
      TEUCHOS_TEST_FOR_EXCEPTION(
        distribution_map.getRawPtr () == NULL, std::invalid_argument,
        "Amesos2::diff_type_get_copy::apply: distribution_map is null.");

      const size_type vals_length = v.size ();
      Teuchos::Array<mv_scalar_t> vals_tmp (vals_length);

      mv->get1dCopy (vals_tmp (), ldx, distribution_map, distribution);
      for (size_type i = 0; i < vals_length; ++i) {
        v[i] = Teuchos::as<S> (vals_tmp[i]);
      }
    }

void MueLuTpetraPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::uninitializePrec(
  PreconditionerBase<Scalar> *prec,
  Teuchos::RCP<const LinearOpSourceBase<Scalar> > *fwdOp,
  ESupportSolveUse *supportSolveUse
  ) const
{
  // Check precondition

  TEUCHOS_ASSERT(prec);

  // Retrieve concrete preconditioner object

  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec =
    Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));

  if (fwdOp) {
    // TODO: Implement properly instead of returning default value
    *fwdOp = Teuchos::null;
  }

  if (supportSolveUse) {
    // TODO: Implement properly instead of returning default value
    *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
  }

  defaultPrec->uninitialize();
}

void performNodalMeshReduction(
    stk::mesh::Part &samplePart,
    stk::mesh::BulkData& bulkData)
{
  const stk::mesh::EntityRank nodeEntityRank(0);
  const stk::mesh::MetaData &metaData = stk::mesh::MetaData::get(bulkData);

  std::vector<stk::mesh::Entity *> sampleNodes;
  stk::mesh::get_selected_entities(samplePart, bulkData.buckets(nodeEntityRank), sampleNodes);

  const stk::mesh::Selector locallyOwned = stk::mesh::MetaData::get(bulkData).locally_owned_part();

  std::vector<stk::mesh::Entity *> relatedEntities;
  typedef boost::indirect_iterator<std::vector<stk::mesh::Entity *>::const_iterator> EntityIterator;
  for (EntityIterator it(sampleNodes.begin()), it_end(sampleNodes.end()); it != it_end; ++it) {
    const stk::mesh::PairIterRelation relations = it->relations();
    typedef stk::mesh::PairIterRelation::first_type RelationIterator;
    for (RelationIterator rel_it = relations.first, rel_it_end = relations.second; rel_it != rel_it_end; ++rel_it) {
      const Teuchos::Ptr<stk::mesh::Entity> relatedEntity(rel_it->entity());
      if (Teuchos::nonnull(relatedEntity) && locallyOwned(*relatedEntity)) {
        relatedEntities.push_back(relatedEntity.get());
      }
    }
  }
  std::sort(relatedEntities.begin(), relatedEntities.end(), stk::mesh::EntityLess());
  relatedEntities.erase(
      std::unique(relatedEntities.begin(), relatedEntities.end(), stk::mesh::EntityEqual()),
      relatedEntities.end());

  std::vector<stk::mesh::Entity *> sampleClosure;
  stk::mesh::find_closure(bulkData, relatedEntities, sampleClosure);

  // Keep only the closure, remove the rest, by decreasing entityRanks
  {
    const stk::mesh::Selector ownedOrShared = metaData.locally_owned_part() | metaData.globally_shared_part();
    typedef boost::indirect_iterator<std::vector<stk::mesh::Entity *>::const_iterator> EntityIterator;
    EntityIterator allKeepersEnd(sampleClosure.end());
    const EntityIterator allKeepersBegin(sampleClosure.begin());
    for (stk::mesh::EntityRank candidateRankCount = metaData.entity_rank_count(); candidateRankCount > 0; --candidateRankCount) {
      const stk::mesh::EntityRank candidateRank = candidateRankCount - 1;
      const EntityIterator keepersBegin = std::lower_bound(allKeepersBegin, allKeepersEnd,
                                                           stk::mesh::EntityKey(candidateRank, 0),
                                                           stk::mesh::EntityLess());
      const EntityIterator keepersEnd = allKeepersEnd;
      std::vector<stk::mesh::Entity *> candidates;
      stk::mesh::get_selected_entities(ownedOrShared, bulkData.buckets(candidateRank), candidates);
      {
        BulkModification modification(bulkData);
        std::set_difference(candidates.begin(), candidates.end(),
                            keepersBegin.base(), keepersEnd.base(),
                            EntityDestructor(modification),
                            stk::mesh::EntityLess());
      }
      allKeepersEnd = keepersBegin;
    }
  }
}

Teuchos::RCP<const Teuchos::Comm<Teuchos::Ordinal> > getComm(const Thyra::VectorSpaceBase<Scalar> &space) {
  typedef Thyra::SpmdVectorSpaceBase<Scalar> SVSB;
  const Teuchos::Ptr<const SVSB> space_downcasted = Teuchos::ptr_dynamic_cast<const SVSB>(Teuchos::ptrFromRef(space));
  if (Teuchos::nonnull(space_downcasted)) {
    return space_downcasted->getComm();
  } else {
    return Teuchos::null;
  }
}

  int
  MUMPS<Matrix,Vector>::solve_impl(
                         const Teuchos::Ptr<MultiVecAdapter<Vector> >  X,
                         const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
  {
    
    typedef FunctionMap<MUMPS,scalar_type> function_map;
    
    using Teuchos::as;
    
    const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
    const size_t nrhs = X->getGlobalNumVectors();
    
    const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
    
    xvals_.resize(val_store_size);
    bvals_.resize(val_store_size);
    
    #ifdef HAVE_AMESOS2_TIMERS
    Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
    Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
    #endif
    Util::get_1d_copy_helper<MultiVecAdapter<Vector>,
      slu_type>::do_get(B, bvals_(),as<size_t>(ld_rhs),
                        ROOTED);
 
    
    int ierr = 0; // returned error code
    mumps_par.nrhs = nrhs;
    mumps_par.lrhs = mumps_par.n;
    mumps_par.job = 3;

    if ( this->root_ ) 
      {                           
        mumps_par.rhs = bvals_.getRawPtr();
      }
    
    #ifdef HAVE_AMESOS2_TIMERS
    Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
    #endif
    
    function_map::mumps_c(&(mumps_par));
    MUMPS_ERROR();
    
    
    #ifdef HAVE_AMESOS2_TIMERS
    Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
    #endif
      
    Util::put_1d_data_helper<
    MultiVecAdapter<Vector>,slu_type>::do_put(X, bvals_(),
                                              as<size_t>(ld_rhs),
                                              ROOTED);
    
    return(ierr);
  }//end solve()

void LibmeshAdjacencies::getLibmeshAdjacencies<libMesh::Elem,libMesh::Node>(
    const Teuchos::Ptr<libMesh::Elem>& entity,
    Teuchos::Array<Teuchos::Ptr<libMesh::Node> >& adjacent_entities
    ) const
{
    int num_nodes = entity->n_nodes();
    adjacent_entities.resize( num_nodes );
    for ( int n = 0; n < num_nodes; ++n )
    {
    	adjacent_entities[n] = Teuchos::ptr( entity->get_node(n) );
    }
}

// Convert a Epetra_MultiVector with assumed block structure dictated by the
// vector space into a Thyra::MultiVectorBase object.
// const Teuchos::RCP<const Thyra::MultiVectorBase<double> > blockEpetraToThyra(const Epetra_MultiVector & e,const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & vs)
void blockEpetraToThyra(const Epetra_MultiVector & epetraX,const Teuchos::Ptr<Thyra::MultiVectorBase<double> > & thyraX) 
{
   TEUCHOS_ASSERT(thyraX->range()->dim()==epetraX.GlobalLength());

   // extract local information from the Epetra_MultiVector
   int leadingDim=0,numVectors=0,localDim=0;
   double * epetraData=0;
   epetraX.ExtractView(&epetraData,&leadingDim);

   numVectors = epetraX.NumVectors();

   blockEpetraToThyra(numVectors,epetraData,leadingDim,thyraX.ptr(),localDim);   

   TEUCHOS_ASSERT(localDim==epetraX.MyLength());
}

    void same_type_get_copy<MV>::apply(const Teuchos::Ptr<const MV>& mv,
				       const Teuchos::ArrayView<typename MV::scalar_t>& v,
				       const size_t ldx,
				       Teuchos::Ptr<const Tpetra::Map<typename MV::local_ordinal_t, typename MV::global_ordinal_t, typename MV::node_t> > distribution_map )
    {
      mv->get1dCopy(v, ldx, distribution_map);
    }

void initializePrec(
  const PreconditionerFactoryBase<Scalar> &precFactory,
  const Teuchos::RCP<const LinearOpBase<Scalar> > &fwdOp,
  const Teuchos::Ptr<PreconditionerBase<Scalar> > &prec,
  const ESupportSolveUse supportSolveUse = SUPPORT_SOLVE_UNSPECIFIED
  )
{
  precFactory.initializePrec(defaultLinearOpSource(fwdOp), prec.get(),
    supportSolveUse);
}

void PointValues_Evaluator<EvalT,TRAITST>::initialize(const Teuchos::RCP<const panzer::PointRule> & pointRule,
                                                      const Teuchos::Ptr<const ArrayT> & userArray,
                                                      // const Teuchos::Ptr<const Kokkos::DynRankView<double,PHX::Device> > & userArray,
                                                      const Teuchos::RCP<const panzer::PureBasis> & pureBasis)
{
  basis = pureBasis;

  if(userArray!=Teuchos::null && basis==Teuchos::null) 
    useBasisValuesRefArray = false;
  else if(userArray==Teuchos::null && basis!=Teuchos::null) 
    useBasisValuesRefArray = true;
  else {
    // this is a conflicting request, throw an exception
    TEUCHOS_ASSERT(false);
  }

  panzer::MDFieldArrayFactory af(pointRule->getName()+"_");
       
  // copy user array data
  if(userArray!=Teuchos::null) {
    TEUCHOS_ASSERT(userArray->rank()==2);
    refPointArray = Kokkos::DynRankView<double,PHX::Device>("refPointArray",userArray->dimension(0),userArray->dimension(1));
    // TEUCHOS_ASSERT(refPointArray.size()==userArray->size());
    for(int i=0;i<userArray->extent_int(0);i++)
      for(int j=0;j<userArray->extent_int(1);j++)
        refPointArray(i,j) = (*userArray)(i,j); 
  }

  // setup all fields to be evaluated and constructed
  pointValues.setupArrays(pointRule,af);

  // the field manager will allocate all of these field
  this->addEvaluatedField(pointValues.coords_ref);
  this->addEvaluatedField(pointValues.node_coordinates);
  this->addEvaluatedField(pointValues.jac);
  this->addEvaluatedField(pointValues.jac_inv);
  this->addEvaluatedField(pointValues.jac_det);
  this->addEvaluatedField(pointValues.point_coords);

  std::string n = "PointValues_Evaluator: " + pointRule->getName();
  this->setName(n);
}

//---------------------------------------------------------------------------//
// Get the coordinates of the entity nodes in canonical order.
void MoabHelpers::getEntityNodeCoordinates(
    const moab::EntityHandle& moab_entity,
    const Teuchos::Ptr<moab::ParallelComm>& moab_mesh,
    Teuchos::Array<double>& coordinates )
{
    const moab::EntityHandle* entity_nodes;
    int num_nodes = 0;
    std::vector<moab::EntityHandle> storage;
    DTK_CHECK_ERROR_CODE(
	moab_mesh->get_moab()->get_connectivity( moab_entity,
						 entity_nodes,
						 num_nodes,
						 false,
						 &storage )
	);

    coordinates.resize( 3 * num_nodes );
    DTK_CHECK_ERROR_CODE(
	moab_mesh->get_moab()->get_coords( entity_nodes,
					   num_nodes,
					   coordinates.getRawPtr() )
	);
}

void DiscreteBoundaryOperator<ValueType>::applyImpl(
    const Thyra::EOpTransp M_trans,
    const Thyra::MultiVectorBase<ValueType> &X_in,
    const Teuchos::Ptr<Thyra::MultiVectorBase<ValueType>> &Y_inout,
    const ValueType alpha, const ValueType beta) const {
  typedef Thyra::Ordinal Ordinal;

  // Note: the name is VERY misleading: these asserts don't disappear in
  // release runs, and in case of failure throw exceptions rather than
  // abort.
  TEUCHOS_ASSERT(this->opSupported(M_trans));
  TEUCHOS_ASSERT(X_in.range()->isCompatible(*this->domain()));
  TEUCHOS_ASSERT(Y_inout->range()->isCompatible(*this->range()));
  TEUCHOS_ASSERT(Y_inout->domain()->isCompatible(*X_in.domain()));

  const Ordinal colCount = X_in.domain()->dim();

  // Loop over the input columns

  for (Ordinal col = 0; col < colCount; ++col) {
    // Get access the the elements of X_in's and Y_inout's column #col
    Thyra::ConstDetachedSpmdVectorView<ValueType> xVec(X_in.col(col));
    Thyra::DetachedSpmdVectorView<ValueType> yVec(Y_inout->col(col));
    const Teuchos::ArrayRCP<const ValueType> xArray(xVec.sv().values());
    const Teuchos::ArrayRCP<ValueType> yArray(yVec.sv().values());

    // Wrap the Trilinos array in an Armadillo vector. const_cast is used
    // because it's more natural to have a const arma::Col<ValueType> array
    // than an arma::Col<const ValueType> one.
    const arma::Col<ValueType> xCol(const_cast<ValueType *>(xArray.get()),
                                    xArray.size(), false /* copy_aux_mem */);
    arma::Col<ValueType> yCol(yArray.get(), yArray.size(), false);

    applyBuiltInImpl(static_cast<TranspositionMode>(M_trans), xCol, yCol, alpha,
                     beta);
  }
}

//---------------------------------------------------------------------------//
// Given an entity, get the ids of its support locations
void LibmeshNodalShapeFunction::entitySupportIds( 
    const DataTransferKit::Entity& entity,
    Teuchos::Array<DataTransferKit::SupportId>& support_ids ) const
{
    // Node case.
    if ( 0 == entity.topologicalDimension() )
    {
	DTK_CHECK( extractGeom<libMesh::Node>(entity)->valid_id() );
	support_ids.assign( 1, extractGeom<libMesh::Node>(entity)->id() );
    }

    // Element case.
    else
    {
	Teuchos::Ptr<libMesh::Elem> elem = extractGeom<libMesh::Elem>(entity);
	int num_nodes = elem->n_nodes();
	support_ids.resize( num_nodes );
	for ( int n = 0; n < num_nodes; ++n )
	{
	    DTK_CHECK( elem->get_node(n)->valid_id() );
	    support_ids[n] = elem->get_node(n)->id();
	}
    }
}

void MultiVecAdapter<Epetra_MultiVector>::get1dCopy(
  const Teuchos::ArrayView<MultiVecAdapter<Epetra_MultiVector>::scalar_t>& av,
  size_t lda,
  Teuchos::Ptr<
    const Tpetra::Map<MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
                      MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
                      MultiVecAdapter<Epetra_MultiVector>::node_t> > distribution_map ) const
{
  using Teuchos::rcpFromPtr;
  using Teuchos::as;

  const size_t num_vecs = getGlobalNumVectors();
  
#ifdef HAVE_AMESOS2_DEBUG
  const size_t requested_vector_length = distribution_map->getNodeNumElements();
  TEUCHOS_TEST_FOR_EXCEPTION( lda < requested_vector_length,
		      std::invalid_argument,
		      "Given stride is not large enough for local vector length" );
  TEUCHOS_TEST_FOR_EXCEPTION( as<size_t>(av.size()) < (num_vecs-1) * lda + requested_vector_length,
		      std::invalid_argument,
		      "MultiVector storage not large enough given leading dimension "
		      "and number of vectors" );
#endif

  // Optimization for ROOTED
  if ( num_vecs == 1 && mv_->Comm().MyPID() == 0 && mv_->Comm().NumProc() == 1 ) {
	  mv_->ExtractCopy(av.getRawPtr(), lda);
  }
  else {
  Epetra_Map e_dist_map
    = *Util::tpetra_map_to_epetra_map<local_ordinal_t,
                                      global_ordinal_t,
                                      global_size_t,
                                      node_t>(*distribution_map);

  multivec_t redist_mv(e_dist_map, as<int>(num_vecs));
  const Epetra_Import importer(e_dist_map, *mv_map_); // Note, target/source order is reversed in Tpetra
  redist_mv.Import(*mv_, importer, Insert);

  // Finally, do copy
  redist_mv.ExtractCopy(av.getRawPtr(), lda);
  }

}

void ObserverImpl::observeSolution (
  double stamp, const Epetra_Vector& nonOverlappedSolution,
  const Teuchos::Ptr<const Epetra_Vector>& nonOverlappedSolutionDot)
{
  // If solution == "Steady" or "Continuation", we need to update the solution
  // from the initial guess prior to writing it out, or we will not get the
  // proper state of things like "Stress" in the Exodus file.
  {
    // Evaluate state field manager
    if(nonOverlappedSolutionDot != Teuchos::null)
      app_->evaluateStateFieldManager(stamp, nonOverlappedSolutionDot.get(),
                                      NULL, nonOverlappedSolution);
    else
      app_->evaluateStateFieldManager(stamp, NULL, NULL, nonOverlappedSolution);

    // Renames the New state as the Old state in preparation for the next step
    app_->getStateMgr().updateStates();

#ifdef ALBANY_PERIDIGM
#if defined(ALBANY_EPETRA)
    const Teuchos::RCP<LCM::PeridigmManager>&
      peridigmManager = LCM::PeridigmManager::self();
    if (Teuchos::nonnull(peridigmManager)) {
      peridigmManager->writePeridigmSubModel(stamp);
      peridigmManager->updateState();
    }
#endif
#endif
  }

  //! update distributed parameters in the mesh
  Teuchos::RCP<DistParamLib> distParamLib = app_->getDistParamLib();
  distParamLib->scatter();
  DistParamLib::const_iterator it;
  Teuchos::RCP<const Epetra_Comm> comm = app_->getEpetraComm();
  for(it = distParamLib->begin(); it != distParamLib->end(); ++it) {
    app_->getDiscretization()->setFieldT(*it->second->overlapped_vector(), it->second->name(),
                                        /*overlapped*/ true);
  }

  StatelessObserverImpl::observeSolution(stamp, nonOverlappedSolution,
                                         nonOverlappedSolutionDot);
}

void
MultiVecAdapter<
MultiVector<Scalar,
            LocalOrdinal,
            GlobalOrdinal,
            Node> >::get1dCopy(const Teuchos::ArrayView<scalar_t>& av,
                               size_t lda,
                               Teuchos::Ptr<
                               const Tpetra::Map<LocalOrdinal,
                               GlobalOrdinal,
                               Node> > distribution_map ) const
{
    using Teuchos::rcpFromPtr;
    using Teuchos::as;

    size_t num_vecs = getGlobalNumVectors();

#ifdef HAVE_AMESOS2_DEBUG
    size_t requested_vector_length = distribution_map->getNodeNumElements();
    TEUCHOS_TEST_FOR_EXCEPTION( lda < requested_vector_length,
                                std::invalid_argument,
                                "Given stride is not large enough for local vector length" );
    TEUCHOS_TEST_FOR_EXCEPTION( as<size_t>(av.size()) < as<size_t>((num_vecs-1) * lda + requested_vector_length),
                                std::invalid_argument,
                                "MultiVector storage not large enough given leading dimension "
                                "and number of vectors" );
#endif

    multivec_t redist_mv(rcpFromPtr(distribution_map), num_vecs);

    typedef Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node> import_type;
    import_type importer (this->getMap (), rcpFromPtr (distribution_map));
    redist_mv.doImport (*mv_, importer, Tpetra::REPLACE);

    // do copy
    redist_mv.get1dCopy (av, lda);
}

void MueLuTpetraPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::initializePrec(
  const Teuchos::RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
  PreconditionerBase<Scalar> *prec,
  const ESupportSolveUse supportSolveUse
  ) const
{
  // Check precondition

  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
  TEUCHOS_ASSERT(prec);

  Teuchos::Time totalTimer(""), timer("");
  totalTimer.start(true);

  const RCP<Teuchos::FancyOStream> out = this->getOStream();
  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  Teuchos::OSTab tab(out);
  if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
    *out << "\nEntering Thyra::MueLuTpetraPreconditionerFactory::initializePrec(...) ...\n";
  }

  // Retrieve wrapped concrete Tpetra matrix from FwdOp

  const Teuchos::RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc->getOp();
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));

  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyraTpetraLinOp;
  const Teuchos::RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = Teuchos::rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraTpetraFwdOp));

  typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinOp;
  const Teuchos::RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));

  typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraCrsMat;
  const Teuchos::RCP<const TpetraCrsMat> tpetraFwdCrsMat = Teuchos::rcp_dynamic_cast<const TpetraCrsMat>(tpetraFwdOp);
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdCrsMat));

  // Retrieve concrete preconditioner object

  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec =
    Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));

  if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
    *out << "\nCreating a new MueLu::TpetraOperator object...\n";
  }
  timer.start(true);

  // Workaround since MueLu interface does not accept const matrix as input
  const Teuchos::RCP<TpetraCrsMat> tpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<TpetraCrsMat>(tpetraFwdCrsMat);

  // Create and compute the initial preconditioner
  
  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> MV;
  Teuchos::RCP<MV> coords, null_space;
  if(paramList_->isType<Teuchos::RCP<MV> >("Coordinates")) coords = paramList_->get<Teuchos::RCP<MV> >("Coordinates");
  if(paramList_->isType<Teuchos::RCP<MV> >("Null Space")) null_space = paramList_->get<Teuchos::RCP<MV> >("Null Space");

  typedef MueLu::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> MueLuOperator;
  const Teuchos::RCP<MueLuOperator> mueluPrecOp = MueLu::CreateTpetraPreconditioner(tpetraFwdCrsMatNonConst, *paramList_,coords,null_space);

  timer.stop();
  if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
    Teuchos::OSTab(out).o() << "> Creation time = " << timer.totalElapsedTime() << " sec\n";
  }

  const Teuchos::RCP<LinearOpBase<Scalar> > thyraPrecOp = Thyra::createLinearOp(Teuchos::RCP<TpetraLinOp>(mueluPrecOp));
  defaultPrec->initializeUnspecified(thyraPrecOp);

  totalTimer.stop();
  if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
    *out << "\nTotal time in Thyra::MueLuTpetraPreconditionerFactory::initializePrec(...) = " << totalTimer.totalElapsedTime() << " sec\n";
  }

  if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
    *out << "\nLeaving Thyra::MueLuTpetraPreconditionerFactory::initializePrec(...) ...\n";
  }
}

void panzer::EquationSet_DefaultImpl<EvalT>::
buildAndRegisterDOFProjectionsToIPEvaluators(PHX::FieldManager<panzer::Traits>& fm,
                                             const panzer::FieldLayoutLibrary& fl,
                                             const Teuchos::RCP<panzer::IntegrationRule>& ir,
                                             const Teuchos::Ptr<const panzer::LinearObjFactory<panzer::Traits> > & lof,
                                             const Teuchos::ParameterList& user_data) const
{
  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;

  Teuchos::RCP<const panzer::UniqueGlobalIndexerBase> globalIndexer;
  if(lof!=Teuchos::null) 
    globalIndexer = lof->getUniqueGlobalIndexerBase();
  
  // DOFs: Scalar value @ basis --> Scalar value @ IP 
  for (DescriptorIterator dof_iter = m_provided_dofs_desc.begin(); dof_iter != m_provided_dofs_desc.end(); ++dof_iter) {

    ParameterList p;
    p.set("Name", dof_iter->first);
    p.set("Basis", fl.lookupLayout(dof_iter->first));
    p.set("IR", ir);

    if(globalIndexer!=Teuchos::null) {
      // build the offsets for this field
      int fieldNum = globalIndexer->getFieldNum(dof_iter->first);
      RCP<const std::vector<int> > offsets = 
          rcp(new std::vector<int>(globalIndexer->getGIDFieldOffsets(m_block_id,fieldNum)));
      p.set("Jacobian Offsets Vector", offsets);
    }
    // else default to the slow DOF call
    
    RCP< PHX::Evaluator<panzer::Traits> > op = 
      rcp(new panzer::DOF<EvalT,panzer::Traits>(p));
    
    this->template registerEvaluator<EvalT>(fm, op);
  }

  // Gradients of DOFs: Scalar value @ basis --> Vector value @ IP

  for(typename std::map<std::string,DOFDescriptor>::const_iterator itr=m_provided_dofs_desc.begin();
      itr!=m_provided_dofs_desc.end();++itr) {
    
    if(itr->second.basis->supportsGrad()) {
      
      // is gradient required for this variable
      if(!itr->second.grad.first) 
        continue; // its not required, quit the loop

      const std::string dof_name =      itr->first;
      const std::string dof_grad_name = itr->second.grad.second;

      ParameterList p;
      p.set("Name", dof_name);
      p.set("Gradient Name", dof_grad_name);
      p.set("Basis", fl.lookupLayout(dof_name)); 
      p.set("IR", ir);
      
      RCP< PHX::Evaluator<panzer::Traits> > op = 
        rcp(new panzer::DOFGradient<EvalT,panzer::Traits>(p));

      this->template registerEvaluator<EvalT>(fm, op);
    }
  }

  // Curl of DOFs: Vector value @ basis --> Vector value @ IP (3D) or Scalar value @ IP (2D)

  for(typename std::map<std::string,DOFDescriptor>::const_iterator itr=m_provided_dofs_desc.begin();
      itr!=m_provided_dofs_desc.end();++itr) {
    
    if(itr->second.basis->supportsCurl()) {

      // is curl required for this variable
      if(!itr->second.curl.first) 
        continue; // its not required, quit the loop

      const std::string dof_name =      itr->first;
      const std::string dof_curl_name = itr->second.curl.second;

      ParameterList p;
      p.set("Name", dof_name);
      p.set("Curl Name", dof_curl_name);
      p.set("Basis", fl.lookupLayout(dof_name)); 
      p.set("IR", ir);

      // this will help accelerate the DOFCurl evaluator when Jacobians are needed
      if(globalIndexer!=Teuchos::null) {
        // build the offsets for this field
        int fieldNum = globalIndexer->getFieldNum(dof_name);
        RCP<const std::vector<int> > offsets = 
            rcp(new std::vector<int>(globalIndexer->getGIDFieldOffsets(m_block_id,fieldNum)));
        p.set("Jacobian Offsets Vector", offsets);
      }
      // else default to the slow DOF call
    
      
      RCP< PHX::Evaluator<panzer::Traits> > op = 
        rcp(new panzer::DOFCurl<EvalT,panzer::Traits>(p));

      this->template registerEvaluator<EvalT>(fm, op);
//.........这里部分代码省略.........

    MPISession::MPISession (Teuchos::Ptr<int> argc, Teuchos::Ptr<char**> argv) 
    {
#ifdef EPETRA_MPI 
      MPI_Init (argc.getRawPtr(), argv.getRawPtr());
#endif // EPETRA_MPI 
    }