Patch to NewField/Field.h (2/3)

[Dave Nystrom]

  // They are not really meant to be called by users.

  // Model-initializer constructor. Used by SubFieldView and 
  // View1Implementation above and by MakeFieldReturn in FieldCreateLeaf.h.

  template<class GT2, class T2, class ET2, class Initializer>
  Field(const Field<GT2, T2, ET2> &model, const Initializer &i)
  : fieldEngine_m(model.fieldEngine(), i)
    { }

  // Centering-Model constructor. This version differs from the one above
  // in that it constructs a field with uninitialized engines with a new
  // centering.  The new field gets its mesh and domain information from
  // the model field.

  template<class Centering, class GT2, class T2, class ET2>
  Field(const Centering &baseInit, const Field<GT2, T2, ET2> &model)
    : fieldEngine_m(baseInit, model.fieldEngine())
  { }


  //---------------------------------------------------------------------------
  // Empty destructor is fine for us.
  
  ~Field() { }


  //---------------------------------------------------------------------------
  // Accessors.

  inline const Engine_t &engine() const
    {
      return fieldEngine_m.engine();
    }

  inline const FieldEngine_t &fieldEngine() const
    {
      return fieldEngine_m;
    }
    
  inline FieldEngine_t &fieldEngine()
    {
      return fieldEngine_m;
    }
    
  inline int numSubFields() const
    {
      return fieldEngine_m.numSubFields();
    }
        
  inline const Domain_t physicalCellDomain() const
    {
      return fieldEngine_m.physicalCellDomain();
    }
        
  inline Domain_t totalCellDomain() const
    {
      return fieldEngine_m.totalCellDomain();
    }

  Domain_t physicalDomain(int iSubfield) const
    {
      return fieldEngine_m.physicalDomain(iSubfield);
    }

  Domain_t totalDomain(int iSubfield) const
    {
      return fieldEngine_m.totalDomain(iSubfield);
    }

  Domain_t physicalDomain() const
    {
      return fieldEngine_m.physicalDomain();
    }

  Domain_t totalDomain() const
    {
      return fieldEngine_m.totalDomain();
    }

  Domain_t domain() const
    {
      return fieldEngine_m.physicalDomain();
    }

  Layout_t layout() const
    {
      return fieldEngine_m.engine().layout();
    }

  //---------------------------------------------------------------------------
  // Instruct the field to make its own copy of its data.
  // Recursively call ourself with subfield views of this field. When we're
  // through, tell the fieldEngine to make a distinct copy of itself.

  void makeOwnCopy()
    {
      if (numSubFields() == 0)
        {
          // Make a distinct copy of the fieldEngine.
          
          fieldEngine_m.makeOwnCopy(*this);
        }
      else
        {
          for (int i = 0; i < numSubFields(); i++)
            (*this)[i].makeOwnCopy();
        }
    }
      
  
  //---------------------------------------------------------------------------
  // Sub-field view creation function.
  // A field consists of (potentially) several sub-fields. This function
  // returns a view of one of these.

  inline typename SubFieldView<This_t>::Type_t
  operator[](int iSubfield) const
    {
      typedef SubFieldView<This_t> Ret_t;
      return Ret_t::make(*this, iSubfield);
    }


  //---------------------------------------------------------------------------
  // View-creation operations. These operator() and read() functions take 
  // zero or more sub-domains, which combine to form a domain with 
  // dimensionality identical to the rank of the field. The zero argument 
  // version returns a view of the physical domain and the 'All'-suffixed
  // versions return a view of the total domain.

  inline typename View1<This_t, Domain_t>::ReadType_t 
  read() const
    {
      typedef View1<This_t, Domain_t> Ret_t;
      return Ret_t::makeRead(*this, physicalDomain());
    }

  inline typename View1<This_t, Domain_t>::ReadType_t 
  readAll() const
    {
      typedef View1<This_t, Domain_t> Ret_t;
      return Ret_t::makeRead(*this, totalDomain());
    }

  template<class Sub1> 
  inline typename View1<This_t, Sub1>::ReadType_t 
  read(const Sub1 &s1) const
    {
      typedef View1<This_t, Sub1> Ret_t;
      return Ret_t::makeRead(*this, s1);
    }

  template<class Sub1, class Sub2> 
  inline typename View2<This_t, Sub1, Sub2>::ReadType_t 
  read(const Sub1 &s1, const Sub2 &s2) const
    {
      typedef View2<This_t, Sub1, Sub2> Ret_t;
      return Ret_t::makeRead(*this, s1, s2);
    }

  template<class Sub1, class Sub2, class Sub3> 
  inline typename View3<This_t, Sub1, Sub2, Sub3>::ReadType_t 
  read(const Sub1 &s1, const Sub2 &s2, const Sub3 &s3) const
    {
      typedef View3<This_t, Sub1, Sub2, Sub3> Ret_t;
      return Ret_t::makeRead(*this, s1, s2, s3);
    }

  inline typename View1<This_t, Domain_t>::Type_t 
  operator()() const
    {
      typedef View1<This_t, Domain_t> Ret_t;
      return Ret_t::make(*this, physicalDomain());
    }

  inline typename View1<This_t, Domain_t>::Type_t 
  all() const
    {
      typedef View1<This_t, Domain_t> Ret_t;
      return Ret_t::make(*this, totalDomain());
    }

  template<class Sub1> 
  inline typename View1<This_t, Sub1>::Type_t 
  operator()(const Sub1 &s1) const
    {
      typedef View1<This_t, Sub1> Ret_t;
      return Ret_t::make(*this, s1);
    }

  template<class Sub1, class Sub2> 
  inline typename View2<This_t, Sub1, Sub2>::Type_t 
  operator()(const Sub1 &s1, const Sub2 &s2) const
    {
      typedef View2<This_t, Sub1, Sub2> Ret_t;
      return Ret_t::make(*this, s1, s2);
    }

  template<class Sub1, class Sub2, class Sub3> 
  inline typename View3<This_t, Sub1, Sub2, Sub3>::Type_t 
  operator()(const Sub1 &s1, const Sub2 &s2, const Sub3 &s3) const
    {
      typedef View3<This_t, Sub1, Sub2, Sub3> Ret_t;
      return Ret_t::make(*this, s1, s2, s3);
    }


  //---------------------------------------------------------------------------
  // Component-forwarding functions. These work quite similar to the ones from
  // Array except we produce a Field with the same GeometryTag.

  inline typename ComponentView<Loc<1>, This_t>::Type_t
  comp(const int &i1) const
  {
    return ComponentView<Loc<1>, This_t>::make(*this, Loc<1>(i1));
  }

  inline typename ComponentView<Loc<2>, This_t>::Type_t
  comp(int i1, int i2) const
  {
    return ComponentView<Loc<2>, This_t>::make(*this, Loc<2>(i1, i2));
  }

  template<class Components>
  inline typename ComponentView<Components, This_t>::Type_t
  comp(const Components &loc) const
  {
    return ComponentView<Components, This_t>::make(*this, loc);
  }


  //---------------------------------------------------------------------------
  // Patch accessor functions returns the i'th patch.

  inline typename Patch<This_t>::Type_t
  patchLocal(EnginePatch::PatchID_t i) const
    {
      return Patch<This_t>::make(*this, i);
    }

  inline int
  numPatchesLocal() const
  {
    return engineFunctor(engine(), EngineNumPatches());
  }

  //---------------------------------------------------------------------------
  // Copy assignment operators. We pack this assignment expression into a
  // PETE binary expression tree node and then use this to construct an
  // array with an expression engine. We then pass this on to an evaluator,
  // which handles the computation. The first three versions handle assigning
  // Arrays and ConstArrays to Arrays and the fourth one handles assigning
  // scalars.

  This_t &operator=(const This_t &rhs)
    {
      assign(*this, rhs, OpAssign());
      return *this;
    }

  const This_t &operator=(const This_t &rhs) const
    {
      return assign(*this, rhs, OpAssign());
    }

  template<class T1>
  const This_t &operator=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpAssign());
    }


  //---------------------------------------------------------------------------
  // Op-assignment operators. 

  // Addition.

  template<class T1>
  const This_t &operator+=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpAddAssign());
    }

  // Subtraction.

  template<class T1>
  const This_t &operator-=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpSubtractAssign());
    }

  // Multiplication.

  template<class T1>
  const This_t &operator*=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpMultiplyAssign());
    }

  // Division.

  template<class T1>
  const This_t &operator/=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpDivideAssign());
    }

  // Modulus.

  template<class T1>
  const This_t &operator%=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpModAssign());
    }

  // Bitwise-Or.

  template<class T1>
  const This_t &operator|=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpBitwiseOrAssign());
    }

  // Bitwise-And.

  template<class T1>
  const This_t &operator&=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpBitwiseAndAssign());
    }

  // Bitwise-Xor.

  template<class T1>
  const This_t &operator^=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpBitwiseXorAssign());
    }

  // Left shift.

  template<class T1>
  const This_t &operator<<=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpLeftShiftAssign());
    }

  // Right shift.

  template<class T1>
  const This_t &operator>>=(const T1 &rhs) const
    {
      return assign(*this, rhs, OpRightShiftAssign());
    }


  //---------------------------------------------------------------------------
  // Automatic update support.
  
  template<class Category>         
  void addUpdater(const Category &bc)
  {
    if (numSubFields() == 0)
      {
        fieldEngine_m.updaters().addUpdater(createUpdater(*this, bc));
      }
    else
      {
        for (int i = 0; i < numSubFields(); i++)
          (*this)[i].addUpdater(bc);
      }
  }

  template<class UpdaterIntitializer>  
  void addUpdaters(const UpdaterIntitializer &init)
    {
      init(*this);
    }
  
  void removeUpdaters()
    {
      if (numSubFields() == 0)
        fieldEngine_m.updaters().erase();
      else
        {
          for (int i = 0; i < numSubFields(); i++)
            (*this)[i].removeUpdaters();
        }
    }
  
  void update(bool makeDirty = false) const
    {
      if (numSubFields() == 0)
        {
          if (makeDirty)
            fieldEngine_m.updaters().setDirty();
          fieldEngine_m.updaters().notifyPreRead();
        }
      else
        {
          for (int i = 0; i < numSubFields(); i++)
            (*this)[i].update(makeDirty);
        }
    }
  
  void setDirty() const
    {
      if (numSubFields() == 0)
        {
          fieldEngine_m.updaters().setDirty();
        }
      else
        {
          for (int i = 0; i < numSubFields(); i++)
            (*this)[i].setDirty();
        }
    }
  
  void clearDirty() const
    {
      if (numSubFields() == 0)
        {
          fieldEngine_m.updaters().clearDirty();
        }
      else
        {
          for (int i = 0; i < numSubFields(); i++)
            (*this)[i].clearDirty();
        }
    }


private:

  FieldEngine_t fieldEngine_m;

};


//----------------------------------------------------------------------
// Set up a little traits class that distinguishes between OpAssign and
// other assignment operators that read the LHS.
//----------------------------------------------------------------------

template<class Op>
struct AssignOpReadWriteTraits
{
  enum { readLHS = true };
};

template<>
struct AssignOpReadWriteTraits<OpAssign>
{
  enum { readLHS = false };
};


//----------------------------------------------------------------------
// Apply the ConformTag to the leaves of the tree.
// Check to see if a given Field conforms.
//----------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, int Dim>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, ConformTag<Dim> >
{
  typedef bool Type_t;
  static Type_t apply1(const Interval<Dim> &d, 
    const ConformTag<Dim> &ct)
    {
      return conforms(d, ct);
    }
  template<int Dim2>
  static Type_t apply1(const Interval<Dim2> &d, 
    const ConformTag<Dim> &ct)
    {
      return false;
    }
  static Type_t apply(const Field<GeometryTag, T, EngineTag> &f,
    const ConformTag<Dim> &ct)
    {
      return apply1(f.physicalDomain(), ct);
    }
};


//----------------------------------------------------------------------
// This specialization of LeafFunctor is used to pass the 
// DataObjectRequest functor down into the FieldEngine. The default 
// behavior, given in the functor below, is to just pass it on to the 
// fieldEngine's engine.
//----------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, class RequestType>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>,
  DataObjectRequest<RequestType> >
{
  typedef Field<GeometryTag, T, EngineTag> Subject_t;
  typedef typename Subject_t::FieldEngine_t FieldEngine_t;
  typedef LeafFunctor<FieldEngine_t, DataObjectRequest<RequestType> > 
    LeafFunctor_t;
  typedef typename LeafFunctor_t::Type_t Type_t;
  enum { dataObject = LeafFunctor_t::dataObject };
  
  inline static
  Type_t apply(const Subject_t &f,
               const DataObjectRequest<RequestType> &functor)
    {
      return LeafFunctor_t::apply(field.fieldEngine(), functor);
    }
};

template<class GeometryTag, class T, class EngineTag, class RequestType>
struct LeafFunctor<FieldEngine<GeometryTag, T, EngineTag>,
  DataObjectRequest<RequestType> >
{
  typedef typename FieldEngine<GeometryTag, T, EngineTag>::Engine_t 
    Engine_t;
  enum { dataObject = Engine_t::dataObject };
  typedef typename DataObjectRequest<RequestType>::Type_t Type_t;
  inline static
  Type_t apply(const FieldEngine<GeometryTag, T, EngineTag> &f,
                   const DataObjectRequest<RequestType> &functor)
    {
      return DataObjectApply<dataObject>::apply(f.engine(), functor);
    }
};


//-----------------------------------------------------------------------------
// This specialization of LeafFunctor is used to get the domain type or the
// (zero-based) domain itself from a Field. Used only by Expression-Engine.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, DomainFunctorTag>
{
  typedef typename Field<GeometryTag, T, EngineTag>::Domain_t Type_t;

  inline static Type_t apply(const Field<GeometryTag, T, EngineTag> &f, 
    const DomainFunctorTag &)
    {
      // Return zero-based domain.
      
      return f.physicalDomain() - f.physicalDomain().firsts();
    }
};


//-----------------------------------------------------------------------------
// This specialization of LeafFunctor is used to pass the ExpressionApply
// functor
// down into the FieldEngine. The default behavior, given in the functor
// below, is to just pass it on to the fieldEngine's engine.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, class Tag>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, ExpressionApply<Tag> >
{
  typedef Field<GeometryTag, T, EngineTag> Subject_t;
  typedef typename Subject_t::FieldEngine_t FieldEngine_t;
  typedef LeafFunctor<FieldEngine_t, ExpressionApply<Tag> > LeafFunctor_t;
  typedef int Type_t;

  inline static
  Type_t apply(const Subject_t &field, 
               const ExpressionApply<Tag> &tag)
    {
      return LeafFunctor_t::apply(field.fieldEngine(), tag);
    }
};

template<class GeometryTag, class T, class EngineTag, class Tag>
struct LeafFunctor<FieldEngine<GeometryTag, T, EngineTag>,
  ExpressionApply<Tag> >
{
  typedef FieldEngine<GeometryTag, T, EngineTag> Subject_t;
  typedef typename Subject_t::Base_t Base_t;
  typedef LeafFunctor<Base_t, ExpressionApply<Tag> > LeafFunctor_t;
  typedef int Type_t;

  inline static
  Type_t apply(const Subject_t &fieldEngine, 
                 const ExpressionApply<Tag> &tag)
  {
    return LeafFunctor_t::apply(fieldEngine, tag);
  }
};

template<class GeometryTag, class T, class EngineTag, class Tag>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, EngineView<Tag> >
{
  typedef Field<GeometryTag, T, EngineTag> Subject_t;
  typedef typename Subject_t::Engine_t Engine_t;
  typedef typename LeafFunctor<Engine_t, EngineView<Tag> >::Type_t NewEngine_t;
  typedef typename NewEngine_t::Tag_t NewEngineTag_t;

  // Don't bother computing NewGeometry tag yet.
  // For now all EngineView operations are equivalent to Interval views.

  typedef Field<GeometryTag, T, NewEngineTag_t> Type_t;

  inline static
  Type_t apply(const Subject_t &field,
               const EngineView<Tag> &tag)
  {
    return Type_t(field, tag);
  }
};


//-----------------------------------------------------------------------------
// Handle general engine functor requests.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, class Tag>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, EngineFunctorTag<Tag> >
{
  typedef typename Field<GeometryTag,T,EngineTag>::Engine_t Engine_t;
  typedef typename EngineFunctor<Engine_t,Tag>::Type_t Type_t;
  inline static
  Type_t apply(const Field<GeometryTag, T, EngineTag> &field,
               const EngineFunctorTag<Tag> &tag)
  {
    return EngineFunctor<Engine_t,Tag>::apply(field.engine(), tag.tag());
  }
};


//---------------------------------------------------------------------------
// A specialization of EngineFunctor for field.
//---------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, class Tag>
struct EngineFunctor<Field<GeometryTag, T, EngineTag>, Tag>
{
  typedef typename Field<GeometryTag, T, EngineTag>::Engine_t Engine_t;
  typedef typename EngineFunctor<Engine_t, Tag>::Type_t Type_t;

  inline static 
  Type_t apply(const Field<GeometryTag, T, EngineTag> &field,
                   const Tag &tag)
    {
      return engineFunctor(field.engine(), tag);
    }
};


//-----------------------------------------------------------------------------
// This version of LeafFunctor is used by Expression-Engines to used to 
// evaluate a Field using indices. 
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, int Dim>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, EvalLeaf<Dim> >
{
  typedef typename Field<GeometryTag, T, EngineTag>::Element_t Type_t;
  inline static
  Type_t apply(const Field<GeometryTag, T, EngineTag> &f, 
    const EvalLeaf<Dim> &t) 
    {
      return t.eval(f.engine());
    }
};


//-----------------------------------------------------------------------------
// These leaf functor specializations are used to notify a field or expression
// that it is going to be read and, therefore, needs to update itself. 
//
// The first LeafFunctor represents default behavior, which is to do nothing. 
// The second handles fields other than those with expression-engines by simply
// calling update(). The third passes the tag to the leaves.
//
// Fields with engines that store internal fields AND don't copy those
// fields' updaters to its list must provide a specialization to get the 
// internal fields to update.
//
// NOTE: we don't use the ExpressionApply machinery here because this really
// operate on the engines.
//
//-----------------------------------------------------------------------------

struct PerformUpdateTag
{
  POOMA_PURIFY_CONSTRUCTORS(PerformUpdateTag)
};

template<class Node>
struct LeafFunctor<Node, PerformUpdateTag>
{
  typedef int Type_t;

  inline static
  Type_t apply(const Node &, const PerformUpdateTag &)
    {
      return 0;
    }
};


template<class GeometryTag, class T, class EngineTag>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, 
  PerformUpdateTag>
{
  typedef Field<GeometryTag, T, EngineTag> Subject_t;
  typedef int Type_t;

  inline static
  Type_t apply(const Subject_t &f, const PerformUpdateTag &)
    {
      f.update();
      return 0;
    }
};

template<class GeometryTag, class T, class Expr>
struct LeafFunctor<Field<GeometryTag, T, ExpressionTag<Expr> >, 
  PerformUpdateTag>
{
  typedef Field<GeometryTag, T, ExpressionTag<Expr> > Subject_t;
  typedef int Type_t;

  inline static
  Type_t apply(const Subject_t &f, const PerformUpdateTag &tag)
    {
      forEach(f.engine().expression(), tag, NullCombine());
      return 0;
    }
};


//-----------------------------------------------------------------------------
// This version of LeafFunctor is used to determine the type resulting from a
// sub-field view. 
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, SubFieldViewFunctorTag>
{
  typedef Field<GeometryTag, T, EngineTag> Type_t;
};

template<class T>
struct LeafFunctor<Scalar<T>, SubFieldViewFunctorTag>
{
  typedef Scalar<T> Type_t;
};


//-----------------------------------------------------------------------------
// This specialization of LeafFunctor is used to apply a view (subsetting) 
// operation to a Field. The domain will always be zero-based since this
// is used only by Expression-Engine. This is why we add the firsts() to
// the domain.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, class Domain>
struct LeafFunctor<Field<GeometryTag, T, EngineTag>, ViewFunctorTag<Domain> >
{
  typedef typename View1<Field<GeometryTag, T, EngineTag>, Domain>::Type_t 
    Type_t;
};


//-----------------------------------------------------------------------------
// Overload the << operator to print a Field to a stream.  This
// uses the 'PrintField' class to perform the formatting of the data.
// It will create a default printer, print out the field with it, and
// return the stream object.
//-----------------------------------------------------------------------------

template <class GeometryTag, class T, class EngineTag>
std::ostream &operator<<(std::ostream &o, 
  const Field<GeometryTag, T, EngineTag> &cf)
{
  Pooma::blockAndEvaluate();
  PrintField().print(o, cf);
  return o;
}

template <class GeometryTag, class T, class EngineTag>
std::fstream &operator<<(std::fstream &f, 
  const Field<GeometryTag, T, EngineTag> &cf)
{
  Pooma::blockAndEvaluate();
  PrintField().print(f, cf);
  return f;
}


//-----------------------------------------------------------------------------
// assign() function for Field assign-op array.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, 
 int Dim2, class T2, class EngineTag2, class Op>
const Field<GeometryTag, T, EngineTag> &
assign(const Field<GeometryTag, T, EngineTag> &lhs, 
       const Array<Dim2, T2, EngineTag2> &rhs, const Op &op)
{
  int nsf = lhs.numSubFields();
  
  if (nsf != 0)
    {
      // We need to take sub-field views in order to get at all of the
      // internal fields.
      
      for (int i = 0; i < nsf; i++)
        assign(lhs[i], rhs, op);
    }
  else
    {
      if (AssignOpReadWriteTraits<Op>::readLHS)
        lhs.update();
     
      // Evaluate. 

      Evaluator<MainEvaluatorTag>().evaluate(lhs, op, rhs);

      // Having done the evaluation, we need to notify the LHS
      // that we've just written.
  
      lhs.setDirty();
    }
        
  return lhs;
}


//-----------------------------------------------------------------------------
// assign() function for Field assign-op Field.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag,
  class GeometryTag2, class T2, class EngineTag2, class Op>
const Field<GeometryTag, T, EngineTag> &
assign(const Field<GeometryTag, T, EngineTag> &lhs,
       const Field<GeometryTag2, T2, EngineTag2> &rhs,
       const Op &op)
{
  int nsf = lhs.numSubFields();
  
  if (nsf != 0)
    {
      // We need to take sub-field views in order to get at all of the
      // internal fields.
      
      for (int i = 0; i < nsf; i++)
        assign(lhs[i], rhs[i], op);
    }
  else
    {
      // We need to notify the RHS and, when not doing pure assignment, the
      // LHS that we are getting ready to read.
      
      forEach(rhs, PerformUpdateTag(), NullCombine());
      
      if (AssignOpReadWriteTraits<Op>::readLHS)
        lhs.update();
  
      // Evaluate.

      Evaluator<MainEvaluatorTag>().evaluate(lhs, op, rhs);

      // Having done the evaluation, we need to notify the LHS
      // that we've just written.
  
      lhs.setDirty();
    }
        
  return lhs;
}


//-----------------------------------------------------------------------------
// assign() function for Field assign-op scalar.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag, class T1, class Op>
const Field<GeometryTag, T, EngineTag> &
assign(const Field<GeometryTag, T, EngineTag> &lhs, const T1 &rhs,
       const Op &op)
{
  int nsf = lhs.numSubFields();
  
  if (nsf != 0)
    {
      // We need to take sub-field views in order to get at all of the
      // internal fields.
      
      for (int i = 0; i < nsf; i++)
        assign(lhs[i], rhs, op);
    }
  else
    {
      if (AssignOpReadWriteTraits<Op>::readLHS)
        lhs.update();

      // Make an array out of the scalar.

      typedef Field<GeometryTag, T, EngineTag> LHS_t;
      Array<LHS_t::dimensions, T1, ConstantFunction>  
rhsExpr(lhs.physicalDomain());
      rhsExpr.engine().setConstant(rhs);
     
      // Evaluate. 

      Evaluator<MainEvaluatorTag>().evaluate(lhs, op, rhsExpr);

      // Having done the evaluation, we need to notify the LHS
      // that we've just written.
  
      lhs.setDirty();
    }
        
  return lhs;
}


//-----------------------------------------------------------------------------
// assign() function for Array assign-op Field.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag,
  int Dim2, class T2, class EngineTag2, class Op>
const Array<Dim2, T2, EngineTag2> &
assign(const Array<Dim2, T2, EngineTag2> &lhs,
       const Field<GeometryTag, T, EngineTag> &rhs, const Op &op)
{
  forEach(rhs, PerformUpdateTag(), NullCombine());

  Evaluator<MainEvaluatorTag>().evaluate(lhs, op, rhs);

  return lhs;
}


//-----------------------------------------------------------------------------
// Compute whether or not a Field is currently compressed.
//
// This is only a sensible thing to do if there are no subfields, hence the
// assertion.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class EngineTag>
inline bool compressed(const Field<GeometryTag, T, EngineTag> &f)
{
  PAssert(f.numSubFields() == 0);
  return compressed(f.engine());
}


//-----------------------------------------------------------------------------
// Compute the fraction of the total domain that is currently compressed.
//
// This is only a sensible thing to do if there are no subfields, hence the
// assertion.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class LTag>
inline double
compressedFraction(
  const Field<GeometryTag, T, MultiPatch<LTag,CompressibleBrick> > &f)
{
  PAssert(f.numSubFields() == 0);
  return compressedFraction(f.engine());
}


//-----------------------------------------------------------------------------
// (Try to) compress all the patches of the Field. Only need to do work with
// multipatch engines.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class LTag>
void
compress(Field<GeometryTag, T, MultiPatch<LTag,CompressibleBrick> > &)
{
  if (numSubFields() == 0)
    {
      compress(f.engine());
    }
  else
    {
      for (int i = 0; i < numSubFields(); i++)
        compress(f[i]);
    }
}


//-----------------------------------------------------------------------------
// Manually uncompress all the patches of the Field. Only need to do work with
// multipatch engines.
//-----------------------------------------------------------------------------

template<class GeometryTag, class T, class LTag>
void
uncompress(Field<GeometryTag, T, MultiPatch<LTag,CompressibleBrick> > &f)
{