diff --git a/include/GooseFEM/Vector.h b/include/GooseFEM/Vector.h
index e0f081a..06a2fde 100644
--- a/include/GooseFEM/Vector.h
+++ b/include/GooseFEM/Vector.h
@@ -1,337 +1,345 @@
 /**
 Methods to switch between storage types based on a mesh.
 
 \file Vector.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_VECTOR_H
 #define GOOSEFEM_VECTOR_H
 
 #include "config.h"
 
 namespace GooseFEM {
 
 /**
 Class to switch between storage types. In particular:
 
 -   "dofval": DOF values [#ndof].
 -   "nodevec": nodal vectors [#nnode, #ndim].
 -   "elemvec": nodal vectors stored per element [#nelem, #nne, #ndim].
 */
 class Vector {
 public:
 
     Vector() = default;
 
     /**
     Constructor.
 
     \param conn connectivity [#nelem, #nne].
     \param dofs DOFs per node [#nnode, #ndim].
     */
     template <class S, class T>
     Vector(const S& conn, const T& dofs);
 
     /**
     \return  Number of elements.
     */
     size_t nelem() const;
 
     /**
     \return  Number of nodes per element.
     */
     size_t nne() const;
 
     /**
     \return  Number of nodes.
     */
     size_t nnode() const;
 
     /**
     \return  Number of dimensions.
     */
     size_t ndim() const;
 
     /**
     \return  Number of DOFs.
     */
     size_t ndof() const;
 
     /**
     \return Connectivity (nodes per element) [#nelem, #nne].
     */
     xt::xtensor<size_t, 2> conn() const;
 
     /**
     \return DOFs per node [#nnode, #ndim]
     */
     xt::xtensor<size_t, 2> dofs() const;
 
     /**
     Copy "nodevec" to another "nodevec".
 
     \param nodevec_src input [#nnode, #ndim]
     \param nodevec_dest input [#nnode, #ndim]
     \return nodevec output [#nnode, #ndim]
     */
     template <class T>
     T Copy(const T& nodevec_src, const T& nodevec_dest) const;
 
     /**
     Copy "nodevec" to another "nodevec".
 
     \param nodevec_src input [#nnode, #ndim]
     \param nodevec_dest output [#nnode, #ndim]
     */
     template <class T>
     void copy(const T& nodevec_src, T& nodevec_dest) const;
 
     /**
     Convert "nodevec" or "elemvec" to "dofval" (overwrite entries that occur more than once).
 
     \param arg nodevec [#nnode, #ndim] or elemvec [#nelem, #nne, #ndim]
     \return dofval [#ndof]
     */
     template <class T>
     xt::xtensor<double, 1> AsDofs(const T& arg) const;
 
     /**
     Convert "nodevec" or "elemvec" to "dofval" (overwrite entries that occur more than once).
 
     \param arg nodevec [#nnode, #ndim] or elemvec [#nelem, #nne, #ndim]
     \param dofval (output) [#ndof]
     */
     template <class T, class R>
-    typename std::enable_if_t<!xt::has_fixed_rank_t<T>::value>
-    asDofs(const T& arg, R& dofval) const;
-
-    /** copydoc asDofs(const T&, R&) const */
-    template <class T, class R>
-    typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 1>>
-    asDofs(const T& arg, R& dofval) const;
-
-    /** copydoc asDofs(const T&, R&) const */
-    template <class T, class R>
-    typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 2>>
-    asDofs(const T& arg, R& dofval) const;
-
-    /** copydoc asDofs(const T&, R&) const */
-    template <class T, class R>
-    typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 3>>
-    asDofs(const T& arg, R& dofval) const;
+    void asDofs(const T& arg, R& dofval) const;
 
     /**
     Convert "dofval" or "elemvec" to "nodevec" (overwrite entries that occur more than once).
 
     \param arg dofval [#ndof] or elemvec [#nelem, #nne, #ndim]
     \return nodevec output [#nnode, #ndim]
     */
     template <class T>
     xt::xtensor<double, 2> AsNode(const T& arg) const;
 
     /**
     Convert "dofval" or "elemvec" to "nodevec" (overwrite entries that occur more than once).
 
     \param arg dofval [#ndof] or elemvec [#nelem, #nne, #ndim]
     \param nodevec output [#nnode, #ndim]
     */
     template <class T, class R>
     void asNode(const T& arg, R& nodevec) const;
 
     /**
     Convert "dofval" or "nodevec" to "elemvec" (overwrite entries that occur more than once).
 
     \param arg dofval [#ndof] or nodevec [#nnode, #ndim].
     \param elemvec output [#nelem, #nne, #ndim].
     */
     template <class T, class R>
     void asElement(const T& arg, R& elemvec) const;
 
     /**
     Convert "dofval" or "nodevec" to "elemvec" (overwrite entries that occur more than once).
 
     \param arg dofval [#ndof] or nodevec [#nnode, #ndim].
     \return elemvec output [#nelem, #nne, #ndim].
     */
     template <class T>
     xt::xtensor<double, 3> AsElement(const T& arg) const;
 
     /**
     Assemble "nodevec" or "elemvec" to "dofval" (adds entries that occur more that once).
 
     \param arg nodevec [#nnode, #ndim] or elemvec [#nelem, #nne, #ndim]
     \return dofval output [#ndof]
     */
     template <class T>
     xt::xtensor<double, 1> AssembleDofs(const T& arg) const;
 
     /**
     Assemble "nodevec" or "elemvec" to "dofval" (adds entries that occur more that once).
 
     \param arg nodevec [#nnode, #ndim] or elemvec [#nelem, #nne, #ndim]
     \param dofval output [#ndof]
     */
     template <class T, class R>
     void assembleDofs(const T& arg, R& dofval) const;
 
     /**
     Assemble "elemvec" to "nodevec" (adds entries that occur more that once.
 
     \param arg elemvec [#nelem, #nne, #ndim]
     \return nodevec output [#nnode, #ndim]
     */
     template <class T>
     xt::xtensor<double, 2> AssembleNode(const T& arg) const;
 
     /**
     Assemble "elemvec" to "nodevec" (adds entries that occur more that once.
 
     \param arg elemvec [#nelem, #nne, #ndim]
     \param nodevec output [#nnode, #ndim]
     */
     template <class T, class R>
     void assembleNode(const T& arg, R& nodevec) const;
 
     /**
     Shape of "dofval".
 
     \return [#ndof]
     */
     std::array<size_t, 1> shape_dofval() const;
 
     /**
     Shape of "nodevec".
 
     \return [#nnode, #ndim]
     */
     std::array<size_t, 2> shape_nodevec() const;
 
     /**
     Shape of "elemvec".
 
     \return [#nelem, #nne, #ndim]
     */
     std::array<size_t, 3> shape_elemvec() const;
 
     /**
     Shape of "elemmat".
 
     \return [#nelem, #nne * #ndim, #nne * #ndim]
     */
     std::array<size_t, 3> shape_elemmat() const;
 
     /**
     Allocated "dofval".
 
     \return [#ndof]
     */
     xt::xtensor<double, 1> allocate_dofval() const;
 
     /**
     Allocated and initialised "dofval".
 
     \param val value to which to initialise.
     \return [#ndof]
     */
     xt::xtensor<double, 1> allocate_dofval(double val) const;
 
     /**
     Allocated "nodevec".
 
     \return [#nnode, #ndim]
     */
     xt::xtensor<double, 2> allocate_nodevec() const;
 
     /**
     Allocated and initialised "nodevec".
 
     \param val value to which to initialise.
     \return [#nnode, #ndim]
     */
     xt::xtensor<double, 2> allocate_nodevec(double val) const;
 
     /**
     Allocated "elemvec".
 
     \return [#nelem, #nne, #ndim]
     */
     xt::xtensor<double, 3> allocate_elemvec() const;
 
     /**
     Allocated and initialised "elemvec".
 
     \param val value to which to initialise.
     \return [#nelem, #nne, #ndim]
     */
     xt::xtensor<double, 3> allocate_elemvec(double val) const;
 
     /**
     Allocated "elemmat".
 
     \return [#nelem, #nne * #ndim, #nne * #ndim]
     */
     xt::xtensor<double, 3> allocate_elemmat() const;
 
     /**
     Allocated and initialised "elemmat".
 
     \param val value to which to initialise.
     \return [#nelem, #nne * #ndim, #nne * #ndim]
     */
     xt::xtensor<double, 3> allocate_elemmat(double val) const;
 
-protected:
+private:
+
+    /** Distribution of \copydoc asDofs(const T&, R&) const */
+    template <class T, class R>
+    typename std::enable_if_t<!xt::has_fixed_rank_t<T>::value>
+    asDofs_impl(const T& arg, R& dofval) const;
+
+    /** Distribution of \copydoc asDofs(const T&, R&) const */
+    template <class T, class R>
+    typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 1>>
+    asDofs_impl(const T& arg, R& dofval) const;
+
+    /** Distribution of \copydoc asDofs(const T&, R&) const */
+    template <class T, class R>
+    typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 2>>
+    asDofs_impl(const T& arg, R& dofval) const;
+
+    /** Distribution of \copydoc asDofs(const T&, R&) const */
+    template <class T, class R>
+    typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 3>>
+    asDofs_impl(const T& arg, R& dofval) const;
+
+    /** Implementation for 'nodevec' input of \copydoc asDofs(const T&, R&) const */
+    template <class T, class R>
+    void asDofs_impl_dofval(const T& arg, R& dofval) const;
 
     /** Implementation for 'nodevec' input of \copydoc asDofs(const T&, R&) const */
     template <class T, class R>
-    void asDofs_nodevec(const T& arg, R& dofval) const;
+    void asDofs_impl_nodevec(const T& arg, R& dofval) const;
 
     /** Implementation for 'elemvec' input of \copydoc asDofs(const T&, R&) const */
     template <class T, class R>
-    void asDofs_elemvec(const T& arg, R& dofval) const;
+    void asDofs_impl_elemvec(const T& arg, R& dofval) const;
 
     /** Implementation for 'dofval' input of \copydoc asNode(const T&, R&) const */
     template <class T, class R>
     void asNode_dofval(const T& arg, R& nodevec) const;
 
     /** Implementation for 'elemvec' input of \copydoc asNode(const T&, R&) const */
     template <class T, class R>
     void asNode_elemvec(const T& arg, R& nodevec) const;
 
     /** Implementation for 'dofval' input of \copydoc asElement(const T&, R&) const */
     template <class T, class R>
     void asElement_dofval(const T& arg, R& elemvec) const;
 
     /** Implementation for 'nodevec' input of \copydoc asElement(const T&, R&) const */
     template <class T, class R>
     void asElement_nodevec(const T& arg, R& elemvec) const;
 
     /** Implementation for 'nodevec' input of \copydoc assembleDofs(const T&, R&) const */
     template <class T, class R>
     void assembleDofs_nodevec(const T& arg, R& dofval) const;
 
     /** Implementation for 'elemvec' input of \copydoc assembleDofs(const T&, R&) const */
     template <class T, class R>
     void assembleDofs_elemvec(const T& arg, R& dofval) const;
 
     /** Implementation for 'elemvec' input of \copydoc assembleNode(const T&, R&) const */
     template <class T, class R>
     void assembleNode_elemvec(const T& arg, R& nodevec) const;
 
 protected:
     xt::xtensor<size_t, 2> m_conn; ///< See conn()
     xt::xtensor<size_t, 2> m_dofs; ///< See dofs()
     size_t m_nelem; ///< See #nelem
     size_t m_nne;   ///< See #nne
     size_t m_nnode; ///< See #nnode
     size_t m_ndim;  ///< See #ndim
     size_t m_ndof;  ///< See #ndof
 };
 
 } // namespace GooseFEM
 
 #include "Vector.hpp"
 
 #endif
diff --git a/include/GooseFEM/Vector.hpp b/include/GooseFEM/Vector.hpp
index 8cb70c6..bf36686 100644
--- a/include/GooseFEM/Vector.hpp
+++ b/include/GooseFEM/Vector.hpp
@@ -1,462 +1,481 @@
 /**
 Implementation of Vector.h
 
 \file Vector.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_VECTOR_HPP
 #define GOOSEFEM_VECTOR_HPP
 
 #include "Vector.h"
 
 namespace GooseFEM {
 
 template <class S, class T>
 inline Vector::Vector(const S& conn, const T& dofs)
 {
     GOOSEFEM_ASSERT(conn.dimension() == 2);
     GOOSEFEM_ASSERT(dofs.dimension() == 2);
 
     m_conn = conn;
     m_dofs = dofs;
 
     m_nelem = m_conn.shape(0);
     m_nne = m_conn.shape(1);
     m_nnode = m_dofs.shape(0);
     m_ndim = m_dofs.shape(1);
     m_ndof = xt::amax(m_dofs)() + 1;
 
     GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
     GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
 }
 
 inline size_t Vector::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t Vector::nne() const
 {
     return m_nne;
 }
 
 inline size_t Vector::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t Vector::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t Vector::ndof() const
 {
     return m_ndof;
 }
 
 inline xt::xtensor<size_t, 2> Vector::conn() const
 {
     return m_conn;
 }
 
 inline xt::xtensor<size_t, 2> Vector::dofs() const
 {
     return m_dofs;
 }
 
 template <class T>
 inline void Vector::copy(const T& nodevec_src, T& nodevec_dest) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_src, this->shape_nodevec()));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec_dest, this->shape_nodevec()));
 
     xt::noalias(nodevec_dest) = nodevec_src;
 }
 
+template <class T, class R>
+void Vector::asDofs(const T& arg, R& dofval) const
+{
+    this->asDofs_impl(arg, dofval);
+}
+
+// asDofs : distribution
+
 template <class T, class R>
 inline
 typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 1>>
-Vector::asDofs(const T& arg, R& dofval) const
+Vector::asDofs_impl(const T& arg, R& dofval) const
 {
     dofval = arg;
 }
 
 template <class T, class R>
 inline
 typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 2>>
-Vector::asDofs(const T& arg, R& dofval) const
+Vector::asDofs_impl(const T& arg, R& dofval) const
 {
-    this->asDofs_nodevec(arg, dofval);
+    this->asDofs_impl_nodevec(arg, dofval);
 }
 
 template <class T, class R>
 inline
 typename std::enable_if_t<std::equal<xt::get_rank<T>::value, 3>>
-Vector::asDofs(const T& arg, R& dofval) const
+Vector::asDofs_impl(const T& arg, R& dofval) const
 {
-    this->asDofs_elemvec(arg, dofval);
+    this->asDofs_impl_elemvec(arg, dofval);
 }
 
 template <class T, class R>
 inline
 typename std::enable_if_t<!xt::has_fixed_rank_t<T>::value>
-Vector::asDofs(const T& arg, R& dofval) const
+Vector::asDofs_impl(const T& arg, R& dofval) const
 {
     if (arg.dimension() == 1) {
         dofval = arg;
     }
     else if (arg.dimension() == 2) {
-        this->asDofs_nodevec(arg, dofval);
+        this->asDofs_impl_nodevec(arg, dofval);
     }
     else if (arg.dimension() == 3) {
-        this->asDofs_elemvec(arg, dofval);
+        this->asDofs_impl_elemvec(arg, dofval);
     }
     else {
         throw std::runtime_error("Vector::asDofs unknown dimension first argument");
     }
 }
 
 template <class T, class R>
 inline void Vector::asNode(const T& arg, R& nodevec) const
 {
     if (arg.dimension() == 1) {
         this->asNode_dofval(arg, nodevec);
     }
     else if (arg.dimension() == 3) {
         this->asNode_elemvec(arg, nodevec);
     }
     else if (arg.dimension() == 2) {
         nodevec = arg;
     }
     else {
         throw std::runtime_error("Vector::asNode unknown dimension first argument");
     }
 }
 
 template <class T, class R>
 inline void Vector::asElement(const T& arg, R& elemvec) const
 {
     if (arg.dimension() == 1) {
         this->asElement_dofval(arg, elemvec);
     }
     else if (arg.dimension() == 2) {
         this->asElement_nodevec(arg, elemvec);
     }
     else if (arg.dimension() == 3) {
         elemvec = arg;
     }
     else {
         throw std::runtime_error("Vector::asElement unknown dimension first argument");
     }
 }
 
 template <class T, class R>
 inline void Vector::assembleDofs(const T& arg, R& dofval) const
 {
     if (arg.dimension() == 2) {
         this->assembleDofs_nodevec(arg, dofval);
     }
     else if (arg.dimension() == 3) {
         this->assembleDofs_elemvec(arg, dofval);
     }
     else if (arg.dimension() == 1) {
         dofval = arg;
     }
     else {
         throw std::runtime_error("Vector::assembleDofs unknown dimension first argument");
     }
 }
 
 template <class T, class R>
 inline void Vector::assembleNode(const T& arg, R& nodevec) const
 {
     if (arg.dimension() == 3) {
         this->assembleNode_elemvec(arg, nodevec);
     }
     else if (arg.dimension() == 2) {
         nodevec = arg;
     }
     else {
         throw std::runtime_error("Vector::assembleNode unknown dimension first argument");
     }
 }
 
+// asDofs : implementation
+
 template <class T, class R>
-inline void Vector::asDofs_nodevec(const T& nodevec, R& dofval) const
+inline void Vector::asDofs_impl_dofval(const T& arg, R& dofval) const
 {
-    GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec()));
+    GOOSEFEM_ASSERT(xt::has_shape(arg, this->shape_dofval()));
+    GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval()));
+
+    dofval = arg;
+}
+
+template <class T, class R>
+inline void Vector::asDofs_impl_nodevec(const T& arg, R& dofval) const
+{
+    GOOSEFEM_ASSERT(xt::has_shape(arg, this->shape_nodevec()));
     GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval()));
 
     dofval.fill(0.0);
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
-            dofval(m_dofs(m, i)) = nodevec(m, i);
+            dofval(m_dofs(m, i)) = arg(m, i);
         }
     }
 }
 
 template <class T, class R>
-inline void Vector::asDofs_elemvec(const T& elemvec, R& dofval) const
+inline void Vector::asDofs_impl_elemvec(const T& arg, R& dofval) const
 {
-    GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
+    GOOSEFEM_ASSERT(xt::has_shape(arg, this->shape_elemvec()));
     GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval()));
 
     dofval.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
-                dofval(m_dofs(m_conn(e, m), i)) = elemvec(e, m, i);
+                dofval(m_dofs(m_conn(e, m), i)) = arg(e, m, i);
             }
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::asNode_dofval(const T& dofval, R& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval()));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec()));
 
     #pragma omp parallel for
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             nodevec(m, i) = dofval(m_dofs(m, i));
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::asNode_elemvec(const T& elemvec, R& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec()));
 
     nodevec.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 nodevec(m_conn(e, m), i) = elemvec(e, m, i);
             }
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::asElement_dofval(const T& dofval, R& elemvec) const
 {
     GOOSEFEM_ASSERT(dofval.size() == m_ndof);
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = dofval(m_dofs(m_conn(e, m), i));
             }
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::asElement_nodevec(const T& nodevec,R& elemvec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec()));
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 elemvec(e, m, i) = nodevec(m_conn(e, m), i);
             }
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::assembleDofs_nodevec(const T& nodevec, R& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec()));
     GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval()));
 
     dofval.fill(0.0);
 
     for (size_t m = 0; m < m_nnode; ++m) {
         for (size_t i = 0; i < m_ndim; ++i) {
             dofval(m_dofs(m, i)) += nodevec(m, i);
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::assembleDofs_elemvec(const T& elemvec, R& dofval) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
     GOOSEFEM_ASSERT(xt::has_shape(dofval, this->shape_dofval()));
 
     dofval.fill(0.0);
 
     for (size_t e = 0; e < m_nelem; ++e) {
         for (size_t m = 0; m < m_nne; ++m) {
             for (size_t i = 0; i < m_ndim; ++i) {
                 dofval(m_dofs(m_conn(e, m), i)) += elemvec(e, m, i);
             }
         }
     }
 }
 
 template <class T, class R>
 inline void Vector::assembleNode_elemvec(const T& elemvec, R& nodevec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
     GOOSEFEM_ASSERT(xt::has_shape(nodevec, this->shape_nodevec()));
 
     xt::xtensor<double, 1> dofval = this->AssembleDofs(elemvec);
     this->asNode(dofval, nodevec);
 }
 
 template <class T>
 inline xt::xtensor<double, 1> Vector::AsDofs(const T& arg) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>(this->shape_dofval());
     this->asDofs(arg, dofval);
     return dofval;
 }
 
 template <class T>
 inline xt::xtensor<double, 2> Vector::AsNode(const T& arg) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>(this->shape_nodevec());
     this->asNode(arg, nodevec);
     return nodevec;
 }
 
 template <class T>
 inline xt::xtensor<double, 3> Vector::AsElement(const T& arg) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>(this->shape_elemvec());
     this->asElement(arg, elemvec);
     return elemvec;
 }
 
 template <class T>
 inline xt::xtensor<double, 1> Vector::AssembleDofs(const T& arg) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>(this->shape_dofval());
     this->assembleDofs(arg, dofval);
     return dofval;
 }
 
 template <class T>
 inline xt::xtensor<double, 2> Vector::AssembleNode(const T& arg) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>(this->shape_nodevec());
     this->assembleNode(arg, nodevec);
     return nodevec;
 }
 
 template <class T>
 inline T Vector::Copy(const T& nodevec_src, const T& nodevec_dest) const
 {
     T ret = T::from_shape(nodevec_dest.shape());
     this->copy(nodevec_src, ret);
     return ret;
 }
 
 inline std::array<size_t, 1> Vector::shape_dofval() const
 {
     std::array<size_t, 1> shape;
     shape[0] = m_ndof;
     return shape;
 }
 
 inline std::array<size_t, 2> Vector::shape_nodevec() const
 {
     std::array<size_t, 2> shape;
     shape[0] = m_nnode;
     shape[1] = m_ndim;
     return shape;
 }
 
 inline std::array<size_t, 3> Vector::shape_elemvec() const
 {
     std::array<size_t, 3> shape;
     shape[0] = m_nelem;
     shape[1] = m_nne;
     shape[2] = m_ndim;
     return shape;
 }
 
 inline std::array<size_t, 3> Vector::shape_elemmat() const
 {
     std::array<size_t, 3> shape;
     shape[0] = m_nelem;
     shape[1] = m_nne * m_ndim;
     shape[2] = m_nne * m_ndim;
     return shape;
 }
 
 inline xt::xtensor<double, 1> Vector::allocate_dofval() const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>(this->shape_dofval());
     return dofval;
 }
 
 inline xt::xtensor<double, 2> Vector::allocate_nodevec() const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>(this->shape_nodevec());
     return nodevec;
 }
 
 inline xt::xtensor<double, 3> Vector::allocate_elemvec() const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>(this->shape_elemvec());
     return elemvec;
 }
 
 inline xt::xtensor<double, 3> Vector::allocate_elemmat() const
 {
     xt::xtensor<double, 3> elemmat = xt::empty<double>(this->shape_elemmat());
     return elemmat;
 }
 
 inline xt::xtensor<double, 1> Vector::allocate_dofval(double val) const
 {
     xt::xtensor<double, 1> dofval = xt::empty<double>(this->shape_dofval());
     dofval.fill(val);
     return dofval;
 }
 
 inline xt::xtensor<double, 2> Vector::allocate_nodevec(double val) const
 {
     xt::xtensor<double, 2> nodevec = xt::empty<double>(this->shape_nodevec());
     nodevec.fill(val);
     return nodevec;
 }
 
 inline xt::xtensor<double, 3> Vector::allocate_elemvec(double val) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>(this->shape_elemvec());
     elemvec.fill(val);
     return elemvec;
 }
 
 inline xt::xtensor<double, 3> Vector::allocate_elemmat(double val) const
 {
     xt::xtensor<double, 3> elemmat = xt::empty<double>(this->shape_elemmat());
     elemmat.fill(val);
     return elemmat;
 }
 
 } // namespace GooseFEM
 
 #endif