用于切片的C++表达式模板



我正在重写一个使用表达式模板的c++FEM库。在以前的版本中,我能够做到,例如

BilinearForm s(mesh, Vh, Vh, [&gp, &mesh](const TrialFunction& u, const TestFunction& v) -> double
{
return mesh.integrate(gp*dot(grad(u), grad(v)) );
});

它将自动评估Vh中每个试验和测试函数的表达式,并向我返回稀疏矩阵。

它的代码相当沉重和繁琐,所以我想从这篇文章中获得一些灵感来重写它。对于在单个元素上定义的函数来说,这是相当简单的。

template<typename E1, typename E2>
class ElementWiseScalarProd : public ElementWiseScalarExpression< ElementWiseScalarProd<E1, E2> >
{
public:
ElementWiseScalarProd(const E1& lhs, const E2& rhs) : m_lhs(lhs), m_rhs(rhs) {}
public:
inline double operator[] (const size_t k) const { return m_lhs[k] * m_rhs[k]; }
public:
inline bool containsTrial() const { return m_lhs.containsTrial() or m_rhs.containsTrial(); }
inline bool containsTest()  const { return m_lhs.containsTest() or m_rhs.containsTest(); }
public:
inline const Element* getElement() const { assert(m_lhs.getElement() == m_rhs.getElement()); return m_lhs.getElement(); }
private:
const E1& m_lhs;
const E2& m_rhs;
};

然而,当我想乘以在整个网格上定义的函数时,事情会变得有点棘手。我的运算符的返回类型变为数组或数组的切片。

template<typename E1, typename E2>
class FunctionProd : public FunctionExpression< FunctionProd<E1, E2> >
{
public:
typedef ElementWiseScalarProd<E1::ElementWiseType, E2::ElementWiseType> ElementWiseType;
public:
inline const ElementWiseType operator[] (const size_t e) const { return ElementWiseType(m_lhs[e], m_rhs[e]); }
public:
inline const Mesh* getMesh() const { assert(m_lhs.getMesh() == m_rhs.getMesh()); return m_lhs.getMesh(); }
private:
const E1& m_lhs;
const E2& m_rhs;
};

似乎,如果我要将在元素上定义的函数和在整个网格上定义的功能相乘,我的FunctionProd::operator[]应该返回一个引用,但这意味着我需要存储它创建的对象,不是吗?有办法绕过这一点吗?

提前感谢

编辑:这里回答了类似的问题在表达式模板中嵌套子表达式

您需要考虑内部表达式的值类别。

你可以用一个特征来区分左值和右值,以及一堆完美的转发。

template <typename T>
struct expression_holder {
using type = T;
};
template <typename T>
struct expression_holder<const T> : expression_holder<T> {
};
template <typename T>
struct expression_holder<T &> {
using type = const T &;
};
template <typename T>
struct expression_holder<T &&> {
using type = T;
};
template <typename T>
using expression_holder_t = typename expression_holder<T>::type;
template<typename E1, typename E2>
class ElementWiseScalarProd : public ElementWiseScalarExpression< ElementWiseScalarSum<E1, E2> >
{
public:
ElementWiseScalarProd(E1&& lhs, E2&& rhs) : m_lhs(std::forward<E1>(lhs)), m_rhs(std::forward<E2>(rhs)) {}
...
private:
expression_holder_t<E1> m_lhs;
expression_holder_t<E2> m_rhs;
};
template<typename E1, typename E2> 
ElementWiseScalarProd<E1, E2> scalar_prod(E1 && lhs, E2 && rhs) {
return ElementWiseScalarProd<E1, E2>(std::forward<E1>(lhs), std::forward<E2>(rhs));
}

最新更新