/** \brief Constructor.

      @param[in]     parlist is a parameter list specifying inputs

      parlist should contain sublists "SOL"->"Risk Measure"->"Convex Combination Risk Measure" and
      within the "Convex Combination Risk Measure" sublist should have the following parameters
      \li "Convex Combination Parameters" (greater than 0 and sum to 1)
      \li Sublists labeled 1 to n with risk measure definitions.
  */
  ConvexCombinationRiskMeasure(Teuchos::ParameterList &parlist)
    : RiskMeasure<Real>(), size_(0), firstReset_(true) {
    Teuchos::ParameterList &list
      = parlist.sublist("SOL").sublist("Risk Measure").sublist("Convex Combination Risk Measure");
    // Get convex combination parameters
    Teuchos::Array<Real> lambda
      = Teuchos::getArrayFromStringParameter<Real>(list,"Convex Combination Parameters");
    lambda_ = lambda.toVector();
    size_ = lambda_.size();
    // Build risk measures
    risk_.clear(); risk_.resize(size_,Teuchos::null);
    parlist_.clear(); parlist_.resize(size_);
    for (uint i = 0; i < size_; ++i) {
      std::ostringstream convert;
      convert << i;
      std::string si = convert.str();
      Teuchos::ParameterList &ilist = list.sublist(si);
      std::string name = ilist.get<std::string>("Name");
      parlist_[i].sublist("SOL").sublist("Risk Measure").set("Name",name);
      parlist_[i].sublist("SOL").sublist("Risk Measure").sublist(name) = ilist;
      risk_[i] = RiskMeasureFactory<Real>(parlist_[i]);
    }
    // Check inputs
    checkInputs();
  }

Adapter *
Adapter::initialize(const Teuchos::RCP<Teuchos::ParameterList>& catalystParams)
{
  // Validate parameters against list for this specific class
  catalystParams->validateParameters(*getValidAdapterParameters(),0);

  if (Adapter::instance)
    delete Adapter::instance;
  Adapter::instance = new Adapter();

  // Register our Grid class with Catalyst so that it can be used in a pipeline.
  if (vtkClientServerInterpreterInitializer *intInit =
      vtkClientServerInterpreterInitializer::GetInitializer()) {
    if (vtkClientServerInterpreter *interp = intInit->GetGlobalInterpreter()) {
      interp->AddNewInstanceFunction("Grid", Private::MakeGrid);
    }
  }

  // Load pipeline file
  Teuchos::Array<std::string> files =
      catalystParams->get<Teuchos::Array<std::string> >("Pipeline Files");
  typedef Teuchos::Array<std::string>::const_iterator FileIterT;
  for (FileIterT it = files.begin(), itEnd = files.end(); it != itEnd; ++it)
    Adapter::instance->addPythonScriptPipeline(*it);

  return Adapter::instance;
}

Teuchos::Array<int> getMyBlockLIDs(
    const Teuchos::ParameterList &blockingParams,
    const Albany::AbstractDiscretization &disc)
{
  Teuchos::Array<int> result;

  const std::string nodeSetLabel = blockingParams.get<std::string>("Node Set");
  const int dofRank = blockingParams.get<int>("Dof");

  const Albany::NodeSetList &nodeSets = disc.getNodeSets();
  const Albany::NodeSetList::const_iterator it = nodeSets.find(nodeSetLabel);
  TEUCHOS_TEST_FOR_EXCEPT(it == nodeSets.end());
  {
    typedef Albany::NodeSetList::mapped_type NodeSetEntryList;
    const NodeSetEntryList &nodeEntries = it->second;

    for (NodeSetEntryList::const_iterator jt = nodeEntries.begin(); jt != nodeEntries.end(); ++jt) {
      typedef NodeSetEntryList::value_type NodeEntryList;
      const NodeEntryList &entries = *jt;
      result.push_back(entries[dofRank]);
    }
  }

  return result;
}

void FunctionTests::checkFunctionsAgree(FunctionPtr f1, FunctionPtr f2, BasisCachePtr basisCache) {
  ASSERT_EQ(f1->rank(), f2->rank())
    << "f1->rank() " << f1->rank() << " != f2->rank() " << f2->rank() << endl;

  int rank = f1->rank();
  int numCells = basisCache->getPhysicalCubaturePoints().dimension(0);
  int numPoints = basisCache->getPhysicalCubaturePoints().dimension(1);
  int spaceDim = basisCache->getPhysicalCubaturePoints().dimension(2);
  Teuchos::Array<int> dim;
  dim.append(numCells);
  dim.append(numPoints);
  for (int i=0; i<rank; i++) {
    dim.append(spaceDim);
  }
  FieldContainer<double> f1Values(dim);
  FieldContainer<double> f2Values(dim);
  f1->values(f1Values,basisCache);
  f2->values(f2Values,basisCache);
  
  double tol = 1e-14;
  double maxDiff;
  EXPECT_TRUE(fcsAgree(f1Values,f2Values,tol,maxDiff))
    << "Test failed: f1 and f2 disagree; maxDiff " << maxDiff << ".\n"
    << "f1Values: \n" << f1Values
    << "f2Values: \n" << f2Values;
}

// ============================================================================
Teuchos::RCP<Recti::Domain::Abstract>
Recti::Domain::Factory::
buildPolygon() const
{
    const Teuchos::ParameterList & pointsList =
            pList_.sublist("Vertices");
    
    Teuchos::Array<Point> vertices;
    
    Teuchos::ParameterList::ConstIterator k;
    for ( k=pointsList.begin(); k!=pointsList.end(); ++k )
    {
        Teuchos::ParameterEntry e = pointsList.entry(k);
        TEUCHOS_ASSERT( e.isList() );
        
        std::string label = pointsList.name(k);
        const Teuchos::ParameterList & coordinates =
                pointsList.sublist(label);
        
        Point X = Teuchos::tuple( coordinates.get<double>("x"),
                                  coordinates.get<double>("y"),
                                  0.0
                                );
        vertices.push_back( X );
    }
    
    return Teuchos::rcp( new Polygon( vertices ) );
}

  Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
  createMeshMap (const LO& blockSize, const Tpetra::Map<LO,GO,Node>& pointMap)
  {
    typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
    typedef Tpetra::Map<LO,GO,Node> map_type;

    //calculate mesh GIDs
    Teuchos::ArrayView<const GO> pointGids = pointMap.getNodeElementList();
    Teuchos::Array<GO> meshGids;
    GO indexBase = pointMap.getIndexBase();

    // Use hash table to track whether we've encountered this GID previously.  This will happen
    // when striding through the point DOFs in a block.  It should not happen otherwise.
    // I don't use sort/make unique because I don't want to change the ordering.
    meshGids.reserve(pointGids.size());
    Tpetra::Details::HashTable<GO,int> hashTable(pointGids.size());
    for (int i=0; i<pointGids.size(); ++i) {
      GO meshGid = (pointGids[i]-indexBase) / blockSize + indexBase;
      if (hashTable.get(meshGid) == -1) {
        hashTable.add(meshGid,1);   //(key,value)
        meshGids.push_back(meshGid);
      }
    }

    Teuchos::RCP<const map_type> meshMap = Teuchos::rcp( new map_type(TOT::invalid(), meshGids(), 0, pointMap.getComm()) );
    return meshMap;

  }

Teuchos::RCP<const Tpetra_Map>
Albany::SolutionResponseFunction::
buildCulledMapT(const Tpetra_Map& x_mapT,
	       const Teuchos::Array<int>& keepDOF) const
{
  int numKeepDOF = std::accumulate(keepDOF.begin(), keepDOF.end(), 0);
  int Neqns = keepDOF.size();
  int N = x_mapT.getNodeNumElements(); // x_mapT is map for solution vector

  TEUCHOS_ASSERT( !(N % Neqns) ); // Assume that all the equations for
                                  // a given node are on the assigned
                                  // processor. I.e. need to ensure
                                  // that N is exactly Neqns-divisible

  int nnodes = N / Neqns;          // number of fem nodes
  int N_new = nnodes * numKeepDOF; // length of local x_new

  Teuchos::ArrayView<const GO> gidsT = x_mapT.getNodeElementList();
  Teuchos::Array<GO> gids_new(N_new);
  int idx = 0;
  for ( int inode = 0; inode < N/Neqns ; ++inode) // For every node
    for ( int ieqn = 0; ieqn < Neqns; ++ieqn )  // Check every dof on the node
      if ( keepDOF[ieqn] == 1 )  // then want to keep this dof
	gids_new[idx++] = gidsT[(inode*Neqns)+ieqn];
  // end cull

  Teuchos::RCP<const Tpetra_Map> x_new_mapT = Tpetra::createNonContigMapWithNode<LO, GO, KokkosNode> (gids_new, x_mapT.getComm(), x_mapT.getNode());

  return x_new_mapT;

}

//---------------------------------------------------------------------------//
// Get the parameteric center of an entity.
void MoabEntityLocalMap::parametricCenter(
    const Entity& entity,
    Teuchos::Array<double>& center ) const
{
    moab::EntityType moab_type = d_moab_mesh->get_moab()->type_from_handle(
	MoabHelpers::extractEntity(entity) );

    switch( moab_type )
    {
	case moab::MBTRI:
	    center.assign( 2, 1.0 / 3.0 );
	    break;
	case moab::MBQUAD:
	    center.assign( 2, 0.0 );
	    break;
	case moab::MBTET:
	    center.assign( 3, 1.0 / 6.0 );
	    break;
	case moab::MBHEX:
	    center.assign( 3, 0.0 );
	    break;
	default:
	    center.resize( 0 );
	    break;
    }
}

// =============================================================================
Teuchos::RCP<Tpetra::Map<ORD> >
VIO::TpetraMesh::Mesh::
getElemsToNodesMap_() const
{
  // create list of elements that need to be accessible from this process


  // Make sure that *all entries that belong to any of the elements in
  // this core are accessible.
  // First mark all the nodes that need to be accessible:
  TEUCHOS_ASSERT( !elems_.is_null() );
  TEUCHOS_ASSERT( !nodes_.is_null() );
  int numNodes = nodes_.size();
  Teuchos::Array<bool> mustBeAccessible( numNodes );
  for ( int k=0; k<elems_.size(); k++ )
      for ( int l=0; l<elems_[k].size(); l++ )
          mustBeAccessible[ elems_[k][l] ] = true;
  // now create the list
  Teuchos::Array<ORD> entryList;
  for ( int k=0; k<numNodes; k++ )
      if ( mustBeAccessible[k] )
          entryList.append( k );

  Teuchos::RCP<Tpetra::Map<ORD> > map =
      Teuchos::rcp( new Tpetra::Map<ORD>( Teuchos::OrdinalTraits<ORD>::invalid(), entryList(), 0, comm_ )
                  );

  return map;
}

/*!
 * \brief Find the set of entities a point neighbors.
 */
void CoarseLocalSearch::search( const Teuchos::ArrayView<const double>& point,
				const Teuchos::ParameterList& parameters,
				Teuchos::Array<Entity>& neighbors ) const
{
    // Find the leaf of nearest neighbors.
    int num_neighbors = 100;
    if ( parameters.isParameter("Coarse Local Search kNN") )
    {
	num_neighbors = parameters.get<int>("Coarse Local Search kNN");
    }
    num_neighbors = 
	std::min( num_neighbors, Teuchos::as<int>(d_entity_map.size()) );
    Teuchos::Array<unsigned> local_neighbors = 
	d_tree->nnSearch( point, num_neighbors );

    // Extract the neighbors.
    neighbors.resize( local_neighbors.size() );
    Teuchos::Array<unsigned>::const_iterator local_it;
    Teuchos::Array<Entity>::iterator entity_it;
    for ( local_it = local_neighbors.begin(),
	 entity_it = neighbors.begin();
	  local_it != local_neighbors.end();
	  ++local_it, ++entity_it )
    {
	DTK_CHECK( d_entity_map.count(*local_it) );
	*entity_it = d_entity_map.find(*local_it)->second;
    }
}

Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > create_tridiag_graph(LocalOrdinal num_rows_per_proc)
{
  Teuchos::RCP<const Teuchos::Comm<int> > comm = getDefaultComm();

  const global_size_t INVALID = Teuchos::OrdinalTraits<global_size_t>::invalid();
  const LocalOrdinal indexBase = 0;

  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowmap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(INVALID, num_rows_per_proc, indexBase, comm));

  Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > crsgraph = Teuchos::rcp(new Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>(rowmap, 0));

  Teuchos::Array<GlobalOrdinal> cols;
 for(LocalOrdinal l_row = 0; (size_t) l_row<rowmap->getNodeNumElements(); l_row++) {
    GlobalOrdinal g_row = rowmap->getGlobalElement(l_row);
    if (g_row == rowmap->getMinAllGlobalIndex() ||
        g_row == rowmap->getMaxAllGlobalIndex()) cols.resize(2);
    else cols.resize(3);

    size_t coloffset = 0;
    if (g_row > rowmap->getMinAllGlobalIndex()) cols[coloffset++] = g_row-1;
    cols[coloffset++] = g_row;
    if (g_row < rowmap->getMaxAllGlobalIndex()) cols[coloffset++] = g_row+1;

    crsgraph->insertGlobalIndices(g_row, cols());
  }

  crsgraph->fillComplete();

  return crsgraph;
}

void FELIX::Hydrology::constructNeumannEvaluators (const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
{
    // Note: we only enter this function if sidesets are defined in the mesh file
    // i.e. meshSpecs.ssNames.size() > 0

    Albany::BCUtils<Albany::NeumannTraits> nbcUtils;

    // Check to make sure that Neumann BCs are given in the input file
    if (!nbcUtils.haveBCSpecified(this->params))
    {
        return;
    }

    // Construct BC evaluators for all side sets and names
    // Note that the string index sets up the equation offset, so ordering is important

    std::vector<std::string> neumannNames(neq + 1);
    Teuchos::Array<Teuchos::Array<int> > offsets;
    offsets.resize(1);

    neumannNames[0] = "Hydraulic Potential";
    neumannNames[1] = "all";
    offsets[0].resize(1);
    offsets[0][0] = 0;

    // Construct BC evaluators for all possible names of conditions
    std::vector<std::string> condNames(1);
    condNames[0] = "neumann";

    nfm.resize(1); // FELIX problem only has one element block

    nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, neumannNames, dof_names, false, 0,
                                            condNames, offsets, dl,
                                            this->params, this->paramLib);
}

 void values(FieldContainer<double> &values, BasisCachePtr sliceBasisCache) {
   vector<GlobalIndexType> sliceCellIDs = sliceBasisCache->cellIDs();
   
   Teuchos::Array<int> dim;
   values.dimensions(dim);
   dim[0] = 1; // one cell
   Teuchos::Array<int> offset(dim.size());
   
   for (int cellOrdinal = 0; cellOrdinal < sliceCellIDs.size(); cellOrdinal++) {
     offset[0] = cellOrdinal;
     int enumeration = values.getEnumeration(offset);
     FieldContainer<double>valuesForCell(dim,&values[enumeration]);
     GlobalIndexType sliceCellID = sliceCellIDs[cellOrdinal];
     int numPoints = sliceBasisCache->getPhysicalCubaturePoints().dimension(1);
     int spaceDim = sliceBasisCache->getPhysicalCubaturePoints().dimension(2);
     FieldContainer<double> spaceTimePhysicalPoints(1,numPoints,spaceDim+1);
     for (int ptOrdinal=0; ptOrdinal<numPoints; ptOrdinal++) {
       for (int d=0; d<spaceDim; d++) {
         spaceTimePhysicalPoints(0,ptOrdinal,d) = sliceBasisCache->getPhysicalCubaturePoints()(cellOrdinal,ptOrdinal,d);
       }
       spaceTimePhysicalPoints(0,ptOrdinal,spaceDim) = _t;
     }
     
     GlobalIndexType cellID = _cellIDMap[sliceCellID];
     BasisCachePtr spaceTimeBasisCache = BasisCache::basisCacheForCell(_spaceTimeMesh, cellID);
     
     FieldContainer<double> spaceTimeRefPoints(1,numPoints,spaceDim+1);
     CamelliaCellTools::mapToReferenceFrame(spaceTimeRefPoints, spaceTimePhysicalPoints, _spaceTimeMesh, cellID);
     spaceTimeRefPoints.resize(numPoints,spaceDim+1);
     spaceTimeBasisCache->setRefCellPoints(spaceTimeRefPoints);
     _spaceTimeFunction->values(valuesForCell, spaceTimeBasisCache);
   }
 }

void BilinearFormUtility::weightCellBasisValues(FieldContainer<double> &basisValues, const FieldContainer<double> &weights, int offset) {
  // weights are (numCells, offset+numFields)
  // basisValues are (numCells, numFields, ...)
  int numCells = basisValues.dimension(0);
  int numFields = basisValues.dimension(1);
  
  Teuchos::Array<int> dimensions;
  basisValues.dimensions(dimensions);

  int numAffectedValues = 1;
  for (int dimIndex=2; dimIndex<dimensions.size(); dimIndex++) {
    numAffectedValues *= dimensions[dimIndex];
  }
  
  Teuchos::Array<int> index(dimensions.size(),0);
  
  for (int cellIndex=0; cellIndex < numCells; cellIndex++) {
    index[0] = cellIndex;
    for (int fieldIndex=0; fieldIndex < numFields; fieldIndex++) {
      index[1] = fieldIndex;
      int enumIndex = basisValues.getEnumeration(index);
      for (int valIndex=enumIndex; valIndex < numAffectedValues + enumIndex; valIndex++) {
        basisValues[valIndex] *= weights(cellIndex,fieldIndex+offset);
      }
    }
  }
  
}

void
Stokhos::PecosOneDOrthogPolyBasis<ordinal_type,value_type>::
getQuadPoints(ordinal_type quad_order,
	      Teuchos::Array<value_type>& quad_points,
	      Teuchos::Array<value_type>& quad_weights,
	      Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
{
  // This assumes Gaussian quadrature
  ordinal_type num_points = 
    static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
  const Pecos::RealArray& gp = pecosPoly->collocation_points(num_points);
  const Pecos::RealArray& gw = pecosPoly->type1_collocation_weights(num_points); 
  quad_points.resize(num_points);
  quad_weights.resize(num_points);
  for (ordinal_type i=0; i<num_points; i++) {
    quad_points[i] = gp[i];
    quad_weights[i] = gw[i];
  }
  
  // Evalute basis at gauss points
  quad_values.resize(num_points);
  for (ordinal_type i=0; i<num_points; i++) {
    quad_values[i].resize(p+1);
    evaluateBases(quad_points[i], quad_values[i]);
  }
}

 RiskVector( Teuchos::ParameterList &parlist,
       const Teuchos::RCP<Vector<Real> > &vec,
       const Real stat = 1. )
   : vec_(vec), augmented_(false), nStat_(0) {
   stat_.clear();
   dual_vec1_ = vec->dual().clone();
   std::string type = parlist.sublist("SOL").sublist("Risk Measure").get("Name","CVaR");
   if ( type == "CVaR"                       ||
        type == "HMCR"                       ||
        type == "KL Divergence"              ||
        type == "Moreau-Yosida CVaR"         ||
        type == "Log-Exponential Quadrangle" ||
        type == "Log-Quantile Quadrangle"    ||
        type == "Quantile-Based Quadrangle"  ||
        type == "Truncated Mean Quadrangle" ) {
     augmented_ = true;
     nStat_     = 1;
     stat_.resize(nStat_,stat);
   }
   else if ( type == "Mixed-Quantile Quadrangle" ) {
     Teuchos::ParameterList &list
       = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mixed-Quantile Quadrangle");
     Teuchos::Array<Real> prob
       = Teuchos::getArrayFromStringParameter<Real>(list,"Probability Array");
     augmented_ = true;
     nStat_     = prob.size();
     stat_.resize(nStat_,stat);
   }
   else if ( type == "Quantile-Radius Quadrangle" ) {
     augmented_ = true;
     nStat_     = 2;
     stat_.resize(nStat_,stat);
   }
 }

void BCManager::modifyPrimalSystem(
    std::string const& val, int offset, std::string const& set)
{
  // should we fill in BC info?
  bool fillRes = (res != Teuchos::null);
  bool fillJac = (jac != Teuchos::null);
  if ((!fillRes) && (!fillJac)) return;

  // get views of the solution and residual vectors
  ArrayRCP<const ST> x_const_view;
  ArrayRCP<ST> f_nonconst_view;
  if (fillRes) {
    x_const_view = sol->get1dView();
    f_nonconst_view = res->get1dViewNonConst();
  }

  // set up some data for replacing jacobian values
  Teuchos::Array<LO> index(1);
  Teuchos::Array<ST> value(1);
  value[0] = 1.0;
  size_t numEntries;
  Teuchos::Array<ST> matrixEntries;
  Teuchos::Array<LO> matrixIndices;

  // loop over all of the nodes in this node set
  std::vector<Albany::GOALNode> nodes = ns[set];
  for (int i=0; i < nodes.size(); ++i)
  {
    Albany::GOALNode node = nodes[i];
    int lunk = disc->getDOF(node.lid, offset);

    // if the node is higher order, we set the value of the DBC to be 0
    // if it is not, we parse the expression from the input file to get the value
    // note: this assumes that bcs are either constant or linear in space
    // anything else would require a linear solve to find coefficients v
    double v = 0.0;
    if (!node.higherOrder)
      v = parseExpression(val, node.coord[0], node.coord[1], node.coord[2], t);

    // modify the residual if necessary
    if (fillRes)
      f_nonconst_view[lunk] = x_const_view[lunk] - v;

    // modify the jacobian if necessary
    if (fillJac)
    {
      index[0] = lunk;
      numEntries = jac->getNumEntriesInLocalRow(lunk);
      matrixEntries.resize(numEntries);
      matrixIndices.resize(numEntries);
      jac->getLocalRowCopy(lunk, matrixIndices(), matrixEntries(), numEntries);
      for (int i=0; i < numEntries; ++i)
        matrixEntries[i] = 0.0;
      jac->replaceLocalValues(lunk, matrixIndices(), matrixEntries());
      jac->replaceLocalValues(lunk, index(), value());
    }
  }
}

FELIX::StokesFO::
StokesFO( const Teuchos::RCP<Teuchos::ParameterList>& params_,
             const Teuchos::RCP<ParamLib>& paramLib_,
             const int numDim_) :
  Albany::AbstractProblem(params_, paramLib_, numDim_),
  numDim(numDim_)
{
  //Set # of PDEs per node based on the Equation Set.  
  //Equation Set is FELIX by default (2 dofs / node -- usual FELIX Stokes FO).
  std::string eqnSet = params_->sublist("Equation Set").get<std::string>("Type", "FELIX"); 
  if (eqnSet == "FELIX") 
    neq = 2; //FELIX FO Stokes system is a system of 2 PDEs
  else if (eqnSet == "Poisson" || eqnSet == "FELIX X-Z") 
    neq = 1; //1 PDE/node for Poisson or FELIX X-Z physics


  // Set the num PDEs for the null space object to pass to ML
  this->rigidBodyModes->setNumPDEs(neq);

  // the following function returns the problem information required for setting the rigid body modes (RBMs) for elasticity problems
  //written by IK, Feb. 2012
  //Check if we want to give ML RBMs (from parameterlist)
  int numRBMs = params_->get<int>("Number RBMs for ML", 0); 
  bool setRBMs = false;   
  if (numRBMs > 0) {
    setRBMs = true; 
    int numScalar = 0;
    if (numRBMs == 2 || numRBMs == 3) 
      rigidBodyModes->setParameters(neq, numDim, numScalar, numRBMs, setRBMs);
    else
      TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"The specified number of RBMs " 
                                     << numRBMs << " is not valid!  Valid values are 0, 2 and 3.");
  }

  // Need to allocate a fields in mesh database
  if (params->isParameter("RequiredFields"))
  {
    // Need to allocate a fields in mesh database
    Teuchos::Array<std::string> req = params->get<Teuchos::Array<std::string> > ("Required Fields");
    for (int i(0); i<req.size(); ++i)
      this->requirements.push_back(req[i]);
  }
  else
  {
    this->requirements.push_back("surface_height");
#ifdef CISM_HAS_FELIX
    this->requirements.push_back("xgrad_surface_height"); //ds/dx which can be passed from CISM
    this->requirements.push_back("ygrad_surface_height"); //ds/dy which can be passed from CISM
#endif
    this->requirements.push_back("temperature");
    this->requirements.push_back("basal_friction");
    this->requirements.push_back("thickness");
    this->requirements.push_back("flow_factor");
    this->requirements.push_back("surface_velocity");
    this->requirements.push_back("surface_velocity_rms");
  }
}

 void interpolate(const std::vector<double>& positions,
   std::vector<double>& results) const 
   {
     Teuchos::Array<double> in(positions.size());
     for (int i=0; i<in.size(); i++) in[i] = positions[i];
     Teuchos::Array<double> out;
     interpolate(in, out);
     for (int i=0; i<out.size(); i++) results[i] = out[i];
   }

//Neumann BCs
void
Albany::PNPProblem::constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
{

   // Note: we only enter this function if sidesets are defined in the mesh file
   // i.e. meshSpecs.ssNames.size() > 0
   
    Albany::BCUtils<Albany::NeumannTraits> nbcUtils;

   // Check to make sure that Neumann BCs are given in the input file

   if(!nbcUtils.haveBCSpecified(this->params)) {
      return;
   }

   // Construct BC evaluators for all side sets and names
   // Note that the string index sets up the equation offset, so ordering is important
   //
   // Currently we aren't exactly doing this right.  I think to do this
   // correctly we need different neumann evaluators for each DOF (velocity,
   // pressure, temperature, flux) since velocity is a vector and the 
   // others are scalars.  The dof_names stuff is only used
   // for robin conditions, so at this point, as long as we don't enable
   // robin conditions, this should work.
   
   std::vector<std::string> nbcNames;
   Teuchos::RCP< Teuchos::Array<std::string> > dof_names =
     Teuchos::rcp(new Teuchos::Array<std::string>);
   // TODO: arbitraty numSpecies
   Teuchos::Array<Teuchos::Array<int> > offsets;
   int idx = 0;
     nbcNames.push_back("C1");
     offsets.push_back(Teuchos::Array<int>(1,idx++));
     if (numSpecies>=2) {
       nbcNames.push_back("C2");
       offsets.push_back(Teuchos::Array<int>(1,idx++));
     }
     if (numSpecies==3) {
       nbcNames.push_back("C3");
       offsets.push_back(Teuchos::Array<int>(1,idx++));
     }
     nbcNames.push_back("Phi");
     offsets.push_back(Teuchos::Array<int>(1,idx++));
     dof_names->push_back("Concentration");
     dof_names->push_back("Potential");

   // Construct BC evaluators for all possible names of conditions
   // Should only specify flux vector components (dudx, dudy, dudz), or dudn, not both
   std::vector<std::string> condNames; //dudx, dudy, dudz, dudn, basal 


   nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, nbcNames,
                                           Teuchos::arcp(dof_names),
                                           true, 0, condNames, offsets, dl,
                                           this->params, this->paramLib);
}