diff --git a/include/GooseFEM/Element.h b/include/GooseFEM/Element.h
index 7fbb436..a5f8f29 100644
--- a/include/GooseFEM/Element.h
+++ b/include/GooseFEM/Element.h
@@ -1,645 +1,668 @@
 /**
 Convenience methods for integration point data.
 
 \file Element.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_ELEMENT_H
 #define GOOSEFEM_ELEMENT_H
 
 #include "config.h"
 #include "Allocate.h"
 
 namespace GooseFEM {
 
 /**
 Element quadrature and interpolation.
 */
 namespace Element {
 
 /**
 Convert nodal vector with ("nodevec", shape:``[nnode, ndim]``) to nodal vector stored per element
 ("elemvec", shape: ``[nelem, nne, ndim]``).
 
 \param conn Connectivity.
 \param nodevec "nodevec".
 \return "elemvec".
 */
 inline xt::xtensor<double, 3> asElementVector(
     const xt::xtensor<size_t, 2>& conn, const xt::xtensor<double, 2>& nodevec);
 
 /**
 Assemble nodal vector stored per element ("elemvec", shape ``[nelem, nne, ndim]``) to nodal vector
 ("nodevec", shape ``[nnode, ndim]``).
 
 \param conn Connectivity.
 \param elemvec "elemvec".
 \return "nodevec".
 */
 inline xt::xtensor<double, 2> assembleNodeVector(
     const xt::xtensor<size_t, 2>& conn, const xt::xtensor<double, 3>& elemvec);
 
 /**
 Check that DOFs leave no holes.
 
 \param dofs DOFs ("nodevec")
 \return ``true`` if there are no holds.
 */
 template <class E>
 inline bool isSequential(const E& dofs);
 
 
 /**
 Check that all of the matrices stored per elemmat (shape: ``[nelem, nne * ndim, nne * ndim]``)
 are diagonal.
 
 \param elemmat Element-vectors ("elemmat")
 \return ``true`` if all element matrices are diagonal.
 */
 bool isDiagonal(const xt::xtensor<double, 3>& elemmat);
 
 /**
 Base quadrature-class.
 This class does not have a specific element-type in mind, it is used mostly internally
 to derive from such that common methods do not have to be reimplementation.
 
 \tparam ne Number of nodes per element.
 \tparam nd Number of dimensions for node vectors.
 \tparam td Number of dimensions for integration point tensors.
 */
 template <size_t ne, size_t nd, size_t td>
 class QuadratureBase {
 public:
     QuadratureBase() = default;
 
     /**
     Constructor
     */
     QuadratureBase(size_t nelem, size_t nip);
 
     /**
     Number of elements.
 
     \return Scalar.
     */
     size_t nelem() const;
 
     /**
     Number of nodes per element.
 
     \return Scalar.
     */
     size_t nne() const;
 
     /**
     Number of dimensions for node vectors.
 
     \return Scalar.
     */
     size_t ndim() const;
 
     /**
     Number of dimensions for integration point tensors.
 
     \return Scalar.
     */
     size_t tdim() const;
 
     /**
     Number of integration points.
 
     \return Scalar.
     */
     size_t nip() const;
 
     /**
     Convert "qscalar" to "qtensor" of certain rank.
     Fully allocated output passed as reference, use AsTensor to allocate and return data.
 
     \param qscalar A "qscalar".
     \param qtensor A "qtensor".
     */
     template <size_t rank = 0, class T>
     void asTensor(const xt::xtensor<T, 2>& qscalar, xt::xtensor<T, 2 + rank>& qtensor) const;
 
     /**
     Convert "qscalar" to "qtensor" of certain rank.
 
     \param "qscalar".
     \return "qtensor".
     */
     template <size_t rank = 0, class T>
     xt::xtensor<T, 2 + rank> AsTensor(const xt::xtensor<T, 2>& qscalar) const;
 
     /**
     Convert "qscalar" to "qtensor" of certain rank.
 
     \param rank Tensor rank.
     \param qscalar A "qscalar".
     \return "qtensor".
     */
     template <class T>
     xt::xarray<T> AsTensor(size_t rank, const xt::xtensor<T, 2>& qscalar) const;
 
     /**
     Get the shape of an "elemvec".
 
     \returns [#nelem, #nne, #ndim].
     */
     std::array<size_t, 3> shape_elemvec() const;
 
     /**
     Get the shape of an "elemmat".
 
     \returns [#nelem, #nne * #ndim, #nne * #ndim].
     */
     std::array<size_t, 3> shape_elemmat() const;
 
     /**
     Get the shape of a "qtensor" of a certain rank
     (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).
     Default: rank = 0, a.k.a. scalar.
 
     \returns [#nelem, #nip, #tdim, ...].
     */
     template <size_t rank = 0>
     std::array<size_t, 2 + rank> shape_qtensor() const;
 
     /**
     Get the shape of a "qtensor" of a certain rank
     (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).
 
     \param rank The tensor rank.
     \returns [#nelem, #nip, #tdim, ...].
     */
     std::vector<size_t> shape_qtensor(size_t rank) const;
 
     /**
     Get the shape of a "qscalar" (a "qtensor" of rank 0)
 
     \returns [#nelem, #nip].
     */
     std::vector<size_t> shape_qscalar() const;
 
     /**
     Get an allocated `xt::xtensor` to store a "elemvec".
     Note: the container is not (zero-)initialised.
 
     \returns `xt::xarray` container of the correct shape.
     */
     template <class T>
     xt::xtensor<T, 3> allocate_elemvec() const;
 
     /**
     Get an allocated and initialised `xt::xarray` to store a "elemvec".
 
     \param val The value to which to initialise all items.
     \returns `xt::xtensor` container of the correct shape.
     */
     template <class T>
     xt::xtensor<T, 3> allocate_elemvec(T val) const;
 
     /**
     Get an allocated `xt::xtensor` to store a "elemmat".
     Note: the container is not (zero-)initialised.
 
     \returns `xt::xarray` container of the correct shape.
     */
     template <class T>
     xt::xtensor<T, 3> allocate_elemmat() const;
 
     /**
     Get an allocated and initialised `xt::xarray` to store a "elemmat".
 
     \param val The value to which to initialise all items.
     \returns `xt::xtensor` container of the correct shape.
     */
     template <class T>
     xt::xtensor<T, 3> allocate_elemmat(T val) const;
 
     /**
     Get an allocated `xt::xtensor` to store a "qtensor" of a certain rank
     (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).
     Default: rank = 0, a.k.a. scalar.
     Note: the container is not (zero-)initialised.
 
     \returns [#nelem, #nip].
     */
     template <size_t rank = 0, class T>
     xt::xtensor<T, 2 + rank> allocate_qtensor() const;
 
     /**
     Get an allocated and initialised `xt::xtensor` to store a "qtensor" of a certain rank
     (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).
     Default: rank = 0, a.k.a. scalar.
 
     \param val The value to which to initialise all items.
     \returns `xt::xtensor` container of the correct shape (and rank).
     */
     template <size_t rank = 0, class T>
     xt::xtensor<T, 2 + rank> allocate_qtensor(T val) const;
 
     /**
     Get an allocated `xt::xarray` to store a "qtensor" of a certain rank
     (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).
     Note: the container is not (zero-)initialised.
 
     \param rank The tensor rank.
     \returns `xt::xarray` container of the correct shape.
     */
     template <class T>
     xt::xarray<T> allocate_qtensor(size_t rank) const;
 
     /**
     Get an allocated and initialised `xt::xarray` to store a "qtensor" of a certain rank
     (0 = scalar, 1, vector, 2 = 2nd-order tensor, etc.).
 
     \param rank The tensor rank.
     \param val The value to which to initialise all items.
     \returns `xt::xtensor` container of the correct shape (and rank).
     */
     template <class T>
     xt::xarray<T> allocate_qtensor(size_t rank, T val) const;
 
     /**
     Get an allocated `xt::xtensor` to store a "qscalar" (a "qtensor" of rank 0).
     Note: the container is not (zero-)initialised.
 
     \returns `xt::xarray` container of the correct shape.
     */
     template <class T>
     xt::xtensor<T, 2> allocate_qscalar() const;
 
     /**
     Get an allocated and initialised `xt::xarray` to store a "qscalar" (a "qtensor" of rank 0).
 
     \param val The value to which to initialise all items.
     \returns `xt::xtensor` container of the correct shape (and rank).
     */
     template <class T>
     xt::xtensor<T, 2> allocate_qscalar(T val) const;
 
     /**
     \cond
     */
     template <size_t rank = 0>
     [[ deprecated ]]
     std::array<size_t, 2 + rank> ShapeQtensor() const;
 
     [[ deprecated ]]
     std::vector<size_t> ShapeQtensor(size_t rank) const;
 
     template <size_t rank = 0, class T>
     [[ deprecated ]]
     xt::xtensor<T, 2 + rank> AllocateQtensor() const;
 
     [[ deprecated ]]
     std::vector<size_t> ShapeQscalar() const;
 
     template <size_t rank = 0, class T>
     [[ deprecated ]]
     xt::xtensor<T, 2 + rank> AllocateQtensor(T val) const;
 
     template <class T>
     [[ deprecated ]]
     xt::xarray<T> AllocateQtensor(size_t rank) const;
 
     template <class T>
     [[ deprecated ]]
     xt::xarray<T> AllocateQtensor(size_t rank, T val) const;
 
     template <class T>
     [[ deprecated ]]
     xt::xtensor<T, 2> AllocateQscalar() const;
 
     template <class T>
     [[ deprecated ]]
     xt::xtensor<T, 2> AllocateQscalar(T val) const;
     /**
     \endcond
     */
 
 protected:
     /**
     Wrapper of constructor, for derived classes.
     */
     void initQuadratureBase(size_t nelem, size_t nip);
 
 protected:
     size_t m_nelem; ///< Number of elements.
     size_t m_nip; ///< Number of integration points per element.
     constexpr static size_t m_nne = ne; ///< Number of nodes per element.
     constexpr static size_t m_ndim = nd; ///< Number of dimensions for nodal vectors.
     constexpr static size_t m_tdim = td; ///< Number of dimensions for integration point tensors.
 };
 
 /**
 Interpolation and quadrature for a generic element in Cartesian coordinates.
 
 Naming convention:
 -    ``elemmat``:  matrices stored per element, [#nelem, #nne * #ndim, #nne * #ndim]
 -    ``elemvec``:  nodal vectors stored per element, [#nelem, #nne, #ndim]
 -    ``qtensor``:  integration point tensor, [#nelem, #nip, #tdim, #tdim]
 -    ``qscalar``:  integration point scalar, [#nelem, #nip]
 */
 template <size_t ne, size_t nd, size_t td>
 class QuadratureBaseCartesian : public QuadratureBase<ne, nd, td> {
 public:
 
     QuadratureBaseCartesian() = default;
 
     virtual ~QuadratureBaseCartesian(){};
 
     /**
     Constructor with custom integration.
     The following is pre-computed during construction:
     -   the shape functions,
     -   the shape function gradients (in local and global) coordinates,
     -   the integration points volumes.
     They can be reused without any cost.
     They only have to be recomputed when the nodal position changes
     (note that they are assumed to be constant under a small-strain assumption).
     In that case use update_x() to update the nodal positions and
     to recompute the above listed quantities.
 
     \param x nodal coordinates (``elemvec``).
     \param xi Integration point coordinates (local coordinates) [#nip].
     \param w Integration point weights [#nip].
     \param N Shape functions in local coordinates [#nip, #nne].
     \param dNdxi Shape function gradient w.r.t. local coordinates [#nip, #nne, #ndim].
     */
     QuadratureBaseCartesian(
         const xt::xtensor<double, 3>& x,
         const xt::xtensor<double, 2>& xi,
         const xt::xtensor<double, 1>& w,
         const xt::xtensor<double, 2>& N,
         const xt::xtensor<double, 3>& dNdxi);
 
     /**
     Update the nodal positions.
     This recomputes:
     -   the shape functions,
     -   the shape function gradients (in local and global) coordinates,
     -   the integration points volumes.
     Under the small deformations assumption this function should not be called.
 
     \param x nodal coordinates (``elemvec``). Shape should match the earlier definition.
     */
     void update_x(const xt::xtensor<double, 3>& x);
 
     /**
     Get the shape function gradients (in global coordinates).
 
     \return ``gradN`` stored per element, per integration point [#nelem, #nip, #nne, #ndim].
     */
     virtual xt::xtensor<double, 4> GradN() const;
 
     /**
     Get the integration volume.
 
     \return volume stored per element, per integration point [#nelem, #nip].
     */
     xt::xtensor<double, 2> dV() const;
 
     /**
     \cond
     */
     [[ deprecated ]]
     xt::xtensor<double, 3> Interp_N_vector(const xt::xtensor<double, 3>& elemvec) const;
 
     [[ deprecated ]]
     virtual void interp_N_vector(
         const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 3>& qvector) const;
     /**
     \endcond
     */
 
     /**
     Interpolate element vector and evaluate at quadrature points.
 
     \f$ \vec{u}(\vec{x}_q) = N_i^e(\vec{x}) \vec{u}_i^e
 
     \param elemvec nodal vector stored per element [#nelem, #nne, #ndim].
     \return qvector [#nelem, #nip, #ndim].
     */
     xt::xtensor<double, 3> InterpQuad_vector(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Same as InterpQuad_vector(), but writing to preallocated return.
 
     \param elemvec nodal vector stored per element [#nelem, #nne, #ndim].
     \param qvector [#nelem, #nip, #ndim].
     */
     virtual void interpQuad_vector(
         const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 3>& qvector) const;
 
     /**
     Element-by-element: dyadic product of the shape function gradients and a nodal vector.
     Typical input: nodal displacements. Typical output: quadrature point strains.
     Within one element::
 
         for e in range(nelem):
             for q in range(nip):
                 for m in range(nne):
                     qtensor(e, q, i, j) += dNdx(e, q, m, i) * elemvec(e, m, j)
 
     Note that the functions and their gradients are precomputed upon construction,
     or updated when calling update_x().
 
     \param elemvec [#nelem, #nne, #ndim]
     \return qtensor [#nelem, #nip, #tdim, #tdim]
     */
     xt::xtensor<double, 4> GradN_vector(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Same as GradN_vector(), but writing to preallocated return.
 
     \param elemvec [#nelem, #nne, #ndim]
     \param qtensor overwritten [#nelem, #nip, #tdim, #tdim]
     */
     virtual void gradN_vector(
         const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 4>& qtensor) const;
 
     /**
     The transposed output of GradN_vector().
     Within one element::
 
         for e in range(nelem):
             for q in range(nip):
                 for m in range(nne):
                     qtensor(e, q, j, i) += dNdx(e, q, m, i) * elemvec(e, m, j)
 
     \param elemvec [#nelem, #nne, #ndim]
     \return qtensor [#nelem, #nip, #tdim, #tdim]
     */
     xt::xtensor<double, 4> GradN_vector_T(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Same as GradN_vector_T(), but writing to preallocated return.
 
     \param elemvec [#nelem, #nne, #ndim]
     \param qtensor overwritten [#nelem, #nip, #tdim, #tdim]
     */
     virtual void gradN_vector_T(
         const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 4>& qtensor) const;
 
     /**
     The symmetric output of GradN_vector().
     Without one element::
 
         for e in range(nelem):
             for q in range(nip):
                 for m in range(nne):
                     qtensor(e, q, i, j) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j)
                     qtensor(e, q, j, i) += 0.5 * dNdx(e, q, m, i) * elemvec(e, m, j)
 
     \param elemvec [#nelem, #nne, #ndim]
     \return qtensor [#nelem, #nip, #tdim, #tdim]
     */
     xt::xtensor<double, 4> SymGradN_vector(const xt::xtensor<double, 3>& elemvec) const;
 
     /**
     Same as SymGradN_vector(), but writing to preallocated return.
 
     \param elemvec [#nelem, #nne, #ndim]
     \param qtensor overwritten [#nelem, #nip, #tdim, #tdim]
     */
     virtual void symGradN_vector(
         const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 4>& qtensor) const;
 
+    /**
+    Element-by-element: integral of a continuous vector-field.
+
+    \f$ \vec{f}_i^e = \int N_i^e(\vec{x}) \vec{f}(\vec{x}) d\Omega_e \f$
+
+    integrated by
+
+    \f$ \vec{f}_i^e = \sum\limits_q N_i^e(\vec{x}_q) \vec{f}(\vec{x}_q) \f$
+
+    \param qvector [#nelem, #nip. #ndim]
+    \return elemvec [#nelem, #nne. #ndim]
+    */
+    xt::xtensor<double, 3> Int_N_vector_dV(const xt::xtensor<double, 3>& qvector) const;
+
+    /**
+    Same as Int_N_vector_dV(), but writing to preallocated return.
+
+    \param qvector [#nelem, #nip. #ndim]
+    \return elemvec overwritten [#nelem, #nne. #ndim]
+    */
+    void int_N_vector_dV(
+        const xt::xtensor<double, 3>& qvector, xt::xtensor<double, 3>& elemvec) const;
+
     /**
     Element-by-element: integral of the scalar product of the shape function with a scalar.
     Within one one element::
 
         for e in range(nelem):
             for q in range(nip):
                 for m in range(nne):
                     for n in range(nne):
                         elemmat(e, m * ndim + i, n * ndim + i) +=
                             N(e, q, m) * qscalar(e, q) * N(e, q, n) * dV(e, q)
 
     with ``i`` a tensor dimension.
     Note that the functions and their gradients are precomputed upon construction,
     or updated when calling update_x().
 
     \param qscalar [#nelem, #nip]
     \return elemmat [#nelem, #nne * #ndim, #nne * #ndim]
     */
     xt::xtensor<double, 3> Int_N_scalar_NT_dV(const xt::xtensor<double, 2>& qscalar) const;
 
     /**
     Same as Int_N_scalar_NT_dV(), but writing to preallocated return.
 
     \param qscalar [#nelem, #nip]
     \param elemmat overwritten [#nelem, #nne * #ndim, #nne * #ndim]
     */
     virtual void int_N_scalar_NT_dV(
         const xt::xtensor<double, 2>& qscalar, xt::xtensor<double, 3>& elemmat) const;
 
     /**
     Element-by-element: integral of the dot product of the shape function gradients with
     a second order tensor. Typical input: stress. Typical output: nodal force.
     Within one one element::
 
         for e in range(nelem):
             for q in range(nip):
                 for m in range(nne):
                     elemvec(e, m, j) += dNdx(e, q, m, i) * qtensor(e, q, i, j) * dV(e, q)
 
     with ``i`` and ``j`` tensor dimensions.
     Note that the functions and their gradients are precomputed upon construction,
     or updated when calling update_x().
 
     \param qtensor [#nelem, #nip, #ndim, #ndim]
     \return elemvec [#nelem, #nne. #ndim]
     */
     xt::xtensor<double, 3> Int_gradN_dot_tensor2_dV(const xt::xtensor<double, 4>& qtensor) const;
 
     /**
     Same as Int_gradN_dot_tensor2_dV(), but writing to preallocated return.
 
     \param qtensor [#nelem, #nip, #ndim, #ndim]
     \param elemvec overwritten [#nelem, #nne. #ndim]
     */
     virtual void int_gradN_dot_tensor2_dV(
         const xt::xtensor<double, 4>& qtensor, xt::xtensor<double, 3>& elemvec) const;
 
     // Integral of the dot product
     // elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV
 
     /**
     Element-by-element: integral of the dot products of the shape function gradients with
     a fourth order tensor. Typical input: stiffness tensor. Typical output: stiffness matrix.
     Within one one element::
 
         for e in range(nelem):
             for q in range(nip):
                 for m in range(nne):
                     for n in range(nne):
                         elemmat(e, m * ndim + j, n * ndim + k) +=
                             dNdx(e, q, m, i) * qtensor(e, q, i, j, k, l) * dNdx(e, q, n, l) * dV(e, q)
 
     with ``i``, ``j``, ``k``, and ``l`` tensor dimensions.
     Note that the functions and their gradients are precomputed upon construction,
     or updated when calling update_x().
 
     \param qtensor [#nelem, #nip, #ndim, #ndim, #ndim, #ndim]
     \return elemmat [#nelem, #nne * #ndim, #nne * #ndim]
     */
     xt::xtensor<double, 3> Int_gradN_dot_tensor4_dot_gradNT_dV(
         const xt::xtensor<double, 6>& qtensor) const;
 
     /**
     Same as Int_gradN_dot_tensor4_dot_gradNT_dV(), but writing to preallocated return.
 
     \param qtensor [#nelem, #nip, #ndim, #ndim, #ndim, #ndim]
     \param elemmat overwritten [#nelem, #nne * #ndim, #nne * #ndim]
     */
     virtual void int_gradN_dot_tensor4_dot_gradNT_dV(
         const xt::xtensor<double, 6>& qtensor, xt::xtensor<double, 3>& elemmat) const;
 
 protected:
     /**
     Constructor alias.
     \param x nodal coordinates (``elemvec``).
     \param xi Integration point coordinates (local coordinates) [#nip].
     \param w Integration point weights [#nip].
     \param N Shape functions in the integration points [#nip, #nne].
     \param dNdxi Shape function gradient w.r.t. local coordinates [#nip, #nne, #ndim].
     */
     void initQuadratureBaseCartesian(
         const xt::xtensor<double, 3>& x,
         const xt::xtensor<double, 2>& xi,
         const xt::xtensor<double, 1>& w,
         const xt::xtensor<double, 2>& N,
         const xt::xtensor<double, 3>& dNdxi);
 
     /**
     Compute m_vol and m_dNdx based on current m_x.
     */
     virtual void compute_dN();
 
 protected:
     using QuadratureBase<ne, nd, td>::m_nelem;
     using QuadratureBase<ne, nd, td>::m_nip;
     using QuadratureBase<ne, nd, td>::m_nne;
     using QuadratureBase<ne, nd, td>::m_ndim;
     using QuadratureBase<ne, nd, td>::m_tdim;
 
     xt::xtensor<double, 3> m_x;    ///< nodal positions stored per element [#nelem, #nne, #ndim]
     xt::xtensor<double, 1> m_w;    ///< weight of each integration point [nip]
     xt::xtensor<double, 2> m_xi;   ///< local coordinate of each integration point [#nip, #ndim]
     xt::xtensor<double, 2> m_N;    ///< shape functions [#nip, #nne]
     xt::xtensor<double, 3> m_dNxi; ///< shape function grad. wrt local  coor. [#nip, #nne, #ndim]
     xt::xtensor<double, 4> m_dNx;  ///< shape function grad. wrt global coor. [#nelem, #nip, #nne, #ndim]
     xt::xtensor<double, 2> m_vol;  ///< integration point volume [#nelem, #nip]
 };
 
 } // namespace Element
 } // namespace GooseFEM
 
 #include "Element.hpp"
 
 #endif
diff --git a/include/GooseFEM/Element.hpp b/include/GooseFEM/Element.hpp
index 1cd4a40..2a3abf5 100644
--- a/include/GooseFEM/Element.hpp
+++ b/include/GooseFEM/Element.hpp
@@ -1,794 +1,834 @@
 /**
 Implementation of Element.h
 
 \file Element.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_ELEMENT_HPP
 #define GOOSEFEM_ELEMENT_HPP
 
 #include "Element.h"
 #include "detail.hpp"
 
 namespace GooseFEM {
 namespace Element {
 
 inline xt::xtensor<double, 3> asElementVector(
     const xt::xtensor<size_t, 2>& conn, const xt::xtensor<double, 2>& nodevec)
 {
     size_t nelem = conn.shape(0);
     size_t nne = conn.shape(1);
     size_t ndim = nodevec.shape(1);
 
     xt::xtensor<double, 3> elemvec = xt::empty<double>({nelem, nne, ndim});
 
     #pragma omp parallel for
     for (size_t e = 0; e < nelem; ++e) {
         for (size_t m = 0; m < nne; ++m) {
             for (size_t i = 0; i < ndim; ++i) {
                 elemvec(e, m, i) = nodevec(conn(e, m), i);
             }
         }
     }
 
     return elemvec;
 }
 
 inline xt::xtensor<double, 2> assembleNodeVector(
     const xt::xtensor<size_t, 2>& conn, const xt::xtensor<double, 3>& elemvec)
 {
     size_t nelem = conn.shape(0);
     size_t nne = conn.shape(1);
     size_t ndim = elemvec.shape(2);
     size_t nnode = xt::amax(conn)() + 1;
 
     GOOSEFEM_ASSERT(elemvec.shape(0) == nelem);
     GOOSEFEM_ASSERT(elemvec.shape(1) == nne);
 
     xt::xtensor<double, 2> nodevec = xt::zeros<double>({nnode, ndim});
 
     for (size_t e = 0; e < nelem; ++e) {
         for (size_t m = 0; m < nne; ++m) {
             for (size_t i = 0; i < ndim; ++i) {
                 nodevec(conn(e, m), i) += elemvec(e, m, i);
             }
         }
     }
 
     return nodevec;
 }
 
 template <class E>
 inline bool isSequential(const E& dofs)
 {
     size_t ndof = xt::amax(dofs)() + 1;
 
     xt::xtensor<int, 1> exists = xt::zeros<int>({ndof});
 
     for (auto& i : dofs) {
         exists[i]++;
     }
 
     for (auto& i : dofs) {
         if (exists[i] == 0) {
             return false;
         }
     }
 
     return true;
 }
 
 inline bool isDiagonal(const xt::xtensor<double, 3>& elemmat)
 {
     GOOSEFEM_ASSERT(elemmat.shape(1) == elemmat.shape(2));
 
     size_t nelem = elemmat.shape(0);
     size_t N = elemmat.shape(1);
 
     double eps = std::numeric_limits<double>::epsilon();
 
     #pragma omp parallel for
     for (size_t e = 0; e < nelem; ++e) {
         for (size_t i = 0; i < N; ++i) {
             for (size_t j = 0; j < N; ++j) {
                 if (i != j) {
                     if (std::abs(elemmat(e, i, j)) > eps) {
                         return false;
                     }
                 }
             }
         }
     }
 
     return true;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline QuadratureBase<ne, nd, td>::QuadratureBase(size_t nelem, size_t nip)
 {
     this->initQuadratureBase(nelem, nip);
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBase<ne, nd, td>::initQuadratureBase(size_t nelem, size_t nip)
 {
     m_nelem = nelem;
     m_nip = nip;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline size_t QuadratureBase<ne, nd, td>::nelem() const
 {
     return m_nelem;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline size_t QuadratureBase<ne, nd, td>::nne() const
 {
     return m_nne;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline size_t QuadratureBase<ne, nd, td>::ndim() const
 {
     return m_ndim;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline size_t QuadratureBase<ne, nd, td>::tdim() const
 {
     return m_tdim;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline size_t QuadratureBase<ne, nd, td>::nip() const
 {
     return m_nip;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank, class T>
 inline void QuadratureBase<ne, nd, td>::asTensor(
     const xt::xtensor<T, 2>& arg,
     xt::xtensor<T, 2 + rank>& ret) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(arg, {m_nelem, m_nne}));
     GooseFEM::asTensor<2, rank>(arg, ret);
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank, class T>
 inline xt::xtensor<T, 2 + rank> QuadratureBase<ne, nd, td>::AsTensor(
     const xt::xtensor<T, 2>& qscalar) const
 {
     return GooseFEM::AsTensor<2, rank>(qscalar, m_tdim);
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xarray<T> QuadratureBase<ne, nd, td>::AsTensor(
     size_t rank,
     const xt::xtensor<T, 2>& qscalar) const
 {
     return GooseFEM::AsTensor(rank, qscalar, m_tdim);
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline std::array<size_t, 3> QuadratureBase<ne, nd, td>::shape_elemvec() const
 {
     return std::array<size_t, 3>{m_nelem, m_nne, m_ndim};
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline std::array<size_t, 3> QuadratureBase<ne, nd, td>::shape_elemmat() const
 {
     return std::array<size_t, 3>{m_nelem, m_nne * m_ndim, m_nne * m_ndim};
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank>
 inline std::array<size_t, 2 + rank> QuadratureBase<ne, nd, td>::shape_qtensor() const
 {
     std::array<size_t, 2 + rank> shape;
     shape[0] = m_nelem;
     shape[1] = m_nip;
     std::fill(shape.begin() + 2, shape.end(), td);
     return shape;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline std::vector<size_t> QuadratureBase<ne, nd, td>::shape_qtensor(size_t rank) const
 {
     std::vector<size_t> shape(2 + rank);
     shape[0] = m_nelem;
     shape[1] = m_nip;
     std::fill(shape.begin() + 2, shape.end(), td);
     return shape;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline std::vector<size_t> QuadratureBase<ne, nd, td>::shape_qscalar() const
 {
     std::vector<size_t> shape(2);
     shape[0] = m_nelem;
     shape[1] = m_nip;
     return shape;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 3> QuadratureBase<ne, nd, td>::allocate_elemvec() const
 {
     xt::xtensor<T, 3> ret = xt::empty<T>(this->shape_elemvec());
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 3> QuadratureBase<ne, nd, td>::allocate_elemvec(T val) const
 {
     xt::xtensor<T, 3> ret = xt::empty<T>(this->shape_elemvec());
     ret.fill(val);
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 3> QuadratureBase<ne, nd, td>::allocate_elemmat() const
 {
     xt::xtensor<T, 3> ret = xt::empty<T>(this->shape_elemmat());
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 3> QuadratureBase<ne, nd, td>::allocate_elemmat(T val) const
 {
     xt::xtensor<T, 3> ret = xt::empty<T>(this->shape_elemmat());
     ret.fill(val);
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank, class T>
 inline xt::xtensor<T, 2 + rank> QuadratureBase<ne, nd, td>::allocate_qtensor() const
 {
     xt::xtensor<T, 2 + rank> ret = xt::empty<T>(this->shape_qtensor<rank>());
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank, class T>
 inline xt::xtensor<T, 2 + rank> QuadratureBase<ne, nd, td>::allocate_qtensor(T val) const
 {
     xt::xtensor<T, 2 + rank> ret = xt::empty<T>(this->shape_qtensor<rank>());
     ret.fill(val);
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xarray<T> QuadratureBase<ne, nd, td>::allocate_qtensor(size_t rank) const
 {
     xt::xarray<T> ret = xt::empty<T>(this->shape_qtensor(rank));
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xarray<T> QuadratureBase<ne, nd, td>::allocate_qtensor(size_t rank, T val) const
 {
     xt::xarray<T> ret = xt::empty<T>(this->shape_qtensor(rank));
     ret.fill(val);
     return ret;
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 2> QuadratureBase<ne, nd, td>::allocate_qscalar() const
 {
     return this->allocate_qtensor<0, T>();
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 2> QuadratureBase<ne, nd, td>::allocate_qscalar(T val) const
 {
     return this->allocate_qtensor<0, T>(val);
 }
 
 /**
 \cond
 */
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank>
 inline std::array<size_t, 2 + rank> QuadratureBase<ne, nd, td>::ShapeQtensor() const
 {
     GOOSEFEM_WARNING("Deprecation warning: ShapeQtensor<rank> -> shape_qtensor<rank>");
     return this->shape_qtensor<rank>();
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline std::vector<size_t> QuadratureBase<ne, nd, td>::ShapeQtensor(size_t rank) const
 {
     GOOSEFEM_WARNING("Deprecation warning: ShapeQtensor(rank) -> shape_qtensor(rank)");
     return this->shape_qtensor(rank);
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline std::vector<size_t> QuadratureBase<ne, nd, td>::ShapeQscalar() const
 {
     GOOSEFEM_WARNING("Deprecation warning: ShapeQscalar -> shape_qscalar");
     return this->shape_qscalar();
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank, class T>
 inline xt::xtensor<T, 2 + rank> QuadratureBase<ne, nd, td>::AllocateQtensor() const
 {
     GOOSEFEM_WARNING("Deprecation warning: AllocateQtensor<rank, T> -> allocate_qtensor<rank, T>");
     return this->allocate_qtensor<rank, T>();
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <size_t rank, class T>
 inline xt::xtensor<T, 2 + rank> QuadratureBase<ne, nd, td>::AllocateQtensor(T val) const
 {
     GOOSEFEM_WARNING("Deprecation warning: AllocateQtensor<rank, T> -> allocate_qtensor<rank, T>");
     return this->allocate_qtensor<rank, T>(val);
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xarray<T> QuadratureBase<ne, nd, td>::AllocateQtensor(size_t rank) const
 {
     GOOSEFEM_WARNING_PYTHON("Deprecation warning: use np.empty(this.shape_qtensor(rank))")
     GOOSEFEM_WARNING("Deprecation warning: AllocateQtensor(rank) -> allocate_qtensor(rank)");
     return this->allocate_qtensor<T>(rank);
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xarray<T> QuadratureBase<ne, nd, td>::AllocateQtensor(size_t rank, T val) const
 {
     GOOSEFEM_WARNING_PYTHON("Deprecation warning: use val * np.ones(this.shape_qtensor(rank))")
     GOOSEFEM_WARNING("Deprecation warning: AllocateQtensor(rank) -> allocate_qtensor(rank)");
     return this->allocate_qtensor<T>(rank, val);
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 2> QuadratureBase<ne, nd, td>::AllocateQscalar() const
 {
     GOOSEFEM_WARNING_PYTHON("Deprecation warning: use np.empty(this.shape_qscalar())")
     GOOSEFEM_WARNING("Deprecation warning: AllocateQscalar -> allocate_qscalar");
     return this->allocate_qtensor<0, T>();
 }
 
 template <size_t ne, size_t nd, size_t td>
 template <class T>
 inline xt::xtensor<T, 2> QuadratureBase<ne, nd, td>::AllocateQscalar(T val) const
 {
     GOOSEFEM_WARNING_PYTHON("Deprecation warning: use np.empty(this.shape_qscalar())")
     GOOSEFEM_WARNING("Deprecation warning: AllocateQscalar -> allocate_qscalar");
     return this->allocate_qtensor<0, T>(val);
 }
 
 /**
 \endcond
 */
 
 template <size_t ne, size_t nd, size_t td>
 inline QuadratureBaseCartesian<ne, nd, td>::QuadratureBaseCartesian(
     const xt::xtensor<double, 3>& x,
     const xt::xtensor<double, 2>& xi,
     const xt::xtensor<double, 1>& w,
     const xt::xtensor<double, 2>& N,
     const xt::xtensor<double, 3>& dNdxi)
 {
     this->initQuadratureBaseCartesian(x, xi, w, N, dNdxi);
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::initQuadratureBaseCartesian(
     const xt::xtensor<double, 3>& x,
     const xt::xtensor<double, 2>& xi,
     const xt::xtensor<double, 1>& w,
     const xt::xtensor<double, 2>& N,
     const xt::xtensor<double, 3>& dNdxi)
 {
     m_x = x;
     m_w = w;
     m_xi = xi;
     m_N = N;
     m_dNxi = dNdxi;
 
     this->initQuadratureBase(m_x.shape(0), m_w.size());
 
     GOOSEFEM_ASSERT(m_x.shape(1) == m_nne);
     GOOSEFEM_ASSERT(m_x.shape(2) == m_ndim);
     GOOSEFEM_ASSERT(xt::has_shape(m_xi, {m_nip, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(m_w, {m_nip}));
     GOOSEFEM_ASSERT(xt::has_shape(m_N, {m_nip, m_nne}));
     GOOSEFEM_ASSERT(xt::has_shape(m_dNxi, {m_nip, m_nne, m_ndim}));
 
     m_dNx = xt::empty<double>({m_nelem, m_nip, m_nne, m_ndim});
     m_vol = xt::empty<double>(this->shape_qscalar());
 
     this->compute_dN();
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::compute_dN()
 {
     #pragma omp parallel
     {
         xt::xtensor<double, 2> J = xt::empty<double>({m_ndim, m_ndim});
         xt::xtensor<double, 2> Jinv = xt::empty<double>({m_ndim, m_ndim});
         m_dNx.fill(0.0);
 
         #pragma omp for
         for (size_t e = 0; e < m_nelem; ++e) {
 
             auto x = xt::adapt(&m_x(e, 0, 0), xt::xshape<m_nne, m_ndim>());
 
             for (size_t q = 0; q < m_nip; ++q) {
 
                 auto dNxi = xt::adapt(&m_dNxi(q, 0, 0), xt::xshape<m_nne, m_ndim>());
                 auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<m_nne, m_ndim>());
 
                 J.fill(0.0);
 
                 for (size_t m = 0; m < m_nne; ++m) {
                     for (size_t i = 0; i < m_ndim; ++i) {
                         for (size_t j = 0; j < m_ndim; ++j) {
                             J(i, j) += dNxi(m, i) * x(m, j);
                         }
                     }
                 }
 
                 double Jdet = detail::tensor<m_ndim>::inv(J, Jinv);
 
                 for (size_t m = 0; m < m_nne; ++m) {
                     for (size_t i = 0; i < m_ndim; ++i) {
                         for (size_t j = 0; j < m_ndim; ++j) {
                             dNx(m, i) += Jinv(i, j) * dNxi(m, i);
                         }
                     }
                 }
 
                 m_vol(e, q) = m_w(q) * Jdet;
             }
         }
     }
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 4> QuadratureBaseCartesian<ne, nd, td>::GradN() const
 {
     return m_dNx;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 2> QuadratureBaseCartesian<ne, nd, td>::dV() const
 {
     return m_vol;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::update_x(const xt::xtensor<double, 3>& x)
 {
     GOOSEFEM_ASSERT(x.shape() == m_x.shape());
     xt::noalias(m_x) = x;
     this->compute_dN();
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 3> QuadratureBaseCartesian<ne, nd, td>::InterpQuad_vector(
     const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 3> qvector = xt::empty<double>({m_nelem, m_nip, m_ndim});
     this->interpQuad_vector(elemvec, qvector);
     return qvector;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::interpQuad_vector(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 3>& qvector) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, this->shape_elemvec()));
     GOOSEFEM_ASSERT(xt::has_shape(qvector, {m_nelem, m_nip, m_ndim}));
 
     qvector.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<m_nne, m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto N = xt::adapt(&m_N(q, 0), xt::xshape<m_nne>());
             auto ui = xt::adapt(&qvector(e, q, 0), xt::xshape<m_ndim>());
 
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t i = 0; i < m_ndim; ++i) {
                     ui(i) += N(m) * u(m, i);
                 }
             }
         }
     }
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 3> QuadratureBaseCartesian<ne, nd, td>::Interp_N_vector(
     const xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_WARNING("Deprecation warning: Interp_N_vector -> InterpQuad_vector");
     GOOSEFEM_WARNING_PYTHON("Deprecation warning: Interp_N_vector -> InterpQuad_vector")
     return this->InterpQuad_vector(elemvec);
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::interp_N_vector(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 3>& qvector) const
 {
     GOOSEFEM_WARNING("Deprecation warning: interp_N_vector -> interpQuad_vector");
     GOOSEFEM_WARNING_PYTHON("Deprecation warning: interp_N_vector -> interpQuad_vector")
     this->interpQuad_vector(elemvec, qvector);
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 4> QuadratureBaseCartesian<ne, nd, td>::GradN_vector(
     const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 4> qtensor = xt::empty<double>({m_nelem, m_nip, m_tdim, m_tdim});
     this->gradN_vector(elemvec, qtensor);
     return qtensor;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::gradN_vector(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 4>& qtensor) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim}));
 
     qtensor.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<m_nne, m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<m_nne, m_ndim>());
             auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<m_tdim, m_tdim>());
 
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t i = 0; i < m_ndim; ++i) {
                     for (size_t j = 0; j < m_ndim; ++j) {
                         gradu(i, j) += dNx(m, i) * u(m, j);
                     }
                 }
             }
         }
     }
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 4> QuadratureBaseCartesian<ne, nd, td>::GradN_vector_T(
     const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 4> qtensor = xt::empty<double>({m_nelem, m_nip, m_tdim, m_tdim});
     this->gradN_vector_T(elemvec, qtensor);
     return qtensor;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::gradN_vector_T(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 4>& qtensor) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim}));
 
     qtensor.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<m_nne, m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<m_nne, m_ndim>());
             auto gradu = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<m_tdim, m_tdim>());
 
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t i = 0; i < m_ndim; ++i) {
                     for (size_t j = 0; j < m_ndim; ++j) {
                         gradu(j, i) += dNx(m, i) * u(m, j);
                     }
                 }
             }
         }
     }
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 4> QuadratureBaseCartesian<ne, nd, td>::SymGradN_vector(
     const xt::xtensor<double, 3>& elemvec) const
 {
     xt::xtensor<double, 4> qtensor = xt::empty<double>({m_nelem, m_nip, m_tdim, m_tdim});
     this->symGradN_vector(elemvec, qtensor);
     return qtensor;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::symGradN_vector(
     const xt::xtensor<double, 3>& elemvec, xt::xtensor<double, 4>& qtensor) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
     GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim}));
 
     qtensor.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto u = xt::adapt(&elemvec(e, 0, 0), xt::xshape<m_nne, m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<m_nne, m_ndim>());
             auto eps = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<m_tdim, m_tdim>());
 
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t i = 0; i < m_ndim; ++i) {
                     for (size_t j = 0; j < m_ndim; ++j) {
                         eps(i, j) += 0.5 * dNx(m, i) * u(m, j);
                         eps(j, i) += 0.5 * dNx(m, i) * u(m, j);
                     }
                 }
             }
         }
     }
 }
 
+template <size_t ne, size_t nd, size_t td>
+inline xt::xtensor<double, 3> QuadratureBaseCartesian<ne, nd, td>::Int_N_vector_dV(
+    const xt::xtensor<double, 3>& qvector) const
+{
+    size_t n = qvector.shape(2);
+    xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, n});
+    this->int_N_vector_dV(qvector, elemvec);
+    return elemvec;
+}
+
+template <size_t ne, size_t nd, size_t td>
+inline void QuadratureBaseCartesian<ne, nd, td>::int_N_vector_dV(
+    const xt::xtensor<double, 3>& qvector, xt::xtensor<double, 3>& elemvec) const
+{
+    size_t n = qvector.shape(2);
+    GOOSEFEM_ASSERT(xt::has_shape(qvector, {m_nelem, m_nip, n}));
+    GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, n}));
+
+    elemvec.fill(0.0);
+
+    #pragma omp parallel for
+    for (size_t e = 0; e < m_nelem; ++e) {
+
+        auto f = &elemvec(e, 0, 0);
+
+        for (size_t q = 0; q < m_nip; ++q) {
+
+            auto N = &m_N(q, 0);
+            auto t = &qvector(e, q, 0);
+            auto& vol = m_vol(e, q);
+
+            for (size_t m = 0; m < m_nne; ++m) {
+                for (size_t i = 0; i < n; ++i) {
+                    f[m * n + i] += N[m] * t[i] * vol;
+                }
+            }
+        }
+    }
+}
+
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 3> QuadratureBaseCartesian<ne, nd, td>::Int_N_scalar_NT_dV(
     const xt::xtensor<double, 2>& qscalar) const
 {
     xt::xtensor<double, 3> elemmat = xt::empty<double>({m_nelem, m_nne * m_ndim, m_nne * m_ndim});
     this->int_N_scalar_NT_dV(qscalar, elemmat);
     return elemmat;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::int_N_scalar_NT_dV(
     const xt::xtensor<double, 2>& qscalar, xt::xtensor<double, 3>& elemmat) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(qscalar, {m_nelem, m_nip}));
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
 
     elemmat.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto M = xt::adapt(&elemmat(e, 0, 0), xt::xshape<m_nne * m_ndim, m_nne * m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto N = xt::adapt(&m_N(q, 0), xt::xshape<m_nne>());
             auto& vol = m_vol(e, q);
             auto& rho = qscalar(e, q);
 
             // M(m * ndim + i, n * ndim + i) += N(m) * scalar * N(n) * dV
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t n = 0; n < m_nne; ++n) {
                     for (size_t i = 0; i < m_ndim; ++i) {
                         M(m * m_ndim + i, n * m_ndim + i) += N(m) * rho * N(n) * vol;
                     }
                 }
             }
         }
     }
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 3> QuadratureBaseCartesian<ne, nd, td>::Int_gradN_dot_tensor2_dV(
     const xt::xtensor<double, 4>& qtensor) const
 {
     xt::xtensor<double, 3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
     this->int_gradN_dot_tensor2_dV(qtensor, elemvec);
     return elemvec;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::int_gradN_dot_tensor2_dV(
     const xt::xtensor<double, 4>& qtensor, xt::xtensor<double, 3>& elemvec) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim}));
     GOOSEFEM_ASSERT(xt::has_shape(elemvec, {m_nelem, m_nne, m_ndim}));
 
     elemvec.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto f = xt::adapt(&elemvec(e, 0, 0), xt::xshape<m_nne, m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<m_nne, m_ndim>());
             auto sig = xt::adapt(&qtensor(e, q, 0, 0), xt::xshape<m_tdim, m_tdim>());
             auto& vol = m_vol(e, q);
 
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t i = 0; i < m_ndim; ++i) {
                     for (size_t j = 0; j < m_ndim; ++j) {
                         f(m, j) += dNx(m, i) * sig(i, j) * vol;
                     }
                 }
             }
         }
     }
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline xt::xtensor<double, 3> QuadratureBaseCartesian<ne, nd, td>::Int_gradN_dot_tensor4_dot_gradNT_dV(
     const xt::xtensor<double, 6>& qtensor) const
 {
     xt::xtensor<double, 3> elemmat = xt::empty<double>({m_nelem, m_ndim * m_nne, m_ndim * m_nne});
     this->int_gradN_dot_tensor4_dot_gradNT_dV(qtensor, elemmat);
     return elemmat;
 }
 
 template <size_t ne, size_t nd, size_t td>
 inline void QuadratureBaseCartesian<ne, nd, td>::int_gradN_dot_tensor4_dot_gradNT_dV(
     const xt::xtensor<double, 6>& qtensor, xt::xtensor<double, 3>& elemmat) const
 {
     GOOSEFEM_ASSERT(xt::has_shape(qtensor, {m_nelem, m_nip, m_tdim, m_tdim, m_tdim, m_tdim}));
     GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
 
     elemmat.fill(0.0);
 
     #pragma omp parallel for
     for (size_t e = 0; e < m_nelem; ++e) {
 
         auto K = xt::adapt(&elemmat(e, 0, 0), xt::xshape<m_nne * m_ndim, m_nne * m_ndim>());
 
         for (size_t q = 0; q < m_nip; ++q) {
 
             auto dNx = xt::adapt(&m_dNx(e, q, 0, 0), xt::xshape<m_nne, m_ndim>());
             auto C = xt::adapt(&qtensor(e, q, 0, 0, 0, 0), xt::xshape<m_tdim, m_tdim, m_tdim, m_tdim>());
             auto& vol = m_vol(e, q);
 
             for (size_t m = 0; m < m_nne; ++m) {
                 for (size_t n = 0; n < m_nne; ++n) {
                     for (size_t i = 0; i < m_ndim; ++i) {
                         for (size_t j = 0; j < m_ndim; ++j) {
                             for (size_t k = 0; k < m_ndim; ++k) {
                                 for (size_t l = 0; l < m_ndim; ++l) {
                                     K(m * m_ndim + j, n * m_ndim + k) +=
                                         dNx(m, i) * C(i, j, k, l) * dNx(n, l) * vol;
                                 }
                             }
                         }
                     }
                 }
             }
         }
     }
 }
 
 
 } // namespace Element
 } // namespace GooseFEM
 
 #endif
diff --git a/python/ElementQuad4.hpp b/python/ElementQuad4.hpp
index 1ae023e..11d9ced 100644
--- a/python/ElementQuad4.hpp
+++ b/python/ElementQuad4.hpp
@@ -1,185 +1,190 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
 
 namespace py = pybind11;
 
 void init_ElementQuad4(py::module& m)
 {
 
     py::class_<GooseFEM::Element::Quad4::Quadrature>(m, "Quadrature")
 
         .def(py::init<const xt::xtensor<double, 3>&>(), "Quadrature", py::arg("x"))
 
         .def(py::init<
                 const xt::xtensor<double, 3>&,
                 const xt::xtensor<double, 2>&,
                 const xt::xtensor<double, 1>&>(),
             "Quadrature",
             py::arg("x"),
             py::arg("xi"),
             py::arg("w"))
 
         .def("update_x",
             &GooseFEM::Element::Quad4::Quadrature::update_x,
             "Update the nodal positions")
 
         .def("dV", &GooseFEM::Element::Quad4::Quadrature::dV, "Integration point volume (qscalar)")
 
         .def("Interp_N_vector",
             &GooseFEM::Element::Quad4::Quadrature::Interp_N_vector,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::Interp_N_vector`.",
             py::arg("elemvec"))
 
         .def("InterpQuad_vector",
             &GooseFEM::Element::Quad4::Quadrature::InterpQuad_vector,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::Interp_N_vector`.",
             py::arg("elemvec"))
 
         .def("GradN_vector",
             py::overload_cast<const xt::xtensor<double, 3>&>(
                 &GooseFEM::Element::Quad4::Quadrature::GradN_vector, py::const_),
             "Dyadic product, returns 'qtensor'",
             py::arg("elemvec"))
 
         .def("GradN_vector_T",
             py::overload_cast<const xt::xtensor<double, 3>&>(
                 &GooseFEM::Element::Quad4::Quadrature::GradN_vector_T, py::const_),
             "Dyadic product, returns 'qtensor'",
             py::arg("elemvec"))
 
         .def("SymGradN_vector",
             py::overload_cast<const xt::xtensor<double, 3>&>(
                 &GooseFEM::Element::Quad4::Quadrature::SymGradN_vector, py::const_),
             "Dyadic product, returns 'qtensor'",
             py::arg("elemvec"))
 
+        .def("Int_N_vector_dV",
+            &GooseFEM::Element::Quad4::Quadrature::Int_N_vector_dV,
+            "See :cpp:class:`GooseFEM::Quad4::Quadrature::Int_N_vector_dV`.",
+            py::arg("qvector"))
+
         .def("Int_N_scalar_NT_dV",
             py::overload_cast<const xt::xtensor<double, 2>&>(
                 &GooseFEM::Element::Quad4::Quadrature::Int_N_scalar_NT_dV, py::const_),
             "Integration, returns 'elemmat'",
             py::arg("qscalar"))
 
         .def("Int_gradN_dot_tensor2_dV",
             py::overload_cast<const xt::xtensor<double, 4>&>(
                 &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor2_dV, py::const_),
             "Integration, returns 'elemvec'",
             py::arg("qtensor"))
 
         .def("Int_gradN_dot_tensor4_dot_gradNT_dV",
             py::overload_cast<const xt::xtensor<double, 6>&>(
                 &GooseFEM::Element::Quad4::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV,
                 py::const_),
             "Integration, returns 'elemvec'",
             py::arg("qtensor"))
 
         // Derived from QuadratureBase
 
         .def("nelem", &GooseFEM::Element::Quad4::Quadrature::nelem, "Number of elements")
 
         .def("nne", &GooseFEM::Element::Quad4::Quadrature::nne, "Number of nodes per element")
 
         .def("ndim", &GooseFEM::Element::Quad4::Quadrature::ndim, "Number of dimensions")
 
         .def("nip", &GooseFEM::Element::Quad4::Quadrature::nip, "Number of integration points")
 
         .def("AsTensor",
             (xt::xarray<double>(GooseFEM::Element::Quad4::Quadrature::*)(
                 size_t, const xt::xtensor<double, 2>&) const)
                 &GooseFEM::Element::Quad4::Quadrature::AsTensor<double>,
             "Convert 'qscalar' to 'qtensor' of certain rank")
 
         .def("shape_elemvec",
             &GooseFEM::Element::Quad4::Quadrature::shape_elemvec,
             "Shape of 'elemvec'")
 
         .def("shape_elemmat",
             &GooseFEM::Element::Quad4::Quadrature::shape_elemmat,
             "Shape of 'elemmat'")
 
         .def("shape_qtensor",
             (std::vector<size_t>(GooseFEM::Element::Quad4::Quadrature::*)(size_t) const)
                 &GooseFEM::Element::Quad4::Quadrature::shape_qtensor,
             "Shape of 'qtensor'",
             py::arg("rank"))
 
         .def("shape_qscalar",
             &GooseFEM::Element::Quad4::Quadrature::shape_qscalar,
             "Shape of 'qscalar'")
 
         // Deprecated
 
         .def("ShapeQtensor",
             (std::vector<size_t>(GooseFEM::Element::Quad4::Quadrature::*)(size_t) const)
                 &GooseFEM::Element::Quad4::Quadrature::shape_qtensor,
             "Shape of 'qtensor'",
             py::arg("rank"))
 
         .def("ShapeQscalar",
             &GooseFEM::Element::Quad4::Quadrature::shape_qscalar,
             "Shape of 'qscalar'")
 
         .def("AllocateQtensor",
             (xt::xarray<double>(GooseFEM::Element::Quad4::Quadrature::*)(size_t) const)
                 &GooseFEM::Element::Quad4::Quadrature::allocate_qtensor<double>,
             "Allocate 'qtensor'",
             py::arg("rank"))
 
         .def("AllocateQtensor",
             (xt::xarray<double>(GooseFEM::Element::Quad4::Quadrature::*)(size_t, double) const)
                 &GooseFEM::Element::Quad4::Quadrature::allocate_qtensor<double>,
             "Allocate 'qtensor'",
             py::arg("rank"),
             py::arg("val"))
 
         .def("AllocateQscalar",
             py::overload_cast<>(
                 &GooseFEM::Element::Quad4::Quadrature::allocate_qscalar<double>, py::const_),
             "Allocate 'qscalar'")
 
         .def("AllocateQscalar",
             py::overload_cast<double>(
                 &GooseFEM::Element::Quad4::Quadrature::allocate_qscalar<double>, py::const_),
             "Allocate 'qscalar'",
             py::arg("val"))
 
         .def("__repr__", [](const GooseFEM::Element::Quad4::Quadrature&) {
             return "<GooseFEM.Element.Quad4.Quadrature>";
         });
 }
 
 void init_ElementQuad4Gauss(py::module& m)
 {
 
     m.def("nip", &GooseFEM::Element::Quad4::Gauss::nip, "Return number of integration point");
 
     m.def("xi", &GooseFEM::Element::Quad4::Gauss::xi, "Return integration point coordinates");
 
     m.def("w", &GooseFEM::Element::Quad4::Gauss::w, "Return integration point weights");
 }
 
 void init_ElementQuad4Nodal(py::module& m)
 {
 
     m.def("nip", &GooseFEM::Element::Quad4::Nodal::nip, "Return number of integration point");
 
     m.def("xi", &GooseFEM::Element::Quad4::Nodal::xi, "Return integration point coordinates");
 
     m.def("w", &GooseFEM::Element::Quad4::Nodal::w, "Return integration point weights");
 }
 
 void init_ElementQuad4MidPoint(py::module& m)
 {
 
     m.def("nip", &GooseFEM::Element::Quad4::MidPoint::nip, "Return number of integration point");
 
     m.def("xi", &GooseFEM::Element::Quad4::MidPoint::xi, "Return integration point coordinates");
 
     m.def("w", &GooseFEM::Element::Quad4::MidPoint::w, "Return integration point weights");
 }
diff --git a/test/basic/ElementQuad4.cpp b/test/basic/ElementQuad4.cpp
index a99b790..82005d8 100644
--- a/test/basic/ElementQuad4.cpp
+++ b/test/basic/ElementQuad4.cpp
@@ -1,203 +1,220 @@
 
 #include <catch2/catch.hpp>
 #include <xtensor/xrandom.hpp>
 #include <xtensor/xmath.hpp>
 #include <GooseFEM/GooseFEM.h>
 
 TEST_CASE("GooseFEM::ElementQuad4", "ElementQuad4.h")
 {
     SECTION("GradN")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
         GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor()));
         auto dNdx = quad.GradN();
         REQUIRE(xt::has_shape(dNdx, {mesh.nelem(), quad.nip(), mesh.nne(), mesh.ndim()}));
     }
 
     SECTION("dV")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
         GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor()));
         auto dV = quad.dV();
         REQUIRE(xt::has_shape(dV, quad.shape_qtensor<0>()));
         REQUIRE(xt::has_shape(dV, quad.shape_qscalar()));
         REQUIRE(xt::allclose(dV, 0.5 * 0.5));
     }
 
     SECTION("InterpQuad_vector")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
         GooseFEM::Element::Quad4::Quadrature quad(vector.AsElement(mesh.coor()));
 
         auto u = vector.allocate_nodevec(1.0);
         auto ue = vector.AsElement(u);
         auto uq = quad.InterpQuad_vector(ue);
 
         REQUIRE(xt::allclose(uq, 1.0));
     }
 
     SECTION("GradN_vector")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
         GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor()));
 
         size_t td = mesh.ndim();
         xt::xtensor<double, 2> F = xt::zeros<double>({td, td});
         xt::xtensor<double, 2> EPS = xt::zeros<double>({td, td});
 
         F(0, 1) = 0.1;
 
         auto coor = mesh.coor();
         auto disp = xt::zeros_like(coor);
 
         for (size_t n = 0; n < mesh.nnode(); ++n) {
             for (size_t i = 0; i < mesh.ndim(); ++i) {
                 for (size_t j = 0; j < mesh.ndim(); ++j) {
                     disp(n, i) += F(i, j) * coor(n, j);
                 }
             }
         }
 
         auto f = quad.GradN_vector(vec.AsElement(disp));
 
         REQUIRE(xt::has_shape(f, {mesh.nelem(), quad.nip(), td, td}));
         REQUIRE(xt::has_shape(f, quad.shape_qtensor<2>()));
 
         for (size_t e = 0; e < mesh.nelem(); ++e) {
             for (size_t q = 0; q < quad.nip(); ++q) {
                 REQUIRE(xt::allclose(xt::view(f, e, q), xt::transpose(F)));
             }
         }
     }
 
     SECTION("GradN_vector_T")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
         GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor()));
 
         size_t td = mesh.ndim();
         xt::xtensor<double, 2> F = xt::zeros<double>({td, td});
         xt::xtensor<double, 2> EPS = xt::zeros<double>({td, td});
 
         F(0, 1) = 0.1;
 
         auto coor = mesh.coor();
         auto disp = xt::zeros_like(coor);
 
         for (size_t n = 0; n < mesh.nnode(); ++n) {
             for (size_t i = 0; i < mesh.ndim(); ++i) {
                 for (size_t j = 0; j < mesh.ndim(); ++j) {
                     disp(n, i) += F(i, j) * coor(n, j);
                 }
             }
         }
 
         auto f = quad.GradN_vector_T(vec.AsElement(disp));
 
         REQUIRE(xt::has_shape(f, {mesh.nelem(), quad.nip(), td, td}));
         REQUIRE(xt::has_shape(f, quad.shape_qtensor<2>()));
 
         for (size_t e = 0; e < mesh.nelem(); ++e) {
             for (size_t q = 0; q < quad.nip(); ++q) {
                 REQUIRE(xt::allclose(xt::view(f, e, q), F));
             }
         }
     }
 
     SECTION("SymGradN_vector")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
         GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor()));
 
         size_t td = mesh.ndim();
         xt::xtensor<double, 2> F = xt::zeros<double>({td, td});
         xt::xtensor<double, 2> EPS = xt::zeros<double>({td, td});
 
         F(0, 1) = 0.1;
         EPS(0, 1) = 0.05;
         EPS(1, 0) = 0.05;
 
         auto coor = mesh.coor();
         auto disp = xt::zeros_like(coor);
 
         for (size_t n = 0; n < mesh.nnode(); ++n) {
             for (size_t i = 0; i < mesh.ndim(); ++i) {
                 for (size_t j = 0; j < mesh.ndim(); ++j) {
                     disp(n, i) += F(i, j) * coor(n, j);
                 }
             }
         }
 
         auto eps = quad.SymGradN_vector(vec.AsElement(disp));
         auto dV = quad.AsTensor<2>(quad.dV());
         auto epsbar = xt::average(eps, dV, {0, 1});
 
         REQUIRE(xt::has_shape(eps, {mesh.nelem(), quad.nip(), td, td}));
         REQUIRE(xt::has_shape(eps, quad.shape_qtensor<2>()));
         REQUIRE(xt::has_shape(dV, {mesh.nelem(), quad.nip(), td, td}));
         REQUIRE(xt::has_shape(dV, quad.shape_qtensor<2>()));
         REQUIRE(xt::has_shape(epsbar, {td, td}));
 
         for (size_t e = 0; e < mesh.nelem(); ++e) {
             for (size_t q = 0; q < quad.nip(); ++q) {
                 REQUIRE(xt::allclose(xt::view(eps, e, q), EPS));
             }
         }
 
         REQUIRE(xt::allclose(epsbar, EPS));
     }
 
+    SECTION("Int_N_vector_dV")
+    {
+        GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
+        GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
+        GooseFEM::Element::Quad4::Quadrature quad(vector.AsElement(mesh.coor()));
+
+        xt::xtensor<double, 3> qvector = xt::empty<double>({mesh.nelem(), quad.nip(), mesh.ndim()});
+        xt::view(qvector, xt::all(), xt::all(), 0) = 2.0;
+        xt::view(qvector, xt::all(), xt::all(), 1) = 5.0;
+
+        auto elemvec = quad.Int_N_vector_dV(qvector);
+        auto nodevec = vector.AssembleNode(elemvec);
+
+        REQUIRE(xt::allclose(xt::view(nodevec, xt::all(), 0), 2.0));
+        REQUIRE(xt::allclose(xt::view(nodevec, xt::all(), 1), 5.0));
+    }
+
     SECTION("Int_N_scalar_NT_dV")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic());
         GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic());
         GooseFEM::Element::Quad4::Quadrature quad(
             vec.AsElement(mesh.coor()),
             GooseFEM::Element::Quad4::Nodal::xi(),
             GooseFEM::Element::Quad4::Nodal::w());
 
         xt::xtensor<double, 2> rho = xt::ones<double>({mesh.nelem(), quad.nip()});
 
         mat.assemble(quad.Int_N_scalar_NT_dV(rho));
         auto M = mat.Todiagonal();
 
         REQUIRE(M.size() == vec.ndof());
         REQUIRE(xt::allclose(M, 1.0));
     }
 
     SECTION("Int_gradN_dot_tensor2_dV")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic());
         GooseFEM::Element::Quad4::Quadrature quad(vec.AsElement(mesh.coor()));
 
         size_t td = mesh.ndim();
         xt::xtensor<double, 2> F = xt::zeros<double>({td, td});
 
         F(0, 1) = 0.1;
 
         auto coor = mesh.coor();
         xt::xtensor<double, 2> disp = xt::zeros<double>(coor.shape());
 
         for (size_t n = 0; n < mesh.nnode(); ++n) {
             for (size_t i = 0; i < mesh.ndim(); ++i) {
                 for (size_t j = 0; j < mesh.ndim(); ++j) {
                     disp(n, i) += F(i, j) * coor(n, j);
                 }
             }
         }
 
         auto eps = quad.SymGradN_vector(vec.AsElement(disp));
         auto Fi = vec.AssembleDofs(quad.Int_gradN_dot_tensor2_dV(eps));
 
         REQUIRE(Fi.size() == vec.ndof());
         REQUIRE(xt::allclose(Fi, 0.0));
     }
 }