function_integral_info.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/assembly/base.hpp>
#include <kernel/analytic/function.hpp>
#include <kernel/util/dist.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
    template<typename DataType_, typename ValueType_, typename GradientType_, typename HessianType_>
    class FunctionIntegralInfo
    {
    public:
      typedef DataType_ DataType;
      typedef ValueType_ ValueType;
      typedef GradientType_ GradientType;
      typedef HessianType_ HessianType;
 
      int max_der;
 
      ValueType_ value;
 
      GradientType_ grad;
 
      HessianType_ hess;
 
      DataType_ norm_h0_sqr;
 
      DataType_ norm_h1_sqr;
 
      DataType_ norm_h2_sqr;
 
      DataType_ norm_l1;
 
      DataType_ norm_lmax;
 
      ValueType_ norm_h0_sqr_comp;
      ValueType_ norm_h1_sqr_comp;
      ValueType_ norm_h2_sqr_comp;
      ValueType_ norm_l1_comp;
      ValueType_ norm_lmax_comp;
 
      DataType divergence_l2_sqr;
 
      DataType vorticity_l2_sqr;
 
    private:
      static void _push_lmax(DataType_& a, const DataType_& b)
      {
        a = Math::max(a, b);
      }
 
      template<int n_, int sn_>
      static void _push_lmax(Tiny::Vector<DataType_, n_, sn_>& a, const Tiny::Vector<DataType_, n_, sn_>& b)
      {
        for(int i(0); i < n_; ++i)
          a[i] = Math::max(a[i], b[i]);
      }
 
      static void _write_to(DataType_* ptr, std::size_t& k, const DataType_& x)
      {
        ptr[k] = x;
        ++k;
      }
 
      template<int n_, int sn_>
      static void _write_to(DataType_* ptr, std::size_t& k, const Tiny::Vector<DataType_, n_, sn_>& x)
      {
        for(int i(0); i < n_; ++i, ++k)
          ptr[k] = x[i];
      }
 
      template<int m_, int n_, int sm_, int sn_>
      static void _write_to(DataType_* ptr, std::size_t& k, const Tiny::Matrix<DataType_, m_, n_, sm_, sn_>& x)
      {
        for(int i(0); i < m_; ++i)
          _write_to(ptr, k, x[i]);
      }
 
      template<int l_, int m_, int n_, int sl_, int sm_, int sn_>
      static void _write_to(DataType_* ptr, std::size_t& k, const Tiny::Tensor3<DataType_, l_, m_, n_, sl_, sm_, sn_>& x)
      {
        for(int i(0); i < l_; ++i)
          _write_to(ptr, k, x[i]);
      }
 
      static void _read_from(DataType_* ptr, std::size_t& k, DataType_& x)
      {
        x = ptr[k];
        ++k;
      }
 
      template<int n_, int sn_>
      static void _read_from(DataType_* ptr, std::size_t& k, Tiny::Vector<DataType_, n_, sn_>& x)
      {
        for(int i(0); i < n_; ++i, ++k)
          x[i] = ptr[k];
      }
 
      template<int m_, int n_, int sm_, int sn_>
      static void _read_from(DataType_* ptr, std::size_t& k, Tiny::Matrix<DataType_, m_, n_, sm_, sn_>& x)
      {
        for(int i(0); i < m_; ++i)
          _read_from(ptr, k, x[i]);
      }
 
      template<int l_, int m_, int n_, int sl_, int sm_, int sn_>
      static void _read_from(DataType_* ptr, std::size_t& k, Tiny::Tensor3<DataType_, l_, m_, n_, sl_, sm_, sn_>& x)
      {
        for(int i(0); i < l_; ++i)
          _read_from(ptr, k, x[i]);
      }
 
      static String _print_val(const DataType_& x, int prec)
      {
        return stringify_fp_sci(x, prec);
      }
 
      template<int n_, int sn_>
      static String _print_val(const Tiny::Vector<DataType_, n_, sn_>& x, int prec)
      {
        String s = "[";
        for(int i(0); i < n_; ++i)
          (s += " ") += stringify_fp_sci(x[i], prec);
        s += " ]";
        return s;
      }
 
      static String _print_val_sqrt(const DataType_& x, int prec)
      {
        return stringify_fp_sci(Math::sqrt(x), prec);
      }
 
      template<int n_, int sn_>
      static String _print_val_sqrt(const Tiny::Vector<DataType_, n_, sn_>& x, int prec)
      {
        String s = "[";
        for(int i(0); i < n_; ++i)
          (s += " ") += stringify_fp_sci(Math::sqrt(x[i]), prec);
        s += " ]";
        return s;
      }
 
      static String _print_norm(const DataType_& x, const DataType_&, int prec)
      {
        return _print_val(x, prec);
      }
 
      template<int n_, int sn_>
      static String _print_norm(const DataType_& x, const Tiny::Vector<DataType_, n_, sn_>& v, int prec)
      {
        return _print_val(x, prec) + " " + _print_val(v, prec);
      }
 
      static String _print_norm_sqrt(const DataType_& x, const DataType_&, int prec)
      {
        return _print_val(Math::sqrt(x), prec);
      }
 
      template<int n_, int sn_>
      static String _print_norm_sqrt(const DataType_& x, const Tiny::Vector<DataType_, n_, sn_>& v, int prec)
      {
        return _print_val(Math::sqrt(x), prec) + " " + _print_val_sqrt(v, prec);
      }
 
      template<int m_, int n_, int sm_, int sn_>
      static DataType _calc_vorticity(const Tiny::Matrix<DataType, m_, n_, sm_, sn_>&)
      {
        // R^m -> R^n with m != n ==> nothing to do here
        return DataType(0);
      }
 
      template<int sm_, int sn_>
      static DataType _calc_vorticity(const Tiny::Matrix<DataType, 2, 2, sm_, sn_>& g)
      {
        // 2D: (dx u2 - dy u1)^2
        // J_ij = d/d_j f_i ==> (J_01 - J_10)^2
        return Math::sqr(g[1][0] - g[0][1]);
      }
 
      template<int sm_, int sn_>
      static DataType _calc_vorticity(const Tiny::Matrix<DataType, 3, 3, sm_, sn_>& g)
      {
        // 3D: (dy u3 - dz u2)^2 + (dz u1 - dx u3)^2 + (dx u2 - dy u1)
        return
          Math::sqr(g[2][1] - g[1][2]) + // dy u3 - dz u2 = J_21 - J_12
          Math::sqr(g[0][2] - g[2][0]) + // dz u1 - dx u3 = J_02 - J_20
          Math::sqr(g[1][0] - g[0][1]);  // dx u2 - dy u1 = J_10 - J_01
      }
 
 
    public:
      FunctionIntegralInfo() :
        max_der(0),
        value(DataType(0)),
        grad(DataType(0)),
        hess(DataType(0)),
        norm_h0_sqr(DataType(0)),
        norm_h1_sqr(DataType(0)),
        norm_h2_sqr(DataType(0)),
        norm_l1(DataType(0)),
        norm_lmax(DataType(0)),
        norm_h0_sqr_comp(DataType(0)),
        norm_h1_sqr_comp(DataType(0)),
        norm_h2_sqr_comp(DataType(0)),
        norm_l1_comp(DataType(0)),
        norm_lmax_comp(DataType(0)),
        divergence_l2_sqr(DataType(0)),
        vorticity_l2_sqr(DataType(0))
      {
      }
 
      void format()
      {
        max_der = 0;
        value = ValueType(DataType(0));
        grad = GradientType(DataType(0));
        hess = HessianType(DataType(0));
        norm_h0_sqr = DataType(0);
        norm_h1_sqr = DataType(0);
        norm_h2_sqr = DataType(0);
        norm_l1 = DataType(0);
        norm_lmax = DataType(0);
        norm_h0_sqr_comp = ValueType(DataType(0));
        norm_h1_sqr_comp = ValueType(DataType(0));
        norm_h2_sqr_comp = ValueType(DataType(0));
        norm_l1_comp = DataType(0);
        norm_lmax_comp = DataType(0);
        divergence_l2_sqr = DataType(0);
        vorticity_l2_sqr = DataType(0);
      }
 
      void push(const FunctionIntegralInfo& other)
      {
        value += other.value;
        grad  += other.grad;
        hess  += other.hess;
        norm_h0_sqr += other.norm_h0_sqr;
        norm_h1_sqr += other.norm_h1_sqr;
        norm_h2_sqr += other.norm_h2_sqr;
        norm_l1 += other.norm_l1;
        _push_lmax(norm_lmax, other.norm_lmax);
        norm_h0_sqr_comp += other.norm_h0_sqr_comp;
        norm_h1_sqr_comp += other.norm_h1_sqr_comp;
        norm_h2_sqr_comp += other.norm_h2_sqr_comp;
        norm_l1_comp += other.norm_l1_comp;
        _push_lmax(norm_lmax_comp, other.norm_lmax_comp);
        divergence_l2_sqr += other.divergence_l2_sqr;
        vorticity_l2_sqr += other.vorticity_l2_sqr;
      }
 
      void add_value(DataType omega, const DataType& v)
      {
        value += omega * v;
        const DataType av = Math::abs(v);
        const DataType sv = Math::sqr(v);
        norm_h0_sqr += omega * sv;
        norm_l1 += omega * av;
        norm_lmax = Math::max(norm_lmax, av);
      }
 
      template<int n_, int sn_>
      void add_value(DataType omega, const Tiny::Vector<DataType, n_, sn_>& v)
      {
        value += omega * v;
        for(int i(0); i < n_; ++i)
        {
          const DataType av = Math::abs(v[i]);
          const DataType sv = Math::sqr(v[i]);
          norm_h0_sqr += omega * sv;
          norm_l1 += omega * av;
          norm_lmax = Math::max(norm_lmax, av);
          norm_h0_sqr_comp[i] += omega * sv;
          norm_l1_comp[i] += omega * av;
          norm_lmax_comp[i] = Math::max(norm_lmax_comp[i], av);
        }
      }
 
      template<int n_, int sn_>
      void add_grad(DataType omega, const Tiny::Vector<DataType, n_, sn_>& g)
      {
        grad += omega * g;
        norm_h1_sqr += omega * g.norm_euclid_sqr();
      }
 
      template<int m_, int n_, int sm_, int sn_>
      void add_grad(DataType omega, const Tiny::Matrix<DataType, m_, n_, sm_, sn_>& g)
      {
        grad += omega * g;
        for(int i(0); i < m_; ++i)
        {
          const DataType sv = g[i].norm_euclid_sqr();
          norm_h1_sqr += omega * sv;
          norm_h1_sqr_comp[i] += omega * sv;
        }
 
        // add divergence component
        divergence_l2_sqr += omega * Math::sqr(g.trace());
 
        // add vorticity component (if possible)
        vorticity_l2_sqr += omega * _calc_vorticity(g);
      }
 
      template<int m_, int n_, int sm_, int sn_>
      void add_hess(DataType omega, const Tiny::Matrix<DataType, m_, n_, sm_, sn_>& h)
      {
        hess += omega * h;
        norm_h2_sqr += omega * h.norm_hessian_sqr();
      }
 
      template<int l_, int m_, int n_, int sl_, int sm_, int sn_>
      void add_hess(DataType omega, const Tiny::Tensor3<DataType, l_, m_, n_, sl_, sm_, sn_>& h)
      {
        hess += omega * h;
        for(int i(0); i < l_; ++i)
        {
          const DataType sv = h[i].norm_hessian_sqr();
          norm_h2_sqr += omega * sv;
          norm_h2_sqr_comp[i] += omega * sv;
        }
      }
 
      void synchronize(const Dist::Comm& comm)
      {
        // nothing to do here?
        if(comm.size() <= 1)
          return;
 
        static constexpr std::size_t max_len = std::size_t(7) +
          (6u*sizeof(ValueType) + sizeof(GradientType) + sizeof(HessianType)) / sizeof(DataType);
        DataType v[max_len];
 
        // write all stuff to array
        std::size_t k = 0u;
        _write_to(v, k, value);
        _write_to(v, k, grad);
        _write_to(v, k, hess);
        _write_to(v, k, norm_h0_sqr);
        _write_to(v, k, norm_h1_sqr);
        _write_to(v, k, norm_h2_sqr);
        _write_to(v, k, norm_l1);
        _write_to(v, k, norm_h0_sqr_comp);
        _write_to(v, k, norm_h1_sqr_comp);
        _write_to(v, k, norm_h2_sqr_comp);
        _write_to(v, k, norm_l1_comp);
        _write_to(v, k, divergence_l2_sqr);
        _write_to(v, k, vorticity_l2_sqr);
        XASSERT(k <= max_len);
 
        // synchronize
        comm.allreduce(v, v, k, Dist::op_sum);
 
        // read stuff form array
        k = std::size_t(0u);
        _read_from(v, k, value);
        _read_from(v, k, grad);
        _read_from(v, k, hess);
        _read_from(v, k, norm_h0_sqr);
        _read_from(v, k, norm_h1_sqr);
        _read_from(v, k, norm_h2_sqr);
        _read_from(v, k, norm_l1);
        _read_from(v, k, norm_h0_sqr_comp);
        _read_from(v, k, norm_h1_sqr_comp);
        _read_from(v, k, norm_h2_sqr_comp);
        _read_from(v, k, norm_l1_comp);
        _read_from(v, k, divergence_l2_sqr);
        _read_from(v, k, vorticity_l2_sqr);
 
        // write lmax values
        k = std::size_t(0u);
        _write_to(v, k, norm_lmax);
        _write_to(v, k, norm_lmax_comp);
 
        // synchronize
        comm.allreduce(v, v, k, Dist::op_max);
 
        // read lmax values
        k = std::size_t(0u);
        _read_from(v, k, norm_lmax);
        _read_from(v, k, norm_lmax_comp);
      }
 
      String print_norms(int precision = 0, std::size_t pad_size = 15u, char pad_char = '.') const
      {
        String s;
        s += String("H0-Norm").pad_back(pad_size, pad_char) + ": "
          + _print_norm_sqrt(norm_h0_sqr, norm_h0_sqr_comp, precision) + "\n";
        if(max_der >= 1)
          s += String("H1-Norm").pad_back(pad_size, pad_char) + ": "
            + _print_norm_sqrt(norm_h1_sqr, norm_h1_sqr_comp, precision) + "\n";
        if(max_der >= 2)
          s += String("H2-Norm").pad_back(pad_size, pad_char) + ": "
            + _print_norm_sqrt(norm_h2_sqr, norm_h2_sqr_comp, precision) + "\n";
        s += String("L1-Norm").pad_back(pad_size, pad_char) + ": "
          + _print_norm(norm_l1, norm_l1_comp, precision) + "\n";
        s += String("Lmax-Norm").pad_back(pad_size, pad_char) + ": "
          + _print_norm(norm_lmax, norm_lmax_comp, precision);
        return s;
      }
 
      String print_field_info(int precision = 0, std::size_t pad_size = 15u, char pad_char = '.') const
      {
        String s;
        s += String("Kinetic Energy").pad_back(pad_size, pad_char) + ": "
          + _print_val(norm_h0_sqr * DataType(0.5), precision) + "\n";
        if(max_der >= 1)
        {
          s += String("Vorticity").pad_back(pad_size, pad_char) + ": "
            + _print_val_sqrt(vorticity_l2_sqr, precision) + "\n";
          s += String("Divergence").pad_back(pad_size, pad_char) + ": "
            + _print_val_sqrt(divergence_l2_sqr, precision) + "\n";
        }
        return s;
      }
    }; // class FunctionIntegralInfo<...>
 
    template<typename FuncIntInfo_, typename CellVector_>
    class FunctionCellIntegralInfo
    {
      public:
      FuncIntInfo_ integral_info;
      CellVector_ vec;
 
      FunctionCellIntegralInfo() = default;
 
      FunctionCellIntegralInfo(const FunctionCellIntegralInfo& other) :
       integral_info(),
       vec(other.vec.clone(LAFEM::CloneMode::Weak))
      {
        integral_info.push(other.integral_info);
      }
 
      FunctionCellIntegralInfo& operator=(const FunctionCellIntegralInfo& other)
      {
        if(&other == this)
          return *this;
        integral_info.format();
        integral_info.push(other.integral_info);
        vec.clone(other.vec, LAFEM::CloneMode::Weak);
      }
 
      explicit FunctionCellIntegralInfo(const FuncIntInfo_& info, const CellVector_& in_vec) :
        integral_info(),
        vec(in_vec.clone(LAFEM::CloneMode::Weak))
      {
        integral_info.push(info);
      }
 
      explicit FunctionCellIntegralInfo(const FuncIntInfo_& info, CellVector_&& in_vec) :
        integral_info(),
        vec(std::move(in_vec))
      {
        integral_info.push(info);
      }
 
    }; // class FunctionCellIntegralInfo<...>
 
    template<typename DataType_, typename Function_>
    struct AnalyticFunctionIntegral
    {
      typedef FunctionIntegralInfo<DataType_,
        typename Analytic::EvalTraits<DataType_, Function_>::ValueType,
        typename Analytic::EvalTraits<DataType_, Function_>::GradientType,
        typename Analytic::EvalTraits<DataType_, Function_>::HessianType> Type;
    };
 
    template<typename Vector_, typename Space_>
    struct DiscreteFunctionIntegral;
 
    template<typename DT_, typename IT_, typename Space_>
    struct DiscreteFunctionIntegral<LAFEM::DenseVector<DT_, IT_>, Space_>
    {
      typedef FunctionIntegralInfo<DT_, DT_,
        Tiny::Vector<DT_, Space_::shape_dim>,
        Tiny::Matrix<DT_, Space_::shape_dim, Space_::shape_dim>> Type;
    };
 
    template<typename DT_, typename IT_, int bs_, typename Space_>
    struct DiscreteFunctionIntegral<LAFEM::DenseVectorBlocked<DT_, IT_, bs_>, Space_>
    {
      typedef FunctionIntegralInfo<DT_, Tiny::Vector<DT_, bs_>,
        Tiny::Matrix<DT_, bs_, Space_::shape_dim>,
        Tiny::Tensor3<DT_, bs_, Space_::shape_dim, Space_::shape_dim>> Type;
    };
 
    namespace Intern
    {
      template<int max_der_>
      struct AnaFunIntJobHelper;
 
      template<>
      struct AnaFunIntJobHelper<0>
      {
        template<typename FID_, typename DT_, typename FEV_, typename IPT_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt)
        {
          fid.add_value(omega, fev.value(ipt));
        }
      };
 
      template<>
      struct AnaFunIntJobHelper<1>
      {
        template<typename FID_, typename DT_, typename FEV_, typename IPT_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt)
        {
          fid.add_value(omega, fev.value(ipt));
          fid.add_grad(omega, fev.gradient(ipt));
        }
      };
 
      template<>
      struct AnaFunIntJobHelper<2>
      {
        template<typename FID_, typename DT_, typename FEV_, typename IPT_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt)
        {
          fid.add_value(omega, fev.value(ipt));
          fid.add_grad(omega, fev.gradient(ipt));
          fid.add_hess(omega, fev.hessian(ipt));
        }
      };
 
      template<int max_der_>
      struct DiscFunIntJobHelper;
 
      template<>
      struct DiscFunIntJobHelper<0>
      {
        template<typename FID_, typename DT_, typename SED_, typename VT_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, const SED_& sed,
          const Tiny::Vector<VT_, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v(DT_(0));
          for(int i(0); i < n; ++i)
          {
            v += lvec[i] * sed.phi[i].value;
          }
          fid.add_value(omega, v);
        }
      };
 
      template<>
      struct DiscFunIntJobHelper<1>
      {
        // version for scalar functions
        template<typename FID_, typename DT_, typename SED_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, const SED_& sed,
          const Tiny::Vector<DT_, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v(DT_(0));
          typename FID_::GradientType g(DT_(0));
          for(int i(0); i < n; ++i)
          {
            v += lvec[i] * sed.phi[i].value;
            g += lvec[i] * sed.phi[i].grad;
          }
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
        }
        // version for vector fields
        template<typename FID_, typename DT_, typename SED_, int d_, int sd_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, const SED_& sed,
          const Tiny::Vector<Tiny::Vector<DT_, d_, sd_>, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v(DT_(0));
          typename FID_::GradientType g(DT_(0));
          for(int i(0); i < n; ++i)
          {
            v.axpy(sed.phi[i].value, lvec[i]);
            g.add_outer_product(lvec[i], sed.phi[i].grad);
          }
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
        }
      };
 
      template<>
      struct DiscFunIntJobHelper<2>
      {
        // version for scalar functions
        template<typename FID_, typename DT_, typename SED_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, const SED_& sed,
          const Tiny::Vector<DT_, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v(DT_(0));
          typename FID_::GradientType g(DT_(0));
          typename FID_::HessianType  h(DT_(0));
          for(int i(0); i < n; ++i)
          {
            v += lvec[i] * sed.phi[i].value;
            g += lvec[i] * sed.phi[i].grad;
            h += lvec[i] * sed.phi[i].hess;
          }
 
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
          fid.add_hess(omega, h);
        }
        // version for vector fields
        template<typename FID_, typename DT_, typename SED_, int d_, int sd_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, const SED_& sed,
          const Tiny::Vector<Tiny::Vector<DT_, d_, sd_>, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType v(DT_(0));
          typename FID_::GradientType g(DT_(0));
          typename FID_::HessianType h(DT_(0));
          for(int i(0); i < n; ++i)
          {
            v.axpy(sed.phi[i].value, lvec[i]);
            g.add_outer_product(lvec[i], sed.phi[i].grad);
            h.add_vec_mat_outer_product(lvec[i], sed.phi[i].hess);
          }
 
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
          fid.add_hess(omega, h);
        }
      };
 
      template<int max_der_>
      struct ErrFunIntJobHelper;
 
      template<>
      struct ErrFunIntJobHelper<0>
      {
        // version for scalar functions
        template<typename FID_, typename DT_, typename FEV_, typename IPT_, typename SED_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt, const SED_& sed,
          const Tiny::Vector<DT_, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v = fev.value(ipt);
          for(int i(0); i < n; ++i)
          {
            v -= lvec[i] * sed.phi[i].value;
          }
 
          fid.add_value(omega, v);
        }
 
        // version for vector fields
        template<typename FID_, typename DT_, typename FEV_, typename IPT_, typename SED_, int d_, int sd_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt, const SED_& sed,
          const Tiny::Vector<Tiny::Vector<DT_, d_, sd_>, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v = fev.value(ipt);
          for(int i(0); i < n; ++i)
          {
            v.axpy(-sed.phi[i].value, lvec[i]);
          }
 
          fid.add_value(omega, v);
        }
      };
 
      template<>
      struct ErrFunIntJobHelper<1>
      {
        // version for scalar functions
        template<typename FID_, typename DT_, typename FEV_, typename IPT_, typename SED_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt, const SED_& sed,
          const Tiny::Vector<DT_, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v = fev.value(ipt);
          typename FID_::GradientType g = fev.gradient(ipt);
          for(int i(0); i < n; ++i)
          {
            v -= lvec[i] * sed.phi[i].value;
            g -= lvec[i] * sed.phi[i].grad;
          }
 
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
        }
 
        // version for vector fields
        template<typename FID_, typename DT_, typename FEV_, typename IPT_, typename SED_, int d_, int sd_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt, const SED_& sed,
          const Tiny::Vector<Tiny::Vector<DT_, d_, sd_>, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v = fev.value(ipt);
          typename FID_::GradientType g = fev.gradient(ipt);
          for(int i(0); i < n; ++i)
          {
            v.axpy(-sed.phi[i].value, lvec[i]);
            g.add_outer_product(lvec[i], sed.phi[i].grad, -DT_(1));
          }
 
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
        }
      };
 
      template<>
      struct ErrFunIntJobHelper<2>
      {
        // version for scalar functions
        template<typename FID_, typename DT_, typename FEV_, typename IPT_, typename SED_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt, const SED_& sed,
          const Tiny::Vector<DT_, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v = fev.value(ipt);
          typename FID_::GradientType g = fev.gradient(ipt);
          typename FID_::HessianType  h = fev.hessian(ipt);
          for(int i(0); i < n; ++i)
          {
            v -= lvec[i] * sed.phi[i].value;
            g -= lvec[i] * sed.phi[i].grad;
            h -= lvec[i] * sed.phi[i].hess;
          }
 
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
          fid.add_hess(omega, h);
        }
 
        // version for vector fields
        template<typename FID_, typename DT_, typename FEV_, typename IPT_, typename SED_, int d_, int sd_, int n_, int sn_>
        static void work(FID_& fid, const DT_ omega, FEV_& fev, const IPT_& ipt, const SED_& sed,
          const Tiny::Vector<Tiny::Vector<DT_, d_, sd_>, n_, sn_>& lvec, int n)
        {
          typename FID_::ValueType    v = fev.value(ipt);
          typename FID_::GradientType g = fev.gradient(ipt);
          typename FID_::HessianType  h = fev.hessian(ipt);
          for(int i(0); i < n; ++i)
          {
            v.axpy(-sed.phi[i].value, lvec[i]);
            g.add_outer_product(lvec[i], sed.phi[i].grad, -DT_(1));
            h.add_vec_mat_outer_product(lvec[i], sed.phi[i].hess, -DT_(1));
          }
 
          fid.add_value(omega, v);
          fid.add_grad(omega, g);
          fid.add_hess(omega, h);
        }
      };
 
      template<typename Fun_, typename Vec_>
      struct ErrCompatHelper
      {
        // invalid combination of function and vector
        static constexpr bool valid = false;
      };
 
      template<typename Fun_, typename DT_, typename IT_>
      struct ErrCompatHelper<Fun_, LAFEM::DenseVector<DT_, IT_>>
      {
        // scalar vector, so function must also be scalar
        static constexpr bool valid = Fun_::ImageType::is_scalar;
      };
 
      template<typename Fun_, typename DT_, typename IT_, int dim_>
      struct ErrCompatHelper<Fun_, LAFEM::DenseVectorBlocked<DT_, IT_, dim_>>
      {
        typedef typename Fun_::ImageType ImageType;
        // blocked vector, so function must be a vector field of same dimension
        static constexpr bool valid = ImageType::is_vector && (ImageType::image_dim == dim_);
      };
    } // namespace Intern
  } // namespace Assembly
} // namespace FEAT