Mesh.h

Go to the documentation of this file.

 
#ifndef GOOSEFEM_MESH_H
#define GOOSEFEM_MESH_H
 
#include "ElementQuad4.h"
#include "MatrixDiagonal.h"
#include "Vector.h"
#include "assertions.h"
#include "config.h"
 
namespace GooseFEM {
 
namespace Mesh {
 
template <class D>
inline std::vector<std::vector<size_t>> nodaltyings(const D& dofs);
 
enum class ElementType {
    Unknown, 
    Quad4, 
    Hex8, 
    Tri3 
};
 
template <class S, class T>
inline ElementType defaultElementType(const S& coor, const T& conn)
{
    GOOSEFEM_ASSERT(coor.dimension() == 2);
    GOOSEFEM_ASSERT(conn.dimension() == 2);
 
    if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) {
        return ElementType::Tri3;
    }
    if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) {
        return ElementType::Quad4;
    }
    if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) {
        return ElementType::Hex8;
    }
 
    throw std::runtime_error("Element-type not implemented");
}
 
namespace detail {
 
template <class T, class R>
inline T renum(const T& arg, const R& mapping)
{
    T ret = T::from_shape(arg.shape());
 
    auto jt = ret.begin();
 
    for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) {
        *jt = mapping(*it);
    }
 
    return ret;
}
 
} // namespace detail
 
inline array_type::tensor<size_t, 2> dofs(size_t nnode, size_t ndim)
{
    return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
}
 
class Renumber {
public:
    Renumber() = default;
 
    template <class T>
    Renumber(const T& dofs)
    {
        size_t n = xt::amax(dofs)() + 1;
        size_t i = 0;
 
        array_type::tensor<size_t, 1> unique = xt::unique(dofs);
 
        m_renum = xt::empty<size_t>({n});
 
        for (auto& j : unique) {
            m_renum(j) = i;
            ++i;
        }
    }
 
    template <class T>
    T apply(const T& list) const
    {
        return detail::renum(list, m_renum);
    }
 
    const array_type::tensor<size_t, 1>& index() const
    {
        return m_renum;
    }
 
private:
    array_type::tensor<size_t, 1> m_renum;
};
 
template <class T>
inline T renumber(const T& dofs)
{
    return Renumber(dofs).apply(dofs);
}
 
template <class D>
class RegularBase {
public:
    using derived_type = D;
 
    auto nelem() const
    {
        return derived_cast().m_nelem;
    }
 
    auto nnode() const
    {
        return derived_cast().m_nnode;
    }
 
    auto nne() const
    {
        return derived_cast().m_nne;
    }
 
    auto ndim() const
    {
        return derived_cast().m_ndim;
    }
 
    auto nelx() const
    {
        return derived_cast().nelx_impl();
    }
 
    auto nely() const
    {
        return derived_cast().nely_impl();
    }
 
    auto h() const
    {
        return derived_cast().m_h;
    }
 
    auto getElementType() const
    {
        return derived_cast().getElementType_impl();
    }
 
    auto coor() const
    {
        return derived_cast().coor_impl();
    }
 
    auto conn() const
    {
        return derived_cast().conn_impl();
    }
 
    auto dofs() const
    {
        return GooseFEM::Mesh::dofs(this->nnode(), this->ndim());
    }
 
    auto dofsPeriodic() const
    {
        array_type::tensor<size_t, 2> ret = this->dofs();
        array_type::tensor<size_t, 2> nodePer = this->nodesPeriodic();
        array_type::tensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
        array_type::tensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
 
        for (size_t j = 0; j < this->ndim(); ++j) {
            xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
        }
 
        return GooseFEM::Mesh::renumber(ret);
    }
 
    auto nodesPeriodic() const
    {
        return derived_cast().nodesPeriodic_impl();
    }
 
    auto nodesOrigin() const
    {
        return derived_cast().nodesOrigin_impl();
    }
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
};
 
template <class D>
class RegularBase2d : public RegularBase<D> {
public:
    using derived_type = D;
 
    auto nodesBottomEdge() const
    {
        return derived_cast().nodesBottomEdge_impl();
    }
 
    auto nodesTopEdge() const
    {
        return derived_cast().nodesTopEdge_impl();
    }
 
    auto nodesLeftEdge() const
    {
        return derived_cast().nodesLeftEdge_impl();
    }
 
    auto nodesRightEdge() const
    {
        return derived_cast().nodesRightEdge_impl();
    }
 
    auto nodesBottomOpenEdge() const
    {
        return derived_cast().nodesBottomOpenEdge_impl();
    }
 
    auto nodesTopOpenEdge() const
    {
        return derived_cast().nodesTopOpenEdge_impl();
    }
 
    auto nodesLeftOpenEdge() const
    {
        return derived_cast().nodesLeftOpenEdge_impl();
    }
 
    auto nodesRightOpenEdge() const
    {
        return derived_cast().nodesRightOpenEdge_impl();
    }
 
    auto nodesBottomLeftCorner() const
    {
        return derived_cast().nodesBottomLeftCorner_impl();
    }
 
    auto nodesBottomRightCorner() const
    {
        return derived_cast().nodesBottomRightCorner_impl();
    }
 
    auto nodesTopLeftCorner() const
    {
        return derived_cast().nodesTopLeftCorner_impl();
    }
 
    auto nodesTopRightCorner() const
    {
        return derived_cast().nodesTopRightCorner_impl();
    }
 
    auto nodesLeftBottomCorner() const
    {
        return derived_cast().nodesBottomLeftCorner_impl();
    }
 
    auto nodesLeftTopCorner() const
    {
        return derived_cast().nodesTopLeftCorner_impl();
    }
 
    auto nodesRightBottomCorner() const
    {
        return derived_cast().nodesBottomRightCorner_impl();
    }
 
    auto nodesRightTopCorner() const
    {
        return derived_cast().nodesTopRightCorner_impl();
    }
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
    friend class RegularBase<D>;
 
    array_type::tensor<size_t, 2> nodesPeriodic_impl() const
    {
        array_type::tensor<size_t, 1> bot = derived_cast().nodesBottomOpenEdge_impl();
        array_type::tensor<size_t, 1> top = derived_cast().nodesTopOpenEdge_impl();
        array_type::tensor<size_t, 1> lft = derived_cast().nodesLeftOpenEdge_impl();
        array_type::tensor<size_t, 1> rgt = derived_cast().nodesRightOpenEdge_impl();
        std::array<size_t, 2> shape = {bot.size() + lft.size() + size_t(3), size_t(2)};
        auto ret = array_type::tensor<size_t, 2>::from_shape(shape);
 
        ret(0, 0) = derived_cast().nodesBottomLeftCorner_impl();
        ret(0, 1) = derived_cast().nodesBottomRightCorner_impl();
 
        ret(1, 0) = derived_cast().nodesBottomLeftCorner_impl();
        ret(1, 1) = derived_cast().nodesTopRightCorner_impl();
 
        ret(2, 0) = derived_cast().nodesBottomLeftCorner_impl();
        ret(2, 1) = derived_cast().nodesTopLeftCorner_impl();
 
        size_t i = 3;
 
        xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
        xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
 
        i += bot.size();
 
        xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
        xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
 
        return ret;
    }
 
    auto nodesOrigin_impl() const
    {
        return derived_cast().nodesBottomLeftCorner_impl();
    }
};
 
template <class D>
class RegularBase3d : public RegularBase<D> {
public:
    using derived_type = D;
 
    auto nelz() const
    {
        return derived_cast().nelz_impl();
    }
 
    auto nodesBottom() const
    {
        return derived_cast().nodesBottom_impl();
    }
 
    auto nodesTop() const
    {
        return derived_cast().nodesTop_impl();
    }
 
    auto nodesLeft() const
    {
        return derived_cast().nodesLeft_impl();
    }
 
    auto nodesRight() const
    {
        return derived_cast().nodesRight_impl();
    }
 
    auto nodesFront() const
    {
        return derived_cast().nodesFront_impl();
    }
 
    auto nodesBack() const
    {
        return derived_cast().nodesBack_impl();
    }
 
    auto nodesFrontBottomEdge() const
    {
        return derived_cast().nodesFrontBottomEdge_impl();
    }
 
    auto nodesFrontTopEdge() const
    {
        return derived_cast().nodesFrontTopEdge_impl();
    }
 
    auto nodesFrontLeftEdge() const
    {
        return derived_cast().nodesFrontLeftEdge_impl();
    }
 
    auto nodesFrontRightEdge() const
    {
        return derived_cast().nodesFrontRightEdge_impl();
    }
 
    auto nodesBackBottomEdge() const
    {
        return derived_cast().nodesBackBottomEdge_impl();
    }
 
    auto nodesBackTopEdge() const
    {
        return derived_cast().nodesBackTopEdge_impl();
    }
 
    auto nodesBackLeftEdge() const
    {
        return derived_cast().nodesBackLeftEdge_impl();
    }
 
    auto nodesBackRightEdge() const
    {
        return derived_cast().nodesBackRightEdge_impl();
    }
 
    auto nodesBottomLeftEdge() const
    {
        return derived_cast().nodesBottomLeftEdge_impl();
    }
 
    auto nodesBottomRightEdge() const
    {
        return derived_cast().nodesBottomRightEdge_impl();
    }
 
    auto nodesTopLeftEdge() const
    {
        return derived_cast().nodesTopLeftEdge_impl();
    }
 
    auto nodesTopRightEdge() const
    {
        return derived_cast().nodesTopRightEdge_impl();
    }
 
    auto nodesBottomFrontEdge() const
    {
        return derived_cast().nodesFrontBottomEdge_impl();
    }
 
    auto nodesBottomBackEdge() const
    {
        return derived_cast().nodesBackBottomEdge_impl();
    }
 
    auto nodesTopFrontEdge() const
    {
        return derived_cast().nodesFrontTopEdge_impl();
    }
 
    auto nodesTopBackEdge() const
    {
        return derived_cast().nodesBackTopEdge_impl();
    }
 
    auto nodesLeftBottomEdge() const
    {
        return derived_cast().nodesBottomLeftEdge_impl();
    }
 
    auto nodesLeftFrontEdge() const
    {
        return derived_cast().nodesFrontLeftEdge_impl();
    }
 
    auto nodesLeftBackEdge() const
    {
        return derived_cast().nodesBackLeftEdge_impl();
    }
 
    auto nodesLeftTopEdge() const
    {
        return derived_cast().nodesTopLeftEdge_impl();
    }
 
    auto nodesRightBottomEdge() const
    {
        return derived_cast().nodesBottomRightEdge_impl();
    }
 
    auto nodesRightTopEdge() const
    {
        return derived_cast().nodesTopRightEdge_impl();
    }
 
    auto nodesRightFrontEdge() const
    {
        return derived_cast().nodesFrontRightEdge_impl();
    }
 
    auto nodesRightBackEdge() const
    {
        return derived_cast().nodesBackRightEdge_impl();
    }
 
    auto nodesFrontFace() const
    {
        return derived_cast().nodesFrontFace_impl();
    }
 
    auto nodesBackFace() const
    {
        return derived_cast().nodesBackFace_impl();
    }
 
    auto nodesLeftFace() const
    {
        return derived_cast().nodesLeftFace_impl();
    }
 
    auto nodesRightFace() const
    {
        return derived_cast().nodesRightFace_impl();
    }
 
    auto nodesBottomFace() const
    {
        return derived_cast().nodesBottomFace_impl();
    }
 
    auto nodesTopFace() const
    {
        return derived_cast().nodesTopFace_impl();
    }
 
    auto nodesFrontBottomOpenEdge() const
    {
        return derived_cast().nodesFrontBottomOpenEdge_impl();
    }
 
    auto nodesFrontTopOpenEdge() const
    {
        return derived_cast().nodesFrontTopOpenEdge_impl();
    }
 
    auto nodesFrontLeftOpenEdge() const
    {
        return derived_cast().nodesFrontLeftOpenEdge_impl();
    }
 
    auto nodesFrontRightOpenEdge() const
    {
        return derived_cast().nodesFrontRightOpenEdge_impl();
    }
 
    auto nodesBackBottomOpenEdge() const
    {
        return derived_cast().nodesBackBottomOpenEdge_impl();
    }
 
    auto nodesBackTopOpenEdge() const
    {
        return derived_cast().nodesBackTopOpenEdge_impl();
    }
 
    auto nodesBackLeftOpenEdge() const
    {
        return derived_cast().nodesBackLeftOpenEdge_impl();
    }
 
    auto nodesBackRightOpenEdge() const
    {
        return derived_cast().nodesBackRightOpenEdge_impl();
    }
 
    auto nodesBottomLeftOpenEdge() const
    {
        return derived_cast().nodesBottomLeftOpenEdge_impl();
    }
 
    auto nodesBottomRightOpenEdge() const
    {
        return derived_cast().nodesBottomRightOpenEdge_impl();
    }
 
    auto nodesTopLeftOpenEdge() const
    {
        return derived_cast().nodesTopLeftOpenEdge_impl();
    }
 
    auto nodesTopRightOpenEdge() const
    {
        return derived_cast().nodesTopRightOpenEdge_impl();
    }
 
    auto nodesBottomFrontOpenEdge() const
    {
        return derived_cast().nodesFrontBottomOpenEdge_impl();
    }
 
    auto nodesBottomBackOpenEdge() const
    {
        return derived_cast().nodesBackBottomOpenEdge_impl();
    }
 
    auto nodesTopFrontOpenEdge() const
    {
        return derived_cast().nodesFrontTopOpenEdge_impl();
    }
 
    auto nodesTopBackOpenEdge() const
    {
        return derived_cast().nodesBackTopOpenEdge_impl();
    }
 
    auto nodesLeftBottomOpenEdge() const
    {
        return derived_cast().nodesBottomLeftOpenEdge_impl();
    }
 
    auto nodesLeftFrontOpenEdge() const
    {
        return derived_cast().nodesFrontLeftOpenEdge_impl();
    }
 
    auto nodesLeftBackOpenEdge() const
    {
        return derived_cast().nodesBackLeftOpenEdge_impl();
    }
 
    auto nodesLeftTopOpenEdge() const
    {
        return derived_cast().nodesTopLeftOpenEdge_impl();
    }
 
    auto nodesRightBottomOpenEdge() const
    {
        return derived_cast().nodesBottomRightOpenEdge_impl();
    }
 
    auto nodesRightTopOpenEdge() const
    {
        return derived_cast().nodesTopRightOpenEdge_impl();
    }
 
    auto nodesRightFrontOpenEdge() const
    {
        return derived_cast().nodesFrontRightOpenEdge_impl();
    }
 
    auto nodesRightBackOpenEdge() const
    {
        return derived_cast().nodesBackRightOpenEdge_impl();
    }
 
    auto nodesFrontBottomLeftCorner() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
 
    auto nodesFrontBottomRightCorner() const
    {
        return derived_cast().nodesFrontBottomRightCorner_impl();
    }
 
    auto nodesFrontTopLeftCorner() const
    {
        return derived_cast().nodesFrontTopLeftCorner_impl();
    }
 
    auto nodesFrontTopRightCorner() const
    {
        return derived_cast().nodesFrontTopRightCorner_impl();
    }
 
    auto nodesBackBottomLeftCorner() const
    {
        return derived_cast().nodesBackBottomLeftCorner_impl();
    }
 
    auto nodesBackBottomRightCorner() const
    {
        return derived_cast().nodesBackBottomRightCorner_impl();
    }
 
    auto nodesBackTopLeftCorner() const
    {
        return derived_cast().nodesBackTopLeftCorner_impl();
    }
 
    auto nodesBackTopRightCorner() const
    {
        return derived_cast().nodesBackTopRightCorner_impl();
    }
 
    auto nodesFrontLeftBottomCorner() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
 
    auto nodesBottomFrontLeftCorner() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
 
    auto nodesBottomLeftFrontCorner() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
 
    auto nodesLeftFrontBottomCorner() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
 
    auto nodesLeftBottomFrontCorner() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
 
    auto nodesFrontRightBottomCorner() const
    {
        return derived_cast().nodesFrontBottomRightCorner_impl();
    }
 
    auto nodesBottomFrontRightCorner() const
    {
        return derived_cast().nodesFrontBottomRightCorner_impl();
    }
 
    auto nodesBottomRightFrontCorner() const
    {
        return derived_cast().nodesFrontBottomRightCorner_impl();
    }
 
    auto nodesRightFrontBottomCorner() const
    {
        return derived_cast().nodesFrontBottomRightCorner_impl();
    }
 
    auto nodesRightBottomFrontCorner() const
    {
        return derived_cast().nodesFrontBottomRightCorner_impl();
    }
 
    auto nodesFrontLeftTopCorner() const
    {
        return derived_cast().nodesFrontTopLeftCorner_impl();
    }
 
    auto nodesTopFrontLeftCorner() const
    {
        return derived_cast().nodesFrontTopLeftCorner_impl();
    }
 
    auto nodesTopLeftFrontCorner() const
    {
        return derived_cast().nodesFrontTopLeftCorner_impl();
    }
 
    auto nodesLeftFrontTopCorner() const
    {
        return derived_cast().nodesFrontTopLeftCorner_impl();
    }
 
    auto nodesLeftTopFrontCorner() const
    {
        return derived_cast().nodesFrontTopLeftCorner_impl();
    }
 
    auto nodesFrontRightTopCorner() const
    {
        return derived_cast().nodesFrontTopRightCorner_impl();
    }
 
    auto nodesTopFrontRightCorner() const
    {
        return derived_cast().nodesFrontTopRightCorner_impl();
    }
 
    auto nodesTopRightFrontCorner() const
    {
        return derived_cast().nodesFrontTopRightCorner_impl();
    }
 
    auto nodesRightFrontTopCorner() const
    {
        return derived_cast().nodesFrontTopRightCorner_impl();
    }
 
    auto nodesRightTopFrontCorner() const
    {
        return derived_cast().nodesFrontTopRightCorner_impl();
    }
 
    auto nodesBackLeftBottomCorner() const
    {
        return derived_cast().nodesBackBottomLeftCorner_impl();
    }
 
    auto nodesBottomBackLeftCorner() const
    {
        return derived_cast().nodesBackBottomLeftCorner_impl();
    }
 
    auto nodesBottomLeftBackCorner() const
    {
        return derived_cast().nodesBackBottomLeftCorner_impl();
    }
 
    auto nodesLeftBackBottomCorner() const
    {
        return derived_cast().nodesBackBottomLeftCorner_impl();
    }
 
    auto nodesLeftBottomBackCorner() const
    {
        return derived_cast().nodesBackBottomLeftCorner_impl();
    }
 
    auto nodesBackRightBottomCorner() const
    {
        return derived_cast().nodesBackBottomRightCorner_impl();
    }
 
    auto nodesBottomBackRightCorner() const
    {
        return derived_cast().nodesBackBottomRightCorner_impl();
    }
 
    auto nodesBottomRightBackCorner() const
    {
        return derived_cast().nodesBackBottomRightCorner_impl();
    }
 
    auto nodesRightBackBottomCorner() const
    {
        return derived_cast().nodesBackBottomRightCorner_impl();
    }
 
    auto nodesRightBottomBackCorner() const
    {
        return derived_cast().nodesBackBottomRightCorner_impl();
    }
 
    auto nodesBackLeftTopCorner() const
    {
        return derived_cast().nodesBackTopLeftCorner_impl();
    }
 
    auto nodesTopBackLeftCorner() const
    {
        return derived_cast().nodesBackTopLeftCorner_impl();
    }
 
    auto nodesTopLeftBackCorner() const
    {
        return derived_cast().nodesBackTopLeftCorner_impl();
    }
 
    auto nodesLeftBackTopCorner() const
    {
        return derived_cast().nodesBackTopLeftCorner_impl();
    }
 
    auto nodesLeftTopBackCorner() const
    {
        return derived_cast().nodesBackTopLeftCorner_impl();
    }
 
    auto nodesBackRightTopCorner() const
    {
        return derived_cast().nodesBackTopRightCorner_impl();
    }
 
    auto nodesTopBackRightCorner() const
    {
        return derived_cast().nodesBackTopRightCorner_impl();
    }
 
    auto nodesTopRightBackCorner() const
    {
        return derived_cast().nodesBackTopRightCorner_impl();
    }
 
    auto nodesRightBackTopCorner() const
    {
        return derived_cast().nodesBackTopRightCorner_impl();
    }
 
    auto nodesRightTopBackCorner() const
    {
        return derived_cast().nodesBackTopRightCorner_impl();
    }
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
    friend class RegularBase<D>;
 
    array_type::tensor<size_t, 2> nodesPeriodic_impl() const
    {
        array_type::tensor<size_t, 1> fro = derived_cast().nodesFrontFace_impl();
        array_type::tensor<size_t, 1> bck = derived_cast().nodesBackFace_impl();
        array_type::tensor<size_t, 1> lft = derived_cast().nodesLeftFace_impl();
        array_type::tensor<size_t, 1> rgt = derived_cast().nodesRightFace_impl();
        array_type::tensor<size_t, 1> bot = derived_cast().nodesBottomFace_impl();
        array_type::tensor<size_t, 1> top = derived_cast().nodesTopFace_impl();
 
        array_type::tensor<size_t, 1> froBot = derived_cast().nodesFrontBottomOpenEdge_impl();
        array_type::tensor<size_t, 1> froTop = derived_cast().nodesFrontTopOpenEdge_impl();
        array_type::tensor<size_t, 1> froLft = derived_cast().nodesFrontLeftOpenEdge_impl();
        array_type::tensor<size_t, 1> froRgt = derived_cast().nodesFrontRightOpenEdge_impl();
        array_type::tensor<size_t, 1> bckBot = derived_cast().nodesBackBottomOpenEdge_impl();
        array_type::tensor<size_t, 1> bckTop = derived_cast().nodesBackTopOpenEdge_impl();
        array_type::tensor<size_t, 1> bckLft = derived_cast().nodesBackLeftOpenEdge_impl();
        array_type::tensor<size_t, 1> bckRgt = derived_cast().nodesBackRightOpenEdge_impl();
        array_type::tensor<size_t, 1> botLft = derived_cast().nodesBottomLeftOpenEdge_impl();
        array_type::tensor<size_t, 1> botRgt = derived_cast().nodesBottomRightOpenEdge_impl();
        array_type::tensor<size_t, 1> topLft = derived_cast().nodesTopLeftOpenEdge_impl();
        array_type::tensor<size_t, 1> topRgt = derived_cast().nodesTopRightOpenEdge_impl();
 
        size_t tface = fro.size() + lft.size() + bot.size();
        size_t tedge = 3 * froBot.size() + 3 * froLft.size() + 3 * botLft.size();
        size_t tnode = 7;
        array_type::tensor<size_t, 2> ret =
            xt::empty<size_t>({tface + tedge + tnode, std::size_t(2)});
 
        size_t i = 0;
 
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesFrontBottomRightCorner_impl();
        ++i;
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesBackBottomRightCorner_impl();
        ++i;
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesBackBottomLeftCorner_impl();
        ++i;
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesFrontTopLeftCorner_impl();
        ++i;
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesFrontTopRightCorner_impl();
        ++i;
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesBackTopRightCorner_impl();
        ++i;
        ret(i, 0) = derived_cast().nodesFrontBottomLeftCorner_impl();
        ret(i, 1) = derived_cast().nodesBackTopLeftCorner_impl();
        ++i;
 
        for (size_t j = 0; j < froBot.size(); ++j) {
            ret(i, 0) = froBot(j);
            ret(i, 1) = bckBot(j);
            ++i;
        }
        for (size_t j = 0; j < froBot.size(); ++j) {
            ret(i, 0) = froBot(j);
            ret(i, 1) = bckTop(j);
            ++i;
        }
        for (size_t j = 0; j < froBot.size(); ++j) {
            ret(i, 0) = froBot(j);
            ret(i, 1) = froTop(j);
            ++i;
        }
        for (size_t j = 0; j < botLft.size(); ++j) {
            ret(i, 0) = botLft(j);
            ret(i, 1) = botRgt(j);
            ++i;
        }
        for (size_t j = 0; j < botLft.size(); ++j) {
            ret(i, 0) = botLft(j);
            ret(i, 1) = topRgt(j);
            ++i;
        }
        for (size_t j = 0; j < botLft.size(); ++j) {
            ret(i, 0) = botLft(j);
            ret(i, 1) = topLft(j);
            ++i;
        }
        for (size_t j = 0; j < froLft.size(); ++j) {
            ret(i, 0) = froLft(j);
            ret(i, 1) = froRgt(j);
            ++i;
        }
        for (size_t j = 0; j < froLft.size(); ++j) {
            ret(i, 0) = froLft(j);
            ret(i, 1) = bckRgt(j);
            ++i;
        }
        for (size_t j = 0; j < froLft.size(); ++j) {
            ret(i, 0) = froLft(j);
            ret(i, 1) = bckLft(j);
            ++i;
        }
 
        for (size_t j = 0; j < fro.size(); ++j) {
            ret(i, 0) = fro(j);
            ret(i, 1) = bck(j);
            ++i;
        }
        for (size_t j = 0; j < lft.size(); ++j) {
            ret(i, 0) = lft(j);
            ret(i, 1) = rgt(j);
            ++i;
        }
        for (size_t j = 0; j < bot.size(); ++j) {
            ret(i, 0) = bot(j);
            ret(i, 1) = top(j);
            ++i;
        }
 
        return ret;
    }
 
    auto nodesOrigin_impl() const
    {
        return derived_cast().nodesFrontBottomLeftCorner_impl();
    }
};
 
template <class S, class T>
inline array_type::tensor<size_t, 2>
overlapping(const S& coor_a, const T& coor_b, double rtol = 1e-5, double atol = 1e-8)
{
    GOOSEFEM_ASSERT(coor_a.dimension() == 2);
    GOOSEFEM_ASSERT(coor_b.dimension() == 2);
    GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
 
    std::vector<size_t> ret_a;
    std::vector<size_t> ret_b;
 
    for (size_t i = 0; i < coor_a.shape(0); ++i) {
 
        auto idx = xt::flatten_indices(xt::argwhere(
            xt::prod(xt::isclose(coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1)));
 
        for (auto& j : idx) {
            ret_a.push_back(i);
            ret_b.push_back(j);
        }
    }
 
    array_type::tensor<size_t, 2> ret = xt::empty<size_t>({size_t(2), ret_a.size()});
    for (size_t i = 0; i < ret_a.size(); ++i) {
        ret(0, i) = ret_a[i];
        ret(1, i) = ret_b[i];
    }
 
    return ret;
}
 
class ManualStitch {
public:
    ManualStitch() = default;
 
    template <class CA, class EA, class NA, class CB, class EB, class NB>
    ManualStitch(
        const CA& coor_a,
        const EA& conn_a,
        const NA& overlapping_nodes_a,
        const CB& coor_b,
        const EB& conn_b,
        const NB& overlapping_nodes_b,
        bool check_position = true,
        double rtol = 1e-5,
        double atol = 1e-8)
    {
        UNUSED(rtol);
        UNUSED(atol);
 
        GOOSEFEM_ASSERT(coor_a.dimension() == 2);
        GOOSEFEM_ASSERT(conn_a.dimension() == 2);
        GOOSEFEM_ASSERT(overlapping_nodes_a.dimension() == 1);
        GOOSEFEM_ASSERT(coor_b.dimension() == 2);
        GOOSEFEM_ASSERT(conn_b.dimension() == 2);
        GOOSEFEM_ASSERT(overlapping_nodes_b.dimension() == 1);
        GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape()));
        GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
        GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1));
 
        if (check_position) {
            GOOSEFEM_CHECK(xt::allclose(
                xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()),
                xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()),
                rtol,
                atol));
        }
 
        size_t nnda = coor_a.shape(0);
        size_t nndb = coor_b.shape(0);
        size_t ndim = coor_a.shape(1);
        size_t nelim = overlapping_nodes_a.size();
 
        size_t nela = conn_a.shape(0);
        size_t nelb = conn_b.shape(0);
        size_t nne = conn_a.shape(1);
 
        m_nel_a = nela;
        m_nel_b = nelb;
        m_nnd_a = nnda;
 
        array_type::tensor<size_t, 1> keep_b =
            xt::setdiff1d(xt::arange<size_t>(nndb), overlapping_nodes_b);
 
        m_map_b = xt::empty<size_t>({nndb});
        xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a;
        xt::view(m_map_b, xt::keep(keep_b)) = xt::arange<size_t>(keep_b.size()) + nnda;
 
        m_conn = xt::empty<size_t>({nela + nelb, nne});
        xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a;
        xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b);
 
        m_coor = xt::empty<size_t>({nnda + nndb - nelim, ndim});
        xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a;
        xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) =
            xt::view(coor_b, xt::keep(keep_b), xt::all());
    }
 
    size_t nmesh() const
    {
        return 2;
    }
 
    size_t nelem() const
    {
        return m_conn.shape(0);
    }
 
    size_t nnode() const
    {
        return m_coor.shape(0);
    }
 
    size_t nne() const
    {
        return m_conn.shape(1);
    }
 
    size_t ndim() const
    {
        return m_coor.shape(1);
    }
 
    const array_type::tensor<double, 2>& coor() const
    {
        return m_coor;
    }
 
    const array_type::tensor<size_t, 2>& conn() const
    {
        return m_conn;
    }
 
    array_type::tensor<size_t, 2> dofs() const
    {
        size_t nnode = this->nnode();
        size_t ndim = this->ndim();
        return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
    }
 
    std::vector<array_type::tensor<size_t, 1>> nodemap() const
    {
        std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
        for (size_t i = 0; i < this->nmesh(); ++i) {
            ret[i] = this->nodemap(i);
        }
        return ret;
    }
 
    std::vector<array_type::tensor<size_t, 1>> elemmap() const
    {
        std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
        for (size_t i = 0; i < this->nmesh(); ++i) {
            ret[i] = this->elemmap(i);
        }
        return ret;
    }
 
    array_type::tensor<size_t, 1> nodemap(size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index <= 1);
 
        if (mesh_index == 0) {
            return xt::arange<size_t>(m_nnd_a);
        }
 
        return m_map_b;
    }
 
    array_type::tensor<size_t, 1> elemmap(size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index <= 1);
 
        if (mesh_index == 0) {
            return xt::arange<size_t>(m_nel_a);
        }
 
        return xt::arange<size_t>(m_nel_b) + m_nel_a;
    }
 
    template <class T>
    T nodeset(const T& set, size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index <= 1);
 
        if (mesh_index == 0) {
            GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a);
            return set;
        }
 
        GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size());
        return detail::renum(set, m_map_b);
    }
 
    template <class T>
    T elemset(const T& set, size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index <= 1);
 
        if (mesh_index == 0) {
            GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a);
            return set;
        }
 
        GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b);
        return set + m_nel_a;
    }
 
private:
    array_type::tensor<double, 2> m_coor;
    array_type::tensor<size_t, 2> m_conn;
    array_type::tensor<size_t, 1> m_map_b;
    size_t m_nnd_a;
    size_t m_nel_a;
    size_t m_nel_b;
};
 
class Stitch {
public:
    Stitch(double rtol = 1e-5, double atol = 1e-8)
    {
        m_rtol = rtol;
        m_atol = atol;
    }
 
    template <class C, class E>
    void push_back(const C& coor, const E& conn)
    {
        GOOSEFEM_ASSERT(coor.dimension() == 2);
        GOOSEFEM_ASSERT(conn.dimension() == 2);
 
        if (m_map.size() == 0) {
            m_coor = coor;
            m_conn = conn;
            m_map.push_back(xt::eval(xt::arange<size_t>(coor.shape(0))));
            m_nel.push_back(conn.shape(0));
            m_el_offset.push_back(0);
            return;
        }
 
        auto overlap = overlapping(m_coor, coor, m_rtol, m_atol);
        size_t index = m_map.size();
 
        ManualStitch stitch(
            m_coor,
            m_conn,
            xt::eval(xt::view(overlap, 0, xt::all())),
            coor,
            conn,
            xt::eval(xt::view(overlap, 1, xt::all())),
            false);
 
        m_coor = stitch.coor();
        m_conn = stitch.conn();
        m_map.push_back(stitch.nodemap(1));
        m_nel.push_back(conn.shape(0));
        m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
    }
 
    size_t nmesh() const
    {
        return m_map.size();
    }
 
    size_t nelem() const
    {
        return m_conn.shape(0);
    }
 
    size_t nnode() const
    {
        return m_coor.shape(0);
    }
 
    size_t nne() const
    {
        return m_conn.shape(1);
    }
 
    size_t ndim() const
    {
        return m_coor.shape(1);
    }
 
    const array_type::tensor<double, 2>& coor() const
    {
        return m_coor;
    }
 
    const array_type::tensor<size_t, 2>& conn() const
    {
        return m_conn;
    }
 
    array_type::tensor<size_t, 2> dofs() const
    {
        size_t nnode = this->nnode();
        size_t ndim = this->ndim();
        return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
    }
 
    std::vector<array_type::tensor<size_t, 1>> nodemap() const
    {
        std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
        for (size_t i = 0; i < this->nmesh(); ++i) {
            ret[i] = this->nodemap(i);
        }
        return ret;
    }
 
    std::vector<array_type::tensor<size_t, 1>> elemmap() const
    {
        std::vector<array_type::tensor<size_t, 1>> ret(this->nmesh());
        for (size_t i = 0; i < this->nmesh(); ++i) {
            ret[i] = this->elemmap(i);
        }
        return ret;
    }
 
    array_type::tensor<size_t, 1> nodemap(size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index < m_map.size());
        return m_map[mesh_index];
    }
 
    array_type::tensor<size_t, 1> elemmap(size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index < m_map.size());
        return xt::arange<size_t>(m_nel[mesh_index]) + m_el_offset[mesh_index];
    }
 
    template <class T>
    T nodeset(const T& set, size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index < m_map.size());
        GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size());
        return detail::renum(set, m_map[mesh_index]);
    }
 
    template <class T>
    T elemset(const T& set, size_t mesh_index) const
    {
        GOOSEFEM_ASSERT(mesh_index < m_map.size());
        GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]);
        return set + m_el_offset[mesh_index];
    }
 
    template <class T>
    T nodeset(const std::vector<T>& set) const
    {
        GOOSEFEM_ASSERT(set.size() == m_map.size());
 
        size_t n = 0;
 
        for (size_t i = 0; i < set.size(); ++i) {
            n += set[i].size();
        }
 
        array_type::tensor<size_t, 1> ret = xt::empty<size_t>({n});
 
        n = 0;
 
        for (size_t i = 0; i < set.size(); ++i) {
            xt::view(ret, xt::range(n, n + set[i].size())) = this->nodeset(set[i], i);
            n += set[i].size();
        }
 
        return xt::unique(ret);
    }
 
    template <class T>
    T nodeset(std::initializer_list<T> set) const
    {
        return this->nodeset(std::vector<T>(set));
    }
 
    template <class T>
    T elemset(const std::vector<T>& set) const
    {
        GOOSEFEM_ASSERT(set.size() == m_map.size());
 
        size_t n = 0;
 
        for (size_t i = 0; i < set.size(); ++i) {
            n += set[i].size();
        }
 
        array_type::tensor<size_t, 1> ret = xt::empty<size_t>({n});
 
        n = 0;
 
        for (size_t i = 0; i < set.size(); ++i) {
            xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i);
            n += set[i].size();
        }
 
        return ret;
    }
 
    template <class T>
    T elemset(std::initializer_list<T> set) const
    {
        return this->elemset(std::vector<T>(set));
    }
 
protected:
    array_type::tensor<double, 2> m_coor; 
    array_type::tensor<size_t, 2> m_conn; 
    std::vector<array_type::tensor<size_t, 1>> m_map; 
    std::vector<size_t> m_nel; 
    std::vector<size_t> m_el_offset; 
    double m_rtol; 
    double m_atol; 
};
 
class Vstack : public Stitch {
public:
    Vstack(bool check_overlap = true, double rtol = 1e-5, double atol = 1e-8)
    {
        m_check_overlap = check_overlap;
        m_rtol = rtol;
        m_atol = atol;
    }
 
    template <class C, class E, class N>
    void push_back(const C& coor, const E& conn, const N& nodes_bot, const N& nodes_top)
    {
        if (m_map.size() == 0) {
            m_coor = coor;
            m_conn = conn;
            m_map.push_back(xt::eval(xt::arange<size_t>(coor.shape(0))));
            m_nel.push_back(conn.shape(0));
            m_el_offset.push_back(0);
            m_nodes_bot.push_back(nodes_bot);
            m_nodes_top.push_back(nodes_top);
            return;
        }
 
        GOOSEFEM_ASSERT(nodes_bot.size() == m_nodes_top.back().size());
 
        size_t index = m_map.size();
 
        double shift = xt::amax(xt::view(m_coor, xt::all(), 1))();
        auto x = coor;
        xt::view(x, xt::all(), 1) += shift;
 
        ManualStitch stitch(
            m_coor,
            m_conn,
            m_nodes_top.back(),
            x,
            conn,
            nodes_bot,
            m_check_overlap,
            m_rtol,
            m_atol);
 
        m_nodes_bot.push_back(stitch.nodeset(nodes_bot, 1));
        m_nodes_top.push_back(stitch.nodeset(nodes_top, 1));
 
        m_coor = stitch.coor();
        m_conn = stitch.conn();
        m_map.push_back(stitch.nodemap(1));
        m_nel.push_back(conn.shape(0));
        m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
    }
 
private:
    std::vector<array_type::tensor<size_t, 1>>
        m_nodes_bot; 
    std::vector<array_type::tensor<size_t, 1>>
        m_nodes_top; 
    bool m_check_overlap; 
};
 
class Reorder {
public:
    Reorder() = default;
 
    template <class T>
    Reorder(const std::initializer_list<T> args)
    {
        size_t n = 0;
        size_t i = 0;
 
        for (auto& arg : args) {
            if (arg.size() == 0) {
                continue;
            }
            n = std::max(n, xt::amax(arg)() + 1);
        }
 
#ifdef GOOSEFEM_ENABLE_ASSERT
        for (auto& arg : args) {
            GOOSEFEM_ASSERT(is_unique(arg));
        }
#endif
 
        m_renum = xt::empty<size_t>({n});
 
        for (auto& arg : args) {
            for (auto& j : arg) {
                m_renum(j) = i;
                ++i;
            }
        }
    }
 
    template <class T>
    T apply(const T& list) const
    {
        T ret = T::from_shape(list.shape());
 
        auto jt = ret.begin();
 
        for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
            *jt = m_renum(*it);
        }
 
        return ret;
    }
 
    const array_type::tensor<size_t, 1>& index() const
    {
        return m_renum;
    }
 
private:
    array_type::tensor<size_t, 1> m_renum;
};
 
template <class E>
inline array_type::tensor<size_t, 1> coordination(const E& conn)
{
    GOOSEFEM_ASSERT(conn.dimension() == 2);
 
    size_t nnode = xt::amax(conn)() + 1;
 
    array_type::tensor<size_t, 1> N = xt::zeros<size_t>({nnode});
 
    for (auto it = conn.begin(); it != conn.end(); ++it) {
        N(*it) += 1;
    }
 
    return N;
}
 
template <class E>
inline std::vector<std::vector<size_t>> elem2node(const E& conn, bool sorted = true)
{
    auto N = coordination(conn);
    auto nnode = N.size();
 
    std::vector<std::vector<size_t>> ret(nnode);
    for (size_t i = 0; i < nnode; ++i) {
        ret[i].reserve(N(i));
    }
 
    for (size_t e = 0; e < conn.shape(0); ++e) {
        for (size_t m = 0; m < conn.shape(1); ++m) {
            ret[conn(e, m)].push_back(e);
        }
    }
 
    if (sorted) {
        for (auto& row : ret) {
            std::sort(row.begin(), row.end());
        }
    }
 
    return ret;
}
 
template <class E, class D>
inline std::vector<std::vector<size_t>> elem2node(const E& conn, const D& dofs, bool sorted = true)
{
    size_t nnode = dofs.shape(0);
    auto ret = elem2node(conn, sorted);
    auto nties = nodaltyings(dofs);
 
    for (size_t m = 0; m < nnode; ++m) {
        if (nties[m].size() <= 1) {
            continue;
        }
        if (m > nties[m][0]) {
            continue;
        }
        size_t n = ret[m].size();
        for (size_t j = 1; j < nties[m].size(); ++j) {
            n += ret[nties[m][j]].size();
        }
        ret[m].reserve(n);
        for (size_t j = 1; j < nties[m].size(); ++j) {
            ret[m].insert(ret[m].end(), ret[nties[m][j]].begin(), ret[nties[m][j]].end());
        }
        if (sorted) {
            std::sort(ret[m].begin(), ret[m].end());
        }
        for (size_t j = 1; j < nties[m].size(); ++j) {
            ret[nties[m][j]] = ret[m];
        }
    }
 
    return ret;
}
 
template <class D>
inline std::vector<std::vector<size_t>> node2dof(const D& dofs, bool sorted = true)
{
    return elem2node(dofs, sorted);
}
 
template <class C, class E>
inline array_type::tensor<double, 2> edgesize(const C& coor, const E& conn, ElementType type)
{
    GOOSEFEM_ASSERT(coor.dimension() == 2);
    GOOSEFEM_ASSERT(conn.dimension() == 2);
    GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
 
    if (type == ElementType::Quad4) {
        GOOSEFEM_ASSERT(coor.shape(1) == 2ul);
        GOOSEFEM_ASSERT(conn.shape(1) == 4ul);
        array_type::tensor<size_t, 1> n0 = xt::view(conn, xt::all(), 0);
        array_type::tensor<size_t, 1> n1 = xt::view(conn, xt::all(), 1);
        array_type::tensor<size_t, 1> n2 = xt::view(conn, xt::all(), 2);
        array_type::tensor<size_t, 1> n3 = xt::view(conn, xt::all(), 3);
        array_type::tensor<double, 1> x0 = xt::view(coor, xt::keep(n0), 0);
        array_type::tensor<double, 1> x1 = xt::view(coor, xt::keep(n1), 0);
        array_type::tensor<double, 1> x2 = xt::view(coor, xt::keep(n2), 0);
        array_type::tensor<double, 1> x3 = xt::view(coor, xt::keep(n3), 0);
        array_type::tensor<double, 1> y0 = xt::view(coor, xt::keep(n0), 1);
        array_type::tensor<double, 1> y1 = xt::view(coor, xt::keep(n1), 1);
        array_type::tensor<double, 1> y2 = xt::view(coor, xt::keep(n2), 1);
        array_type::tensor<double, 1> y3 = xt::view(coor, xt::keep(n3), 1);
        array_type::tensor<double, 2> ret = xt::empty<double>(conn.shape());
        xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0));
        xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0));
        xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0));
        xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0));
        return ret;
    }
 
    throw std::runtime_error("Element-type not implemented");
}
 
template <class C, class E>
inline array_type::tensor<double, 2> edgesize(const C& coor, const E& conn)
{
    return edgesize(coor, conn, defaultElementType(coor, conn));
}
 
template <class C, class E>
inline array_type::tensor<double, 2> centers(const C& coor, const E& conn, ElementType type)
{
    GOOSEFEM_ASSERT(coor.dimension() == 2);
    GOOSEFEM_ASSERT(conn.dimension() == 2);
    GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
    array_type::tensor<double, 2> ret = xt::zeros<double>({conn.shape(0), coor.shape(1)});
 
    if (type == ElementType::Quad4) {
        GOOSEFEM_ASSERT(coor.shape(1) == 2);
        GOOSEFEM_ASSERT(conn.shape(1) == 4);
        for (size_t i = 0; i < 4; ++i) {
            auto n = xt::view(conn, xt::all(), i);
            ret += xt::view(coor, xt::keep(n), xt::all());
        }
        ret /= 4.0;
        return ret;
    }
 
    throw std::runtime_error("Element-type not implemented");
}
 
template <class C, class E>
inline array_type::tensor<double, 2> centers(const C& coor, const E& conn)
{
    return centers(coor, conn, defaultElementType(coor, conn));
}
 
template <class T, class C, class E>
inline array_type::tensor<size_t, 1>
elemmap2nodemap(const T& elem_map, const C& coor, const E& conn, ElementType type)
{
    GOOSEFEM_ASSERT(elem_map.dimension() == 1);
    GOOSEFEM_ASSERT(coor.dimension() == 2);
    GOOSEFEM_ASSERT(conn.dimension() == 2);
    GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
    GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0));
    size_t N = coor.shape(0);
 
    array_type::tensor<size_t, 1> ret = N * xt::ones<size_t>({N});
 
    if (type == ElementType::Quad4) {
        GOOSEFEM_ASSERT(coor.shape(1) == 2);
        GOOSEFEM_ASSERT(conn.shape(1) == 4);
 
        for (size_t i = 0; i < 4; ++i) {
            array_type::tensor<size_t, 1> t = N * xt::ones<size_t>({N});
            auto old_nd = xt::view(conn, xt::all(), i);
            auto new_nd = xt::view(conn, xt::keep(elem_map), i);
            xt::view(t, xt::keep(old_nd)) = new_nd;
            ret = xt::where(xt::equal(ret, N), t, ret);
        }
 
        return ret;
    }
 
    throw std::runtime_error("Element-type not implemented");
}
 
template <class T, class C, class E>
inline array_type::tensor<size_t, 1>
elemmap2nodemap(const T& elem_map, const C& coor, const E& conn)
{
    return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn));
}
 
template <class C, class E>
inline array_type::tensor<double, 2> nodal_mass(const C& coor, const E& conn, ElementType type)
{
    auto dof = dofs(coor.shape(0), coor.shape(1));
    GooseFEM::MatrixDiagonal M(conn, dof);
    GooseFEM::Vector vector(conn, dof);
    array_type::tensor<double, 2> rho = xt::ones<double>(conn.shape());
 
    if (type == ElementType::Quad4) {
        GooseFEM::Element::Quad4::Quadrature quad(
            vector.AsElement(coor),
            GooseFEM::Element::Quad4::Nodal::xi(),
            GooseFEM::Element::Quad4::Nodal::w());
        M.assemble(quad.Int_N_scalar_NT_dV(rho));
    }
    else {
        throw std::runtime_error("Element-type not implemented");
    }
 
    return vector.AsNode(M.data());
}
 
template <class C, class E>
inline array_type::tensor<double, 2> nodal_mass(const C& coor, const E& conn)
{
    return nodal_mass(coor, conn, defaultElementType(coor, conn));
}
 
namespace detail {
 
// todo: remove after upstream fix
template <class T>
array_type::tensor<double, 1> average_axis_0(const T& data, const T& weights)
{
    GOOSEFEM_ASSERT(data.dimension() == 2);
    GOOSEFEM_ASSERT(xt::has_shape(data, weights.shape()));
 
    array_type::tensor<double, 1> ret = xt::zeros<double>({weights.shape(1)});
 
    for (size_t j = 0; j < weights.shape(1); ++j) {
        double norm = 0.0;
        for (size_t i = 0; i < weights.shape(0); ++i) {
            ret(j) += data(i, j) * weights(i, j);
            norm += weights(i, j);
        }
        ret(j) /= norm;
    }
    return ret;
}
 
} // namespace detail
 
template <class C, class E>
inline array_type::tensor<double, 1>
center_of_gravity(const C& coor, const E& conn, ElementType type)
{
    // todo: remove after upstream fix
    return detail::average_axis_0(coor, nodal_mass(coor, conn, type));
    // return xt::average(coor, nodal_mass(coor, conn, type), 0);
}
 
template <class C, class E>
inline array_type::tensor<double, 1> center_of_gravity(const C& coor, const E& conn)
{
    // todo: remove after upstream fix
    return detail::average_axis_0(coor, nodal_mass(coor, conn, defaultElementType(coor, conn)));
    // return xt::average(coor, nodal_mass(coor, conn, defaultElementType(coor, conn)), 0);
}
 
template <class D>
inline std::vector<std::vector<size_t>> nodaltyings(const D& dofs)
{
    size_t nnode = dofs.shape(0);
    size_t ndim = dofs.shape(1);
    auto nodemap = node2dof(dofs);
    std::vector<std::vector<size_t>> ret(nnode);
 
    for (size_t m = 0; m < nnode; ++m) {
        auto r = nodemap[dofs(m, 0)];
        std::sort(r.begin(), r.end());
        ret[m] = r;
#ifdef GOOSEFEM_ENABLE_ASSERT
        for (size_t i = 1; i < ndim; ++i) {
            auto t = nodemap[dofs(m, i)];
            std::sort(t.begin(), t.end());
            GOOSEFEM_ASSERT(std::equal(r.begin(), r.end(), t.begin()));
        }
#endif
    }
 
    return ret;
}
 
} // namespace Mesh
} // namespace GooseFEM
 
#endif