diff --git a/fdaPDE/core.h b/fdaPDE/core.h
index e904448d..1f0852da 100644
--- a/fdaPDE/core.h
+++ b/fdaPDE/core.h
@@ -14,8 +14,8 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef __FDAPDE_CORE_MOUDLE_H__
-#define __FDAPDE_CORE_MOUDLE_H__
+#ifndef __FDAPDE_CORE_MODULE_H__
+#define __FDAPDE_CORE_MODULE_H__
 
 #include "utils.h"
 #include "fields.h"
@@ -26,4 +26,4 @@
 #include "splines.h"
 #include "multithreading.h"
 
-#endif   // __FDAPDE_CORE_MOUDLE_H__
+#endif   // __FDAPDE_CORE_MODULE_H__
diff --git a/fdaPDE/fields.h b/fdaPDE/fields.h
index 63595edf..676af1c3 100644
--- a/fdaPDE/fields.h
+++ b/fdaPDE/fields.h
@@ -14,8 +14,8 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef __FDAPDE_FIELDS_MOUDLE_H__
-#define __FDAPDE_FIELDS_MOUDLE_H__
+#ifndef __FDAPDE_FIELDS_MODULE_H__
+#define __FDAPDE_FIELDS_MODULE_H__
 
 #include "fields/differentiable_field.h"
 #include "fields/dot_product.h"
@@ -27,4 +27,4 @@
 #include "fields/vector_expressions.h"
 #include "fields/vector_field.h"
 
-#endif   // __FDAPDE_FIELDS_MOUDLE_H__
+#endif   // __FDAPDE_FIELDS_MODULE_H__
diff --git a/fdaPDE/finite_elements/basis/lagrangian_element.h b/fdaPDE/finite_elements/basis/lagrangian_element.h
index 8376ae13..559d3615 100644
--- a/fdaPDE/finite_elements/basis/lagrangian_element.h
+++ b/fdaPDE/finite_elements/basis/lagrangian_element.h
@@ -27,7 +27,7 @@ namespace fdapde {
 namespace core {
 
 // A Lagrangian basis of degree R over an M-dimensional element
-template <unsigned int M_, unsigned int R_> class LagrangianElement {
+template <int M_, int R_> class LagrangianElement {
    private:
     std::array<std::array<double, M_>, ct_binomial_coefficient(M_ + R_, R_)> nodes_;   // nodes of the Lagrangian basis
     std::array<MultivariatePolynomial<M_, R_>, ct_binomial_coefficient(M_ + R_, R_)> basis_;
@@ -35,9 +35,9 @@ template <unsigned int M_, unsigned int R_> class LagrangianElement {
     void compute_coefficients_(const std::array<std::array<double, M_>, ct_binomial_coefficient(M_ + R_, R_)>& nodes);
    public:
     // compile time informations
-    static constexpr unsigned int R = R_; // basis order
-    static constexpr unsigned int M = M_; // input space dimension
-    static constexpr unsigned int n_basis = ct_binomial_coefficient(M + R, R);
+    static constexpr int R = R_; // basis order
+    static constexpr int M = M_; // input space dimension
+    static constexpr int n_basis = ct_binomial_coefficient(M + R, R);
     typedef MultivariatePolynomial<M, R> ElementType;
     typedef typename std::array<MultivariatePolynomial<M, R>, n_basis>::const_iterator const_iterator;
 
@@ -59,11 +59,11 @@ template <unsigned int M_, unsigned int R_> class LagrangianElement {
 // implementative details
   
 // compute coefficients via Vandermonde matrix
-template <unsigned int M_, unsigned int R_>
+template <int M_, int R_>
 void LagrangianElement<M_, R_>::compute_coefficients_(
   const std::array<std::array<double, M_>, ct_binomial_coefficient(M_ + R_, R_)>& nodes) {
     // build vandermonde matrix
-    constexpr unsigned int n_basis = ct_binomial_coefficient(M_ + R_, R_);
+    constexpr int n_basis = ct_binomial_coefficient(M_ + R_, R_);
     constexpr std::array<std::array<unsigned, M_>, n_basis> poly_table = MultivariatePolynomial<M, R>::poly_table;
 
     // Vandermonde matrix construction
diff --git a/fdaPDE/finite_elements/basis/multivariate_polynomial.h b/fdaPDE/finite_elements/basis/multivariate_polynomial.h
index 70d7da65..20d8c5c2 100644
--- a/fdaPDE/finite_elements/basis/multivariate_polynomial.h
+++ b/fdaPDE/finite_elements/basis/multivariate_polynomial.h
@@ -47,10 +47,9 @@ namespace core {
 //   2   0   0  | 2
 //
 // ct_poly_exp() computes the above table for any choice of N and R at compile time.
-template <unsigned int N, unsigned int R>
-using PolyTable = std::array<std::array<unsigned, N>, ct_binomial_coefficient(R + N, R)>;
+template <int N, int R> using PolyTable = std::array<std::array<unsigned, N>, ct_binomial_coefficient(R + N, R)>;
 
-template <unsigned int N, unsigned int R> constexpr PolyTable<N, R> ct_poly_exp() {
+template <int N, int R> constexpr PolyTable<N, R> ct_poly_exp() {
     const int monomials = ct_binomial_coefficient(R + N, R);   // number of monomials in polynomial p(x)
     // initialize empty array
     PolyTable<N, R> coeff {};
@@ -98,10 +97,10 @@ template <unsigned int N, unsigned int R> constexpr PolyTable<N, R> ct_poly_exp(
 //     PolyTable                         PolyGradTable
 //
 // ct_grad_exp() computes the PolyGradTable above from a PolyTable for any choice of N and R at compile time.
-template <unsigned int N, unsigned int R>
+template <int N, int R>
 using PolyGradTable = std::array<std::array<std::array<unsigned, N>, ct_binomial_coefficient(R + N, R)>, N>;
 
-template <unsigned int N, unsigned int R> constexpr PolyGradTable<N, R> ct_grad_exp(const PolyTable<N, R> poly_table) {
+template <int N, int R> constexpr PolyGradTable<N, R> ct_grad_exp(const PolyTable<N, R> poly_table) {
     const int monomials = ct_binomial_coefficient(R + N, R);   // number of monomials in polynomial p(x)
     // initialize empty array
     PolyGradTable<N,R> coeff {};
@@ -120,9 +119,9 @@ template <unsigned int N, unsigned int R> constexpr PolyGradTable<N, R> ct_grad_
 }
 
 // recursive template based unfolding of monomial product x1^i1*x2^i2*...*xN^iN.
-template <unsigned int N,   // template recursion loop variable
-	  typename P,       // point where to evaluate the monomial
-	  typename V>       // a row of the polynomial PolyTable (i.e. an array of coefficients [i1 i2 ... iN])
+template <int N,        // template recursion loop variable
+          typename P,   // point where to evaluate the monomial
+          typename V>   // a row of the polynomial PolyTable (i.e. an array of coefficients [i1 i2 ... iN])
 struct MonomialProduct {
     static constexpr double unfold(const P& p, const V& v) {
         return v[N] == 0 ? MonomialProduct<N - 1, P, V>::unfold(p, v) :
@@ -135,11 +134,11 @@ template <typename P, typename V> struct MonomialProduct<0, P, V> {   // base ca
 };
 
 // unfold the sum of monomials at compile time to produce the complete polynomial expression
-template <unsigned int I,   // template recursion loop variable
-	  unsigned int N,   // total number of monomials to unfold
-	  unsigned int M,   // polynomial space dimension
-	  typename P,       // point where to evaluate the polynomial
-	  typename V>       // the whole polynomial PolyTable
+template <int I,        // template recursion loop variable
+          int N,        // total number of monomials to unfold
+          int M,        // polynomial space dimension
+          typename P,   // point where to evaluate the polynomial
+          typename V>   // the whole polynomial PolyTable
 struct MonomialSum {
     static constexpr double unfold(const std::array<double, N>& c, const P& p, const V& v) {
         return (c[I] * MonomialProduct<M - 1, SVector<M>, std::array<unsigned, M>>::unfold(p, v[I])) +
@@ -147,14 +146,14 @@ struct MonomialSum {
     }
 };
 // end of recursion
-template <unsigned int N, unsigned int M, typename P, typename V> struct MonomialSum<0, N, M, P, V> {
+template <int N, int M, typename P, typename V> struct MonomialSum<0, N, M, P, V> {
     static constexpr double unfold(const std::array<double, N>& c, const P& p, const V& v) {
         return c[0] * MonomialProduct<M - 1, SVector<M>, std::array<unsigned, M>>::unfold(p, v[0]);
     }
 };
 
 // functor implementing the derivative of a multivariate N-dimensional polynomial of degree R along a given direction
-template <unsigned int N, unsigned int R>
+template <int N, int R>
 class PolynomialDerivative : public VectorExpr<N, N, PolynomialDerivative<N, R>> {
    private:
     // compile time informations
@@ -184,17 +183,17 @@ class PolynomialDerivative : public VectorExpr<N, N, PolynomialDerivative<N, R>>
 };
 
 // class representing a multivariate polynomial of degree R defined over a space of dimension N
-template <unsigned int N, unsigned int R>
+template <int N, int R>
 class MultivariatePolynomial : public ScalarExpr<N, MultivariatePolynomial<N, R>> {
    private:
-    static const constexpr unsigned MON = ct_binomial_coefficient(R + N, R);
+    static const constexpr int MON = ct_binomial_coefficient(R + N, R);
     std::array<double, MON> coeff_vector_;   // vector of coefficients
     VectorField<N, N, PolynomialDerivative<N, R>> gradient_;
    public:
     // compute this at compile time once, let public access
     static constexpr PolyTable<N, R> poly_table = ct_poly_exp<N, R>();
-    static constexpr unsigned int order = R;                   // polynomial order
-    static constexpr unsigned int input_space_dimension = N;   // input space dimension
+    static constexpr int order = R;                   // polynomial order
+    static constexpr int input_space_dimension = N;   // input space dimension
 
     // constructor
     MultivariatePolynomial() = default;
@@ -202,7 +201,7 @@ class MultivariatePolynomial : public ScalarExpr<N, MultivariatePolynomial<N, R>
         // prepare gradient vector
         std::array<PolynomialDerivative<N, R>, N> gradient;
         // define i-th element of gradient field
-        for (size_t i = 0; i < N; ++i) { gradient[i] = PolynomialDerivative<N, R>(coeff_vector, i); }
+        for (std::size_t i = 0; i < N; ++i) { gradient[i] = PolynomialDerivative<N, R>(coeff_vector, i); }
         gradient_ = VectorField<N, N, PolynomialDerivative<N, R>>(gradient);
     };
     // evaluate polynomial at point
diff --git a/fdaPDE/linear_algebra.h b/fdaPDE/linear_algebra.h
index de0f32b9..f45c5a90 100644
--- a/fdaPDE/linear_algebra.h
+++ b/fdaPDE/linear_algebra.h
@@ -14,12 +14,12 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef __FDAPDE_LINEAR_ALGEBRA_MOUDLE_H__
-#define __FDAPDE_LINEAR_ALGEBRA_MOUDLE_H__
+#ifndef __FDAPDE_LINEAR_ALGEBRA_MODULE_H__
+#define __FDAPDE_LINEAR_ALGEBRA_MODULE_H__
 
 #include "linear_algebra/kronecker_product.h"
 #include "linear_algebra/smw.h"
 #include "linear_algebra/sparse_block_matrix.h"
 #include "linear_algebra/vector_space.h"
 
-#endif   // __FDAPDE_LINEAR_ALGEBRA_MOUDLE_H__
+#endif   // __FDAPDE_LINEAR_ALGEBRA_MODULE_H__
diff --git a/fdaPDE/linear_algebra/vector_space.h b/fdaPDE/linear_algebra/vector_space.h
index 2bfa21c5..95541fa3 100644
--- a/fdaPDE/linear_algebra/vector_space.h
+++ b/fdaPDE/linear_algebra/vector_space.h
@@ -27,7 +27,7 @@ namespace core {
 
 // a template class to perform geometric operations in general vector and affine spaces.
 // M is the vector space dimension, N is the dimension of the embedding space
-template <unsigned int M, unsigned int N> class VectorSpace {
+template <int M, int N> class VectorSpace {
    private:
     std::array<SVector<N>, M> basis_ {};   // the set of vectors generating the vector space
     SVector<N> offset_ {};   // a point throught which the vector space passes. (offset_ \neq 0 \implies affine space)
@@ -51,12 +51,12 @@ template <unsigned int M, unsigned int N> class VectorSpace {
 
 // implementation details
 
-template <unsigned int M, unsigned int N>
+template <int M, int N>
 SVector<N> VectorSpace<M, N>::orthogonal_project(const SVector<N>& u, const SVector<N>& v) const {
     return ((v.dot(u) / u.squaredNorm()) * (u.array())).matrix();
 }
 
-template <unsigned int M, unsigned int N> void VectorSpace<M, N>::orthonormalize() {
+template <int M, int N> void VectorSpace<M, N>::orthonormalize() {
     std::array<SVector<N>, M> orthonormal_basis;
     orthonormal_basis[0] = basis_[0] / basis_[0].norm();
     // implementation of the modified Gram-Schmidt method, see theory for details
@@ -70,7 +70,7 @@ template <unsigned int M, unsigned int N> void VectorSpace<M, N>::orthonormalize
 }
 
 // projects the point x onto the space.
-template <unsigned int M, unsigned int N> SVector<M> VectorSpace<M, N>::project_onto(const SVector<N>& x) {
+template <int M, int N> SVector<M> VectorSpace<M, N>::project_onto(const SVector<N>& x) {
     if constexpr (M == N) {   // you are projecting over the same space!
         return x;
     } else {
@@ -82,7 +82,7 @@ template <unsigned int M, unsigned int N> SVector<M> VectorSpace<M, N>::project_
 }
 
 // project the point x into the space
-template <unsigned int M, unsigned int N> SVector<N> VectorSpace<M, N>::project_into(const SVector<N>& x) {
+template <int M, int N> SVector<N> VectorSpace<M, N>::project_into(const SVector<N>& x) {
     // build the projection operator on the space spanned by the basis
     Eigen::Matrix<double, N, Eigen::Dynamic> A;
     A.resize(N, basis_.size());
@@ -93,12 +93,12 @@ template <unsigned int M, unsigned int N> SVector<N> VectorSpace<M, N>::project_
 }
 
 // euclidean distance from x
-template <unsigned int M, unsigned int N> double VectorSpace<M, N>::distance(const SVector<N>& x) {
+template <int M, int N> double VectorSpace<M, N>::distance(const SVector<N>& x) {
     return (x - project_into(x)).squaredNorm();
 }
 
 // develops the linear combination of basis vectors with respect to the given vector of coefficients.
-template <unsigned int M, unsigned int N>
+template <int M, int N>
 SVector<N> VectorSpace<M, N>::operator()(const std::array<double, M>& coeffs) const {
     SVector<N> result = offset_;
     for (std::size_t i = 0; i < M; ++i) result += coeffs[i] * basis_[i];
diff --git a/fdaPDE/mesh.h b/fdaPDE/mesh.h
index 4e930a1d..682acae4 100644
--- a/fdaPDE/mesh.h
+++ b/fdaPDE/mesh.h
@@ -14,8 +14,8 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef __FDAPDE_MESH_MOUDLE_H__
-#define __FDAPDE_MESH_MOUDLE_H__
+#ifndef __FDAPDE_MESH_MODULE_H__
+#define __FDAPDE_MESH_MODULE_H__
 
 #include "mesh/mesh.h"
 #include "mesh/element.h"
@@ -25,4 +25,4 @@
 #include "mesh/point_location/barycentric_walk.h"
 #include "mesh/point_location/adt.h"
 
-#endif   // __FDAPDE_MESH_MOUDLE_H__
+#endif   // __FDAPDE_MESH_MODULE_H__
diff --git a/fdaPDE/mesh/element.h b/fdaPDE/mesh/element.h
index 08636762..9353f497 100644
--- a/fdaPDE/mesh/element.h
+++ b/fdaPDE/mesh/element.h
@@ -29,24 +29,24 @@ namespace fdapde {
 namespace core {
 
 // M local dimension, N embedding dimension, R order of the element (defaulted to linear finite elements)
-template <unsigned int M, unsigned int N, unsigned int R = 1> class Element;   // forward declaration
+template <int M, int N, int R = 1> class Element;   // forward declaration
 
 // using declarations for manifold specialization
-template <unsigned int R> using SurfaceElement = Element<2, 3, R>;
-template <unsigned int R> using NetworkElement = Element<1, 2, R>;
+template <int R> using SurfaceElement = Element<2, 3, R>;
+template <int R> using NetworkElement = Element<1, 2, R>;
 
 // number of degrees of freedom associated to an M-dimensional element of degree R
-constexpr unsigned int ct_nnodes(const unsigned int M, const unsigned int R) {
+constexpr int ct_nnodes(const int M, const int R) {
     return ct_factorial(M + R) / (ct_factorial(M) * ct_factorial(R));
 }
 // number of vertices of an M-dimensional element
-constexpr unsigned int ct_nvertices(const unsigned int M) { return M + 1; }
+constexpr int ct_nvertices(const int M) { return M + 1; }
 // number of edges of an M-dimensional element
-constexpr unsigned int ct_nedges(const unsigned int M) { return (M * (M + 1)) / 2; }
+constexpr int ct_nedges(const int M) { return (M * (M + 1)) / 2; }
 
 // A single mesh element. This object represents the main **geometrical** abstraction of a physical element.
 // No functional information is carried by instances of this object.
-template <unsigned int M, unsigned int N, unsigned int R> class Element {
+template <int M, int N, int R> class Element {
    private:
     std::size_t ID_;                                         // ID of this element
     std::array<std::size_t, ct_nvertices(M)> node_ids_ {};   // ID of nodes composing the element
@@ -60,11 +60,11 @@ template <unsigned int M, unsigned int N, unsigned int R> class Element {
     SMatrix<N, M> J_;       // [J_]_ij = (coords_[j][i] - coords_[0][i])
     SMatrix<M, N> inv_J_;   // J^{-1} (Penrose pseudo-inverse for manifold elements)
    public:
-    static constexpr unsigned int nodes = ct_nnodes(M, R);
-    static constexpr unsigned int vertices = ct_nvertices(M);
-    static constexpr unsigned int local_dimension = M;
-    static constexpr unsigned int embedding_dimension = N;
-    static constexpr unsigned int order = R;
+    static constexpr int nodes = ct_nnodes(M, R);
+    static constexpr int vertices = ct_nvertices(M);
+    static constexpr int local_dimension = M;
+    static constexpr int embedding_dimension = N;
+    static constexpr int order = R;
 
     // constructor
     Element() = default;
@@ -74,7 +74,7 @@ template <unsigned int M, unsigned int N, unsigned int R> class Element {
     // getters
     const std::array<SVector<N>, ct_nvertices(M)> coords() const { return coords_; }
     const std::vector<int> neighbors() const { return neighbors_; }
-    const unsigned int ID() const { return ID_; }
+    const int ID() const { return ID_; }
     Eigen::Matrix<double, N, M> barycentric_matrix() const { return J_; }
     Eigen::Matrix<double, M, N> inv_barycentric_matrix() const { return inv_J_; }
     std::array<std::size_t, ct_nvertices(M)> node_ids() const { return node_ids_; }
@@ -99,7 +99,7 @@ template <unsigned int M, unsigned int N, unsigned int R> class Element {
 
 // implementation details
 
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 Element<M, N, R>::Element(
   std::size_t ID, const std::array<std::size_t, ct_nvertices(M)>& node_ids,
   const std::array<SVector<N>, ct_nvertices(M)>& coords, const std::vector<int>& neighbors, bool boundary) :
@@ -127,7 +127,7 @@ Element<M, N, R>::Element(
     }
 }
 
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 SVector<M + 1> Element<M, N, R>::to_barycentric_coords(const SVector<N>& x) const {
     // solve: barycenitrc_matrix_*z = (x-ref)
     SVector<M> z = inv_J_ * (x - coords_[0]);
@@ -137,7 +137,7 @@ SVector<M + 1> Element<M, N, R>::to_barycentric_coords(const SVector<N>& x) cons
     return result;
 }
 
-template <unsigned int M, unsigned int N, unsigned int R> SVector<N> Element<M, N, R>::mid_point() const {
+template <int M, int N, int R> SVector<N> Element<M, N, R>::mid_point() const {
     // The center of gravity of an element has all its barycentric coordinates equal to 1/(M+1).
     // The midpoint of an element can be computed by mapping it to a cartesian reference system
     SVector<M> barycentric_mid_point;
@@ -145,7 +145,7 @@ template <unsigned int M, unsigned int N, unsigned int R> SVector<N> Element<M,
     return J_ * barycentric_mid_point + coords_[0];
 }
 
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 std::pair<SVector<N>, SVector<N>> Element<M, N, R>::bounding_box() const {
     // define lower-left and upper-right corner of bounding box
     SVector<N> ll, ur;
@@ -164,14 +164,14 @@ std::pair<SVector<N>, SVector<N>> Element<M, N, R>::bounding_box() const {
 }
 
 // check if a point is contained in the element
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 template <bool is_manifold>
 typename std::enable_if<!is_manifold, bool>::type Element<M, N, R>::contains(const SVector<N>& x) const {
     // A point is inside the element if all its barycentric coordinates are positive
     return (to_barycentric_coords(x).array() >= -10 * std::numeric_limits<double>::epsilon()).all();
 }
 
-template <unsigned int M, unsigned int N, unsigned int R> VectorSpace<M, N> Element<M, N, R>::spanned_space() const {
+template <int M, int N, int R> VectorSpace<M, N> Element<M, N, R>::spanned_space() const {
     // build a basis for the space containing the element
     std::array<SVector<N>, M> basis;
     for (size_t i = 0; i < M; ++i) { basis[i] = (coords_[i + 1] - coords_[0]); }
@@ -179,7 +179,7 @@ template <unsigned int M, unsigned int N, unsigned int R> VectorSpace<M, N> Elem
 }
 
 // specialization for manifold elements of contains() routine.
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 template <bool is_manifold>
 typename std::enable_if<is_manifold, bool>::type Element<M, N, R>::contains(const SVector<N>& x) const {
     // check if the point is contained in the affine space spanned by the element
diff --git a/fdaPDE/mesh/mesh.h b/fdaPDE/mesh/mesh.h
index 69fabe5c..f0305f46 100644
--- a/fdaPDE/mesh/mesh.h
+++ b/fdaPDE/mesh/mesh.h
@@ -33,30 +33,30 @@ namespace fdapde {
 namespace core {
 
 // trait to detect if the mesh is a manifold
-template <unsigned int M, unsigned int N> struct is_manifold {
+template <int M, int N> struct is_manifold {
     static constexpr bool value = (M != N);
 };
 
 // trait to detect if the mesh is a linear network
-template <unsigned int M, unsigned int N> struct is_linear_network {
+template <int M, int N> struct is_linear_network {
     static constexpr bool value = std::conditional<(M == 1 && N == 2), std::true_type, std::false_type>::type::value;
 };
 
 // trait to select a proper neighboring storage structure depending on the type of mesh. In case of linear networks this
 // information is stored as a sparse matrix where entry (i,j) is set to 1 if and only if elements i and j are neighbors
-template <unsigned int M, unsigned int N> struct neighboring_structure {
+template <int M, int N> struct neighboring_structure {
     using type = typename std::conditional<is_linear_network<M, N>::value, SpMatrix<int>, DMatrix<int>>::type;
 };
 
 // access to domain's triangulation, M: tangent space dimension, N: embedding space dimension, R: FEM mesh order
-template <unsigned int M, unsigned int N, unsigned int R = 1> class Mesh {
+template <int M, int N, int R = 1> class Mesh {
    private:
     // coordinates of points costituting the vertices of mesh elements
     DMatrix<double> points_ {};
-    unsigned int num_nodes_ = 0;
+    int num_nodes_ = 0;
     // identifiers of points (as row indexes in points_ matrix) composing each element, by row
     Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> elements_ {};
-    unsigned int num_elements_ = 0;
+    int num_elements_ = 0;
     // vector of binary coefficients such that, boundary_[j] = 1 \iff node j is on boundary
     DMatrix<int> boundary_ {};
     std::size_t dof_;   // degrees of freedom, i.e. the maximmum ID in the dof_table
@@ -79,12 +79,12 @@ template <unsigned int M, unsigned int N, unsigned int R = 1> class Mesh {
          const typename neighboring_structure<M, N>::type& neighbors, const DMatrix<int>& boundary);
 
     // getters
-    const Element<M, N, R>& element(unsigned int ID) const { return cache_[ID]; }
-    Element<M, N, R>& element(unsigned int ID) { return cache_[ID]; }
-    SVector<N> node(unsigned int ID) const { return points_.row(ID); }
+    const Element<M, N, R>& element(int ID) const { return cache_[ID]; }
+    Element<M, N, R>& element(int ID) { return cache_[ID]; }
+    SVector<N> node(int ID) const { return points_.row(ID); }
     bool is_on_boundary(size_t j) const { return boundary_(j) == 1; }
-    unsigned int elements() const { return num_elements_; }
-    unsigned int nodes() const { return num_nodes_; }
+    int elements() const { return num_elements_; }
+    int nodes() const { return num_nodes_; }
     std::array<std::pair<double, double>, N> range() const { return range_; }
 
     // support for the definition of a finite element basis
@@ -154,17 +154,17 @@ template <unsigned int M, unsigned int N, unsigned int R = 1> class Mesh {
 };
 
 // alias exports
-template <unsigned int R = 1> using Mesh2D = Mesh<2, 2, R>;
-template <unsigned int R = 1> using Mesh3D = Mesh<3, 3, R>;
+template <int R = 1> using Mesh2D = Mesh<2, 2, R>;
+template <int R = 1> using Mesh3D = Mesh<3, 3, R>;
 // manifold cases
-template <unsigned int R = 1> using SurfaceMesh = Mesh<2, 3, R>;
-template <unsigned int R = 1> using NetworkMesh = Mesh<1, 2, R>;
+template <int R = 1> using SurfaceMesh = Mesh<2, 3, R>;
+template <int R = 1> using NetworkMesh = Mesh<1, 2, R>;
 
 // implementative details
 
 // builds a node enumeration for the support of a basis of order R. This fills both the elements_ table
 // and recompute the boundary informations. (support only for order 2 basis)
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 void Mesh<M, N, R>::dof_enumerate(const DMatrix<int>& boundary) {
     // algorithm initialization
     int next = num_nodes_;   // next valid ID to assign
@@ -213,7 +213,7 @@ void Mesh<M, N, R>::dof_enumerate(const DMatrix<int>& boundary) {
 }
 
 // construct directly from raw eigen matrix (used from wrappers)
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 Mesh<M, N, R>::Mesh(
   const DMatrix<double>& points, const DMatrix<int>& elements,
   const typename neighboring_structure<M, N>::type& neighbors, const DMatrix<int>& boundary) :
@@ -251,7 +251,7 @@ Mesh<M, N, R>::Mesh(
 }
 
 // construct a mesh from .csv files
-template <unsigned int M, unsigned int N, unsigned int R>
+template <int M, int N, int R>
 Mesh<M, N, R>::Mesh(
   const std::string& points, const std::string& elements, const std::string& neighbors, const std::string& boundary) {
     // open and parse CSV files
@@ -305,7 +305,7 @@ Mesh<M, N, R>::Mesh(
 }
 
 // fill the cache_ data structure with pointers to element objects
-template <unsigned int M, unsigned int N, unsigned int R> void Mesh<M, N, R>::fill_cache() {
+template <int M, int N, int R> void Mesh<M, N, R>::fill_cache() {
     // reserve space for cache
     cache_.reserve(num_elements_);
 
@@ -351,7 +351,7 @@ template <unsigned int M, unsigned int N, unsigned int R> void Mesh<M, N, R>::fi
 }
 
 // produce the matrix of dof coordinates
-template <unsigned int M, unsigned int N, unsigned int R> DMatrix<double> Mesh<M, N, R>::dof_coords() const {
+template <int M, int N, int R> DMatrix<double> Mesh<M, N, R>::dof_coords() const {
     if constexpr (R == 1)
         return points_;   // for order 1 meshes dofs coincide with vertices
     else {
diff --git a/fdaPDE/mesh/point_location/barycentric_walk.h b/fdaPDE/mesh/point_location/barycentric_walk.h
index 956ec224..cb0730c4 100644
--- a/fdaPDE/mesh/point_location/barycentric_walk.h
+++ b/fdaPDE/mesh/point_location/barycentric_walk.h
@@ -26,7 +26,7 @@ namespace fdapde {
 namespace core {
 
 // barycentric walk strategy for point location problem, works only for 2D and 3D *convex* triangualtions
-template <unsigned int M, unsigned int N, unsigned int R> class BarycentricWalk : public PointLocator<M, N, R> {
+template <int M, int N, int R> class BarycentricWalk : public PointLocator<M, N, R> {
     static_assert(M == N, "barycentric walk cannot be applied to manifold domains");
    private:
     const Mesh<M, N, R>& mesh_;
diff --git a/fdaPDE/mesh/point_location/naive_search.h b/fdaPDE/mesh/point_location/naive_search.h
index 03dcffc8..0ecee181 100644
--- a/fdaPDE/mesh/point_location/naive_search.h
+++ b/fdaPDE/mesh/point_location/naive_search.h
@@ -27,7 +27,7 @@ namespace fdapde {
 namespace core {
 
 // naive search strategy for point location problem
-template <unsigned int M, unsigned int N, unsigned int R> class NaiveSearch : public PointLocator<M, N, R> {
+template <int M, int N, int R> class NaiveSearch : public PointLocator<M, N, R> {
    private:
     const Mesh<M, N, R>& mesh_;
    public:
diff --git a/fdaPDE/mesh/point_location/point_locator.h b/fdaPDE/mesh/point_location/point_locator.h
index 0d6fbdb8..4894ff89 100644
--- a/fdaPDE/mesh/point_location/point_locator.h
+++ b/fdaPDE/mesh/point_location/point_locator.h
@@ -27,7 +27,7 @@ namespace fdapde {
 namespace core {
 
 // interface for point locations algorithms
-template <unsigned int M, unsigned int N, unsigned int R> struct PointLocator {
+template <int M, int N, int R> struct PointLocator {
     // solves the point location problem. returns nullptr if p is not found
     virtual const Element<M, N, R>* locate(const SVector<N>& p) const = 0;
 
diff --git a/fdaPDE/mesh/reference_element.h b/fdaPDE/mesh/reference_element.h
index 7f83d00f..64b2171c 100644
--- a/fdaPDE/mesh/reference_element.h
+++ b/fdaPDE/mesh/reference_element.h
@@ -25,8 +25,8 @@ namespace fdapde {
 namespace core {
 
 // Definition of the M-dimensional unit reference simplices of order R.
-template <unsigned int M, unsigned int R> struct ReferenceElement;
-template <unsigned int M, unsigned int R> using point_list = std::array<std::array<double, M>, R>;
+template <int M, int R> struct ReferenceElement;
+template <int M, int R> using point_list = std::array<std::array<double, M>, R>;
 
 template <>   // 1D first order basis
 struct ReferenceElement<1, 1> {
diff --git a/fdaPDE/optimization.h b/fdaPDE/optimization.h
index 5567a3f4..ae53eb12 100644
--- a/fdaPDE/optimization.h
+++ b/fdaPDE/optimization.h
@@ -14,8 +14,8 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef __FDAPDE_OPTIMIZATION_MOUDLE_H__
-#define __FDAPDE_OPTIMIZATION_MOUDLE_H__
+#ifndef __FDAPDE_OPTIMIZATION_MODULE_H__
+#define __FDAPDE_OPTIMIZATION_MODULE_H__
 
 #include "optimization/grid.h"
 #include "optimization/newton.h"
@@ -25,4 +25,4 @@
 #include "optimization/callbacks/backtracking_line_search.h"
 #include "optimization/callbacks/wolfe_line_search.h"
 
-#endif   // __FDAPDE_OPTIMIZATION_MOUDLE_H__
+#endif   // __FDAPDE_OPTIMIZATION_MODULE_H__
diff --git a/fdaPDE/optimization/gradient_descent.h b/fdaPDE/optimization/gradient_descent.h
index 7a186cb3..9387f848 100644
--- a/fdaPDE/optimization/gradient_descent.h
+++ b/fdaPDE/optimization/gradient_descent.h
@@ -25,7 +25,7 @@ namespace fdapde {
 namespace core {
 
 // implementation of the gradient descent method for unconstrained nonlinear optimization
-template <unsigned int N> class GradientDescent {
+template <int N> class GradientDescent {
    private:
     std::size_t max_iter_;   // maximum number of iterations before forced stop
     double tol_;             // tolerance on error before forced stop
diff --git a/fdaPDE/optimization/newton.h b/fdaPDE/optimization/newton.h
index 3676ef29..52cda5ac 100644
--- a/fdaPDE/optimization/newton.h
+++ b/fdaPDE/optimization/newton.h
@@ -25,7 +25,7 @@ namespace fdapde {
 namespace core {
 
 // implementation of the newton method for unconstrained nonlinear optimization
-template <unsigned int N> class Newton {
+template <int N> class Newton {
    private:
     std::size_t max_iter_;   // maximum number of iterations before forced stop
     double tol_;             // tolerance on error before forced stop
diff --git a/fdaPDE/pde/pde.h b/fdaPDE/pde/pde.h
index 61a99046..dcf1febf 100644
--- a/fdaPDE/pde/pde.h
+++ b/fdaPDE/pde/pde.h
@@ -53,9 +53,9 @@ template <typename D,   // triangulated problem domain
 class PDE : public PDEBase {
    public:
     typedef D DomainType;   // triangulated domain
-    static constexpr std::size_t M = DomainType::local_dimension;
-    static constexpr std::size_t N = DomainType::embedding_dimension;
-    static constexpr std::size_t R = DomainType::order;
+    static constexpr int M = DomainType::local_dimension;
+    static constexpr int N = DomainType::embedding_dimension;
+    static constexpr int R = DomainType::order;
     typedef E OperatorType;   // differential operator in its strong-formulation
     static_assert(
       std::is_base_of<DifferentialExpr<OperatorType>, OperatorType>::value, "E is not a valid differential operator");
diff --git a/fdaPDE/splines/basis/spline.h b/fdaPDE/splines/basis/spline.h
index 258da851..6d52ef4d 100644
--- a/fdaPDE/splines/basis/spline.h
+++ b/fdaPDE/splines/basis/spline.h
@@ -32,7 +32,7 @@ namespace core {
 
 // A spline of order R centered in knot u_i.
 // Template parameter M is used to keep track of the order of the spline while developing the Cox-DeBoor recursion.
-template <unsigned int R, unsigned int M = R> class Spline : public ScalarExpr<1, Spline<R>> {
+template <int R, int M = R> class Spline : public ScalarExpr<1, Spline<R>> {
    private:
     DVector<double> knots_;
     std::size_t i_;   // knot index where this basis is centered
@@ -57,7 +57,7 @@ template <unsigned int R, unsigned int M = R> class Spline : public ScalarExpr<1
     // compute derivative of order K as a ScalarExpr
     // d^K/dx^K N_ij(x) = j/(u_i+j - u_i)*[d^{K-1}/dx^{K-1} N_i,j-1(x)]  - j/(u_i+j+1 - u_i+1)*
     // [d^{K-1}/dx^{K-1} N_i+1,j-1(x)]
-    template <unsigned int K> auto derive() const {
+    template <int K> auto derive() const {
         if constexpr (K == 1)   // end of recursion
             return (R * a_) * Spline<R - 1, M>(knots_, i_) - (R * b_) * Spline<R - 1, M>(knots_, i_ + 1);
         else   // exploit Cox-DeBoor recursive formula
@@ -67,7 +67,7 @@ template <unsigned int R, unsigned int M = R> class Spline : public ScalarExpr<1
 };
 
 // partial template specialization for order 0 splines (end of recursion)
-template <unsigned int M> class Spline<0, M> : public ScalarExpr<1, Spline<0>> {
+template <int M> class Spline<0, M> : public ScalarExpr<1, Spline<0>> {
    private:
     DVector<double> knots_;
     std::size_t i_;   // knot index where this basis is centered
diff --git a/fdaPDE/splines/basis/spline_basis.h b/fdaPDE/splines/basis/spline_basis.h
index f4ef8f33..28d5af3b 100644
--- a/fdaPDE/splines/basis/spline_basis.h
+++ b/fdaPDE/splines/basis/spline_basis.h
@@ -24,7 +24,7 @@ namespace fdapde {
 namespace core {
 
 // a spline basis of order R built over a given set of knots
-template <unsigned int R> class SplineBasis {
+template <int R> class SplineBasis {
    private:
     DVector<double> knots_ {};   // vector of knots
     std::vector<Spline<R>> basis_ {};
diff --git a/fdaPDE/splines/spline_assembler.h b/fdaPDE/splines/spline_assembler.h
index a89aa7e4..0090e3e7 100644
--- a/fdaPDE/splines/spline_assembler.h
+++ b/fdaPDE/splines/spline_assembler.h
@@ -45,7 +45,7 @@ template <typename D, typename B, typename I> class Assembler<SPLINE, D, B, I> {
 template <typename D, typename B, typename I>
 template <typename E>
 SpMatrix<double> Assembler<SPLINE, D, B, I>::discretize_operator(const E& op) {
-    constexpr std::size_t R = B::order;
+    constexpr int R = B::order;
     std::size_t M = basis_.size();
     std::vector<fdapde::Triplet<double>> triplet_list;
     SpMatrix<double> discretization_matrix;
diff --git a/fdaPDE/utils.h b/fdaPDE/utils.h
index c3e726a9..7e20d6d1 100644
--- a/fdaPDE/utils.h
+++ b/fdaPDE/utils.h
@@ -14,8 +14,8 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef __FDAPDE_UTILS_MOUDLE_H__
-#define __FDAPDE_UTILS_MOUDLE_H__
+#ifndef __FDAPDE_UTILS_MODULE_H__
+#define __FDAPDE_UTILS_MODULE_H__
 
 #include "utils/assert.h"
 #include "utils/compile_time.h"
@@ -27,4 +27,4 @@
 #include "utils/integration/integrator.h"
 #include "utils/integration/integrator_tables.h"
 
-#endif   // __FDAPDE_UTILS_MOUDLE_H__
+#endif   // __FDAPDE_UTILS_MODULE_H__
diff --git a/fdaPDE/utils/compile_time.h b/fdaPDE/utils/compile_time.h
index 797f5dd6..6dd36d7d 100644
--- a/fdaPDE/utils/compile_time.h
+++ b/fdaPDE/utils/compile_time.h
@@ -23,15 +23,15 @@ namespace fdapde {
 namespace core {
 
 // compile time evaluation of the factorial function
-constexpr std::size_t ct_factorial(const std::size_t n) { return n ? (n * ct_factorial(n - 1)) : 1; }
+constexpr int ct_factorial(const int n) { return n ? (n * ct_factorial(n - 1)) : 1; }
 
 // compile time evaluation of binomial coefficient
-constexpr std::size_t ct_binomial_coefficient(const std::size_t N, const std::size_t M) {
+constexpr int ct_binomial_coefficient(const int N, const int M) {
     return ct_factorial(N) / (ct_factorial(M) * ct_factorial(N - M));
 }
 
 // compile time computation of the sum of elements in array
-template <typename T, std::size_t N>
+template <typename T, int N>
 constexpr typename std::enable_if<std::is_arithmetic<T>::value, T>::type ct_array_sum(std::array<T, N> A) {
     T result = 0;
     for (int j = 0; j < N; ++j) result += A[j];
diff --git a/fdaPDE/utils/integration/integrator.h b/fdaPDE/utils/integration/integrator.h
index 90bac37e..3f9b0396 100644
--- a/fdaPDE/utils/integration/integrator.h
+++ b/fdaPDE/utils/integration/integrator.h
@@ -29,32 +29,32 @@ namespace core {
 
 // A set of utilities to perform numerical integration
 // M: dimension of the domain of integration, R finite element order, K number of quadrature nodes
-template <unsigned int M, unsigned int R, unsigned int K = standard_fem_quadrature_rule<M, R>::K> class Integrator {
+template <int M, int R, int K = standard_fem_quadrature_rule<M, R>::K> class Integrator {
    private:
     IntegratorTable<M, K> integration_table_;
    public:
     Integrator() : integration_table_(IntegratorTable<M, K>()) {};
 
     // integrate a callable F over a mesh element e
-    template <unsigned int N, typename F> double integrate(const Element<M, N, R>& e, const F& f) const;
+    template <int N, typename F> double integrate(const Element<M, N, R>& e, const F& f) const;
     // integrate a callable F over a triangualtion m
-    template <unsigned int N, typename F> double integrate(const Mesh<M, N, R>& m, const F& f) const;
+    template <int N, typename F> double integrate(const Mesh<M, N, R>& m, const F& f) const;
     // computes \int_e [f * \phi] where \phi is a basis function over the *reference element*.
-    template <unsigned int N, typename F, typename B>
+    template <int N, typename F, typename B>
     double integrate(const Element<M, N, R>& e, const F& f, const B& phi) const;
     // integrate the weak form of operator L to produce its (i,j)-th discretization matrix element
-    template <typename L, unsigned int N, typename F> double integrate(const Element<M, N, R>& e, F& f) const;
+    template <typename L, int N, typename F> double integrate(const Element<M, N, R>& e, F& f) const;
 
     // getters
-    template <unsigned int N> DMatrix<double> quadrature_nodes(const Mesh<M, N, R>& m) const;
+    template <int N> DMatrix<double> quadrature_nodes(const Mesh<M, N, R>& m) const;
     std::size_t num_nodes() const { return integration_table_.num_nodes; }
 };
 
 // implementative details
 
 // integration of bilinear form
-template <unsigned int M, unsigned int R, unsigned int K>
-template <typename L, unsigned int N, typename F>
+template <int M, int R, int K>
+template <typename L, int N, typename F>
 double Integrator<M, R, K>::integrate(const Element<M, N, R>& e, F& f) const {
     // apply quadrature rule
     double value = 0;
@@ -73,8 +73,8 @@ double Integrator<M, R, K>::integrate(const Element<M, N, R>& e, F& f) const {
 // perform integration of \int_e [f * \phi] using a basis system defined over the reference element and the change of
 // variables formula: \int_e [f(x) * \phi(x)] = \int_{E} [f(J(X)) * \Phi(X)] |detJ|
 // where J is the affine mapping from the reference element E to the physical element e
-template <unsigned int M, unsigned int R, unsigned int K>
-template <unsigned int N, typename F, typename B>
+template <int M, int R, int K>
+template <int N, typename F, typename B>
 double Integrator<M, R, K>::integrate(const Element<M, N, R>& e, const F& f, const B& Phi) const {
     double value = 0;
     for (size_t iq = 0; iq < integration_table_.num_nodes; ++iq) {
@@ -94,8 +94,8 @@ double Integrator<M, R, K>::integrate(const Element<M, N, R>& e, const F& f, con
 }
 
 // integrate a callable F over a mesh element e. Do not require any particular structure for F
-template <unsigned int M, unsigned int R, unsigned int K>
-template <unsigned int N, typename F>
+template <int M, int R, int K>
+template <int N, typename F>
 double Integrator<M, R, K>::integrate(const Element<M, N, R>& e, const F& f) const {
     double value = 0;
     for (size_t iq = 0; iq < integration_table_.num_nodes; ++iq) {
@@ -115,8 +115,8 @@ double Integrator<M, R, K>::integrate(const Element<M, N, R>& e, const F& f) con
 }
 
 // integrate a callable F over the entire mesh m.
-template <unsigned int M, unsigned int R, unsigned int K>
-template <unsigned int N, typename F>
+template <int M, int R, int K>
+template <int N, typename F>
 double Integrator<M, R, K>::integrate(const Mesh<M, N, R>& m, const F& f) const {
     double value = 0;
     // cycle over all mesh elements
@@ -125,8 +125,8 @@ double Integrator<M, R, K>::integrate(const Mesh<M, N, R>& m, const F& f) const
 }
 
 // returns all quadrature points on the mesh
-template <unsigned int M, unsigned int R, unsigned int K>
-template <unsigned int N>
+template <int M, int R, int K>
+template <int N>
 DMatrix<double> Integrator<M, R, K>::quadrature_nodes(const Mesh<M, N, R>& m) const {
     DMatrix<double> quadrature_nodes;
     quadrature_nodes.resize(m.elements() * integration_table_.num_nodes, N);
diff --git a/fdaPDE/utils/integration/integrator_tables.h b/fdaPDE/utils/integration/integrator_tables.h
index fe5f75f0..67528a3f 100644
--- a/fdaPDE/utils/integration/integrator_tables.h
+++ b/fdaPDE/utils/integration/integrator_tables.h
@@ -13,15 +13,15 @@ namespace core {
 // ** "https://people.sc.fsu.edu/~jburkardt/datasets/datasets.html"
 
 // N: dimension of the integration domain, K: number of nodes of the formula. Args: type of quadrature formula
-template <unsigned int N, unsigned int K, typename... Args> struct IntegratorTable;
+template <int N, int K, typename... Args> struct IntegratorTable;
 // identification tags for quadrature rules
 struct NewtonCotes { };
 struct GaussLegendre { };
 
 // trait for selecting a standard quadrature rule for finite elements, in case K is not supplied
-template <unsigned int N, unsigned int R>
+template <int N, int R>
 struct standard_fem_quadrature_rule {
-    static constexpr unsigned int quadrature(const unsigned int dim, const unsigned int order) {
+    static constexpr int quadrature(const int dim, const int order) {
         switch (dim) {
         case 1:   // 1D elements
             switch (order) {
@@ -53,7 +53,7 @@ struct standard_fem_quadrature_rule {
         }
         return 0;   // error
     }
-    static constexpr unsigned int K = quadrature(N, R);
+    static constexpr int K = quadrature(N, R);
 };
 
 // 1D linear elements (gaussian integration)
diff --git a/fdaPDE/utils/symbols.h b/fdaPDE/utils/symbols.h
index adc9db49..f1a5d008 100644
--- a/fdaPDE/utils/symbols.h
+++ b/fdaPDE/utils/symbols.h
@@ -22,8 +22,8 @@
 #include <memory>
 
 // static structures, allocated on stack at compile time.
-template <unsigned int N> using SVector = Eigen::Matrix<double, N, 1>;
-template <unsigned int N, unsigned int M = N> using SMatrix = Eigen::Matrix<double, N, M>;
+template <int N> using SVector = Eigen::Matrix<double, N, 1>;
+template <int N, int M = N> using SMatrix = Eigen::Matrix<double, N, M>;
 
 // dynamic size, head-appocated, structures.
 template <typename T> using DMatrix    = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
@@ -72,7 +72,7 @@ struct matrix_hash {
 };
 
 // oredering relation for SVector<N>, allows SVector<N> to be keys of std::map
-template <unsigned int N> struct s_vector_compare {
+template <int N> struct s_vector_compare {
     bool operator()(const SVector<N>& lhs, const SVector<N>& rhs) const {
         return std::lexicographical_compare(lhs.data(), lhs.data() + lhs.size(), rhs.data(), rhs.data() + rhs.size());
     };