diff --git a/develop/ElementHex8.cpp b/develop/ElementHex8.cpp
index 957d328..897098f 100644
--- a/develop/ElementHex8.cpp
+++ b/develop/ElementHex8.cpp
@@ -1,157 +1,174 @@
 
 #include "support.h"
 
 using T2  = cppmat::tiny::cartesian::tensor2 <double,3>;
 using T2s = cppmat::tiny::cartesian::tensor2s<double,3>;
 
 // =================================================================================================
 
-TEST_CASE("xGooseFEM::ElementHex8", "ElementHex8.h")
+TEST_CASE("GooseFEM::ElementHex8", "ElementHex8.h")
 {
 
-using T2 = xt::xtensor_fixed<double, xt::xshape<3,3>>;
-
 // =================================================================================================
 
 SECTION( "int_N_scalar_NT_dV" )
 {
   // mesh
-  xGooseFEM::Mesh::Hex8::Regular mesh(3,3,3);
+  GooseFEM::Mesh::Hex8::Regular mesh(3,3,3);
 
   // vector-definition, and a diagonal matrix
-  xGooseFEM::Vector         vec(mesh.conn(), mesh.dofsPeriodic());
-  xGooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector         vec(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic());
 
   // element definition, with nodal quadrature
-  xGooseFEM::Element::Hex8::Quadrature quad(
+  GooseFEM::Element::Hex8::Quadrature quad(
     vec.asElement(mesh.coor()),
-    xGooseFEM::Element::Hex8::Nodal::xi(),
-    xGooseFEM::Element::Hex8::Nodal::w()
+    GooseFEM::Element::Hex8::Nodal::xi(),
+    GooseFEM::Element::Hex8::Nodal::w()
   );
 
   // scalar per quadrature point (e.g. mass-density "rho")
-  xt::xtensor<double,2> rho = xt::ones<double>({mesh.nelem(), quad.nip()});
+  GooseFEM::ArrD rho = GooseFEM::ArrD::Constant({mesh.nelem(), quad.nip()}, 1.);
 
   // evaluate integral and assemble diagonal matrix (e.g. mass matrix)
   mat.assemble(quad.int_N_scalar_NT_dV(rho));
 
   // check matrix
   // - get the matrix
-  xt::xtensor<double,1> M = mat.asDiagonal();
+  GooseFEM::ColD M = mat.asDiagonal();
   // - check the size
   REQUIRE( M.size() == vec.ndof() );
   // - check each component
-  for ( size_t i = 0 ; i < M.size() ; ++i )
+  for ( auto i = 0 ; i < M.size() ; ++i )
     EQ( M(i), 1 );
 }
 
 // =================================================================================================
 
 SECTION( "symGradN_vector" )
 {
   // mesh
-  xGooseFEM::Mesh::Hex8::FineLayer mesh(27,27,27);
+  GooseFEM::Mesh::Hex8::FineLayer mesh(27,27,27);
 
   // vector-definition
-  xGooseFEM::Vector vec(mesh.conn(), mesh.dofs());
+  GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
 
   // element definition, with Gauss quadrature
-  xGooseFEM::Element::Hex8::Quadrature quad( vec.asElement(mesh.coor()) );
+  GooseFEM::Element::Hex8::Quadrature quad( vec.asElement(mesh.coor()) );
 
-  // macroscopic deformation gradient and strain
+  // macroscopic deformation gradient
+  // - allocate
+  T2 F;
   // - zero-initialize
-  T2 F   = xt::zeros<double>({3,3});
-  T2 EPS = xt::zeros<double>({3,3});
+  F.setZero();
   // - set non-zero components
-  F  (0,1) = 0.1;
-  EPS(0,1) = 0.05;
-  EPS(1,0) = 0.05;
+  F(0,1) = 0.1;
+
+  // convert the macroscopic strain tensor
+  T2 EPS = .5 * ( F + F.T() );
 
-  // nodal coordinates and displacement
-  xt::xtensor<double,2> coor = mesh.coor();;
-  xt::xtensor<double,2> disp = xt::zeros<double>(coor.shape());
+  // nodal coordinates
+  GooseFEM::MatD coor = mesh.coor();;
+
+  // nodal displacement
+  // - allocate
+  GooseFEM::MatD disp(mesh.nnode(), mesh.ndim());
+  // - zero-initialize
+  disp.setZero();
 
   // apply macroscopic deformation gradient
   for ( size_t n = 0 ; n < mesh.nnode() ; ++n )
-    for ( size_t i = 0 ; i < F.shape()[0] ; ++i )
-      for ( size_t j = 0 ; j < F.shape()[1] ; ++j )
+    for ( size_t i = 0 ; i < F.ndim() ; ++i )
+      for ( size_t j = 0 ; j < F.ndim() ; ++j )
         disp(n,i) += F(i,j) * coor(n,j);
 
   // compute quadrature point tensors
-  xt::xtensor<double,4> eps = quad.symGradN_vector(vec.asElement(disp));
-
-  // integration point volume
-  xt::xtensor<double,4> dV = eps;
-  quad.dV(dV);
+  GooseFEM::ArrD eps = quad.symGradN_vector(vec.asElement(disp));
 
-  auto epsbar = xt::sum(dV*eps, {0,1}) / xt::sum(dV, {0,1});
+  // compute volume averaged tensor
+  GooseFEM::ArrD epsbar = eps.average(quad.dV(eps.shape(-1)), {0,1});
 
   // check
+  // - temporary tensor, to view the tensors
+  cppmat::view::cartesian::tensor2s<double,3> Eps;
   // - check sizes
-  REQUIRE( eps.shape()[0] == mesh.nelem() );
-  REQUIRE( eps.shape()[1] == quad.nip()   );
-  REQUIRE( eps.shape()[2] == mesh.ndim()   );
-  REQUIRE( eps.shape()[3] == mesh.ndim()   );
+  REQUIRE( eps.shape(0) == mesh.nelem() );
+  REQUIRE( eps.shape(1) == quad.nip()   );
+  REQUIRE( eps.shape(2) == Eps.size()   );
   // - check all components
   for ( size_t e = 0 ; e < mesh.nelem() ; ++e ) {
     for ( size_t k = 0 ; k < quad.nip() ; ++k ) {
-      auto Eps = xt::view(eps, e, k, xt::all(), xt::all());
-      for ( size_t i = 0 ; i < Eps.shape()[0] ; ++i )
-        for ( size_t j = 0 ; j < Eps.shape()[1] ; ++j )
+      Eps.setMap(&eps(e,k));
+      for ( size_t i = 0 ; i < Eps.ndim() ; ++i )
+        for ( size_t j = 0 ; j < Eps.ndim() ; ++j )
           EQ( Eps(i,j), EPS(i,j) );
     }
   }
 
   // check macroscopic tensor
-  for ( size_t i = 0 ; i < epsbar.shape()[0] ; ++i )
-    for ( size_t j = 0 ; j < epsbar.shape()[1] ; ++j )
-      EQ( epsbar(i,j), EPS(i,j) );
+  // - convert to tensor object
+  T2s Epsbar = T2s::Copy(epsbar.begin(), epsbar.end());
+  // - check all components
+  for ( size_t i = 0 ; i < Epsbar.ndim() ; ++i )
+    for ( size_t j = 0 ; j < Epsbar.ndim() ; ++j )
+      EQ( Epsbar(i,j), EPS(i,j) );
 }
 
 // =================================================================================================
 
 SECTION( "symGradN_vector, int_gradN_dot_tensor2s_dV" )
 {
   // mesh
-  xGooseFEM::Mesh::Hex8::FineLayer mesh(27,27,27);
+  GooseFEM::Mesh::Hex8::FineLayer mesh(27,27,27);
 
   // vector-definition
-  xGooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic());
 
   // element definition, with Gauss quadrature
-  xGooseFEM::Element::Hex8::Quadrature quad( vec.asElement(mesh.coor()) );
+  GooseFEM::Element::Hex8::Quadrature quad( vec.asElement(mesh.coor()) );
 
-  // macroscopic deformation gradient and strain
+  // macroscopic deformation gradient
+  // - allocate
+  T2 F;
   // - zero-initialize
-  T2 F = xt::zeros<double>({3,3});
+  F.setZero();
   // - set non-zero components
   F(0,1) = 0.1;
 
-  // nodal coordinates and displacement
-  xt::xtensor<double,2> coor = mesh.coor();;
-  xt::xtensor<double,2> disp = xt::zeros<double>(coor.shape());
+  // nodal coordinates
+  GooseFEM::MatD coor = mesh.coor();;
+
+  // nodal displacement
+  // - allocate
+  GooseFEM::MatD disp(mesh.nnode(), mesh.ndim());
+  // - zero-initialize
+  disp.setZero();
 
   // apply macroscopic deformation gradient
   for ( size_t n = 0 ; n < mesh.nnode() ; ++n )
-    for ( size_t i = 0 ; i < F.shape()[0] ; ++i )
-      for ( size_t j = 0 ; j < F.shape()[1] ; ++j )
+    for ( size_t i = 0 ; i < F.ndim() ; ++i )
+      for ( size_t j = 0 ; j < F.ndim() ; ++j )
         disp(n,i) += F(i,j) * coor(n,j);
 
   // compute quadrature point tensors
-  xt::xtensor<double,4> eps = quad.symGradN_vector(vec.asElement(disp));
+  GooseFEM::ArrD eps = quad.symGradN_vector(vec.asElement(disp));
 
   // nodal force vector (should be zero, as it is only sensitive to periodic fluctuations)
-  xt::xtensor<double,1> Fi = vec.assembleDofs(quad.int_gradN_dot_tensor2_dV(eps));
+  GooseFEM::ColD Fi = vec.assembleDofs(quad.int_gradN_dot_tensor2s_dV(eps));
+
+  for ( size_t i = 0 ; i < vec.ndof() ; ++i )
+    if ( std::abs(Fi(i)) > 1.e-12 )
+      std::cout << i << ", " << Fi(i) << std::endl;
 
   // check
   // - size
   REQUIRE( Fi.size() == vec.ndof() );
   // - check all components
   for ( size_t i = 0 ; i < vec.ndof() ; ++i )
     EQ( Fi(i), 0 );
 
 }
 
 // =================================================================================================
 
 }
diff --git a/develop/ElementQuad4.cpp b/develop/ElementQuad4.cpp
index 5ffbf64..1d997e8 100644
--- a/develop/ElementQuad4.cpp
+++ b/develop/ElementQuad4.cpp
@@ -1,154 +1,171 @@
 
 #include "support.h"
 
+using T2  = cppmat::tiny::cartesian::tensor2 <double,2>;
+using T2s = cppmat::tiny::cartesian::tensor2s<double,2>;
+
 // =================================================================================================
 
-TEST_CASE("xGooseFEM::ElementQuad4", "ElementQuad4.h")
+TEST_CASE("GooseFEM::ElementQuad4", "ElementQuad4.h")
 {
 
-using T2 = xt::xtensor_fixed<double, xt::xshape<2,2>>;
-
 // =================================================================================================
 
 SECTION( "int_N_scalar_NT_dV" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(3,3);
+  GooseFEM::Mesh::Quad4::Regular mesh(3,3);
 
   // vector-definition, and a diagonal matrix
-  xGooseFEM::Vector         vec(mesh.conn(), mesh.dofsPeriodic());
-  xGooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector         vec(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::MatrixDiagonal mat(mesh.conn(), mesh.dofsPeriodic());
 
   // element definition, with nodal quadrature
-  xGooseFEM::Element::Quad4::Quadrature quad(
+  GooseFEM::Element::Quad4::Quadrature quad(
     vec.asElement(mesh.coor()),
-    xGooseFEM::Element::Quad4::Nodal::xi(),
-    xGooseFEM::Element::Quad4::Nodal::w()
+    GooseFEM::Element::Quad4::Nodal::xi(),
+    GooseFEM::Element::Quad4::Nodal::w()
   );
 
   // scalar per quadrature point (e.g. mass-density "rho")
-  xt::xtensor<double,2> rho = xt::ones<double>({mesh.nelem(), quad.nip()});
+  GooseFEM::ArrD rho({mesh.nelem(), quad.nip()});
+  rho.setConstant(1.);
 
   // evaluate integral and assemble diagonal matrix (e.g. mass matrix)
   mat.assemble(quad.int_N_scalar_NT_dV(rho));
 
   // check matrix
   // - get the matrix
-  xt::xtensor<double,1> M = mat.asDiagonal();
+  GooseFEM::ColD M = mat.asDiagonal();
   // - check the size
   REQUIRE( M.size() == vec.ndof() );
   // - check each component
-  for ( size_t i = 0 ; i < M.size() ; ++i )
+  for ( auto i = 0 ; i < M.size() ; ++i )
     EQ( M(i), 1 );
 }
 
 // =================================================================================================
 
 SECTION( "symGradN_vector" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::FineLayer mesh(27,27);
+  GooseFEM::Mesh::Quad4::FineLayer mesh(27,27);
 
   // vector-definition
-  xGooseFEM::Vector vec(mesh.conn(), mesh.dofs());
+  GooseFEM::Vector vec(mesh.conn(), mesh.dofs());
 
   // element definition, with Gauss quadrature
-  xGooseFEM::Element::Quad4::Quadrature quad( vec.asElement(mesh.coor()) );
+  GooseFEM::Element::Quad4::Quadrature quad( vec.asElement(mesh.coor()) );
 
-  // macroscopic deformation gradient and strain
+  // macroscopic deformation gradient
+  // - allocate
+  T2 F;
   // - zero-initialize
-  T2 F   = xt::zeros<double>({2,2});
-  T2 EPS = xt::zeros<double>({2,2});
+  F.setZero();
   // - set non-zero components
-  F  (0,1) = 0.1;
-  EPS(0,1) = 0.05;
-  EPS(1,0) = 0.05;
+  F(0,1) = 0.1;
+
+  // convert the macroscopic strain tensor
+  T2 EPS = .5 * ( F + F.T() );
+
+  // nodal coordinates
+  GooseFEM::MatD coor = mesh.coor();;
 
-  // nodal coordinates and displacement
-  xt::xtensor<double,2> coor = mesh.coor();;
-  xt::xtensor<double,2> disp = xt::zeros<double>(coor.shape());
+  // nodal displacement
+  // - allocate
+  GooseFEM::MatD disp(mesh.nnode(), mesh.ndim());
+  // - zero-initialize
+  disp.setZero();
 
   // apply macroscopic deformation gradient
   for ( size_t n = 0 ; n < mesh.nnode() ; ++n )
-    for ( size_t i = 0 ; i < F.shape()[0] ; ++i )
-      for ( size_t j = 0 ; j < F.shape()[1] ; ++j )
+    for ( size_t i = 0 ; i < F.ndim() ; ++i )
+      for ( size_t j = 0 ; j < F.ndim() ; ++j )
         disp(n,i) += F(i,j) * coor(n,j);
 
   // compute quadrature point tensors
-  xt::xtensor<double,4> eps = quad.symGradN_vector(vec.asElement(disp));
-
-  // integration point volume
-  xt::xtensor<double,4> dV = eps;
-  quad.dV(dV);
+  GooseFEM::ArrD eps = quad.symGradN_vector(vec.asElement(disp));
 
-  auto epsbar = xt::sum(dV*eps, {0,1}) / xt::sum(dV, {0,1});
+  // compute volume averaged tensor
+  GooseFEM::ArrD epsbar = eps.average(quad.dV(eps.shape(-1)), {0,1});
 
   // check
+  // - temporary tensor, to view the tensors
+  cppmat::view::cartesian::tensor2s<double,2> Eps;
   // - check sizes
-  REQUIRE( eps.shape()[0] == mesh.nelem() );
-  REQUIRE( eps.shape()[1] == quad.nip()   );
-  REQUIRE( eps.shape()[2] == mesh.ndim()   );
-  REQUIRE( eps.shape()[3] == mesh.ndim()   );
+  REQUIRE( eps.shape(0) == mesh.nelem() );
+  REQUIRE( eps.shape(1) == quad.nip()   );
+  REQUIRE( eps.shape(2) == Eps.size()   );
   // - check all components
   for ( size_t e = 0 ; e < mesh.nelem() ; ++e ) {
     for ( size_t k = 0 ; k < quad.nip() ; ++k ) {
-      auto Eps = xt::view(eps, e, k, xt::all(), xt::all());
-      for ( size_t i = 0 ; i < Eps.shape()[0] ; ++i )
-        for ( size_t j = 0 ; j < Eps.shape()[1] ; ++j )
+      Eps.setMap(&eps(e,k));
+      for ( size_t i = 0 ; i < Eps.ndim() ; ++i )
+        for ( size_t j = 0 ; j < Eps.ndim() ; ++j )
           EQ( Eps(i,j), EPS(i,j) );
     }
   }
 
   // check macroscopic tensor
-  for ( size_t i = 0 ; i < epsbar.shape()[0] ; ++i )
-    for ( size_t j = 0 ; j < epsbar.shape()[1] ; ++j )
-      EQ( epsbar(i,j), EPS(i,j) );
+  // - convert to tensor object
+  T2s Epsbar = T2s::Copy(epsbar.begin(), epsbar.end());
+  // - check all components
+  for ( size_t i = 0 ; i < Epsbar.ndim() ; ++i )
+    for ( size_t j = 0 ; j < Epsbar.ndim() ; ++j )
+      EQ( Epsbar(i,j), EPS(i,j) );
 }
 
 // =================================================================================================
 
 SECTION( "symGradN_vector, int_gradN_dot_tensor2s_dV" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::FineLayer mesh(27,27);
+  GooseFEM::Mesh::Quad4::FineLayer mesh(27,27);
 
   // vector-definition
-  xGooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector vec(mesh.conn(), mesh.dofsPeriodic());
 
   // element definition, with Gauss quadrature
-  xGooseFEM::Element::Quad4::Quadrature quad( vec.asElement(mesh.coor()) );
+  GooseFEM::Element::Quad4::Quadrature quad( vec.asElement(mesh.coor()) );
 
-  // macroscopic deformation gradient and strain
+  // macroscopic deformation gradient
+  // - allocate
+  T2 F;
   // - zero-initialize
-  T2 F = xt::zeros<double>({2,2});
+  F.setZero();
   // - set non-zero components
   F(0,1) = 0.1;
 
-  // nodal coordinates and displacement
-  xt::xtensor<double,2> coor = mesh.coor();;
-  xt::xtensor<double,2> disp = xt::zeros<double>(coor.shape());
+  // nodal coordinates
+  GooseFEM::MatD coor = mesh.coor();;
+
+  // nodal displacement
+  // - allocate
+  GooseFEM::MatD disp(mesh.nnode(), mesh.ndim());
+  // - zero-initialize
+  disp.setZero();
 
   // apply macroscopic deformation gradient
   for ( size_t n = 0 ; n < mesh.nnode() ; ++n )
-    for ( size_t i = 0 ; i < F.shape()[0] ; ++i )
-      for ( size_t j = 0 ; j < F.shape()[1] ; ++j )
+    for ( size_t i = 0 ; i < F.ndim() ; ++i )
+      for ( size_t j = 0 ; j < F.ndim() ; ++j )
         disp(n,i) += F(i,j) * coor(n,j);
 
   // compute quadrature point tensors
-  xt::xtensor<double,4> eps = quad.symGradN_vector(vec.asElement(disp));
+  GooseFEM::ArrD eps = quad.symGradN_vector(vec.asElement(disp));
 
   // nodal force vector (should be zero, as it is only sensitive to periodic fluctuations)
-  xt::xtensor<double,1> Fi = vec.assembleDofs(quad.int_gradN_dot_tensor2_dV(eps));
+  GooseFEM::ColD Fi = vec.assembleDofs(quad.int_gradN_dot_tensor2s_dV(eps));
 
   // check
   // - size
   REQUIRE( Fi.size() == vec.ndof() );
   // - check all components
   for ( size_t i = 0 ; i < vec.ndof() ; ++i )
     EQ( Fi(i), 0 );
 
 }
 
 // =================================================================================================
 
 }
diff --git a/develop/Iterate.cpp b/develop/Iterate.cpp
index f72505a..06e2229 100644
--- a/develop/Iterate.cpp
+++ b/develop/Iterate.cpp
@@ -1,29 +1,29 @@
 
 #include "support.h"
 
 // =================================================================================================
 
-TEST_CASE("xGooseFEM::Iterate", "Iterate.h")
+TEST_CASE("GooseFEM::Iterate", "Iterate.h")
 {
 
 // =================================================================================================
 
 SECTION( "StopList" )
 {
-  xGooseFEM::Iterate::StopList stop(5);
+  GooseFEM::Iterate::StopList stop(5);
 
   REQUIRE( stop.stop(5.e+0, 1.e-3) == false );
   REQUIRE( stop.stop(5.e+1, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-1, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-2, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-3, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-4, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-4, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-4, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-4, 1.e-3) == false );
   REQUIRE( stop.stop(5.e-4, 1.e-3) == true  );
 }
 
 // =================================================================================================
 
 }
diff --git a/develop/MatrixDiagonal.cpp b/develop/MatrixDiagonal.cpp
index f0209cc..d877738 100644
--- a/develop/MatrixDiagonal.cpp
+++ b/develop/MatrixDiagonal.cpp
@@ -1,82 +1,98 @@
 
 #include "support.h"
 
 // =================================================================================================
 
-TEST_CASE("xGooseFEM::MatrixDiagonal", "MatrixDiagonal.h")
+TEST_CASE("GooseFEM::MatrixDiagonal", "MatrixDiagonal.h")
 {
 
 // =================================================================================================
 
 SECTION( "dot" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(2,2);
+  GooseFEM::Mesh::Quad4::Regular mesh(2,2);
 
   // random matrix and column
-  xt::xtensor<double,1> a = xt::random::rand<double>({mesh.nnode()*mesh.ndim()});
-  xt::xtensor<double,1> b = xt::random::rand<double>({mesh.nnode()*mesh.ndim()});
-  xt::xtensor<double,1> c;
-
-  // compute product
+  GooseFEM::MatD a = GooseFEM::MatD::Random(mesh.nnode()*mesh.ndim(),mesh.nnode()*mesh.ndim());
+  GooseFEM::ColD b = GooseFEM::ColD::Random(mesh.nnode()*mesh.ndim());
+  GooseFEM::ColD c;
+
+  // set diagonal
+  for ( auto i = 0 ; i < a.rows() ; ++i ) {
+    for ( auto j = 0 ; j < a.cols() ; ++j ) {
+      if ( i != j ) a(i,j)  = 0.0;
+      else          a(i,j) += 1.0;
+    }
+  }
+
+  // compute product using Eigen
   c = a * b;
 
-  // convert to xGooseFEM
+  // convert to GooseFEM
   // - allocate
-  xGooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs());
-  xt::xtensor<double,1> C;
+  GooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs());
+  GooseFEM::ColD C;
   // - set
-  for ( size_t i = 0 ; i < a.shape()[0] ; ++i )
-    A(i,i) = a(i);
+  for ( auto i = 0 ; i < a.rows() ; ++i )
+    A(i,i) = a(i,i);
 
   // compute product
-  C = A.dot(b);
+  C = A * b;
 
   // check
   // - size
   REQUIRE( C.size() == c.size() );
   // - components
-  for ( size_t i = 0 ; i < c.shape()[0] ; ++i )
+  for ( auto i = 0 ; i < c.rows() ; ++i )
     EQ( C(i), c(i) );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 SECTION( "solve" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(2,2);
+  GooseFEM::Mesh::Quad4::Regular mesh(2,2);
 
   // random matrix and column
-  xt::xtensor<double,1> a = xt::random::rand<double>({mesh.nnode()*mesh.ndim()});
-  xt::xtensor<double,1> b = xt::random::rand<double>({mesh.nnode()*mesh.ndim()});
-  xt::xtensor<double,1> c;
-
-  // compute product
+  GooseFEM::MatD a = GooseFEM::MatD::Random(mesh.nnode()*mesh.ndim(),mesh.nnode()*mesh.ndim());
+  GooseFEM::ColD b = GooseFEM::ColD::Random(mesh.nnode()*mesh.ndim());
+  GooseFEM::ColD c;
+
+  // set diagonal
+  for ( auto i = 0 ; i < a.rows() ; ++i ) {
+    for ( auto j = 0 ; j < a.cols() ; ++j ) {
+      if ( i != j ) a(i,j)  = 0.0;
+      else          a(i,j) += 1.0;
+    }
+  }
+
+  // compute product using Eigen
   c = a * b;
 
-  // convert to xGooseFEM
+  // convert to GooseFEM
   // - allocate
-  xGooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs());
-  xt::xtensor<double,1> B, C;
+  GooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs());
+  GooseFEM::ColD B, C;
   // - set
-  for ( size_t i = 0 ; i < a.shape()[0] ; ++i )
-    A(i,i) = a(i);
+  for ( auto i = 0 ; i < a.rows() ; ++i )
+    A(i,i) = a(i,i);
 
   // compute product
-  C = A.dot(b);
+  C = A * b;
 
   // solve
   B = A.solve(C);
 
   // check
   // - size
   REQUIRE( B.size() == b.size() );
   // - components
-  for ( size_t i = 0 ; i < b.shape()[0] ; ++i )
+  for ( auto i = 0 ; i < b.rows() ; ++i )
     EQ( B(i), b(i) );
 }
 
-// // =================================================================================================
+// =================================================================================================
 
 }
diff --git a/develop/Vector.cpp b/develop/Vector.cpp
index ba1130e..7079dc7 100644
--- a/develop/Vector.cpp
+++ b/develop/Vector.cpp
@@ -1,178 +1,178 @@
 
 #include "support.h"
 
 // =================================================================================================
 
-TEST_CASE("xGooseFEM::Vector", "Vector.h")
+TEST_CASE("GooseFEM::Vector", "Vector.h")
 {
 
 // =================================================================================================
 
 SECTION( "asDofs - nodevec" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(2,2);
+  GooseFEM::Mesh::Quad4::Regular mesh(2,2);
 
   // vector-definition
-  xGooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
 
   // velocity field
   // - allocate
-  xt::xtensor<double,2> v = xt::empty<double>({mesh.nnode(), std::size_t(2)});
+  GooseFEM::MatD v(mesh.nnode(),2);
   // - set periodic
   v(0,0) = 1.0;  v(0,1) = 0.0;
   v(1,0) = 1.0;  v(1,1) = 0.0;
   v(2,0) = 1.0;  v(2,1) = 0.0;
   v(3,0) = 1.5;  v(3,1) = 0.0;
   v(4,0) = 1.5;  v(4,1) = 0.0;
   v(5,0) = 1.5;  v(5,1) = 0.0;
   v(6,0) = 1.0;  v(6,1) = 0.0;
   v(7,0) = 1.0;  v(7,1) = 0.0;
   v(8,0) = 1.0;  v(8,1) = 0.0;
 
   // convert to DOFs
-  xt::xtensor<double,1> V = vector.asDofs(v);
+  GooseFEM::ColD V = vector.asDofs(v);
 
   // check
   // - size
   REQUIRE( V.size() == mesh.nnodePeriodic() * mesh.ndim() );
   // - individual entries
   EQ( V(0), v(0,0) );
   EQ( V(1), v(0,1) );
   EQ( V(2), v(1,0) );
   EQ( V(3), v(1,1) );
   EQ( V(4), v(3,0) );
   EQ( V(5), v(3,1) );
   EQ( V(6), v(4,0) );
   EQ( V(7), v(4,1) );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 SECTION( "asDofs - elemvec" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(2,2);
+  GooseFEM::Mesh::Quad4::Regular mesh(2,2);
 
   // vector-definition
-  xGooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
 
   // velocity field
   // - allocate
-  xt::xtensor<double,2> v = xt::empty<double>({mesh.nnode(), std::size_t(2)});
+  GooseFEM::MatD v(mesh.nnode(),2);
   // - set periodic
   v(0,0) = 1.0;  v(0,1) = 0.0;
   v(1,0) = 1.0;  v(1,1) = 0.0;
   v(2,0) = 1.0;  v(2,1) = 0.0;
   v(3,0) = 1.5;  v(3,1) = 0.0;
   v(4,0) = 1.5;  v(4,1) = 0.0;
   v(5,0) = 1.5;  v(5,1) = 0.0;
   v(6,0) = 1.0;  v(6,1) = 0.0;
   v(7,0) = 1.0;  v(7,1) = 0.0;
   v(8,0) = 1.0;  v(8,1) = 0.0;
 
 
   // convert to DOFs - element - DOFs
-  xt::xtensor<double,1> V = vector.asDofs(vector.asElement(vector.asDofs(v)));
+  GooseFEM::ColD V = vector.asDofs(vector.asElement(vector.asDofs(v)));
 
   // check
   // - size
   REQUIRE( V.size() == mesh.nnodePeriodic() * mesh.ndim() );
   // - individual entries
   EQ( V(0), v(0,0) );
   EQ( V(1), v(0,1) );
   EQ( V(2), v(1,0) );
   EQ( V(3), v(1,1) );
   EQ( V(4), v(3,0) );
   EQ( V(5), v(3,1) );
   EQ( V(6), v(4,0) );
   EQ( V(7), v(4,1) );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 SECTION( "asDofs - assembleDofs" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(2,2);
+  GooseFEM::Mesh::Quad4::Regular mesh(2,2);
 
   // vector-definition
-  xGooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
 
   // force field
   // - allocate
-  xt::xtensor<double,2> f = xt::empty<double>({mesh.nnode(), std::size_t(2)});
+  GooseFEM::MatD f(mesh.nnode(),2);
   // - set periodic
   f(0,0) = -1.0;  f(0,1) = -1.0;
   f(1,0) =  0.0;  f(1,1) = -1.0;
   f(2,0) =  1.0;  f(2,1) = -1.0;
   f(3,0) = -1.0;  f(3,1) =  0.0;
   f(4,0) =  0.0;  f(4,1) =  0.0;
   f(5,0) =  1.0;  f(5,1) =  0.0;
   f(6,0) = -1.0;  f(6,1) =  1.0;
   f(7,0) =  0.0;  f(7,1) =  1.0;
   f(8,0) =  1.0;  f(8,1) =  1.0;
 
   // assemble as DOFs
-  xt::xtensor<double,1> F = vector.assembleDofs(f);
+  GooseFEM::ColD F = vector.assembleDofs(f);
 
   // check
   // - size
   REQUIRE( F.size() == mesh.nnodePeriodic() * mesh.ndim() );
   // - 'analytical' result
   EQ( F(0), 0 );
   EQ( F(1), 0 );
   EQ( F(2), 0 );
   EQ( F(3), 0 );
   EQ( F(4), 0 );
   EQ( F(5), 0 );
   EQ( F(6), 0 );
   EQ( F(7), 0 );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 SECTION( "asDofs - assembleNode" )
 {
   // mesh
-  xGooseFEM::Mesh::Quad4::Regular mesh(2,2);
+  GooseFEM::Mesh::Quad4::Regular mesh(2,2);
 
   // vector-definition
-  xGooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
+  GooseFEM::Vector vector(mesh.conn(), mesh.dofsPeriodic());
 
   // force field
   // - allocate
-  xt::xtensor<double,2> f = xt::empty<double>({mesh.nnode(), std::size_t(2)});
+  GooseFEM::MatD f(mesh.nnode(),2);
   // - set periodic
   f(0,0) = -1.0;  f(0,1) = -1.0;
   f(1,0) =  0.0;  f(1,1) = -1.0;
   f(2,0) =  1.0;  f(2,1) = -1.0;
   f(3,0) = -1.0;  f(3,1) =  0.0;
   f(4,0) =  0.0;  f(4,1) =  0.0;
   f(5,0) =  1.0;  f(5,1) =  0.0;
   f(6,0) = -1.0;  f(6,1) =  1.0;
   f(7,0) =  0.0;  f(7,1) =  1.0;
   f(8,0) =  1.0;  f(8,1) =  1.0;
 
   // convert to element, assemble as DOFs
-  xt::xtensor<double,1> F = vector.assembleDofs( vector.asElement(f) );
+  GooseFEM::ColD F = vector.assembleDofs( vector.asElement(f) );
 
   // check
   // - size
   REQUIRE( F.size() == mesh.nnodePeriodic() * mesh.ndim() );
   // - 'analytical' result
   EQ( F(0), 0 );
   EQ( F(1), 0 );
   EQ( F(2), 0 );
   EQ( F(3), 0 );
   EQ( F(4), 0 );
   EQ( F(5), 0 );
   EQ( F(6), 0 );
   EQ( F(7), 0 );
 }
 
 // -------------------------------------------------------------------------------------------------
 
 // =================================================================================================
 
 }
diff --git a/develop/support.h b/develop/support.h
index 580c8b8..77234e6 100644
--- a/develop/support.h
+++ b/develop/support.h
@@ -1,24 +1,13 @@
 
 #ifndef SUPPORT_H
 #define SUPPORT_H
 
 #include <catch2/catch.hpp>
-
 #include <cppmat/cppmat.h>
+#include <Eigen/Eigen>
 
-#include <xtensor/xrandom.hpp>
-#include <xtensor/xmath.hpp>
-
-#include "../include/xGooseFEM/GooseFEM.h"
+#include "../include/GooseFEM/GooseFEM.h"
 
 #define EQ(a,b) REQUIRE_THAT( (a), Catch::WithinAbs((b), 1.e-12) );
 
-
-template <class E, class W>
-inline auto Average(E&& e, W&& weights)
-{
-    return xt::sum(e*weights, {0,1}) / xt::sum(weights, {0,1});
-}
-
-
 #endif