ElementData.hpp

Go to the documentation of this file.

// Copyright (C) 2010-2011 von Karman Institute for Fluid Dynamics, Belgium
//
// This software is distributed under the terms of the
// GNU Lesser General Public License version 3 (LGPLv3).
// See doc/lgpl.txt and doc/gpl.txt for the license text.

#ifndef cf3_solver_actions_Proto_ElementData_hpp
#define cf3_solver_actions_Proto_ElementData_hpp

#include <boost/array.hpp>

#include <boost/fusion/algorithm/iteration/for_each.hpp>
#include <boost/fusion/adapted/mpl.hpp>
#include <boost/fusion/mpl.hpp>
#include <boost/fusion/container/vector/convert.hpp>
#include <boost/fusion/container/vector.hpp>
#include <boost/fusion/view/filter_view.hpp>

#include <boost/mpl/assert.hpp>
#include <boost/mpl/for_each.hpp>
#include <boost/mpl/range_c.hpp>
#include <boost/mpl/transform.hpp>
#include <boost/mpl/vector_c.hpp>

#include "common/Component.hpp"
#include "common/FindComponents.hpp"

#include "math/VariablesDescriptor.hpp"
#include "math/LSS/BlockAccumulator.hpp"

#include "mesh/Elements.hpp"
#include "mesh/Field.hpp"
#include "mesh/Mesh.hpp"
#include "mesh/Region.hpp"
#include "mesh/Space.hpp"
#include "mesh/Dictionary.hpp"
#include "mesh/ElementData.hpp"
#include "mesh/Connectivity.hpp"

#include "ElementMatrix.hpp"
#include "ElementOperations.hpp"
#include "ElementTransforms.hpp"
#include "FieldSync.hpp"
#include "Terminals.hpp"

namespace cf3 {
namespace solver {
namespace actions {
namespace Proto {

namespace detail
{  
  template<typename EquationDataT>
  struct GetNbNodes
  {
    template<typename BoolT, typename EqDataT>
    struct apply
    {
    };
    
    template<typename EqDataT>
    struct apply<boost::mpl::false_, EqDataT>
    {
      static const Uint value = boost::remove_pointer
      <
        typename boost::remove_reference
        <
          typename boost::fusion::result_of::front<EquationDataT>::type
        >::type
      >::type::EtypeT::nb_nodes;
    };
    
    template<typename EqDataT>
    struct apply<boost::mpl::true_, EqDataT>
    {
      static const Uint value = 0;
    };
    
    static const Uint value = apply<typename boost::fusion::result_of::empty<EquationDataT>::type, EquationDataT>::value;
  };
  
  template<typename MatchingGrammarT>
  struct HasSubGrammar :
    boost::proto::or_
    <
      boost::proto::when
      <
        MatchingGrammarT,
        boost::mpl::true_()
      >,
      boost::proto::when
      <
        boost::proto::terminal< boost::proto::_ >,
        boost::mpl::false_()
      >,
      boost::proto::when
      <
        boost::proto::nary_expr<boost::proto::_, boost::proto::vararg<boost::proto::_> >,
        boost::proto::fold< boost::proto::_, boost::mpl::false_(), boost::mpl::max< boost::proto::_state, boost::proto::call< HasSubGrammar<MatchingGrammarT> > >() >
      >
    >
  {
  };
  
  template<Uint I>
  struct UsesVar : HasSubGrammar< boost::proto::terminal< Var<boost::mpl::int_<I>, boost::proto::_> > >
  {
  };
  
  template<typename ExprT, typename MatchingGrammarT>
  struct HasSubExpr
  {
    static const bool value = boost::result_of<HasSubGrammar<MatchingGrammarT>(ExprT)>::type::value;
  };
  
  template<Uint I>
  struct MatchImplicitEval :
    boost::proto::or_
    <
      boost::proto::when
      <
        boost::proto::or_
        <
          boost::proto::function< boost::proto::terminal< SFOp<boost::proto::_> >, boost::proto::vararg<boost::proto::_> >,
          boost::proto::terminal< SFOp<boost::proto::_> >,
          boost::proto::function<boost::proto::terminal<NodalValuesTag>, FieldTypes>,
          ElementMatrixGrammar,
          ElementMatrixGrammarIndexed< boost::mpl::int_<0>, boost::mpl::int_<0> >
        >,
        boost::mpl::false_()
      >,
      boost::proto::when
      <
        boost::proto::terminal< Var<boost::mpl::int_<I>, boost::proto::_> >,
        boost::mpl::true_()
      >
    >
  {
  };
  
  template<Uint I>
  struct HasEvalVar :
    boost::proto::or_
    <
      boost::proto::call< MatchImplicitEval<I> >,
      boost::proto::when
      <
        boost::proto::terminal< boost::proto::_ >,
        boost::mpl::false_()
      >,
      boost::proto::when
      <
        boost::proto::nary_expr<boost::proto::_, boost::proto::vararg<boost::proto::_> >,
        boost::proto::fold< boost::proto::_, boost::mpl::false_(), boost::mpl::max< boost::proto::_state, boost::proto::call< HasEvalVar<I> > >() >
      >
    >
  {
  };
  
  
}

template<typename ETYPE>
class GeometricSupport
{
public:
  typedef ETYPE EtypeT ;

  typedef typename EtypeT::NodesT ValueT;

  typedef const ValueT& ValueResultT;

  EIGEN_MAKE_ALIGNED_OPERATOR_NEW

  GeometricSupport(const mesh::Elements& elements) :
    m_coordinates(elements.geometry_fields().coordinates()),
    m_connectivity_array(elements.geometry_space().connectivity().array())
  {
  }

  void set_element(const Uint element_idx)
  {
    m_element_idx = element_idx;
    const mesh::Connectivity::ConstRow row = m_connectivity_array[element_idx];
    std::copy(row.begin(), row.end(), m_connectivity.begin());
    mesh::fill(m_nodes, m_coordinates, m_connectivity);
  }

  ValueResultT nodes() const
  {
    return m_nodes;
  }

  Real volume() const
  {
    return EtypeT::volume(m_nodes);
  }

  const typename EtypeT::CoordsT& coordinates(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    EtypeT::SF::compute_value(mapped_coords, m_sf);
    m_eval_result.noalias() = m_sf * m_nodes;
    return m_eval_result;
  }

  const typename EtypeT::CoordsT& coordinates() const
  {
    return m_eval_result;
  }

  const typename EtypeT::JacobianT& jacobian(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    EtypeT::compute_jacobian(mapped_coords, m_nodes, m_jacobian_matrix);
    return m_jacobian_matrix;
  }

  const typename EtypeT::JacobianT& jacobian() const
  {
    return m_jacobian_matrix;
  }

  const typename EtypeT::JacobianT& jacobian_inverse() const
  {
    return m_jacobian_inverse;
  }

  Real jacobian_determinant(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    m_jacobian_determinant = EtypeT::jacobian_determinant(mapped_coords, m_nodes);
    return m_jacobian_determinant;
  }

  Real jacobian_determinant() const
  {
    return m_jacobian_determinant;
  }

  const typename EtypeT::CoordsT& normal(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    EtypeT::normal(mapped_coords, m_nodes, m_normal_vector);
    return m_normal_vector;
  }

  const typename EtypeT::CoordsT& normal() const
  {
    return m_normal_vector;
  }

  void compute_shape_functions(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    EtypeT::SF::compute_value(mapped_coords, m_sf);
  }

  void compute_jacobian(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    compute_jacobian_dispatch(boost::mpl::bool_<EtypeT::dimension == EtypeT::dimensionality>(), mapped_coords);
  }

  void compute_coordinates() const
  {
    m_eval_result.noalias() = m_sf * m_nodes;
  }

  void compute_normal(const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    compute_normal_dispatch(boost::mpl::bool_<EtypeT::dimension - EtypeT::dimensionality == 1>(), mapped_coords);
  }

  const boost::array<Uint, EtypeT::nb_nodes>& element_connectivity() const
  {
    return m_connectivity;
  }

private:
  void compute_normal_dispatch(boost::mpl::false_, const typename EtypeT::MappedCoordsT&) const
  {
  }

  void compute_normal_dispatch(boost::mpl::true_, const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    EtypeT::normal(mapped_coords, m_nodes, m_normal_vector);
  }

  void compute_jacobian_dispatch(boost::mpl::false_, const typename EtypeT::MappedCoordsT&) const
  {
  }

  void compute_jacobian_dispatch(boost::mpl::true_, const typename EtypeT::MappedCoordsT& mapped_coords) const
  {
    EtypeT::compute_jacobian(mapped_coords, m_nodes, m_jacobian_matrix);
    bool is_invertible;
    m_jacobian_matrix.computeInverseAndDetWithCheck(m_jacobian_inverse, m_jacobian_determinant, is_invertible);
    cf3_assert(is_invertible);
  }

  ValueT m_nodes;

  const common::Table<Real>& m_coordinates;

  const mesh::Connectivity::ArrayT& m_connectivity_array;
  
  boost::array<Uint, EtypeT::nb_nodes> m_connectivity;

  Uint m_element_idx;

private:
  mutable typename EtypeT::SF::ValueT m_sf;
  mutable typename EtypeT::CoordsT m_eval_result;
  mutable typename EtypeT::JacobianT m_jacobian_matrix;
  mutable typename EtypeT::JacobianT m_jacobian_inverse;
  mutable Real m_jacobian_determinant;
  mutable typename EtypeT::CoordsT m_normal_vector;
};

inline mesh::Field& find_field(mesh::Elements& elements, const std::string& tag)
{
  mesh::Mesh& mesh = common::find_parent_component<mesh::Mesh>(elements);
  return common::find_component_recursively_with_tag<mesh::Field>(mesh, tag);
}

template<typename ETYPE, typename SupportEtypeT, Uint Dim, bool IsEquationVar>
class EtypeTVariableData
{
public:
  typedef ETYPE EtypeT;

private:
  template<Uint VarDim, int Dummy=0>
  struct InterpolationImpl
  {
    typedef Eigen::Matrix<Real, 1, VarDim> MatrixT;
    typedef const MatrixT& result_type;

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW

    template<typename NodeValuesT>
    result_type operator()(const typename EtypeT::SF::ValueT& sf, const NodeValuesT& values) const
    {
      stored_result.noalias() = sf * values;
      return stored_result;
    }

    mutable MatrixT stored_result;
  };

  template<int Dummy>
  struct InterpolationImpl<1, Dummy>
  {
    typedef Real result_type;

    template<typename NodeValuesT>
    result_type operator()(const typename EtypeT::SF::ValueT& sf, const NodeValuesT& values) const
    {
      stored_result = sf * values;
      return stored_result;
    }

    mutable Real stored_result;
  };

public:

  typedef GeometricSupport<SupportEtypeT> SupportT;

  typedef Eigen::Matrix<Real, EtypeT::nb_nodes, Dim> ValueT;

  typedef Eigen::Matrix<Real, EtypeT::nb_nodes*Dim, 1> ElementVectorT;

  typedef const ValueT& ValueResultT;

  typedef typename EtypeT::MappedCoordsT MappedCoordsT;

  typedef typename InterpolationImpl<Dim>::result_type EvalT;

  // Specialization of InterpolationImpl should ensure that 1x1 matrices are replaced by Reals
  BOOST_MPL_ASSERT_NOT(( boost::is_same<EvalT, const Eigen::Matrix<Real, 1, 1>&> ));

  typedef typename EtypeT::SF::GradientT GradientT;

  typedef Eigen::Matrix<Real, 1, Dim * EtypeT::nb_nodes> DivergenceLinT;

  static const Uint dimension = Dim;

  static const bool is_equation_variable = IsEquationVar;

  EIGEN_MAKE_ALIGNED_OPERATOR_NEW

  const mesh::Connectivity& get_connectivity(const std::string& tag, mesh::Elements& elements)
  {
    const mesh::Mesh& mesh = common::find_parent_component<mesh::Mesh>(elements);
    Handle<mesh::Dictionary const> dict = common::find_component_ptr_with_tag<mesh::Dictionary>(mesh, tag);
    if(is_null(dict))
      dict = mesh.geometry_fields().handle<mesh::Dictionary>(); // fall back to the geometry if the dict is not found by tag
    return elements.space(*dict).connectivity();
  }

  template<typename VariableT>
  EtypeTVariableData(const VariableT& placeholder, mesh::Elements& elements, const SupportT& support) :
    m_field(find_field(elements, placeholder.field_tag())),
    m_connectivity_array(get_connectivity(placeholder.field_tag(), elements).array()),
    m_support(support),
    offset(m_field.descriptor().offset(placeholder.name())),
    m_need_sync(false)
  {
  }
  
  ~EtypeTVariableData()
  {
    if(common::PE::Comm::instance().is_active())
    {
      const Uint my_sync = m_need_sync ? 1 : 0;
      Uint global_sync = 0;
      common::PE::Comm::instance().all_reduce(common::PE::plus(), &my_sync, 1, &global_sync);
      if(global_sync != 0)
        FieldSynchronizer::instance().insert(m_field, true);
    }
  }

  void set_element(const Uint element_idx)
  {
    m_element_idx = element_idx;
    mesh::fill(m_element_values, m_field, m_connectivity_array[element_idx], offset);
  }
  
  void update_block_connectivity(math::LSS::BlockAccumulator& block_accumulator)
  {
    block_accumulator.neighbour_indices(m_connectivity_array[m_element_idx]);
  }

  const SupportT& support() const
  {
    return m_support;
  }

  ValueResultT value() const
  {
    return m_element_values;
  }

  const ElementVectorT& element_vector() const
  {
    for(int i = 0; i != Dim; ++i)
    {
      m_element_vector.template segment<EtypeT::nb_nodes>(EtypeT::nb_nodes*i) = m_element_values.col(i);
    }
    return m_element_vector;
  }
  
  template<typename NodeValsT>
  void add_nodal_values(const NodeValsT& vals)
  {
    const mesh::Connectivity::ConstRow row = m_connectivity_array[m_element_idx];
    for(Uint i = 0; i != dimension; ++i)
    {
      m_element_values.col(i) += vals.template block<EtypeT::nb_nodes, 1>(i*EtypeT::nb_nodes, 0);
      for(Uint j = 0; j != EtypeT::nb_nodes; ++j)
        m_field[row[j]][offset+i] = m_element_values(j,i);
    }
    
    m_need_sync = true;
  }
  
  template<typename NodeValsT>
  void add_nodal_values_component(const NodeValsT& vals, const Uint component_idx)
  {
    const mesh::Connectivity::ConstRow row = m_connectivity_array[m_element_idx];
    for(Uint i = 0; i != EtypeT::nb_nodes; ++i)
    {
      m_element_values(i, component_idx) += vals[i];
      m_field[row[i]][offset+component_idx] = m_element_values(i, component_idx);
    }
    
    m_need_sync = true;
  }
  
  void compute_values(const MappedCoordsT& mapped_coords) const
  {
    compute_values_dispatch(boost::mpl::bool_<EtypeT::dimension == EtypeT::dimensionality>(), mapped_coords);
  }
  
  EvalT eval(const MappedCoordsT& mapped_coords) const
  {
    EtypeT::SF::compute_value(mapped_coords, m_sf);
    return m_eval(m_sf, m_element_values);
  }

  EvalT eval() const
  {
    return m_eval.stored_result;
  }

  const typename EtypeT::SF::ValueT& shape_function(const MappedCoordsT& mapped_coords) const
  {
    EtypeT::SF::compute_value(mapped_coords, m_sf);
    return m_sf;
  }

  const typename EtypeT::SF::ValueT& shape_function() const
  {
    return m_sf;
  }

  const GradientT& nabla(const MappedCoordsT& mapped_coords) const
  {
    EtypeT::SF::compute_gradient(mapped_coords, m_mapped_gradient_matrix);
    m_gradient.noalias() = m_support.jacobian(mapped_coords).inverse() * m_mapped_gradient_matrix;
    return m_gradient;
  }

  const GradientT& nabla() const
  {
    return m_gradient;
  }

private:
  void compute_values_dispatch(boost::mpl::false_, const MappedCoordsT& mapped_coords) const
  {
    EtypeT::SF::compute_value(mapped_coords, m_sf);
    m_eval(m_sf, m_element_values);
  }

  void compute_values_dispatch(boost::mpl::true_, const MappedCoordsT& mapped_coords) const
  {
    compute_values_dispatch(boost::mpl::false_(), mapped_coords);
    EtypeT::SF::compute_gradient(mapped_coords, m_mapped_gradient_matrix);
    m_gradient.noalias() = m_support.jacobian_inverse() * m_mapped_gradient_matrix;
  }

  ValueT m_element_values;

  mutable ElementVectorT m_element_vector;

  mesh::Field& m_field;

  const mesh::Connectivity::ArrayT& m_connectivity_array;

  const SupportT& m_support;

  Uint m_element_idx;

  mutable typename EtypeT::SF::ValueT m_sf;
  mutable typename EtypeT::SF::GradientT m_mapped_gradient_matrix;
  mutable GradientT m_gradient;

  InterpolationImpl<Dim> m_eval;
  
  bool m_need_sync;

public:
  const Uint offset;
};

template<typename SupportEtypeT, Uint Dim, bool IsEquationVar>
class EtypeTVariableData<ElementBased<Dim>, SupportEtypeT, Dim, IsEquationVar>
{
public:
  typedef ElementBased<Dim> EtypeT;

  typedef Eigen::Map< Eigen::Matrix<Real, 1, Dim> > ValueResultT;

  typedef GeometricSupport<SupportEtypeT> SupportT;

  static const Uint dimension = Dim;

  static const bool is_equation_variable = IsEquationVar;

  template<typename VariableT>
  EtypeTVariableData(const VariableT& placeholder, mesh::Elements& elements, const SupportT& support) :
    m_field(find_field(elements, placeholder.field_tag())),
    m_support(support),
    m_elements_begin(m_field.dict().space(elements).size() ? m_field.dict().space(elements).connectivity()[0][0] : 0),
    offset(m_field.descriptor().offset(placeholder.name()))
  {
  }

  void set_element(const Uint element_idx)
  {
    m_field_idx = element_idx + m_elements_begin;
  }

  ValueResultT value() const
  {
    return ValueResultT(&m_field[m_field_idx][offset]);
  }

  typedef typename SupportEtypeT::MappedCoordsT MappedCoordsT;
  typedef ValueResultT EvalT;

  EvalT eval(const MappedCoordsT& mapped_coords=MappedCoordsT()) const
  {
    return value();
  }

  // Dummy types for compatibility with higher order elements
  RealMatrix& nabla(RealMatrix mapped_coords = RealMatrix()) const
  {
    cf3_assert(false); // should not be used
    return m_dummy_result;
  }

  const RealMatrix& shape_function(RealMatrix mapped_coords = RealMatrix()) const
  {
    cf3_assert(false); // should not be used
    return m_dummy_result;
  }
  
  void compute_values(const MappedCoordsT& mapped_coords) const
  {
  }

private:
  mesh::Field& m_field;
  const SupportT& m_support;
  const Uint m_elements_begin;
  Uint m_field_idx;
  RealMatrix m_dummy_result; // only there for compilation purposes during the checking of the variable types. Never really used.

public:
  const Uint offset;
};

template<typename SupportEtypeT, bool IsEquationVar>
class EtypeTVariableData<ElementBased<1>, SupportEtypeT, 1, IsEquationVar>
{
public:
  typedef ElementBased<1> EtypeT;

  typedef Real& ValueResultT;

  typedef GeometricSupport<SupportEtypeT> SupportT;

  static const Uint dimension = 1;

  static const bool is_equation_variable = IsEquationVar;

  template<typename VariableT>
  EtypeTVariableData(const VariableT& placeholder, mesh::Elements& elements, const SupportT& support) :
    m_field(find_field(elements, placeholder.field_tag())),
    m_support(support),
    m_elements_begin(m_field.dict().space(elements).size() ? m_field.dict().space(elements).connectivity()[0][0] : 0),
    offset(m_field.descriptor().offset(placeholder.name()))
  {
  }

  void set_element(const Uint element_idx)
  {
    m_field_idx = element_idx + m_elements_begin;
  }

  Real& value() const
  {
    cf3_assert(m_field_idx < m_field.size());
    return m_field[m_field_idx][offset];
  }

  typedef Real EvalT;

  typedef typename SupportEtypeT::MappedCoordsT MappedCoordsT;

  EvalT eval(const MappedCoordsT& mapped_coords=MappedCoordsT()) const
  {
    return value();
  }
  const RealMatrix& nabla(RealMatrix mapped_coords = RealMatrix()) const
  {
    cf3_assert(false); // should not be used
    return m_dummy_result;
  }

  const RealMatrix& shape_function(RealMatrix mapped_coords = RealMatrix()) const
  {
    cf3_assert(false); // should not be used
    return m_dummy_result;
  }
  
  void compute_values(const MappedCoordsT& mapped_coords) const
  {
  }

private:
  mesh::Field& m_field;
  const SupportT& m_support;
  const Uint m_elements_begin;
  Uint m_field_idx;
  RealMatrix m_dummy_result; // only there for compilation purposes during the checking of the variable types. Never really used.

public:
  const Uint offset;
};

struct IsEquationData
{
  template<typename DataT, int Dummy = 0>
  struct apply
  {
    typedef boost::mpl::bool_<boost::remove_pointer<DataT>::type::is_equation_variable> type;
  };

  template<int Dummy>
  struct apply<boost::mpl::void_, Dummy>
  {
    typedef boost::mpl::false_ type;
  };
};

template<typename VariablesT, typename SupportEtypeT, typename ShapeFunctionsT, typename EquationVariablesT, typename MatrixSizesT, typename EMatrixSizeT>
struct MakeVarData
{
  template<typename I>
  struct apply
  {
    typedef typename boost::mpl::at<VariablesT, I>::type VarT;
    typedef typename boost::mpl::at<ShapeFunctionsT, I>::type EtypeT;
    typedef typename boost::mpl::at<EquationVariablesT, I>::type IsEquationVar;
    typedef typename boost::mpl::if_<IsEquationVar, EMatrixSizeT, typename boost::mpl::at<MatrixSizesT, I>::type>::type MatSize;

    template<typename AVarT, typename AnETypeT>
    struct GetEETypeT
    {
      typedef typename boost::mpl::if_c<AnETypeT::order == 0, ElementBased<FieldWidth<VarT, SupportEtypeT>::value>, AnETypeT>::type type;
    };

    template<typename AnETypeT>
    struct GetEETypeT<boost::mpl::void_,  AnETypeT>
    {
      typedef boost::mpl::void_ type;
    };

    typedef typename GetEETypeT<VarT, EtypeT>::type EEtypeT;

    typedef typename boost::mpl::if_
    <
      boost::mpl::is_void_<VarT>,
      boost::mpl::void_,
      EtypeTVariableData<EEtypeT, SupportEtypeT, FieldWidth<VarT, SupportEtypeT>::value, IsEquationVar::value>*
    >::type type;
  };
};

template<typename VariablesT, typename EtypeT, typename EquationVariablesT, typename MatrixSizesT, typename EMatrixSizeT>
struct MakeVarData<VariablesT, EtypeT, EtypeT, EquationVariablesT, MatrixSizesT, EMatrixSizeT>
{
  template<typename I>
  struct apply
  {
    typedef typename boost::mpl::at<VariablesT, I>::type VarT;
    typedef typename boost::mpl::at<EquationVariablesT, I>::type IsEquationVar;
    typedef typename boost::mpl::if_<IsEquationVar, EMatrixSizeT, typename boost::mpl::at<MatrixSizesT, I>::type>::type MatSize;

    typedef typename boost::mpl::if_
    <
      boost::mpl::is_void_<VarT>,
      boost::mpl::void_,
      EtypeTVariableData<EtypeT, EtypeT, FieldWidth<VarT, EtypeT>::value, IsEquationVar::value>*
    >::type type;
  };
};

template<typename IsEqVarT, typename VariableSFT>
struct FilterElementField
{
  typedef typename boost::mpl::if_c<VariableSFT::order == 0, boost::mpl::false_, IsEqVarT>::type type;
};


template<typename IsEqVarT>
struct FilterElementField<IsEqVarT, boost::mpl::void_>
{
  typedef boost::mpl::false_ type;
};

template<typename EquationVariablesT, typename VariablesEtypeTT, typename SupportEtypeT>
struct FilterEquationVars
{
  typedef typename boost::mpl::eval_if
  <
    boost::mpl::is_sequence<VariablesEtypeTT>,
    boost::mpl::transform
    <
      typename boost::mpl::copy<EquationVariablesT, boost::mpl::back_inserter< boost::mpl::vector0<> > >::type,
      VariablesEtypeTT,
      FilterElementField<boost::mpl::_1, boost::mpl::_2>
    >,
    boost::mpl::transform
    <
      typename boost::mpl::copy<EquationVariablesT, boost::mpl::back_inserter< boost::mpl::vector0<> > >::type,
      FilterElementField<boost::mpl::_1, VariablesEtypeTT>
    >
  >::type type;
};

template<typename VariablesT, typename VariablesEtypeTT, typename SupportEtypeT, typename EquationVariablesInT>
class ElementData
{
public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  typedef typename boost::fusion::result_of::size<VariablesT>::type NbVarsT;

  typedef SupportEtypeT SupportShapeFunction;

  static const Uint dimension = SupportShapeFunction::dimension;

  typedef typename MatrixSizePerVar<VariablesT, VariablesEtypeTT, SupportEtypeT>::type MatrixSizesT;

  typedef typename FilterEquationVars<EquationVariablesInT, VariablesEtypeTT, SupportEtypeT>::type EquationVariablesT;

  typedef typename ElementMatrixSize<MatrixSizesT, EquationVariablesT>::type EMatrixSizeT;

  typedef Eigen::Matrix<Real, EMatrixSizeT::value, EMatrixSizeT::value> ElementMatrixT;

  typedef Eigen::Matrix<Real, EMatrixSizeT::value, 1> ElementVectorT;

  typedef typename boost::fusion::result_of::as_vector
  <
    typename boost::mpl::transform
    <
      typename boost::mpl::copy<boost::mpl::range_c<int,0,NbVarsT::value>, boost::mpl::back_inserter< boost::mpl::vector_c<Uint> > >::type, //range from 0 to NbVarsT
      MakeVarData<VariablesT, SupportEtypeT, VariablesEtypeTT, EquationVariablesT, MatrixSizesT, EMatrixSizeT>
    >::type
  >::type VariablesDataT;

  typedef boost::fusion::filter_view< VariablesDataT, IsEquationData > EquationDataT;
  
  static const Uint nb_lss_nodes = detail::GetNbNodes<EquationDataT>::value;

  ElementData(VariablesT& variables, mesh::Elements& elements) :
    m_variables(variables),
    m_elements(elements),
    m_support(elements),
    m_equation_data(m_variables_data)
  {
    boost::mpl::for_each< boost::mpl::range_c<int, 0, NbVarsT::value> >(InitVariablesData(m_variables, m_elements, m_variables_data, m_support));
    for(Uint i = 0; i != CF3_PROTO_MAX_ELEMENT_MATRICES; ++i)
    {
      m_element_matrices[i].setZero();
      m_element_vectors[i].setZero();
    }
    
    init_block_accumulator(typename boost::fusion::result_of::empty<EquationDataT>::type());
  }
  
  // Init LSS block accumulator
  void init_block_accumulator(boost::mpl::false_)
  {
    typedef typename boost::mpl::transform
    <
      typename boost::mpl::copy<VariablesT, boost::mpl::back_inserter< boost::mpl::vector0<> > >::type,
      FieldWidth<boost::mpl::_1, SupportEtypeT>
    >::type NbEqsPerVarT;
    
    block_accumulator.resize(nb_lss_nodes, ElementMatrixSize<NbEqsPerVarT, EquationVariablesT>::type::value);
  }
  
  // No LSS, so nothing to be done
  void init_block_accumulator(boost::mpl::true_)
  {
  }

  ~ElementData()
  {
    boost::mpl::for_each< boost::mpl::range_c<int, 0, NbVarsT::value> >(DeleteVariablesData(m_variables_data));
  }

  void set_element(const Uint element_idx)
  {
    m_element_idx = element_idx;
    m_support.set_element(element_idx);
    boost::mpl::for_each< boost::mpl::range_c<int, 0, NbVarsT::value> >(SetElement(m_variables_data, element_idx));
    boost::fusion::for_each(m_equation_data, FillRhs(m_element_rhs));
    update_blocks(typename boost::fusion::result_of::empty<EquationDataT>::type());
  }

  void update_blocks(boost::mpl::false_)
  {
    boost::fusion::front(m_equation_data)->update_block_connectivity(block_accumulator);
    indices_converted = false;
  }

  void update_blocks(boost::mpl::true_)
  {
  }

  template<typename ExprT>
  void precompute_element_matrices(const typename SupportEtypeT::MappedCoordsT& mapped_coords, const ExprT& e)
  {
    // TODO: Add some granularity in here
    m_support.compute_shape_functions(mapped_coords);
    m_support.compute_coordinates();
    m_support.compute_jacobian(mapped_coords);
    m_support.compute_normal(mapped_coords);
    boost::mpl::for_each< boost::mpl::range_c<int, 0, NbVarsT::value> >(PrecomputeData<ExprT>(m_variables_data, mapped_coords));
  }

  template<typename I>
  struct DataType
  {
    typedef typename boost::remove_pointer
    <
      typename boost::remove_reference
      <
        typename boost::fusion::result_of::at
        <
          VariablesDataT, I
        >::type
      >::type
    >::type type;
  };

  template<typename I>
  struct VariableType
  {
    typedef typename boost::remove_pointer
    <
      typename boost::remove_reference
      <
        typename boost::fusion::result_of::at
        <
          VariablesT, I
        >::type
      >::type
    >::type type;
  };

  typedef GeometricSupport<SupportEtypeT> SupportT;

  template<typename I>
  typename DataType<I>::type& var_data(const I&)
  {
    return *boost::fusion::at<typename IndexType<I>::type>(m_variables_data);
  }

  template<typename I>
  const typename VariableType<I>::type& variable(const I&)
  {
    return boost::fusion::at<I>(m_variables);
  }

  SupportT& support()
  {
    return m_support;
  }

  const SupportT& support() const
  {
    return m_support;
  }

  ElementMatrixT& element_matrix(const int i)
  {
    return m_element_matrices[i];
  }

  ElementVectorT& element_vector(const int i)
  {
    return m_element_vectors[i];
  }

  ElementVectorT& element_rhs()
  {
    return m_element_rhs;
  };

  mutable math::LSS::BlockAccumulator block_accumulator;
  mutable bool indices_converted; // Indicate if the indices in the block accumulator have been converted to LSS indices

private:
  VariablesT& m_variables;

  mesh::Elements& m_elements;

  SupportT m_support;

  VariablesDataT m_variables_data;

  Uint m_element_idx;

  ElementMatrixT m_element_matrices[CF3_PROTO_MAX_ELEMENT_MATRICES];
  ElementVectorT m_element_vectors[CF3_PROTO_MAX_ELEMENT_MATRICES];
  ElementVectorT m_element_rhs;

  const EquationDataT m_equation_data;


  struct InitVariablesData
  {
    InitVariablesData(VariablesT& vars, mesh::Elements& elems, VariablesDataT& vars_data, const SupportT& sup) :
      m_variables(vars),
      m_elements(elems),
      m_variables_data(vars_data),
      m_support(sup)
    {
    }

    template<typename I>
    void operator()(const I&)
    {
      apply(boost::fusion::at<I>(m_variables), boost::fusion::at<I>(m_variables_data));
    }

    void apply(const boost::mpl::void_&, const boost::mpl::void_&)
    {
    }

    template<typename VarT, typename VarDataT>
    void apply(const VarT& v, VarDataT*& d)
    {
      d = new VarDataT(v, m_elements, m_support);
    }

    VariablesT& m_variables;
    mesh::Elements& m_elements;
    VariablesDataT& m_variables_data;
    const SupportT& m_support;
  };

  struct DeleteVariablesData
  {
    DeleteVariablesData(VariablesDataT& vars_data) : variables_data(vars_data)
    {
    }

    template<typename I>
    void operator()(const I&)
    {
      apply(boost::fusion::at<I>(variables_data));
    }

    void apply(const boost::mpl::void_&)
    {
    }

    template<typename T>
    void apply(T*& d)
    {
      delete d;
    }

    VariablesDataT& variables_data;
  };

  struct SetElement
  {
    SetElement(VariablesDataT& vars_data, const Uint elem_idx) :
      variables_data(vars_data),
      element_idx(elem_idx)
    {
    }

    template<typename I>
    void operator()(const I&)
    {
      apply(boost::fusion::at<I>(variables_data));
    }

    void apply(const boost::mpl::void_&)
    {
    }

    template<typename T>
    void apply(T*& d)
    {
      d->set_element(element_idx);
    }

    VariablesDataT& variables_data;
    const Uint element_idx;
  };

  template<typename ExprT>
  struct PrecomputeData
  {
    PrecomputeData(VariablesDataT& vars_data, const typename SupportEtypeT::MappedCoordsT& mapped_coords) :
      m_variables_data(vars_data),
      m_mapped_coords(mapped_coords)
    {
    }

    template<typename I>
    void operator()(const I&)
    {
      apply(typename boost::result_of<detail::UsesVar<I::value>(ExprT)>::type(), I(), boost::fusion::at<I>(m_variables_data));
    }

    template<typename I, typename T>
    void apply(boost::mpl::false_, I, T)
    {
    }

    template<typename I, typename T>
    void apply(boost::mpl::true_, I, T*& d)
    {
//       static const bool needs_eval = boost::result_of<detail::HasEvalVar<I::value>(ExprT)>::type::value;
      d->compute_values(m_mapped_coords);
    }

    template<Uint Dim, bool IsEquationVar>
    void apply(boost::mpl::true_, EtypeTVariableData<ElementBased<Dim>, SupportEtypeT, Dim, IsEquationVar>*&)
    {
    }

    // Variable is not used - do nothing
    template<typename T>
    void apply(boost::mpl::false_, T*& d)
    {
    }

  private:
    VariablesDataT& m_variables_data;
    const typename SupportEtypeT::MappedCoordsT& m_mapped_coords;
  };

  struct FillRhs
  {
    FillRhs(ElementVectorT& elm_v) :
      element_vector(elm_v)
    {
    }

    template<typename DataT>
    void operator()(const DataT* data) const
    {
      typename DataT::ValueResultT v = data->value();
      static const Uint rows = DataT::ValueT::RowsAtCompileTime;
      static const Uint cols = DataT::ValueT::ColsAtCompileTime;
      for(Uint i = 0; i != cols; ++i)
      {
        element_vector.template segment<rows>( (data->offset + i)*rows ) = v.col(i);
      }
    }

    ElementVectorT& element_vector;
  };
};

} // namespace Proto
} // namespace actions
} // namespace solver
} // namespace cf3

#endif // cf3_solver_actions_Proto_ElementData_hpp